2.1.1 Elevation change

Assuming negligible uplift, the change in elevation over time is described by the continuity of mass equation expressed as the divergence of sediment flow (Tucker et al., 2001):

where z is elevation [m], t is time [s], q_s is sediment flow per unit width (vector) [kg m⁻¹ s⁻¹], d_s is the net erosion–deposition rate [kg m⁻² s⁻¹], and ρ_s is sediment mass density [kg m⁻³].

In r.sim.terrain, the net erosion–deposition rate d_s driven by overland flow is estimated at different levels of complexity based on the simulation mode selected by the user. Gravitational diffusion is then applied to the changed topography to simulate the smoothing effects of localized soil transport between rainfall events. The change in elevation due to gravitational diffusion is a function of the diffusion coefficient and the Laplacian of elevation (Thaxton, 2004):

where ε_g is the diffusion coefficient [m² s⁻¹].

The discrete implementation follows Thaxton (2004):

where Δz_s is elevation change [m] caused by net erosion or deposition during time interval Δt₁ (Eq. 2), and Δz_g is the diffusion-driven elevation change [m] during time interval Δt₁ (Eq. 3).

2.1.2 Erosion–deposition regimes

Following experimental observations and qualitative arguments, Foster et al. (1977) proposed that the sum of the ratio of the net erosion–deposition rate d_s to the detachment capacity D_c [kg m⁻² s⁻¹] and the ratio of the sediment flow rate $q_{\mathrm{s}}=\left|{\mathit{q}}_{\mathrm{s}}\right|$ to the sediment transport capacity T_c [kg m⁻¹ s⁻¹] is a conserved quantity (unity):

The net erosion and deposition rate d_s can then be expressed as being proportional to the difference between the sediment transport capacity T_c and the actual sediment flow rate q_s:

This principle is used in several erosion models including the Water Erosion Prediction Project (WEPP) (Flanagan et al., 2013) and SIMWE (Mitas and Mitasova, 1998).

Using this concept, it is possible to identify two limiting erosion–deposition regimes. When T_c≫D_c leading to T_c≫q_s, the erosion regime is detachment capacity limited and net erosion is equal to the detachment capacity:

For this case, the transport capacity of overland flow exceeds the detachment capacity, and thus sediment flow, erosion, and sediment transport are limited by the detachment capacity. Therefore, no deposition occurs. An example of this case is when a strong storm producing intense overland flow over compacted clay soils causes high-capacity flows to transport light clay particles, while the detachment of compacted soils is limited. When D_c≫T_c, sediment flow is at sediment transport capacity q_s=T_c, leading to a transport-capacity-limited regime with deposition reaching its maximum extent for the given water flow. Net erosion–deposition is computed as the divergence of transport capacity multiplied by a unit vector s₀ in the direction of flow:

This case may occur, for example, during a moderate storm with overland flow over sandy soils with high detachment capacity but low transport capacity. For $\mathrm{0}<(D_{\mathrm{c}}/T_{\mathrm{c}})<\mathrm{\infty }$, the spatial pattern of net erosion–deposition is variable and depends on the difference between the sediment transport capacity and the actual sediment flow rate at the given location.

The detachment capacity D_c and the sediment transport capacity T_c are estimated using shear stress and stream power equations, respectively, expressed as power functions of water flow properties and slope angle. The relations between the topographic parameters of well-known empirical equations for erosion modeling, such as the Universal Soil Loss Equation (USLE) and stream power, were presented by Moore and Burch (1986) and used to develop simple, GIS-based models for limiting erosion–deposition cases such as RUSLE3D and USPED (Mitasova and Mitas, 2001). The SIMWE model estimates T_c and D_c using modified equations and parameters developed for the WEPP model (Flanagan et al., 2013; Mitasova et al., 2013).

The simulation modes in r.sim.terrain include (Fig. 2)

• the process-based SIMWE model for steady-state and unsteady shallow overland flow in variable erosion–deposition regimes with d_s computed by solving the shallow water flow and sediment transport continuity equations,

• the RUSLE3D model for detachment-capacity-limited cases with d_s given by Eq. (8), and

• the USPED model for transport-capacity-limited regimes with d_s given by Eq. (9).

The following sections explain the computation of d_s for these three modes in more detail.

2.2.1 Shallow water flow

The SIMWE model simulates shallow overland water flow controlled by spatially variable topographic, soil, land cover, and rainfall parameters using a Green function Monte Carlo path sampling method. The steady-state shallow water flow continuity equation relates the change in water depth across space to source, defined in our case as rainfall excess rate:

