Soil erosion processes

The SPHY model allows to model soil erosion with 6 different soil erosion models, i.e. MUSLE (1), MMF (2), INCA (3), SHETRAN (4), DHSVM (5) and HSPF (6). The MUSLE model is an empirical model, which is forced by accumulated runoff, as generated by the hydrological part of the SPHY model. The implementation of MUSLE in SPHY was part of study in which three different soil erosion model concepts were compared (Eekhout & de Vente, 2020). All other models are process-based models that determine detachment of soil particles separately for raindrop impact and accumulated runoff. Subsequently, these two different detachment processes are summed and sediment taken into transport is determined accounting for immediate deposition. The first of these 5 process-based models that was implemented was the MMF model (Eekhout et al., 2018). The other 4 process-based soil erosion models were part of soil erosion model ensemble, with the aim to assess the uncertainty of process-based soil erosion models in climate change impact studies (Eekhout et al., 2021).

All soil erosion models make use of model parameters related to the infiltration excess surface runoff equation. It is therefore advisable to use this equation, i.e. Infil_excess = 1. Furthermore, most of the process-based soil erosion models use the LAI (leaf area index) from the vegetation module to determine the canopy cover, hence, it is also advisable to use the vegetation module when applying one of the process-based soil erosion models.

MUSLE

The Modified Universal Soil Loss Equation (MUSLE) is a modification of the USLE, where the rainfall erosivity factor is replaced by a runoff factor, and applied at a daily time step. MUSLE is incorporated in various widely used hydrological models, such as SWAT, in which a separate hydrological module is used to calculate runoff. MUSLE is determined as follows (Williams, 1995):

equation 59

$$sed = 11.8(Q_{surf} \cdot q_{peak} \cdot A)^{0.56} K \cdot LS \cdot C\cdot P \cdot CFRG$$

Where $sed$ is the sediment yield $(kg m^{-2}day^{-1})$, $Q_{surf}$ is the surface runoff depth (mm), $q_{peak}$ is the peak runoff rate $(m^3 s^{-1})$, A is the cell area $(m^2)$, K is the soil erodibility factor $(kg h MJ^{-1} mm^{-1})$, LS is the topographic factor (-), C is the crop and management factor (-), P is the erosion control practice factor (-) and CFRG is the coarse fragment factor (-).

The surface runoff $Q_{surf}$ is determined within SPHY as described in other parts of the manual. The peak runoff rate is determined as follows:

equation 60

$$q_{peak}=\frac{\alpha_{tc}Q_{surf}A}{3.6t_{conc}}$$

Where $\alpha_{tc}$ is the fraction of daily rainfall that occurs during the time of concentration (-) and $t_{conc}$ is the time of concentration (hr). The time of concentration is the amount of time from the beginning of a rainfall event until the entire cell area is contributing to flow at the cell outlet.

The fraction of daily rainfall that occurs during the time of concentration is determined as follows:

equation 61

$$\alpha_{tc}=1-e^{2 t_{conc} ln(1-\alpha_{0.5})}$$

Where $\alpha_{0.5}$ is the fraction of the daily rain falling in the half-hour highest intensity (-), which is obtained from a model parameter of the infiltration excess surface runoff equation and can be determined within the calibration of the hydrological model.

The time of concentration is determined accounting for both channel flow and overland flow:

equation 62

$$t_{conc}=\frac{t_{conc,ch}+t_{conc,ov}}{60}$$

Where $t_{conc,ch}$ is the channel flow time of concentration (minutes) and $t_{conc,ov}$ the overland flow time of concentration (minutes). The channel flow time of concentration ($t_{conc,ch}$) is determined using the Kirpich (1940) method:

equation 63

$$t_{conc,ch}=0.0195L^{0.77}S^{-0.385}$$

Where L is the slope length (m) and S the slope $(m\cdot m^{-1})$.

The overland flow time of concentration $t_{conc,ov}$ is determined using the Kerby (1959) method:

equation 64

$$t_{conc,ov}=1.44(LN)^{0.467}S^{-0.235}$$

Where N is the retardance coefficient (-).

The soil erodibility factor was determined using the equation developed by Wischmeier et al. (1971):

equation 65

$$K=\frac{0.00021M^{1.14}(12-OM)+3.25(c_{soilstr}-2)+2.5(c_{perm}-3)}{100}$$

Where K is the particle-size parameter (-), is the organic matter content (%), $c_{soilstr}$ is the soil structure class (-) and $c_{perm}$ is the profile permeability class (-).

The particle-size parameter is calculated as follows:

equation 66

$$M=(m_{silt}+m_{vfs})(100-m_{c})$$

Where $m_{silt}$ is the silt content (%), $m_{vfs}$ is the very fine sand content (%) and $m_{c}$ is the clay content (%). The profile permeability classes are defined according to the saturated hydraulic conductivity.

