# River Discharge¶

## Flow accumulation¶

Flow accumulation (FA) calculations are core component of landscape evolution models as they are often used as proxy to estimate flow discharge, sediment load, river width, bedrock erosion as well as sediment deposition.

Note

Until recently conventional FA algorithms were **serial** and limited to small spatial problems. With ever growing high resolution digital elevation dataset, new methods based on **parallel** approaches have been proposed over the last decade.

In addition, nearly all of these parallel approaches assume a **single flow direction** (SFD). This assumption makes the emergent flow network highly sensitive to the underlying mesh geometry and most dendritic shape of obtained stream networks is often an artefact of the surface triangulation. To reduce this effect, authors have proposed to consider not only the steepest downhill direction but also to represent other directions appropriately weighted by slope (**multiple flow direction** - MFD). Using MFD algorithms prevent locking of erosion pathways along a single direction and help to route flow over flat regions into multiple branches.

## Single and multiple flow directions¶

`gospl`

allows for both SFD and MFD routing by using an adapted version of the parallel implicit drainage area (IDA) method from Richardson et al. (2014) to unstructured meshes. It consists in writing the FA calculation as a *sparse matrix system of linear equations* and takes full advantage of purpose-built, efficient linear algebra routines available in PETSc.

The river discharge is computed from the calculated FA and the net precipitation rate \(\mathrm{P}\). At node \(\mathrm{i}\), the river discharge (\(\mathrm{q_i}\)) is determined as follows:

where \(\mathrm{b_i}\) is the local volume of water \(\mathrm{\Omega_i P_i}\) where \(\mathrm{\Omega_i}\) is the voronoi area and \(\mathrm{P_i}\) the local precipitation value available for runoff during a given time step. \(\mathrm{N_d}\) is the number of donors with a donor defined as a node that drains into \(\mathrm{i}\) (as an example the donor of vertex 5 in the SFD sketch in the above figure is 1). To find the donors of each node, the method consists in finding their receivers first. Then, the receivers of each donor is saved into a receiver matrix, noting that the nodes, which are local minima, are their own receivers.

The transpose of the matrix is then used to get the donor matrix. When the previous equation is applied to all nodes and considering the MFD case illustrated above, the following relations are obtained:

The choice of weights \(\mathrm{w_{m,n}}\) depends on the number of flow directions that is used. The weights range between zero and one and sum to one for each node:

The number of flow direction paths is user-defined and can vary from 1 (*i.e.* SFD) up to 6 (*i.e.* MFD) depending of the grid neighbourhood complexity. The weights are calculated based on the number of downslope neighbours and are proportional to the slope.

## Linear solver¶

In matrix form the system defined above is equivalent to **W q** = **b** or:

The vector **q** corresponds to the unknown river discharge (volume of water flowing on a given node per year) and the elements of **W** left blank are zeros.

Note

As explained in Richardson et al. (2014), the above system is implicit as the river discharge for a given vertex depends on its neighbours unknown flow discharge. The matrix **W** is sparse and is composed of diagonal terms set to unity (identity matrix) and off-diagonal terms corresponding to at most the immediate neighbours of each vertex (typically lower than 6 in constrained Delaunay triangulation).

In `gospl`

, this matrix is built in parallel using compressed sparse row matrix functionality available from SciPy.

Once the matrix has been constructed, PETSc library is used to solve matrices and vectors across the decomposed domain. The performance of the **IDA** algorithm is strongly dependent on the choice of solver and preconditioner. In `gospl`

, the solution for **q** is obtained using the *Richardson solver* with block Jacobi preconditioning (*bjacobi*). This choice was made based on convergence results but can be changed if better solver and preconditioner combinations are found.

Iterative methods allow for an initial guess to be provided. When this initial guess is close to the solution, the number of iterations required for convergence dramatically decreases. This option is used in `gospl`

by allocation the river discharge solution from previous time step as an initial guess. It allows to decrease the number of iterations of the IDA solver as discharge often exhibits small change between successive time intervals.