# Erosion rate and sediment flux#

## Stream power law#

River incision, associated sediment transport and subsequent deposition are critical elements of landscape evolution models.

In goSPL, a detachment-limited approach is implemented. The detachment-limited hypothesis supposes that erosion is not limited by a transport capacity but instead by the ability of rivers to remove material from the bed.
The sediment erosion rate is expressed using a **stream power** formulation function of river discharge, precipitation and slope (Howard et al., 1994).

Note

The proposed formulation differs from traditional ones in that it relates explicitly incision to precipitation following the approach described in Murphy et al. (2016). This modification allows to account for local climatic control of bedrock river incision induced by chemical weathering. Weathering is known to strongly contribute to erosion and to largely influence sediment load distributed to downstream catchments. Studies have shown that the physical strength of bedrock which varies with the degree of chemical weathering, increases systematically with local rainfall rate. The proposed formulation expresses in a simple mathematical form the following two feedbacks between chemical weathering and erosion:

weathering tends to weaken bedrocks, increasing their erodibility.

enhanced erosion rates exposes fresh rocks at the surface therefore contributing to increased weathering.

The volumetric entrainment flux of sediment per unit bed area \(\mathrm{E}\) is of the following form:

where \(\mathrm{\kappa}\) is the precipitation-independent sediment erodibility parameter, \(\mathrm{P}\) the annual precipitation rate, \(\mathrm{d}\) is a positive exponent, \(\mathrm{Q}=\bar{P}A\) is the water discharge (computed in the River Discharge section with \(\mathrm{A}\) the flow accumulation) and \(\mathrm{S}\) is the river local slope. \(\mathrm{m}\) and \(\mathrm{n}\) are scaling exponents. In goSPL, \(\mathrm{\kappa}\) is user defined and the coefficients \(\mathrm{m}\) and \(\mathrm{n}\) are set to 0.5 and 1 respectively. \(\mathrm{E}\) is in \(\mathrm{m/y}\) and therefore the erodibility dimension is \(\mathrm{m\,y^{-0.5}}\).

The elevation (\(\mathrm{\eta_i}\)) will change due to local river erosion rate and is defined implicitly by:

where \(\mathrm{\lambda_{i,rcv}}\) is the length of the edges connecting the considered vertex to its receiver. Rearranging the above equation gives:

with the coefficient \(\mathrm{K_{f,i|rcv} = \kappa P^d_i \sqrt{Q_i} \Delta t / \lambda_{i,rcv}}\).

## Matrix formalism#

In matrix form the system defined above is equivalent to:

Using the example presented in the figure 1 in section River Discharge and considering the **multiple flow direction** scenario, the matrix system based on the receivers distribution is defined as:

with

The above system is **implicit** and the matrix is **sparse**. The SciPy compressed sparse row matrix functionality is used here again to build \(\mathrm{\boldsymbol\Gamma}\) on each local domain. The SciPy matrix format (*e.g.* csr_matrix) is efficiently loaded as a PETSc Python matrix and the system is then solved using *Richardson solver* with block Jacobi preconditioning (*bjacobi*) using an initial guess for the solution set to vertices elevation.

## Sediment entrainment#

Once the erosion rates have been obtained, the sediment flux moving out at every node \(\mathrm{Q_s^{out}}\) equals the flux of sediment flowing in plus the local erosion rate. \(\mathrm{Q_s^{out}}\) takes the following form:

\(\mathrm{\Omega}\) is the voronoi area of the considered vertex.

The solution of the above equation requires the calculation of the incoming sediment volume from upstream nodes \(\mathrm{Q_s^{in}}\). At node \(\mathrm{i}\), this equation is equivalent to:

where \(\mathrm{e_{i} = E_{i} \Omega_i}\) and \(\mathrm{N_d}\) is the number of donors. Assuming that river sediment concentration is distributed in a similar way as the water discharge we write the following set of equalities for our example:

It is worth noting that in this system, the matrix **W** is the same as the one proposed for the River Discharge and therefore does not have to be built. As for the previous system, this one is solved using the PETSc solver previously defined to find the \(\mathrm{q_{s,i}}\) values implicitly.

## SPL with sediment deposition#

An alternative method to the detachment-limited approach proposed above consists in accounting for the role played by sediment in modulating erosion and deposition rates. It follows the model of Yuan et al, 2019, whereby the deposition flux depends on a deposition coefficient \(G\) and is proportional to the ratio between cell area \(\mathrm{\Omega}\) and water discharge \(\mathrm{Q}=\bar{P}A\).

The approach considers the local balance between erosion and deposition and is based on sediment flux resulting from net upstream erosion.

where \(\mathrm{\lambda_{i,rcv}}\) is the length of the edges connecting the considered vertex to its receiver and \(\mathrm{\Omega_i}\) is the area (voronoi) of the node \(i\).

\(\mathrm{Q_{s_i}}\) is the upstream incoming sediment flux in m3/yr and \(\mathrm{G'}\) is equal to \(\mathrm{G \Omega_i / \bar{P}A}\).

The upstream incoming sediment flux is obtained from the total sediment flux \(\mathrm{Q_{t_i}}\) where:

which gives:

This system of coupled equations is solved implicitly using PETSc by assembling the matrix and vectors using the nested submatrix and subvectors and by using the `fieldsplit`

preconditioner combining two separate preconditioners for the collections of variables.

The `TFQMR`

(transpose-free QMR (quasi minimal residual)) KSP solver is used to solve the coupled system with sub KSPs set to `preonly`

and preconditioner set to `hypre`

. (See PETSC documentation for more details about the solver and preconditoner options and settings).