where q is the water flow per unit width (vector) [m² s⁻¹], h is the depth of overland flow [m], v is the water flow velocity vector [m s⁻¹] whose magnitude is computed with Manning's equation $v=n^{-\mathrm{1}}{h}^{\mathrm{2}/\mathrm{3}}{s}^{\mathrm{1}/\mathrm{2}}$, n is Manning's coefficient [s m${}^{-\mathrm{1}/\mathrm{3}}$], s is slope (unitless), and i_e is the rainfall excess rate [m s⁻¹] (i.e., rainfall intensity − infiltration − vegetation intercept).

An approximation of diffusive wave effects is incorporated by adding a diffusion term proportional to ∇²[h^5∕3]:

where ε_w is a spatially variable diffusion coefficient [m^4∕3 s⁻¹].

The path sampling method solves the continuity equation for h^5∕3 through the accumulation of the evolving source (Mitasova et al., 2004). The solution assumes that water flow velocity is largely controlled by the slope of the terrain and surface roughness and that its change at a given location during the simulated event is negligible. The initial number of particles per grid cell is proportional to the rainfall excess rate i_e (source). The water depth h^5∕3 at time τ during the simulated rainfall event is computed as a function of particle (walkers) density at each grid cell. Particles are routed across the landscape by finding a new position for each walker at time τ+Δτ:

where $\mathit{r}=(x,y)$ is the mth walker position [m], Δτ is the particle routing time step [s], and g is a random vector with Gaussian components with variance Δτ [m].

The mathematical background of the method, including the computation of the temporal evolution of water depth and incorporation of approximate momentum through an increased diffusion rate in the prevailing direction of flow, is presented by Mitas and Mitasova (1998) and Mitasova et al. (2004).

2.2.2 Sediment flow and net erosion–deposition

The SIMWE model simulates the sediment flow over complex topography with spatially variable overland flow, soil, and land cover properties by solving the sediment flow continuity equation using a Green function Monte Carlo path sampling method. Steady-state sediment flow q_s is approximated by the bivariate continuity equation, which relates the change in sediment flow rate to effective sources and sinks:

The sediment flow rate q_s is a function of water flow and sediment concentration (Mitas and Mitasova, 1998):

where ρ_s is sediment mass density in the water column [kg m⁻³], c is sediment concentration [particle m⁻³], and ϱ=ρ_sch is the mass of sediment transported by water per unit area [kg m⁻²].

The sediment flow equation (Eq. 13), like the water flow equation, has been rewritten to include a small diffusion term that is proportional to the mass of water-carried sediment per unit area ∇²ϱ (Mitas and Mitasova, 1998):

where ε_s is the diffusion constant [m² s⁻¹].

On the left-hand side of Eq. (15), the first term describes local diffusion, the second term is drift driven by water flow, and the third term represents a velocity-dependent “potential” acting on the mass of transported sediment. The initial number of particles per grid cell is proportional to the soil detachment capacity D_c (source). The particles are then routed across the landscape by finding a new position for each walker at time τ+Δτ:

while the updated weight is

where $u=\left(D_{\mathrm{c}}/T_{\mathrm{c}}\right)\left|\mathit{v}\right|$.

Sediment flow is computed as the product of weighted particle densities and the water flow velocity (Eq. 14), and the net erosion–deposition rate d_s is computed as the divergence of sediment flow using Eq. (13). See Mitas and Mitasova (1998) and Mitasova et al. (2004) for more details on the Green function Monte Carlo solution and equations for computing D_c and T_c.

This model can simulate erosion regimes from prevailing detachment-limited conditions when T_c≫D_c to prevailing transport-capacity-limited conditions when D_c≫T_c and the erosion–deposition patterns between these conditions. At each landscape evolution time step, the regime can change based on the ratio between the sediment detachment capacity D_c and the sediment transport capacity T_c and the actual sediment flow rate. If the landscape evolution time step is shorter than the time to concentration (i.e., the time for water to reach steady state), then net erosion–deposition is derived from unsteady flow.

2.3.1 Erosivity factor

The erosivity factor R in USLE and RUSLE is the combination of the total energy and peak intensity of a rainfall event, representing the interaction between the detachment of sediment particles and the transport capacity of the flow. It can be calculated as the product of the the kinetic energy of the rainfall event E and its maximum 30 min intensity I₃₀ (Brown and Foster, 1987; Renard et al., 1997; Panagos et al., 2015, 2017). In this model, however, the erosivity factor is derived at each time step as a function of kinetic energy, rainfall depth, rainfall intensity, and time. First, rain energy is derived from rainfall intensity (Brown and Foster, 1987; Yin et al., 2017):

