Erosion and sediment flux#

Stream power law#

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

In goSPL, both detachment and transport-limited approaches are 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. In the transport-limited case, the concentration of sediment in the flow hinders the erosive ability of the flow.

In both cases, the sediment erosion rate is expressed using a stream power formulation function of river discharge, precipitation, slope and upstream sediment flux (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:

  1. weathering tends to weaken bedrocks, increasing their erodibility.

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

Detachment-limited#

Here we illustrate the approach implemented in goSPL based on the detachment-limited case with linear dependency on slope.

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

\[\mathrm{E} = \mathrm{\kappa P^d Q^m S^n}\]

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 by default to 0.5 and 1 respectively (but could also be tuned).

In the detachment-limited case with default values of \(\mathrm{m}\) and \(\mathrm{n}\), the elevation (\(\mathrm{\eta_i}\)) will change due to local river erosion rate and is defined implicitly by:

\[\mathrm{\frac{\eta_i^{t+\Delta t}-\eta_i^t}{\Delta t}} = \mathrm{-\kappa P^d_i \sqrt{Q_i} \frac{\eta_i^{t+\Delta t} - \eta_{rcv}^{t+\Delta t}}{\lambda_{i,rcv}}}\]

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

\[\mathrm{(1+K_{f,i|rcv}) \eta_i^{t+\Delta t} - K_{f,i|rcv} \eta_{i,rcv}^{t+\Delta t}} = \mathrm{\eta_i^t}\]

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:

\[\mathrm{\boldsymbol\Gamma \boldsymbol\eta} \mathrm{^{t+\Delta t}} = \mathrm{\boldsymbol\eta} \mathrm{^{t}}\]

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:

\[\begin{split}\begin{align} \begin{bmatrix} \mathrm{\gamma_{1,1}} & - \mathrm{\gamma_{1,2}} & & - \mathrm{\gamma_{1,4}} & - \mathrm{\gamma_{1,5}} & & & & & \\ & \mathrm{\gamma_{2,2}} & - \mathrm{\gamma_{2,3}} & - \mathrm{\gamma_{2,4}} & & & & & & \\ & & \mathrm{\gamma_{3,3}} & & & & & -\mathrm{\gamma_{3,8}} & -\mathrm{\gamma_{3,9}} & \\ & & - \mathrm{\gamma_{4,3}} & \mathrm{\gamma_{4,4}} & - \mathrm{\gamma_{4,5}} & - \mathrm{\gamma_{4,6}}& & -\mathrm{\gamma_{4,8}} & & \\ & & && \mathrm{\gamma_{5,5}} & - \mathrm{\gamma_{5,6}} & & & & \\ & & & && \mathrm{\gamma_{6,6}} & & -\mathrm{\gamma_{6,8}} & & \\ & & & & & - \mathrm{\gamma_{7,6}} & \mathrm{\gamma_{7,7}} & -\mathrm{\gamma_{7,8}} & & \\ & & & & & & &\mathrm{\gamma_{8,8}} & -\mathrm{\gamma_{8,9}} & \\ & & & & & & & & 1 & \\ & & & & & & - \mathrm{\gamma_{10,7}} & -\mathrm{\gamma_{10,8}} & -\mathrm{\gamma_{10,9}} & \mathrm{\gamma_{10,10}} \end{bmatrix} \begin{bmatrix} \mathrm{\eta_1^{t+\Delta t}} \\ \mathrm{\eta_2^{t+\Delta t}} \\ \mathrm{\eta_3^{t+\Delta t}} \\ \mathrm{\eta_4^{t+\Delta t}} \\ \mathrm{\eta_5^{t+\Delta t}} \\ \mathrm{\eta_6^{t+\Delta t}} \\ \mathrm{\eta_7^{t+\Delta t}} \\ \mathrm{\eta_8^{t+\Delta t}} \\ \mathrm{\eta_9^{t+\Delta t}} \\ \mathrm{\eta_{10}^{t+\Delta t}} \end{bmatrix} = \begin{bmatrix} \mathrm{\eta_1^{t}} \\ \mathrm{\eta_2^{t}} \\ \mathrm{\eta_3^{t}} \\ \mathrm{\eta_4^{t}} \\ \mathrm{\eta_5^{t}} \\ \mathrm{\eta_6^{t}} \\ \mathrm{\eta_7^{t}} \\ \mathrm{\eta_8^{t}} \\ \mathrm{\eta_9^{t}} \\ \mathrm{\eta_{10}^{t}} \end{bmatrix}\label{eq:SPACE8} \end{align}\end{split}\]

with

\[\begin{split}\begin{array}{ccl} \mathrm{\gamma_{i,j}} & = & \mathrm{w_{i,j}K_{f,i|j}} \qquad i \ne j \\ \mathrm{\gamma_{1,1}} & = & 1 + \mathrm{\sum\limits_{j=2,4,5} w_{1,j}K_{f,1|j} } \\ \mathrm{\gamma_{2,2}} & = &1 + \mathrm{\sum\limits_{j=3,4} w_{2,j}K_{f,2|j} } \\ \mathrm{\gamma_{3,3}} & = & 1 + \mathrm{\sum\limits_{j=8,9} w_{3,j}K_{f,3|j} } \\ \mathrm{\gamma_{4,4}} & = & 1 + \mathrm{\sum\limits_{j=3,5,6,8} w_{4,j}K_{f,4|j} }\\ \mathrm{\gamma_{5,5}} & = & 1 + \mathrm{ w_{5,6}K_{f,5|6} } \\ \mathrm{\gamma_{6,6}} & = & 1 + \mathrm{ w_{6,8}K_{f,6|8}} \\ \mathrm{\gamma_{7,7}} & = & 1 + \mathrm{ \sum\limits_{j=6,8} w_{7,j}K_{f,7|j}} \\ \mathrm{\gamma_{8,8}} & = & 1 + \mathrm{w_{8,9}K_{f,8|9}} \\ \mathrm{\gamma_{10,10}} & = & 1 + \mathrm{\sum\limits_{j=7,8,9} w_{10,j}K_{f,10|j}} \end{array}\end{split}\]

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.

../_images/ero.png

Fig. 2 Flow accumulation patterns and associated erosion based on a radially symmetric surface defined with a central, high region and a series of distal low-lying valleys. Resulting topography after 100,000 years of evolution under uniform precipitation for the multiple flow direction algorithms. Patterns of flow accumulation after 20,000 and 50,000 years are presented as well as estimated landscape erosion at the end of the simulation time.#

Transport-limited#

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. Here we illustrate the corresponding system of equation assuming a linear dependency on slope (\(\mathrm{m}\) and \(\mathrm{n}\) are set by default to 0.5 and 1 respectively).

\[\mathrm{\frac{\eta_i^{t+\Delta t}-\eta_i^t}{\Delta t}} = \mathrm{-\kappa P^d_i \sqrt{Q_i} \frac{\eta_i^{t+\Delta t} - \eta_{rcv}^{t+\Delta t}}{\lambda_{i,rcv}}} + \mathrm{G' Q_{s_i} / \Omega_i}\]

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:

\[\mathrm{Q_{t_i}^{t+\Delta t} - \sum_{ups} w_{i,j} Q_{t_u}^{t+\Delta t}}= \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Omega_i}{\Delta t}}\]

which gives:

\[\mathrm{Q_{s_i}} = \mathrm{Q_{t_i}} - \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Omega_i}{\Delta t}}\]

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).