The coarse fragment factor is determined as follows:

equation 67

$$CFRG=e^{-0.053p_{rock}}$$

Where $p_{rock}$ is the rock content in the root zone layer (%).

The topographic factor LS is the expected ratio of soil loss per unit area from a field slope of 22.1 m length with uniform slope of 9%. We applied the following equation (Wischmeier et al., 1971):

equation 68

$$LS = \Big(\frac{L}{22.1}\Big)^{m}(65.41 sin^2(\alpha_{hill})+0.0065)$$

Where m is an exponential term (-) and $\alpha_{hill}$ is the slope angle (°). The exponential term m is calculated as follows:

equation 69

$$m=0.6(1-e^{-35.835S})$$

MMF

The Morgan-Morgan-Finney model (Morgan and Duzant, 2008a) was originally implemented as an annual model, however, Eekhout et al. (2018) included MMF as a daily model in SPHY. In MMF, the detachment by raindrop impact is a function of the highest daily precipitation intensity and the canopy cover, which are obtained from the infiltration excess surface runoff equation and the vegetation module, respectively. Detachment by runoff is a function of the accumulated runoff. Both detachment by raindrop impact and runoff are determined for each of the three textural classes (sand, silt and clay) separately and later aggregated to determine the total detachment. Immediate deposition is a function of the particle fall number, in which the flow velocity is determined with the Manning equation. A detailed explanation is given below.

Detachment by raindrop impact

Detachment by raindrop impact $(F, kg m^{-2})$ is determined for each of the soil texture classes separately and subsequently summed and is calculated as follows:

Equation 70

$$F_i=K_i \frac{\%i}{100}(1-GC)KE \cdot 10^{-3}$$

Where F is the detachment by raindrop impact for textural class $i \space (kg m^{-2})$, $i$ the textural class, K the detachability of the soil by raindrop impact $(gJ^{-1})$, GC the ground cover (-) and KE the kinetic energy of the effective precipitation $(J m^{-2})$. The ground cover GC, expressed as a proportion between zero and unity, protects the soil from detachment and is determined by the proportion of vegetation and rocks covering the surface. The ground cover is set to 1 in case the surface is covered with snow, which is determined from the snow module.

The total kinetic energy of the effective precipitation (KE) is calculated as follows:

Equation 71

$$KE=KE_{LD}+KE_{DT}$$

Where $KE_{LD}$ is the kinetic energy of the leaf drainage $(J m^{-2})$ and $KE_{DT}$ is the kinetic energy of the direct throughfall $(J m^{-2})$.

The kinetic energy of the leaf drainage $KE_{LD}$ is based on (Brandt, 1990):

Equation 72

$$KE_{LD}= \begin{cases} 0 &\text{for } PH<0.15 \\ LD(15.8 \sqrt{PH}-5.87) &\text{for } PH \ge 0.15 \end{cases}$$

Where LD is the leaf drainage (mm) and PH is the plant height (m), specified for each landuse class.

Leaf drainage is determined as:

$$LD=P_{eff}CC$$

Equation 73

Where $P_{eff}$ is the precipitation throughfall (mm) and CC is the canopy cover (fraction, -). The canopy cover is either introduced by a landuse-class specific tabular value or determined by the vegetation module. When the vegetation module is used, the canopy cover is obtained from the LAI (Equation 6), maximized by 1. The effective precipitation from the hydrological model is first corrected for the slope angle, following (Choi et al., 2017):

Equation 74

$$P_{eff}=P_{eff}\space cos \space S$$

Where $P_{eff}$ is the effective precipitation (mm) and S is the slope (°).

The kinetic energy of the direct throughfall is based on a relationship described by Brown and Foster (1987):

Equation 75

$$KE_{DT}=DT \Big(0.29(1-0.72e^{-0.05I})\Big)\cdot 100$$

Where DT is the direct throughfall (mm) and I is the intensity of the erosive precipitation $(mm \space h^{-1})$. The intensity of the erosive precipitation is a model parameter and varies according to geographical location.

Direct throughfall (DT) is calculated as follows:

Equation 76

$$DT=P_{eff}-LD$$

Detachment by runoff

Detachment by runoff (H) is calculated as follows:

equation 77

$$H_i=DR_i \frac{\%i}{100}Q^{1.5}(1-GC)sin^{0.3}S \cdot10^{-3}$$

Where H is the detachment by runoff $(kg \space m^{-2})$, DR the detachability of the soil by runoff $(g \space mm^{-1})$, Q is the volume of accumulated runoff (mm) and S is the slope $(m \space m^{-1})$.

Sediment transported

A proportion of the detached soil is deposited in the cell of its origin as a function of the abundance of vegetation and the surface roughness. The percentage of the detached sediment that is deposited (DEP) is estimated from the relationship obtained by (Tollner et al., 1976) and calculated separately for each texture class:

equation 78