where e_r is unit rain energy [MJ ha⁻¹ mm⁻¹], i_r is rainfall intensity [mm h⁻¹], b is empirical coefficient, i₀ is reference rainfall intensity [mm h⁻¹], and e₀ is reference energy [MJ ha⁻¹ mm⁻1]. The parameters for this equation were derived from observed data published for different regions by Panagos et al. (2017). Then the event-based erosivity index R_e is calculated as the product of unit rain energy, rainfall depth, rainfall intensity, and time:

where R_e is the event-based erosivity index [MJ mm ha⁻¹ h⁻¹], v_r is the rainfall depth [mm] derived from v_r=i_rΔt, and Δt is the change in time [s].

2.3.2 Flow accumulation

The upslope contributing area per unit width a is determined by flow accumulation (the number of grid cells draining into a given grid cell) multiplied by grid cell width (Fig. 3d). Flow accumulation is calculated using a multiple flow direction algorithm (Metz et al., 2009) based on A^T least-cost path searches (Ehlschlaeger, 1989). The multiple flow direction algorithm implemented in GRASS GIS, as the module r.watershed is computationally efficient, does not require sink filling and can navigate nested depressions and other obstacles.

Figure 3Water and sediment flows modeled with spatially variable land cover for Patterson Branch, Fort Bragg, NC: (a) water depth [m] simulated by SIMWE for a 10 min event with 50 mm h⁻¹ in the subwatershed; (b) flow accumulation for RUSLE3D in the subwatershed; (c) erosion and deposition [kg m⁻² s⁻¹] simulated by SIMWE in drainage area 1; and (d) erosion [kg m⁻² s⁻¹] modeled by RUSLE3D in drainage area 1.

Download

2.3.3 3-D topographic factor

The 3-D topographic factor LS_3-D is calculated as a function of the upslope contributing area and the slope (Fig. 3e):

where LS_3-D is the dimensionless topographic factor, a is upslope contributing area per unit width [m], a₀ is the length of the standard USLE plot [22.1 m], β is the angle of the slope [^∘], m is an empirical coefficient, n is an empirical coefficient, and β₀ is the slope of the standard USLE plot [5.14^∘].

The empirical coefficients m and n for the upslope contributing area and the slope can range from 0.2 to 0.6 and 1.0 to 1.3, respectively, with low values representing dominant sheet flow and high values representing dominant rill flow.

2.3.4 Detachment-limited erosion rate

The erosion rate is a function of the event-based erosivity factor, soil erodibility factor, 3-D topographic factor, land cover factor, and prevention measures factor (Fig. 3d):

where E is soil erosion rate (soil loss) [kg m⁻² min⁻¹], R_e is the event-based erosivity factor [MJ mm ha⁻¹ h⁻¹], K is the soil erodibility factor [t ha h ha⁻¹ MJ⁻¹ mm⁻¹], LS_3-D is the dimensionless topographic (length–slope) factor, C is the dimensionless land cover factor, and P is the dimensionless prevention measures factor.

The detachment-limited erosion represented by RUSLE3D leads to the simulated change in elevation:

which is combined with Eq. (3) for gravitational diffusion.

2.4.1 Topographic sediment transport factor

Using the unit stream power concept presented by Moore and Burch (1986), the 3-D topographic factor (Eq. 20) for RUSLE3D is modified to represent the topographic sediment transport factor (LS_T) – the topographic component of overland flow at sediment transport capacity:

where LS_T is the topographic sediment transport factor, a is the upslope contributing area per unit width [m], β is the angle of the slope [^∘], m is an empirical coefficient, and n is an empirical coefficient.

2.4.2 Transport-limited sediment flow and net erosion–deposition

Sediment flow at transport capacity is a function of the event-based rainfall factor, soil erodibility factor, topographic sediment transport factor, land cover factor, and prevention measures factor:

where T is sediment flow at transport capacity [kg m⁻¹ s⁻¹], R_e is the event-based rainfall factor [MJ mm ha⁻¹ h⁻¹], K is the soil erodibility factor [t ha h ha⁻¹ MJ⁻¹ mm⁻¹], C is the dimensionless land cover factor, and P is the dimensionless prevention measures factor.

Net erosion–deposition is estimated as the divergence of sediment flow, assuming that sediment flow is equal to sediment transport capacity:

where d_s is net erosion–deposition [kg m⁻² s⁻¹], α is the aspect of the topography (i.e., the direction of flow) [^∘]. With USPED, the simulated change in elevation Δz_s=d_s is derived from Eq. (2) for landscape evolution and then Eq. (3) for gravitational diffusion.