Dynamic soil#

goSPL could also simulate dynamic soil production and transport. In this case, the implementation tracks a layer of regolith (defined as unconsolidated and potentially mobile sediment, such as soil or alluvium) in combination with sediment entrainment–deposition erosion law.

Note

In this case, the SPACE model of Shobe et al. (2017) is used in place of entrainment–deposition law. As described in the BasicHySa governing equations from Terrainbento (Appendix B20 from Barnhart et al. (2019)).

Following Barnhart et al. (2019)), the governing equations can be summarized as:

\[ \begin{align}\begin{aligned}\mathrm{\eta} = \mathrm{H + \eta_b}\\\mathrm{\frac{\partial H}{ \partial t}} = \mathrm{P_0 e^{- H/H_s} - \kappa_s Q^m S^n (1-e^{-H/H_\star}) - \nabla \cdot {q_{h}} + G' Q_{s_i} / \Omega_i}\\\mathrm{\frac{\partial \eta_b}{ \partial t}} = -\mathrm{P_0 e^{- H/H_s} - \kappa_b Q^m S^n e^{-H/H_\star}}\end{aligned}\end{align} \]

with \(H\) the soil thickness and \(\eta_b\) the bedrock elevation. \(\mathrm{P_0}\) is the soil production maximum rate, \(\mathrm{H_s}\) the soil production decay depth, and \(\mathrm{H_\star}\) the roughness length scale. In addition we differentiate erodibility coefficients \(\kappa_b\) and \(\kappa_s\) for bedrock and soil respectively.

The top soil layer is subject to diffusion \(\mathrm{q_{h}}\) and follows:

\[\mathrm{q_{h}} = \mathrm{-D \left(1 - e^{-H/H_0} \right) \nabla \cdot {\eta}}\]

with \(\mathrm{H_0}\) the soil transport decay depth for non-linear diffusion and \(\mathrm{D}\) the coefficient of diffusion.

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{Q_s^{out}} = \mathrm{Q_s^{in} + E \Omega}\]

\(\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:

\[\mathrm{q_{s,i}} = \mathrm{e_{i} + \sum_{d=1}^{N_d}q_{s,d}}\]

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:

\[\begin{split}\begin{align} \begin{bmatrix} 1 & & & & & & & & & \\ \mathrm{-w_{1,2}} & 1 & & & & & & & & \\ & \mathrm{-w_{2,3}} & 1 & \mathrm{-w_{4,3}} & & & & & & \\ \mathrm{-w_{1,4}} & \mathrm{-w_{2,4}} & & 1 & & & & & & \\ \mathrm{-w_{1,5}} & & & \mathrm{-w_{4,5}} & 1 & & & & & \\ & & & \mathrm{-w_{4,6}} & \mathrm{-w_{5,6}} & 1 & \mathrm{-w_{7,6}} & & & \\ & & & & & & 1 & & & \mathrm{-w_{10,7}}\\ & & \mathrm{-w_{3,8}} & \mathrm{-w_{4,8}} & & \mathrm{-w_{6,8}} & \mathrm{-w_{7,8}} & 1 & & \mathrm{-w_{10,8}} \\ & & \mathrm{-w_{3,9}} & & & & & \mathrm{-w_{8,9}} & 1 & \mathrm{-w_{10,9}} \\ & & & & & & & & & 1 \end{bmatrix} \begin{bmatrix} \mathrm{q_{s,1}} \\ \mathrm{q_{s,2}} \\ \mathrm{q_{s,3}} \\ \mathrm{q_{s,4}} \\ \mathrm{q_{s,5}} \\ \mathrm{q_{s,6}} \\ \mathrm{q_{s,7}} \\ \mathrm{q_{s,8}} \\ \mathrm{q_{s,9}} \\ \mathrm{q_{s,10}} \end{bmatrix} = \begin{bmatrix} \mathrm{e_1} \\ \mathrm{e_2} \\ \mathrm{e_3} \\ \mathrm{e_4} \\ \mathrm{e_5} \\ \mathrm{e_6} \\ \mathrm{e_7} \\ \mathrm{e_8} \\ \mathrm{e_9} \\ \mathrm{e_{10}} \end{bmatrix} \end{align}\end{split}\]

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.