$$DEP=44.1N_{fi}^{0.29}$$

Where $N_f$ is the particle fall number (-), defined as:

Equation 79

$$N_{fi}=\frac{lv_{si}}{vd}$$

Where l is the length of a grid cell (m), $v_s$ the particle fall velocity $(m\space s^{-1})$, $v$ the flow velocity $(m\space s^{-1})$ and d the depth of flow (m).

The particle fall velocities $v_s$ are estimated from:

Equation 80

$$v_s=\frac{\frac{1}{18}\delta^2(\rho_s+\rho)g}{\eta}$$

Where $\delta$ is the diameter of the particle (m), the sediment density $(kg \space m^{-3})$ , the flow density $(kg \space m^{-3})$ (Abrahams et al., 2001), gravitational acceleration $(m \space s^{2})$ and the fluid viscosity $(kg \space m^{-1} \space s^{-1})$.

The flow velocity $v$ from Equation 78 is obtained by the Manning formula:

equation 81

$$v=\frac{1}{n'}d^{\frac{2}{3}}S^{\frac{1}{2}}$$

Where n' is the modified Manning's roughness coefficient $(s \space m^{-\frac{1}{3}})$, which is a combination of the Manning's roughness coefficient for the soil surface and vegetation, defined as (Petryk and Bosmajian, 1975):

Equation 82

$$n'=\sqrt{n^2_{soil}+n^2_{vegetation}}$$

Where $n_{soil}$ is the Manning's roughness coefficient for soil $(kg \space m^{-3})$ and the Manning's roughness coefficient for vegetation $(kg \space m^{-3})$. The Manning’s roughness coefficient for soil can either be defined by bare soil (Figure 4a) or tilled soil (Figure 4b). The Manning’s roughness coefficient for vegetation can either be obtained for regular spaced vegetation (Figure 4c) or irregular spaced vegetation (Figure 4d).

For tilled conditions (Figure 4b) the following equation is applied to obtain the Manning's roughness coefficient for the soil:

Equation 83

$$n_{soil}=exp(-2.1132+0.0349RFR)$$

Where RFR is the surface roughness parameter $(cm \space m^{-1})$.

The Manning's roughness coefficient for regular spaced vegetation (Figure 4c) is obtained from the following equation (Jin et al., 2000):

equation 84

$$n_{vegetation}=\frac{d^{\frac{2}{3}}}{\sqrt{\frac{2g}{D \space NV}}}$$

Where D is the stem diameter (m) and NV the stem density $(stems \space m^{-2})$.

Equations 79, 81 and 84 require a flow depth d, a model parameter that can be used in the model calibration. The value for d should be taken such that it corresponds to a water depth from runoff generated within the cell margins, i.e. without accumulation of flow from upstream located cells.

Figure 4. Surface and vegetation roughness options: (a) bare soil, (b) tilled soil (Equation 68), (c) regular vegetation (Equation 69), and (d) irregular vegetation.

The amount of sediment that is taken into transport is determined from the sum of the detached sediment from raindrop impact ($F_i$; Equation 70) and runoff ($H_i$; equation 77), subtracting the proportion of the sediment that is deposited within the cell of its origin ($DEP_i$; equation 78):

Equation 85

$$G_i=(F_i+H_i) \Big(1-\frac{DEP_i}{100}\Big)$$

Where G is the amount of sediment taken into transport for textural class $i \space (kg \space m^{-2})$. The amount of sediment that is routed to downstream cells is the summation of the individual amounts for clay, silt and sand.

INCA

The Integrated Catchments model for Sediments (INCA-Sed; Lazar et al., 2010) is originally applied in a semi-distributed manner, however, here the model is implemented in a spatially distributed manner. Detachment by raindrop impact is a function of the daily precipitation intensity and the canopy cover, for which the latter is obtained from the vegetation module. For model calibration purposes, we included the ground cover as a model parameter in the detachment by raindrop impact formulation. Detachment by runoff is a function of the sediment transport, the surface runoff and the detachment by raindrop impact. Sediment that is taken into transport is determined from the before mentioned formulations, accounting for sediment storage.

Detachment by raindrop impact

Detachment by raindrop impact $(S_{sp})$ is calculated as follows:

equation 86

$$S_{sp}=(1-C_g)C_{X1}\space p_{Sed}E_{sp}^{(\frac{10}{10-V})}8.64 \cdot 10^{10}$$

Where $S_{sp}$ the detachment by raindrop impact $(kg \space km^{-2})$, $C_g$ is the ground cover (-), $C_{X1}$ a scaling parameter $(s\space m^{-1})$, $p_{Sed}$ the precipitation throughfall (mm), $E_{sp}$ a soil specific erosion potential parameter $(kg \space m^{-2} \space s^{-1})$ and V the vegetation cover (-), here estimated with the canopy cover from the vegetation module, multiplied by 10.

Detachment by runoff

Detachment by runoff $(S_{FL})$ is calculated as follows:

equation 87

$$S_{FL}=\frac{K(S_{TC}-S_{SP})}{S_{TC}+K}$$

Where $S_{FL}$ is the detachment by runoff $(kg\space km^{-2})$, K a function of runoff $(kg\space km^{-2})$ and the transport capacity $(kg\space km^{-2})$.

The function K is calculated as follows:

equation 88

$$K=a_1E_{FL}\Big(\frac{A\space q_{DR}}{L}-a_2\Big)^{a_s}\cdot 86400$$

Where $E_{FL}$ is the soil erosion potential $(kg \space km^{-2}\space s^{-1})$, A the grid cell area $(km^2)$, $q_{DR}$ the routed runoff $(m^3 \space s^{-1} km^{-2})$, L the slope length (km), $a_1$ is the flow erosion scaling factor $(s \space m^{-2})$, $a_2$ the flow erosion direct runoff threshold $(m^2\space s^{-1})$ and $a_3$ the flow erosion non-linear coefficient (-).

Sediment transport capacity$(S_{TC})$ is calculated as follows:

equation 89

$$S_{TC}=a_4\Big(\frac{A \space q_{DR}}{L}-a_5\Big)^{a_6} \cdot 86400$$

Where $a_4$ is the transport capacity scaling factor $(kg \space m^{-2} \space km^{-2})$, $a_5$ the transport capacity direct runoff threshold $(m^2 \space s^{-1})$ and $a_6$ the transport capacity non-linear coefficient (-).

Sediment transported

The amount of sediment that is taken into transport depends on the amount of sediment in the sediment storage. The daily change in sediment storage is calculated as follows:

equation 90

$$\frac{dS_{store}}{dt}=\begin{cases} S_{SP}-S_{TC} &\text{for } S_{store}+S_{SP}>S_{TC} \\ \frac{-K(S_{SP}-S_{TC})}{S_{TC}+K} &\text{for } S_{store}+S_{SP}<S_{TC} \end{cases}$$

Where $S_{store}$ is the sediment storage $(kg \space km^{-2})$, which is subsequently updated following:

equation 91

$$S_{store}=S_{store}+\frac{dS_{store}}{dt}$$

The amount of sediment taken into transport is calculated as follows:

equation 92

$$M_{out}=\begin{cases} S_{TC} &\text{for } S_{store}+S_{SP}>S_{TC} \\ S_{SP}+S_{FL} &\text{for } S_{store}+S_{SP}<S_{TC} \end{cases}$$

Where $M_{out}$ is the mass of sediment taken into transport $(kg \space km^{-2})$.

SHETRAN

The SHETRAN model (Lukey et al., 1995) is a sediment transport model implemented in the Système Hydrologique Européen (SHE) hydrological model. The detachment by raindrop impact formulations are similar to the ones used in DHSVM, with some small differences in the leaf drip formulations. The canopy cover is obtained from the vegetation module. Detachment by runoff is a function of the shear stress and critical shear stress, which are both a function of the water depth. We obtained the water depth using the Manning equation, assuming a triangular-shaped flow profile, with the width-depth ratio as model parameter. For model calibration purposes, we included the ground cover as a model parameter in the detachment by runoff formulation. Immediate deposition of sediment is determined with a sediment transport equation.

Detachment by raindrop impact

Detachment by raindrop impact is determined with the following empirical equation, which is derived from Wicks (1988):

equation 93

$$D_r=k_rF_w(1-C_g-C_r)(M_r+M_d)$$

where $D_r$ is the rate of detachment of soil $(kg\space m^{-2}\space s^{-1})$, $k_r$ the raindrop impact soil erodibility coefficient $(J^{-1})$, $F_w$ the protective effect of ponding (-), $C_g$ the proportion of ground shielded by near ground cover (fraction, -), $C_r$ the proportion of ground shielded by ground level (rock) cover (fraction, -), $M_r$ the momentum squared of raindrops reaching the ground $(kg^2\space s^{-3})$ and $M_d$ the momentum squared of leaf drip reaching the ground $(kg^2\space s^{-3})$.

The original SHETRAN model accounts for the protective effect of ponding on detachment by raindrop impact by model parameter $F_w$(Park et al., 1982). The hydrological model SPHY does not account for ponding, hence, we assume $F_w=1$.

The momentum squared of raindrops reaching the ground $(M_r)$ is based on the formulations by Marshall & Palmer (1948):

equation 94

$$M_r=(1-C_c)a_1I^{b_1}$$

Where $I$ is the rainfall intensity $(mm\space h^{-1})$ and $a_1$ and $b_1$ are coefficients dependent on $I$ and are given in Table 2. The rainfall intensity is obtained from the infiltration excess surface runoff equation and $a_1$ and $b_1$ are determined inside the model code.

Table 2: Values for the empirical coefficients $a_1$ and $b_1$ used to determine the momentum squared of raindrops.

Range for I (mm h-1)	a1	b1
0 - 10	2.6893 ∙ 10-8	1.6896
10 - 50	3.7514 ∙ 10-8	1.5545
50 - 100	6.1192 ∙ 10-8	1.4242
≥ 100	11.737 ∙ 10-8	1.2821

The momentum squared of leaf drip reaching the ground $(M_d)$ is calculated as follows:

equation 95

$$M_d=\frac{\pi}{6}V^2_d\rho^2d^3_lL_d\space DRAINA$$

Where $V_d$ is the leaf drip fall speed $(m\space s^{-1})$, $\rho$ the density of water $(kg \space m^{-3})$, $d_l$ the leaf drip diameter (m), $L_d$ the proportion of drainage that falls as leaf drip (fraction, -) and $DRAINA$ the water drainage rate from canopy$(m \space s^{-1})$. The proportion of drainage that falls as leaf drip $(L_d)$ is assumed to be equal to the canopy cover $(C_c)$. The water drainage rate from canopy $(DRAINA)$ is assumed to be equal to the daily precipitation intensity in $m \space s^{-1}$.

The leaf drip fall speed $(V_d)$ is calculated as follows:

equation 96

$$V_d=\sqrt{\frac{M}{\beta}g\Big(1-e^{\frac{2X\beta}{M}}\Big)}$$

Where M the average mass of leaf drips (kg), $\beta$ the friction constant $(kg \space m^{-1})$, g the acceleration due to gravity $(m \space s^{-2})$ and X is the average leaf drip fall distance (m).

The fraction $\frac{M}{\beta}$ is a function of the leaf drip diameter $d_l$ and two coefficients, $a_2$ and $b_2$.

equation 97

$$\frac{M}{\beta}=a_2+b_2d_l$$

where $a_2$ and $b_2$ are given in Table 3 and are determined inside the model code.

Table 3: Values for the empirical coefficients and used to determine the fraction $\frac{M}{\beta}$.

Range for dl (m)	Range for X (m)	a2	b2
< 0.0033	all	0	2200
≥ 0.0033	< 7.5	1.93	1640
≥ 0.0033	≥ 7.5	5.14	6600

Detachment by runoff

Detachment by runoff is determined using the approach of Ariathurai & Arulanandan (1978):

equation 98

$$D_q=\begin{cases} k_f(1-C_g-C_r)\Big[\frac{\tau}{\tau_{cr}}-1\Big] &\text{for } \tau>\tau_{cr} \\ 0 &\text{for } \tau<\tau_{cr} \end{cases}$$

Where $D_q$ is the rate of detachment of soil per unit area $(kg\space m^{-2}\space s^{-1})$, $k_f$ the overland flow soil erodibility coefficient $(kg\space m^{-2}\space s^{-1})$, $\tau$ the shear stress due to overland flow $(N \space m^{-2})$, $\tau_{cr}$ the critical shear stress for initiation of sediment motion $(N \space m^{-2})$.

The shear stress due to overland flow $(\tau)$ is given by:

equation 99

$$\tau = \rho ghS$$

with $\rho$ the water density $(kg \space m^{-3})$, g the acceleration due to gravity $(m\space s^{-2})$, h the water depth (m) and S the water surface slope in the direction of the flow $(m\space m^{-1})$.

The water depth (h) is determined with the Manning equation. We assumed a triangular shaped profile on which the Manning equation is applied, where the width-to-depth ratio is a model parameter. First the flow area is determined with an algebraic re-arrangement of the Manning equation:

equation 100

$$A=\Bigg[\frac{Q\space n\Big(2\sqrt{\frac{WD^2+1}{WD}}\Big)^{2/3}}{\sqrt{S}}\Bigg]^{3/4}$$

where Q is the discharge $(m^3\space s^{-1})$, n the Manning's coefficient $(s\space m^{-1/3})$ and WD the width-to-depth ratio (-). The discharge (Q) is obtained from the hydrological model and the Manning's coefficient (n) is defined per land use class.

The water depth (h) is calculated as follows:

equation 101

$$h=\sqrt{\frac{A}{WD}}$$

The critical shear stress $\tau_{cr}$ is calculated as follows:

equation 102

$$\tau_{cr}=(\rho_s- \rho)\space g\space D_{50} \space a_3 \space R_*^{b_3}$$

Where $\rho_s$ is the density of sediment particles $(kg \space m^{-3})$, $D_{50}$ the median sediment particle diameter (m), $R_*$ the particle Reynolds number (-), and $a_3$ and $b_3$ are given in Table 4 and are determined inside the model code.

Table 4. Values for the empirical coefficients $a_3$ and $b_3$ used to determine the particle Reynolds number ($R_*$).

Range for R*	a3	b3
0.03 - 1	0.1	-0.3
1 – 6	0.1	-0.62
6 – 30	0.033	0
30 – 135	0.013	0.28
135 – 400	0.03	0.1
> 400	0.056	0

The particle Reynolds number $R_*$ is calculated as follows:

equation 103

$$R_*=max\Big[0.03,\frac{D_{50}\sqrt{\tau/\rho}}{v}\Big]$$

where $v$ is the water viscosity $(m^2 \space s^{-1})$.

Sediment transported

The total sediment taken into transport (G) is calculated as follows:

equation 104

$$G = \begin{cases} D_r+D_q &\text{for } D_r+D_q<TC \\ TC &\text{for } D_r+D_q>TC \end{cases}$$

Where TC the transport capacity $(kg \space m^{-2} \space s^{-1})$, which is calculated as follows:

equation 105

$$TC=G_{tot}\frac{\rho_s}{A_{cell}}$$

Where $G_{tot}$ is the capacity particulate transport rate for overland flow $(m^3 \space s^{-1})$ and $A_{cell}$ the cells area $(m^2)$.

The capacity particulate transport rate for overland flow $(G_{tot})$ is determined with the formulations from Engelund & Hansen (1967):

equation 106

$$G = \begin{cases} \frac{0.05Q^2S^{3/2}}{\sqrt{gh}(\frac{\rho_s}{\rho}-1)^2D_{50}l}&\text{for } h>0 \\ 0 &\text{for } h<0 \end{cases}$$

Where $l$ is the width of the flow (m), Q the water flow rate $(m^3 \space s^{-1})$. The width of the flow $l$ is determined as:

equation 107

$$l=WD\cdot h$$

DHSVM

The distributed hydrology-soil-vegetation model (DHSVM; Doten et al., 2006) simulates hillslope erosion based on detachment energy of raindrops, leaf drip and surface runoff. The detachment by raindrop impact formulations originate from the SHESED model (Wicks & Bathurst, 1996). These formulations require hourly precipitation intensity as input. While the SPHY hydrological model runs at a daily time step, the model includes a sub-daily infiltration formulation. This formulation determines hourly precipitation intensity, which was subsequently used as input for the DHSVM model. Furthermore, the detachment by raindrop impact formulations require the canopy cover as input, which was obtained from the vegetation module.

Detachment by runoff is determined from a detachment coefficient, the settling velocity, and the transport capacity. The detachment coefficient is a function of the soil cohesion, which is determined from the sum of the soil cohesion and root cohesion. The transport capacity is based on the unit stream power approach from the KINEROS model (Woolhiser et al., 1990), which requires the water depth of the flow as input. We obtained the water depth by applying the Manning equation, assuming a triangular-shaped flow profile, with the width-depth ratio as model parameter.

Detachment by raindrop impact

Detachment by raindrop impact is based on the sum of the momentum squared for rain $(M_R)$ and the momentum squared for leaf drip$(M_D)$(Wicks & Bathurst, 1996):

equation 108

$$D_R=k_rF_w(1-C_G)[(1-C_C)M_R+M_D]$$

where $D_R$ is the soil detached by raindrop impact $(kg \space m^{-2} \space s^{-1})$, $k_r$ the raindrop soil erodibility coefficient $(J^{-1})$, $F_w$ the protective effect of ponding (-), $C_G$ the ground cover (fraction, -), $C_C$ the canopy cover (fraction, -), $M_R$ the momentum squared for rain $(kg^2\space s^{-3})$ and $M_D$ the momentum squared for leaf drip $(kg^2\space s^{-3})$.

The original DHSVM model accounts for the protective effect of ponding on detachment by raindrop impact by model parameter $F_w$(Park, 1982). The hydrological model SPHY does not account for ponding, hence, we assume $F_w=1$.

The momentum squared for the rain $(M_R)$ is determined as follows:

equation 109

$$M_R=\alpha I^{\beta}$$

Where $I$ is the rainfall intensity $(mm\space h^{-1})$ and and are empirical coefficients (Wicks, 1988). The rainfall intensity is determined from the infiltration excess surface runoff formulations of the hydrological model. In these formulations, the hourly rainfall is assumed to decrease linearly over time, where the fraction of the daily rainfall that falls in the first hour is a model parameter. Hence, from these assumptions the hourly rainfall intensity was determined as input for the DHSVM model. Values for $\alpha$ and for $\beta$ each rainfall intensity interval are given in Table 5 and are determined inside the model code.

Table 5. Values for the empirical coefficients $\alpha$ and $\beta$ used to determine the momentum squared for the rain.

Range for (mm h-1)	α	β
0 – 10	2.69 ∙ 10-8	1.6896
10 - 50	3.75 ∙ 10-8	1.5545
50 - 100	6.12 ∙ 10-8	1.4242
≥ 100	11.75 ∙ 10-8	1.2821

Momentum squared for leaf drip $(M_D)$ is calculated as follows:

equation 110

$$M_D=\frac{\big(\frac{V\rho \pi D^3}{6}\big)DRIP\%\space DRAIN}{\big(\frac{\pi D^3}{6}\big)}$$

where V is the leaf drip fall velocity $(m \space s^{-1})$, the density of water $(kg \space m^{-3})$, D the leaf drip diameter (m), $DRIP\%$ the proportion of drainage that falls as leaf drip and$DRAIN$ the canopy drainage rate $(m \space s^{-1})$. The proportion of the drainage that falls as leaf drip ($DRIP\%$) is assumed to be equal to the canopy cover $(C_C)$. The canopy drainage rate ($DRAIN$) is assumed to be equal to the daily precipitation intensity in $m \space s^{-1}$.

The leaf drip fall speed V is calculated as follows (Epema & Riezebos, 1983):

equation 111

$$V=\sqrt{\frac{M}{\beta_V}g\bigg(1-e^{-\frac{2X\beta}{M}}\bigg)}$$

where X is the average leaf drip fall distance (m), M the average mass of leaf drips (kg), $\beta_V$ a friction constant $(kg \space m^{-1})$ and the acceleration due to gravity $m \space s^{-2}$.

The fraction $\frac{M}{\beta_V}$ is a function of the leaf drip diameter D and two coefficients, a and b.

equation 112

$$\frac{M}{\beta_V}=a+b \space D$$

where a and b are a function of the drip diameter and fall distance and are given in Table 6 and are determined inside the model code.

Table 6. Values for the empirical coefficients a and b used to determine the fraction$\frac{M}{\beta_V}$.

Drip diameter (m)	Fall distance (m)	a	b
< 0.0033	all	0	2200
≥ 0.0033	< 7.5	1.93	1640
≥ 0.0033	≥ 7.5	5.14	6600

Detachment by runoff

Detachment by runoff is calculated as follows:

equation 113

$$D_{of}=\beta_{de}\space dy\space v_s \space TC$$

where $D_{of}$ is the soil detachment by overland flow $(kg \space m^{-2} \space s^{-1})$, $\beta_{de}$ the detachment efficiency (-), dy the length of a grid cell (m), $v_s$ the settling velocity ($m \space s^{-1}$) and the transport capacity ($m^3$ sediment $m^{-3}$ water).

The detachment efficiency $\beta_{de}$ is calculated as follows:

equation 114

$$\beta_{de}=0.79e^{-0.6\space COH}$$

Where COH is the soil cohesion (kPa), which is determined from the combination of soil cohesion ($COH_S$) and root cohesion ($COH_r$).

The settling velocity ($v_s$) is calculated following the method as used in the KINEROS model (Woolhiser et al., 1990):

equation 115

$$v_s=\sqrt{\frac{(\frac{4}{3}g\frac{\rho_s}{\rho}-1)d_{50}}{C_d}}$$

where $\rho_s$ the sediment density $(kg \space m^{-3})$, $d_{50}$ the median grain size (m) and $C_d$ the drag coefficient, which is a function of the particle Reynolds number:

equation 116

$$C_d=\frac{24}{R_n}+\frac{3}{\sqrt{R_n}}+0.34$$

where $R_n$ is the particle Reynolds number, defined as:

equation 117

$$R_n=\frac{v_{s_0}d_{50}}{\nu}$$

where $v_{s_0}$ is an initial estimate of the settling velocity $(m \space s^{-1})$ and $\nu$ is the kinematic viscosity of water $(m^2 \space s^{-1})$, assumed to be equal to $1 \cdot 10^{-6} \space m^2 \space s^{-1}$. The initial estimate of the settling velocity $v_{s_0}$ is calculated as follows:

equation 118

$$v_{s_0}=\sqrt{\bigg(\frac{4}{3}g\frac{\rho_s}{\rho}-1\bigg)d_{50}}$$

The transport capacity (TC) is determined according to the unit stream power method from the KINEROS model (Woolhiser et al., 1990):

equation 119

$$TC=\frac{0.05}{d_{50}\big(\frac{\rho_s}{\rho}-1\big)^2}\sqrt{\frac{S \space h}{g}}(SP-SP_{cr})$$

where S is the slope $(m \space m^{-1})$, h the water depth (m), SP the stream power $(kg \space m \space s^{-3})$ and $SP_{cr}$ the critical stream power $(kg \space m \space s^{-3})$.

The water depth (h) is determined with the Manning equation. We assumed a triangular shaped profile on which the Manning equation is applied, where the width-to-depth ratio is a model parameter. First the flow area is determined with an algebraic re-arrangement of the Manning equation:

equation 120

$$A=\Bigg[\frac{Q\space n\bigg(2\sqrt{\frac{WD^2+1}{WD}}\bigg)^{2/3}}{\sqrt{S}}\Bigg]^{3/4}$$

where Q is the discharge $(m^3 \space s^{-1})$, n the Manning's coefficient $(s\space m^{-1/3})$ and WD the width-to-depth ratio (-). The discharge (Q) is obtained from the hydrological model and the Manning's coefficient (n) is defined per land use class.

The water depth (h) is calculated as follows:

equation 121

$$h=\sqrt{\frac{A}{WD}}$$

The stream power (SP) is calculated as follows:

equation 122

$$SP=\rho \space g \space Q \space S$$

Sediment transported

The sediment taken into transport is simply the sum of detachment by raindrop impact and runoff:

equation 123

$$sed=D_R+D_{of}$$

Where sed is the sediment taken into transport $(kg \space m^2 \space s^{-1})$.

HSPF

The Hydrological Simulation Program-Fortran (HSPF) model (Bicknell et al., 1993) simulates detachment by raindrop impact with daily precipitation intensity as input. The soil erodibility is based on the USLE K-factor, here estimated using the method proposed by Wischmeier et al. (1971). Detached sediment by raindrop impact is stored in the sediment storage, which decreases as a result of soil crusting, simulated by a reduction parameter. The amount of detached sediment by raindrop impact taken into transport is a function of the sediment storage and the transport capacity. Detachment by runoff is a function of surface runoff and a coefficient for scour of the soil matrix.

Detachment by raindrop impact

The original HSPF model accounts for the surface water storage (SURS), for instance as a result of ponding. Since the hydrological model SPHY does not account for ponding, we assume the surface water storage to be equal to 0. The detachment by raindrop impact, which is called washoff of detached sediment by raindrop impact in Bicknell et al. (1993), is calculated as follows:

equation 124

$$WSSD=\begin{cases} DETS&\text{for } STCAP>DETS \\ STCAP &\text{for } STCAP<DETS \end{cases}$$

Where WSSD is the detachment by raindrop impact $(ton \space acre^{-1})$, DETS is the sediment storage $(ton \space acre^{-1})$ and STCAP is the capacity for removing detached sediment $(ton \space acre^{-1})$.

The sediment storage is calculated as follows:

equation 125

$$DETS=DETS(1-AFFIX)+DET$$

Where DET is the sediment detached from the soil matrix by rainfall $(ton \space acre^{-1})$ and AFFIX is the fraction by which DETS decreases each day as a result of soil compaction (-).

The sediment detached from the soil matrix by rainfall DET is calculated as follows:

equation 126

$$DET=DELT60\space (1-CR)\space SMPF\space KRER \bigg(\frac{RAIN}{DELT60}\bigg)^{JRER}$$

Where DELT60 is the number of hours per interval (-), CR the fraction of the land covered by vegetation (-), SMPF the supporting management practice factor (-), KRER the detachment coefficient dependent on soil properties (-), RAIN the rainfall $(inch \space interval^{-1})$ and JRER the detachment exponent dependent on soil properties (-).

The supporting management practice factor SMPF is assumed to be 1 for all land use classes. The detachment coefficient dependent on soil properties KRER is estimated with the USLE K-factor developed by Wischmeier et al. (1971):

equation 127

$$KRER=\frac{0.00021M^{1.14}(12-OM)+3.25(c_{soilstr}-2)+2.5(c_{perm}-3)}{100}$$

Where KRER is the particle-size parameter (-), OM is the organic matter content (%), $c_{soilstr}$ is the soil structure class (-) and $c_{perm}$ is the profile permeability class (-).

The particle-size parameter is calculated as follows:

equation 128

$$M=(m_{silt}+m_{vfs})(100-m_c)$$

Where $m_{silt}$ is the silt content (%), $m_{vfs}$ is the very fine sand content (%) and $m_{c}$ is the clay content (%).

The capacity for removing detached sediment STCAP is calculated as follows:

equation 129

$$STCAP=DELT60\space KSER \Big(\frac{SURO}{DELT60}\Big)^{JSER}$$

Where KSER the coefficient for transport of detached sediment (-), SURO the surface outflow of water $(inch \space interval^{-1})$ and JSER the exponent for transport of detached sediment (-). The surface water storage SURO is estimated by the (routed) runoff from the hydrological model.

Detachment by runoff

Detachment by runoff, which is called scour of matrix soil in Bicknell et al. (1993), is calculated as follows:

equation 130

$$SCRSD=DELT60 \space KGER \Big( \frac{SURO}{DELT60} \Big)^{JGER}$$

Where SCRSD is the scour of matrix soil $(ton \space acre^{-1})$, KGER is the coefficient for scour of the matrix soil (-) and JGER the exponent for scour of the matrix soil (-).

Sediment transported

The sediment taken into transport is simply the sum of detachment by raindrop impact and runoff:

equation 131

$$SOSED=WSSD+SCRSD$$

Where SOSED is the total removal of soil and sediment from the surface by water $(ton \space acre^{-1})$.