Assessing the large-scale impacts of environmental change using a coupled hydrology and soil erosion model

Assessing the impacts of environmental change on soil erosion and sediment yield at the large catchment scale remains one of the main challenges in soil erosion modelling studies. Here, we present a process-based soil erosion model, based on the integration of the Morgan–Morgan–Finney erosion model in a daily based hydrological model. The model overcomes many of the limitations of previous large-scale soil erosion models, as it includes a more complete representation of crucial processes like surface runoff generation, dynamic vegetation development, and sediment deposition, and runs at the catchment scale with a daily time step. This makes the model especially suited for the evaluation of the impacts of environmental change on soil erosion and sediment yield at regional scales and over decadal periods. The model was successfully applied in a large catchment in southeastern Spain. We demonstrate the model's capacity to perform impact assessments of environmental change scenarios, specifically simulating the scenario impacts of intra- and inter-annual variations in climate, land management, and vegetation development on soil erosion and sediment yield.

1 Introduction

Climate change will likely affect soil erosion and sediment yield at a wide variety of scales (Burt et al., 2016; Li et al., 2006; Nearing et al., 2004). However, assessing the impacts of environmental change on soil erosion and sediment yield at the large catchment scale remains one of the main challenges in soil erosion modelling (Poesen, 2018; de Vente et al., 2013). Most soil erosion and sediment yield models adopt simplified model formulations, are applied at low temporal resolutions, and in many cases only partly represent the impacts of changes in land use or climate conditions. This often leads to unreliable results that do not sufficiently increase process understanding or support decision making (de Vente et al., 2013). To overcome part of these limitations, we present a process-based, large-scale, soil erosion model, coupled to a hydrological model, which accounts for the most relevant factors determining soil erosion by water, including saturated and infiltration excess surface runoff, dynamic vegetation development, and sediment deposition.

Essentially, soil erosion by water occurs due to the impact of raindrops, overland flow, and river flow. Therefore, it is crucial to quantify raindrop impact, overland flow, and possible interactions with vegetation cover. Soil erosion due to the impact of raindrops is a function of the amount and the size of the raindrops that reach the soil surface (Morgan, 2005). Vegetation cover can reduce the impact of raindrops by interception, separating the precipitation into direct throughfall, which has a high impact, and leaf drainage, which has a lower impact (Morgan and Nearing, 2011). Therefore, assessment of soil erosion needs to account for spatial and temporal changes in vegetation cover. Nevertheless, most large-scale soil erosion assessments do not consider dynamic seasonal and inter-annual vegetation development. Previous studies have also shown that cropping patterns, inter-annual precipitation variability, and land use changes can have a significant impact on soil erosion rates (Maetens et al., 2012).

Soil erosion is also a function of runoff (Morgan, 2005). Runoff generation depends on surface and subsurface processes and is a function of precipitation volume and intensity, soil moisture, and soil hydraulic properties (Kirkby, 1988). Process-based models often incorporate a separate hydrological model to simulate surface runoff generation, which then directly forces soil erosion by runoff (Morgan, 2005). Surface runoff may be generated by several distinct processes, of which saturation excess and infiltration excess are the most common (Beven, 2012). Saturation excess surface runoff occurs when the soil water content reaches saturation, while infiltration excess surface runoff occurs when the precipitation intensity exceeds the soils infiltration capacity. However, many large-scale soil erosion models only consider saturation excess surface runoff, disregarding the infiltration excess surface runoff mechanism. Infiltration excess surface runoff is a sub-daily process that is often only implemented in event-scale models (Morgan et al., 1998; Nearing et al., 1989). Infiltration excess surface runoff is an especially important process in areas where the major part of soil erosion takes place during extreme rainfall events with high precipitation intensities (Farnsworth and Milliman, 2003; González-Hidalgo et al., 2007, 2012; López-Bermúdez et al., 2002; Mulligan, 1998).

Assessing catchment sediment yield further requires the evaluation of the sediment transport capacity and sediment deposition. Sediment deposition occurs when the sediment transport capacity of the runoff is exceeded and basically depends on the interaction between the texture of the detached soil material, flow velocity, and the roughness of the surface (Rose et al., 1983). Large particles are more likely to be deposited close to the source, while small particles are more easily brought and maintained into transport. The roughness of the surface is a combination of the soil roughness and vegetation roughness (Chow, 1959). Many studies have highlighted the importance of accounting for sediment transport and sediment deposition, claiming that a large part of the eroded sediment is often deposited close to its source (Walling, 1983; de Vente et al., 2007). However, many large-scale soil erosion models insufficiently consider this process or only simulate soil detachment processes.

From the available soil erosion models, process-based models aim to incorporate the most relevant processes driving soil detachment, sediment transport, and deposition, as described in the previous paragraphs, see Morgan (2005) for an overview of process-based models. Process-based models are often run at small spatial (hillslope to small catchment) and temporal scales (sub-hourly to daily time steps). Most detailed assessments are obtained from event-scale model applications in process-based models, such as WEPP (Nearing et al., 1989) and EUROSEM (Morgan et al., 1998), which require detailed input data, such as (sub-)hourly precipitation and topographic data, and incorporate a large number of calibration parameters (Govers, 2011). While these models present a strong potential to provide increased process understanding, it is often infeasible to obtain all required input data for large catchments. Furthermore, the high uncertainty regarding how input data and model parameters will change under scenarios of environmental change severely limits their application in large-scale assessments.

At larger scales, soil erosion is often assessed using empirical erosion models, see de Vente et al. (2013) for an overview. These models are derived from field studies where soil erosion has been observed under different land use, management, soil, climate, and topographical conditions. The most well known and applied empirical model is the Universal Soil Loss Equation (Wischmeier and Smith, 1978) and its derivatives RULSE (Renard et al., 1997) and MUSLE (Williams, 1995). While the empirical formulations of the USLE were obtained at plot-scale, the model is often applied at much larger scales, sometimes in combination with a sediment transport capacity equation or a sediment delivery ratio to assess sediment yield (Van Rompaey et al., 2001; de Vente et al., 2008). Due to their simplicity, the USLE and its derivatives can be applied with a relatively limited amount of input data. However, their main restriction is the limited number of processes accounted for (e.g. the USLE and RUSLE based models only consider sheet and rill erosion) and the limited potential to evaluate the impacts of changes in climate and land management (Govers, 2011; de Vente et al., 2013). Furthermore, these models are typically applied at annual time steps, largely neglecting intra-annual variation of climate and vegetation conditions.

Figure 1 Overview of the model: (a) representation of a single cell; (b) the hydrological processes, where ET is evapotranspiration; and (c) the soil erosion processes, where P is precipitation, DT is direct throughfall (Eq. 14), LD is leaf drainage (Eq. 13), F is detachment by raindrop impact (Eq. 18), H is detachment by runoff (Eq. 19), D is immediate deposition (Eq. 20), and G is sediment routed to downstream cells (Eq. 27).

Most current soil erosion models have a limited potential for application at larger temporal and spatial scales (i.e. process-based models) or lack sufficient representation of the underlying soil detachment and sediment transport processes and sensitivity to changes in land use or climate (i.e. empirical models), making them of limited use for scenario studies and process understanding. Here, we present a process-based soil erosion model based on the integration of the Morgan–Morgan–Finney erosion model (Morgan and Duzant, 2008) and the spatially distributed hydrological model “Spatial Processes in HYdrology” (Terink et al., 2015). This integrated model overcomes many of the limitations of previous large-scale soil erosion models, as it includes a more complete representation of crucial processes, such as surface runoff generation, dynamic vegetation development, and sediment deposition, and runs at the catchment scale with a daily time step. This makes the model especially suitable for the evaluation of the inter- and intra-annual impacts of environmental change on soil erosion and sediment yield at large spatial and temporal scales. In the next paragraphs we first present the different model components and enhancements as compared to previous models. We then illustrate its functionality and potential for scenario studies by applying the model to the upper Segura catchment in southeastern Spain under present and projected future climate conditions.

2 Model description

2.1 Model overview

The SPHY–MMF model presented here is an integration of the (Modified) Morgan–Morgan–Finney soil erosion model into the SPHY hydrological model (version 2.1). Figure 1 shows the main hydrological and soil erosion processes considered by the model. SPHY is a spatially distributed “leaky-bucket” type model that simulates hydrological processes on a cell-by-cell basis at a daily time step (Terink et al., 2015). The model is written in the Python programming language using the PCRaster dynamic modelling framework (Karssenberg et al., 2010). MMF is a conceptual soil erosion model that is originally applied with an annual time step. Here we present a modification of the model at a daily time step, fully integrated with the SPHY model. MMF receives input from the SPHY model, such as effective precipitation (throughfall), runoff, and canopy cover for calculation of erosion and deposition processes. The model parameters are presented in Table S1 in the Supplement. A video abstract (Eekhout and de Vente, 2018) is available that explains the most important aspects of the model.

2.2 Hydrological model

SPHY simulates most relevant hydrological processes (Fig. 1b), such as interception, evapotranspiration, dynamic evolution of vegetation cover, surface runoff, and lateral and vertical soil moisture flow at a daily time step. The model is described in full detail by Terink et al. (2015); therefore, here we only provide a summary of the processes that are simulated by the model, some hydrological processes that have been changed with respect to the original SPHY model, and a detailed description of the processes that are related to the integration of MMF.

SPHY requires daily precipitation and temperature maps as input. Effective precipitation is determined by subtracting canopy storage and interception from precipitation. Canopy storage is determined from the leaf area index (LAI), which is derived from normalized differenced vegetation index (NDVI) images. Reference evapotranspiration is determined using the Hargreaves equation (Hargreaves and Samani, 1985), which is subsequently multiplied by the crop coefficient to obtain the potential evapotranspiration. The crop coefficient is determined from the NDVI images using a linear relationship. Actual evapotranspiration is obtained by multiplying the potential evapotranspiration by a reduction factor for water deficit or water surplus, which are functions of current soil water content, soil hydraulic properties, and plant-specific water need. Surface runoff is determined by a daily implementation of the Green–Ampt formula (Heber Green and Ampt, 1911) and is a function of infiltration, effective precipitation, and soil hydraulic properties. The soil profile consists of three layers, i.e. root-zone, subzone, and groundwater layer. Water can percolate from the root-zone to the subzone and from the subzone to the groundwater layer. Water travels from the subzone to the root zone through capillary rise. Water drains from the root zone as lateral flow and from the groundwater layer as baseflow. The total runoff is the sum of surface runoff, lateral flow, and baseflow. All soil processes are functions of current water content (in the respective layers) and soil hydraulic properties, i.e. saturated hydraulic conductivity, saturated water content, field capacity, and wilting point. Water is routed using a single flow algorithm (Karssenberg et al., 2010). A flow recession coefficient accounts for flow delay from channel friction. When reservoirs are present, the user can opt to include an advanced routing scheme accounting for reservoir storage and outflow.

2.2.1 Evapotranspiration

The actual evapotranspiration (ET_a) is determined by multiplying the potential evapotranspiration with the reduction parameters for water surplus and water deficit conditions. In the current version of the hydrological model we have changed the reduction parameter for water shortage conditions using the method proposed by Allen et al. (1998):

$$\begin{array}{}\text{(1)}& K_{\mathrm{s}}={\displaystyle \frac{\mathrm{TAW}-D_{\mathrm{r}}}{(\mathrm{1}-p)\mathrm{TAW}}},\end{array}$$

where K_s is the reduction parameter for water shortage (–), TAW is the total available water in the root zone (mm), D_r is the root-zone depletion (mm), and p is the depletion fraction (–). The total available water TAW is defined as

$$\begin{array}{}\text{(2)}& \mathrm{TAW}={\mathit{\theta }}_{\mathrm{FC}}-{\mathit{\theta }}_{\mathrm{WP}},\end{array}$$

where θ_FC is the soil water content at field capacity (mm) and θ_WP is the soil water content at wilting point (mm). The root-zone depletion D_r is defined as

$$\begin{array}{}\text{(3)}& D_{\mathrm{r}}={\mathit{\theta }}_{\mathrm{FC}}-\mathit{\theta },\end{array}$$

where θ is the current soil water content (mm). The depletion fraction p is defined as the fraction of TAW that a crop can extract from the root zone without suffering water stress, which is determined by the following equation:

$$\begin{array}{}\text{(4)}& p=p_{\mathrm{tabular}}+\mathrm{0.04}(\mathrm{5}-{\mathrm{ET}}_{\mathrm{pot}}),\end{array}$$

where p_tabular is a land use specific tabular value of the depletion fraction (–) and ET_pot is the potential evapotranspiration (mm). Values for the land use specific tabular value of the depletion fraction can be obtained from Table 22 in Allen et al. (1998).

Open water evaporation is determined in the reservoir cells. In these cells all soil hydraulic processes are turned off and runoff equals precipitation minus open water evaporation. Reservoir cells cannot dry up, i.e. we assume that there is always water present in the reservoir cells. Open water evaporation is determined as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{ET}}_{\mathrm{open}-\mathrm{water}}={\mathit{kc}}_{\mathrm{open}-\mathrm{water}}{\mathrm{ET}}_{\mathrm{ref}},\end{array}$$

where kc_open−water is the crop coefficient value for open water evaporation (–) and ET_ref is the reference evapotranspiration (mm). We set kc_open−water to a value of 1.2, following Allen et al. (1998). In each time step the open water evaporation is subtracted from the reservoir storage.

2.2.2 Infiltration excess surface runoff

The original SPHY model simulates saturated surface runoff but not infiltration excess surface runoff. Therefore, we have incorporated an infiltration excess equation at a daily time step, based on the Green–Ampt formula (Heber Green and Ampt, 1911). We assumed a constant infiltration rate (f) (mm), which is determined for each cell and each day by

$$\begin{array}{}\text{(6)}& f={\displaystyle \frac{K_{\mathrm{eff}}}{\mathrm{24}}}{\left[\mathrm{1}+{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{sat}}-\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{sat}}}}\right]}^{\mathit{\lambda }},\end{array}$$

where K_eff is the effective hydraulic conductivity (mm), θ_sat is the saturated water content (mm), θ is the actual water content (mm), and λ is a calibration parameter (–). Bouwer (1969) suggested an approximation of K_eff≈0.5 K_sat. We included a calibration parameter (k) in order to be able to change the value of K_eff as a fraction of K_sat (K_eff=kK_sat).

Infiltration excess surface runoff occurs when the precipitation intensity exceeds the infiltration rate f (Beven, 2012). Analysis of hourly precipitation time series for 25 years (1991–2015) from 5 precipitation stations in a large catchment in southeastern Spain showed that, on average, the highest precipitation intensity was recorded in the first hour of a rain storm and decreased linearly until the end of the storm. We assumed a triangular-shaped precipitation intensity p(t) (mm hr⁻¹) according to the following:

$$\begin{array}{}\text{(7)}& p\left(t\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\alpha }}^{\mathrm{2}}Pt+\mathit{\alpha }P,\end{array}$$

where α is the fraction of daily rainfall that occurs in the hour with the highest intensity (–), P is the daily rainfall (mm), and t is an hourly time step. Daily infiltration excess surface runoff Q_surf is determined as

$$\begin{array}{}\text{(8)}& Q_{\mathrm{surf}}=\left\{\begin{array}{rl}{\displaystyle \frac{{\left(\mathit{\alpha }P-f\right)}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}P}}& {\displaystyle }\mathrm{if}\mathit{\alpha }P>f\\ {\displaystyle }\mathrm{0}& {\displaystyle }\mathrm{if}\mathit{\alpha }P\le f\end{array}\right.\end{array}$$

2.2.3 Dynamic vegetation processes

SPHY–MMF contains a dynamic vegetation module that allows for the characterisation of the seasonal and inter-annual differences in vegetation cover and resulting canopy storage, interception, and precipitation throughfall. The latter is subsequently used in both the hydrological and soil erosion model. A time series of NDVI images is used as input for the dynamic vegetation module. NDVI images may only be available for a limited period, e.g. from 2000 to the present for MODIS NDVI images (Didan, 2015). Therefore, in Sect. 3 we present a method to obtain NDVI images for historical and future model assessments. The LAI is determined from the individual NDVI images using the following logarithmic relation, which is valid for vegetation that is evenly distributed over a surface (Sellers et al., 1996):

$$\begin{array}{}\text{(9)}& \mathrm{LAI}={\mathrm{LAI}}_{\max }{\displaystyle \frac{\log (\mathrm{1}-\mathrm{fPAR})}{\log (\mathrm{1}-{\mathrm{fPAR}}_{\max })}},\end{array}$$

where LAI_max is the maximum LAI (–), fPAR is the fraction of the photosynthetically active radiation (–), and fPAR_max is the maximum fPAR (–), which is set to 0.95 (Sellers et al., 1996). The maximum LAI (LAI_max) is vegetation dependent; values for several vegetation types can be found in Sellers et al. (1996). The fraction of the photosynthetically active radiation fPAR is determined as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathrm{fPAR}={\displaystyle \frac{(\mathit{SR}-{\mathit{SR}}_{\min })({\mathrm{fPAR}}_{\max }-{\mathrm{fPAR}}_{\min })}{{\mathit{SR}}_{\max }-{\mathit{SR}}_{\min }}}+{\mathrm{fPAR}}_{\min },\end{array}$$

where SR is a transformation of NDVI (–), SR_min and SR_max are the minimum and maximum SR values (–), respectively, and fPAR_min is the minimum fPAR (–), which is set to 0.001 (Sellers et al., 1996). fPAR is bounded by fPAR_min and fPAR_max. SR is determined as follows:

$$\begin{array}{}\text{(11)}& \mathit{SR}={\displaystyle \frac{\mathrm{1}+\mathrm{NDVI}}{\mathrm{1}-\mathrm{NDVI}}}.\end{array}$$

SR_min and SR_max are determined with Eq. (11), applying an NDVI value corresponding to the 5 % and 98 % quantiles, respectively.

2.3 Soil erosion simulation with a daily based Morgan–Morgan–Finney model

The soil erosion model is based on the Modified MMF model (Morgan and Duzant, 2008). In the current model (Fig. 1c), total soil erosion is calculated from detachment by raindrop impact and detachment by runoff, while sediment yield is calculated by routing detached sediment, considering within cell deposition and sediment transport capacity. Detachment of soil particles from raindrop impact is determined from the total rainfall energy, which is determined for direct throughfall and leaf drainage, respectively. Detachment of soil particles by runoff is determined from the accumulated runoff from the hydrological model. Both soil erosion equations account for the fraction of the soil covered by stones and vegetation or snow and are determined separately for three texture classes (sand, silt, and clay). Within cell deposition is calculated as a function of vegetation and surface roughness. The remainder of the detached sediment is taken into transport and routed through the catchment, taking the transport capacity of the flow and the trapping efficiency of the reservoirs into account. The next paragraphs provide a detailed description of all of these processes.

2.3.1 Estimation of rainfall energy

The total kinetic energy of the effective precipitation (KE, J m⁻²) is used to determine the detachment of soil particles by raindrop impact and is defined as

$$\begin{array}{}\text{(12)}& \mathrm{KE}={\mathrm{KE}}_{\mathrm{DT}}+{\mathrm{KE}}_{\mathrm{LD}},\end{array}$$

where KE_LD is the kinetic energy of the leaf drainage (J m⁻²) and KE_DT is the kinetic energy of the direct throughfall (J m⁻²).

The kinetic energy of the leaf drainage is based on Brandt (1990):

$$\begin{array}{}\text{(13)}& {\mathrm{KE}}_{\mathrm{LD}}=\left\{\begin{array}{rl}{\displaystyle }\mathrm{0}& {\displaystyle }\mathrm{for}\mathrm{PH}<\mathrm{0.15}\\ {\displaystyle }\mathrm{LD}(\mathrm{15.8}{\mathrm{PH}}^{\mathrm{0.5}}-\mathrm{5.87})& {\displaystyle }\mathrm{for}\mathrm{PH}\ge \mathrm{0.15},\end{array}\right.\end{array}$$

where LD is the leaf drainage (mm) and PH is the plant height (m), specified for each land use class.

The kinetic energy of the direct throughfall is based on a relationship described by Marshall and Palmer (1948), which is representative of a wide range of environments (Morgan, 2005):

$$\begin{array}{}\text{(14)}& {\mathrm{KE}}_{\mathrm{DT}}=\mathrm{DT}(\mathrm{8.95}+\mathrm{8.44}{\log }_{\mathrm{10}}I),\end{array}$$

where DT is the direct throughfall (mm) and I is the intensity of the erosive precipitation (mm h⁻¹). The intensity of the erosive precipitation is a model parameter and varies according to geographical location. Morgan and Duzant (2008) proposed 10 mm h⁻¹ for temperate climates, 25 mm h⁻¹ for tropical climates, and 30 mm h⁻¹ for strongly seasonal climates (e.g. Mediterranean, tropical monsoon).

The leaf drainage (LD) (the precipitation that reaches the soil surface as flow or drips from the leaves and stems of the vegetation) and direct throughfall (DT) (the precipitation that reaches the soil surface directly through gaps in the vegetation cover), from Eqs. (13) and (14), are obtained from the effective precipitation, also known as throughfall (P_eff, mm). The effective precipitation from the hydrological model is first corrected for the slope angle, following Choi et al. (2017):

$$\begin{array}{}\text{(15)}& P_{\mathrm{eff}}=P_{\mathrm{eff}}\cos S,\end{array}$$

where P_eff is the effective precipitation (mm) and S is the slope (^∘).

Leaf drainage is determined as

$$\begin{array}{}\text{(16)}& \mathrm{LD}=P_{\mathrm{eff}}\mathrm{CC},\end{array}$$

where CC is the canopy cover (proportion between zero and unity). The canopy cover is either introduced by a land use class specific tabular value or determined by the vegetation module. When the vegetation module is used, the canopy cover is obtained from the LAI (Eq. 9), maximised by 1.

Direct throughfall becomes

$$\begin{array}{}\text{(17)}& \mathrm{DT}=P_{\mathrm{eff}}-\mathrm{LD}.\end{array}$$

2.3.2 Detachment of soil particles

Detachment of soil particles is determined separately for raindrop impact and accumulated runoff and is subsequently summed.

The detachment of soil particles by raindrop impact (F, kg m⁻²) and the detachment of soil particles by runoff (H, kg m⁻²) are determined for each of the soil texture classes separately and subsequently summed. The detachment of soil particles by raindrop impact is calculated as

$$\begin{array}{}\text{(18)}& F_i=K_i{\displaystyle \frac{\mathit{\%}i}{\mathrm{100}}}(\mathrm{1}-\mathrm{GC})\mathrm{KE}\times {\mathrm{10}}^{-\mathrm{3}},\end{array}$$

where K is the detachability of the soil by raindrop impact (g J⁻¹), i is the textural class, with c for clay, z for silt, and s for sand, and GC is the ground cover (–). The detachability of the soil for each texture class is included as a model parameter, for which Quansah (1982) proposed K_c=0.1, K_z=0.5, and K_s=0.3 g J⁻¹. The ground cover, 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 one when the surface is covered with snow, which is determined in the hydrological model.

The detachment of soil particles by runoff (H, kg m⁻²) is calculated as

$$\begin{array}{}\text{(19)}& H_i={\mathrm{DR}}_i{\displaystyle \frac{\mathit{\%}i}{\mathrm{100}}}{Q}^{\mathrm{1.5}}(\mathrm{1}-\mathrm{GC}){\sin }^{\mathrm{0.3}}S\times {\mathrm{10}}^{-\mathrm{3}},\end{array}$$

where DR is the detachability of the soil by runoff (g mm⁻¹) and Q is the volume of accumulated runoff (mm). The detachability of the soil for each texture class is included as a model parameter for which Quansah (1982) proposed DR_c=1.0, DR_z=1.6, and DR_s=1.5 g mm⁻¹.

2.3.3 Immediate deposition of detached particles

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:

$$\begin{array}{}\text{(20)}& {\mathrm{DEP}}_i=\mathrm{44.1}{N}_{{\mathrm{f}}_i}^{\mathrm{0.29}}.\end{array}$$

In the abovementioned equation, N_f is the particle fall number (–) defined as

$$\begin{array}{}\text{(21)}& N_{\mathrm{f}{}_i}={\displaystyle \frac{l{v}_{s_i}}{vd}},\end{array}$$

where l is the length of a grid cell (m), v_s is the particle fall velocity (m s⁻¹), v is the flow velocity (m s⁻¹), and d is the depth of flow (m). Particle fall velocities are estimated from

$$\begin{array}{}\text{(22)}& v_s={\displaystyle \frac{\mathrm{1}/\mathrm{18}{\mathit{\delta }}^{\mathrm{2}}({\mathit{\rho }}_s-\mathit{\rho })g}{\mathit{\eta }}},\end{array}$$

where δ is the diameter of the particle (m), ρ_s is the sediment density (=2650 kg m⁻³), ρ is the flow density (Abrahams et al., 2001), g is gravitational acceleration (taken as 9.81 m s⁻²), and η is the fluid viscosity (Morgan and Duzant, 2008). When Eq. 22 is applied to the three texture sizes of 2 µm for clay, 60 µm for silt, and 200 µm for sand, this gives respective values of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ for clay, $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ for silt, and 0.02 m s⁻¹ for sand.

The flow velocity v from Eq. (21) is obtained by the Manning formula:

$$\begin{array}{}\text{(23)}& v={\displaystyle \frac{\mathrm{1}}{n^{\prime }}}{d}^{\mathrm{2}/\mathrm{3}}{S}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where n^′ is the modified Manning roughness coefficient (s m${}^{-\mathrm{1}/\mathrm{3}}$), which is a combination of the Manning roughness coefficient for the soil surface and vegetation. The Manning roughness coefficient is defined as per Petryk and Bosmajian (1975):

$$\begin{array}{}\text{(24)}& n^{\prime }=\sqrt{n_{\mathrm{soil}}^{\mathrm{2}}+n_{\mathrm{vegetation}}^{\mathrm{2}}}.\end{array}$$

The Manning roughness coefficient for bare soil (n_soil) is set to 0.015 s m${}^{-\mathrm{1}/\mathrm{3}}$, as suggested by Morgan and Duzant (2008) (Fig. 2a). For tilled conditions (Fig. 2b) the following equation is applied:

$$\begin{array}{}\text{(25)}& n_{\mathrm{soil}}=\exp (-\mathrm{2.1132}+\mathrm{0.0349}\mathrm{RFR}),\end{array}$$

where RFR is the surface roughness parameter (cm m⁻¹). Values for RFR are tillage implementation specific and can be obtained from Table 4 in Morgan and Duzant (2008).

Figure 2 Surface and vegetation roughness options: (a) bare soil, (b) tilled soil (Eq. 25), (c) regular vegetation (Eq. 26), and (d) irregular vegetation.

The Manning roughness coefficient for regularly spaced vegetation (n_vegetation) (Fig. 2c) is obtained from the following equation (Jin et al., 2000):

$$\begin{array}{}\text{(26)}& n_{\mathrm{vegetation}}={\displaystyle \frac{d^{\mathrm{2}/\mathrm{3}}}{\sqrt{\frac{\mathrm{2}g}{D\mathit{NV}}}}},\end{array}$$

where D is the stem diameter (m) and NV is the stem density (stems m⁻²). Stem diameter and stem density may be difficult to obtain for certain land use classes with irregularly spaced vegetation (e.g. forest, shrubland); therefore, users may opt to use tabular values for n_vegetation, e.g. from Chow (1959) (Fig. 2d). Preliminary model runs showed that Eq. (26) results in unrealistically high flow velocity values for land use classes where the stem density is very low, such as in orchards. Therefore, in these conditions where the influence of vegetation on flow velocity is negligible, n_vegetation can be set to zero.

Figure 3 Sediment routing: (a) transport of sediment through the catchment (Eq. 28), and (b) trapping efficiency at the reservoirs (Eq. 30).

Equations (21), (23), and (26) require a flow depth (d), which is 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.

2.3.4 Sediment deposition and transport

The amount of sediment that is routed to downstream cells is determined from the sum of the detached sediment from raindrop impact (Eq. 18) and accumulated runoff (Eq. 19), subtracting the proportion of the sediment that is deposited within the cell of its origin (Eq. 20):

$$\begin{array}{}\text{(27)}& G=(F_i+H_i)(\mathrm{1}-({\mathrm{DEP}}_i/\mathrm{100}\left)\right).\end{array}$$

The amount of sediment that is routed to downstream cells is the summation of the individual amounts for clay, silt, and sand.

Sediment is routed using a routing scheme that takes both the transport capacity (TC; ton ha⁻¹) of the accumulated runoff and the trapping efficiency of the reservoirs (TE; –) into account. The transport capacity TC (Fig. 3a) of the accumulated runoff is based on Prosser and Rustomji (2000):

$$\begin{array}{}\text{(28)}& \mathrm{TC}={\mathrm{flow}}_{\mathrm{factor}}{q}_{\mathrm{surf}}^{\mathit{\beta }}{S}^{\mathit{\gamma }},\end{array}$$

where flow_factor is a spatially distributed roughness factor (–), q_surf is accumulated runoff per unit width (m² day⁻¹), S is the local energy gradient (^∘), approximated by the slope, and β and γ are model parameters (–). As suggested by Prosser and Rustomji (2000), γ is set to 1.4 and β is used for model calibration.

The roughness factor (flow_factor) is determined as follows:

$$\begin{array}{}\text{(29)}& {\mathrm{flow}}_{\mathrm{factor}}={\displaystyle \frac{v_{\mathrm{actual}}}{v_{\mathrm{bare}}}},\end{array}$$

where v_actual is the actual flow velocity (m s⁻¹) and v_bare is the flow velocity for bare soil conditions (m s⁻¹). The actual flow velocity (v_actual) is obtained from Eqs. (23)–(26), applying a water depth (d_actual) of 0.25 m, which coincides with deeper rills from Morgan and Duzant (2008). The flow velocity for bare soil conditions (v_bare) is obtained from Eq. 23, applying values for $n^{\prime }=\mathrm{0.015}$ s m${}^{-\mathrm{1}/\mathrm{3}}$ and d_bare=0.005 m (Morgan and Duzant, 2008).

Reservoir sediment trapping efficiency (TE) (Fig. 3b), the percentage of sediment trapped by the reservoir, is calculated according to Brown (1943):

$$\begin{array}{}\text{(30)}& \mathrm{TE}=\mathrm{100}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{0.0021}D\frac{C}{A_{\mathrm{basin}}}}}\right],\end{array}$$

where D is a constant (–) within the range from 0.046 to 1, depending on the reservoir operation characteristics that we set at 0.1, C is the reservoir capacity (m³), and A_basin is the drainage area of the subcatchment (km²).

3 Model application

To illustrate the model performance, its functionality, and its capacity for use in scenario studies, we applied the model to the upper Segura River catchment (2589 km²) under present and future projected climate conditions. The upper Segura catchment is located in the headwaters of the Segura River in southeastern Spain (Fig. 4). The elevation ranges between 411 and 2055 m a.s.l. (Fig. 4). The climate in the catchment is classified as temperate (Cfa and Cfb according to the Köppen–Geiger climate classification, 80 %) and semi-arid (BSk, 20 %). The catchment-average annual precipitation is 570 mm (for the period from 1981 to 2000) and the mean annual temperature is 13.2 ^∘C (for the period from 1981 to 2000).

The main land use types are forest (45 %), shrubland (40 %), cereal fields (7 %), and almond orchards (4 %) (Fig. 4), based on a detailed land use map (MAPAMA, 2010). Agriculture accounts for 14 % of the catchment's surface area. Huerta is defined as small-scale traditional vegetable and/or fruit orchards, which are common in the study area. The main soil classes are Leptosols (38 %), Luvisols (27 %), Cambisols (16 %), and Calcisols (11 %), based on the SoilGrids database (Hengl et al., 2017). There are five reservoirs located in the catchment (Fig. 4b) with a total capacity of 663 Hm³, which are mainly used to store water for irrigation purposes.

Figure 4 Location and characteristics of the upper Segura River catchment: (a) location of the catchment within Europe; (b) the hydrological calibration area (orange), the channels (light blue), the reservoirs (dark blue), and the calibration reservoirs (red dots); (c) digital elevation model (Farr et al., 2007); (d) land use map (MAPAMA, 2010); and (e) soil texture map (Hengl et al., 2017).

3.1 Input data

All input data were prepared at a 200 m grid size. Daily precipitation data were obtained from the SPREAD dataset (Serrano-Notivoli et al., 2017), with a 5 km spatial resolution. Daily temperature data were obtained from the Spain02 dataset (Herrera et al., 2016), with a 0.11^∘ resolution. Precipitation and temperature data were subsequently interpolated on the model grid using bivariate interpolation (Akima, 1996). Soil textural fractions (sand, clay, and silt) and soil organic matter content were obtained from the global SoilGrids dataset (Hengl et al., 2017) at 250 m resolution. The soil hydraulic properties (saturated hydraulic conductivity, saturated water content, field capacity, and wilting point) were obtained by applying pedotransfer functions (Saxton and Rawls, 2006). A digital elevation model was obtained from the SRTM dataset (Farr et al., 2007) at 30 m resolution and was resampled to the model grid using bilinear interpolation (Fig. 4d). The spatially distributed rock fraction map was obtained by applying the empirical formulations from Poesen et al. (1998), which determine rock fraction based on slope gradient.

Table 1 Input parameters for the soil erosion model.

^a T represents tillage, N.E. represents no erosion, N.V. represents no vegetation, ^b day of the year, ^c obtained from Chow (1959), and ^d not applicable.

Both the hydrological and the soil erosion model require land use specific input. We used a detailed land use map (MAPAMA, 2010) that identifies 25 land use classes within the study area. Values for the land use specific tabular value of the depletion fraction (p_tabular) to calculate actual evapotranspiration were obtained from Table 22 in Allen et al. (1998). Values for the maximum LAI (LAI_max) were obtained from Sellers et at. (1996). The soil erosion model requires land use specific input for plant height (PH), stem density (NV), stem diameter (D), ground cover fraction (GC) and, optionally, the Manning roughness coefficient for vegetation (n_vegetation). The user needs to specify whether the land use class is non-erodible (e.g. pavement and water), tilled, or non-vegetated (e.g. bare soil or tilled orchards). We obtained values for each of these parameters through observations from aerial photographs, expert judgement, and as part of the calibration procedure (Table 1). The tillage parameter RFR was set to six, which corresponds to “Cultivator tillage” in Table 4 of Morgan and Duzant (2008). The input parameters change when a crop is harvested; therefore, we varied the input parameters according to the sowing–harvest cycle representing the cropping cycle for horticulture and cereals. Table 1 shows the values of all the land use specific input parameters after calibration.

We applied the dynamic vegetation module to obtain crop coefficients and vegetation cover from NDVI, which we obtained from the 16-day Moderate Resolution Imaging Spectroradiometer data (Didan, 2015) for the period from 2000 to 2012. For model calibration (2001–2010) we used each of the individual NDVI images, after gap-filling (mainly due to cloud cover) with the long-term average 16-day period NDVI for the period from 2000 to 2012.

For the model validation period (1981–2000) no NDVI images of sufficient quality and resolution were available; therefore, we prepared NDVI model input accounting for the intra- and inter-annual variability. The intra-annual variability was obtained from the long-term average 16-day period NDVI for the period from 2000 to 2012. The inter-annual variability was determined based on a log-linear relationship between the annual precipitation sum, annual mean temperature, annual maximum temperature, and yearly averaged NDVI for each of the 25 land use classes for the period from 2000 to 2012:

$$\begin{array}{ll}{\displaystyle }{\mathrm{NDVI}}_{\mathrm{year}}=& {\displaystyle }{\mathit{\beta }}_{\mathrm{0}}+\log \left(P_{\mathrm{year}}\right){\mathit{\beta }}_{\mathrm{1}}+\log \left(P_{\mathrm{year}-\mathrm{1}}\right){\mathit{\beta }}_{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+\log \left(T{\mathrm{avg}}_{\mathrm{year}}\right){\mathit{\beta }}_{\mathrm{3}}+\log \left(T{\mathrm{avg}}_{\mathrm{year}-\mathrm{1}}\right){\mathit{\beta }}_{\mathrm{4}}\\ \text{(31)}& {\displaystyle }& {\displaystyle }+\log \left(T{\max }_{\mathrm{year}}\right){\mathit{\beta }}_{\mathrm{5}}+\log \left(T{\max }_{\mathrm{year}-\mathrm{1}}\right){\mathit{\beta }}_{\mathrm{6}},\end{array}$$

where NDVI is the yearly averaged NDVI, P is the annual precipitation sum, Tavg is the annual mean temperature, Tmax is the annual maximum temperature, and β₀₋₆ are coefficients of the log-linear model. We used the annual climate indices of two years, the current year, and the previous year, to account for the climate lag that may influence the vegetation development. A stepwise model selection procedure was applied for each of the 25 land use classes, selecting the best combination of variables from Eq. (31) with the lowest AIC (Akaike information criterion) in R (version 3.4.0), using the stepAIC algorithm from the MASS package (Venables and Ripley, 2002).

3.2 Model calibration & validation

We calibrated the model for the period from 2001 to 2010. To prevent overfitting and achieve the most realistic model calibration we set most of the potential calibration parameters at literature values and maintained the other parameters within reasonable physical limits of the parameter domain. We used daily discharge time series from the Segura River Basin Agency (Confederación Hidrográfica del Segura) for the Fuensanta Reservoir (Fig. 4b) to determine model performance for the hydrological model. The calibration procedure consisted of two steps. Firstly, we optimized by minimising the difference between the observed and simulated discharge sum at the Fuensanta discharge station. We adjusted the calibration parameter λ from Eq. 6 to obtain a surface runoff ratio between 2 and 10 %, which previous studies reported to be representative for catchments with similar conditions (Descroix et al., 2001; Love et al., 2010; Mekki et al., 2006; Reaney et al., 2014). We used parameters from the dynamic vegetation module and soil hydraulic properties to optimise the difference in discharge sum (percent bias). Secondly, using the parameter set from the first step, we optimized the Nash–Sutcliffe model efficiency (Nash and Sutcliffe, 1970) at the Fuensanta discharge station by calibrating the routing parameter (kx). The calibration resulted in a NSE of 0.47 for the daily discharge data, a NSE of 0.76 for the monthly discharge data, and a percent bias of 2.3 % (Fig. 5a). Model validation for the period from 1981 to 2000 resulted in a NSE of 0.25 for the daily discharge data, a NSE of 0.39 for the monthly discharge data, and a percent bias of −18.7 % (5b). These model efficiency values are similar to previous SPHY model applications (Terink et al., 2015).

Cerdan et al. (2010)Maetens et al. (2012)

Table 2 Calibration and validation of hillslope soil erosion and comparison with literature data (Mg km⁻² yr⁻¹), where NA stands for not available.

Figure 5 Discharge time series for the calibration (a) and validation period (b). The solid line corresponds to the observed time series and the dashed orange line corresponds to the simulated time series.

To calibrate the soil erosion model we first optimized the detached material going into transport (G) for 8 land use classes (aggregated from 25 land use classes), based on literature data (Cerdan et al., 2010; Maetens et al., 2012) (Table 2). The resulting optimized parameter values are presented in Table 1. Next, we optimized percent bias in the prediction of sediment yield at the reservoirs using measured reservoir sediment yield data from four reservoirs (Avendaño-Salas et al., 1997). Annual sediment yield at the reservoirs was determined from the measured reservoir volume at two moments in time and the measured bulk density. Figure 4b shows the locations of the reservoirs used for calibration. The calibration procedure focused on the most sensitive parameter from the sediment transport capacity equation (Eq. 28), i.e. the β parameter. We obtained a percent bias of 0.0 % in the calibration and −19.8 % in the validation (see Fig. 6).

Figure 6 Average yearly sediment yield (SY) at the reservoirs for the calibration (blue) and validation period (red).

3.3 Results

Here we present a selection of model results to illustrate the main capabilities of the SPHY–MMF model. Soil erosion shows an important intra-annual variability due to seasonal changes in climate forcing and vegetation cover (Fig. 7). For crops with little to no ground cover (i.e. tree crops and vineyard), soil erosion follows the precipitation sum, with high values in the winter, spring, and autumn months and low values in the summer months. Some crops show a distinct peak in the vegetation development in the spring (April–May), e.g. huerta and horticulture. While this period has a relatively high precipitation sum, soil erosion decreases as a consequence of the increased vegetation cover indicated by the NDVI in this period.

The temporal variation of the vegetation development of cereals and horticulture shows a slightly distinct pattern from the other land use classes. Both crops show an increase in the spring months (March–May), which indicates the rapid growth of these crops in these months. However, during the summer months (June–August) the NDVI decreases, which coincides with the period when the crops are harvested, followed by the post-harvest period. In the latter period, we assume bare soil conditions for these crops. For both crops this ultimately results in the highest annual erosion rates in the post-harvest period (October).

To illustrate the models capacity to perform scenario studies we evaluated the impacts of the application of sustainable land management and the impacts of a future climate change scenario. First, we evaluated the application of a sustainable land management scenario in cereal and horticulture fields. We assumed that after harvest the cereal and horticulture fields are not tilled until the next sowing. We parameterized this by setting the plant height to zero and leaving all other input parameters unchanged, i.e. stem density, diameter and ground cover, which in the conventional scenario are also set to zero. This will affect immediate deposition of detached particles through Eq. (26). The sustainable land management scenario leads to lower soil erosion estimates for these two land use classes (dashed lines in Fig. 7), as a result of an increase of immediate deposition of detached particles.

Figure 7 Monthly precipitation sum (mm), NDVI (–), and soil erosion (Mg km⁻² yr⁻¹) per land use class for the period from 1981 to 2000. The gray area indicates the period when cereals and horticulture are harvested and model parameters are changed to simulate bare soil conditions. The dashed lines in the lower panel show the soil erosion of cereals and horticulture under the sustainable land management scenario.

Next, we simulated the impacts of a projected climate change scenario, by comparing predicted soil erosion rates and sediment yield under the reference scenario (1981–2000) with a future scenario (2081–2100). We used a future emission scenario from the Representative Concentration Pathways (van Vuuren et al., 2011) that describes a continuous increase of GHG emissions throughout the 21st century, i.e. RCP8.5. For this exercise we used projected climate data for RCP8.5 obtained from a regional climate model (CLMcom MPI-ESM-LR) from the EURO-CORDEX initiative (Jacob et al., 2014), for the period from 2081 to 2100. The climate forcing (precipitation and temperature) was bias-corrected using quantile mapping (Themeßl et al., 2012) and we applied the log-linear relationship (Eq. 31) to construct future NDVI input based on future climate conditions. Figure 8 shows the precipitation and vegetation response under the reference and future scenarios. The annual precipitation sum decreases in the future scenario, with a catchment-averaged decrease of 128 mm. However, heavy precipitation, defined as the 95th percentile of daily precipitation, considering only rainy days (Jacob et al., 2014), increases by 27 % on average in the catchment. The NDVI increases in the western part of the catchment due to increasing temperatures and decreases in the eastern part of the catchment due to a combination of decreasing precipitation and increasing temperatures.

Figure 8 Average annual precipitation sum (mm), heavy precipitation (mm), and NDVI (–) for the reference (1981–2000) and future (2081–2100) scenarios.

In the reference scenario the highest specific sediment yield (SSY) is projected in the river network (Fig. 9), where accumulated runoff causes an increase of soil erosion rates (Eq. 19). In the future climate scenario, the catchment-median SSY increases from 43.3 to 55.2 Mg km⁻² yr⁻¹, an increase of 27.7 %. This shows that the increase in heavy precipitation has a more pronounced impact on soil erosion than the decrease of annual precipitation sum. The increase of heavy precipitation leads to both an increase of detachment by raindrop impact (Eq. 18) and an increase of detachment by runoff (Eq. 19), as a consequence of an increase in surface runoff due to infiltration excess surface runoff. However, as a result of the increased vegetation cover in the western part of the catchment SSY decreases, despite the increased extreme precipitation intensities (Fig. 8).

Reservoir sediment yield (SY) decreases in all five reservoirs between 42.4 % and 59.0 % in the future climate scenario. While it is likely that a decrease of SSY in the western part of the catchment causes a decrease of reservoir SY, it is less obvious why in the eastern part of the catchment an increase in SSY is not reflected in an increase in reservoir SY. The explanation for this lies in the fact that a decrease in precipitation sum causes a decrease of accumulated runoff and, subsequently, a decrease of sediment transport capacity (Eq. 28), increased sediment deposition, and decreased reservoir SY.

Figure 9 Specific sediment yield (Mg km⁻² yr⁻¹) and reservoir sediment yield (Gg yr⁻¹) for the reference (1981–2000) scenario and the change (%) for the future (2081–2100) scenario.

4 Discussion

The SPHY–MMF model, based on integration of the SPHY hydrological model with the MMF soil erosion model, provides an important step toward simulating the regional scale impacts of environmental change on soil erosion and sediment yield. The processes included in the model provide the flexibility and accuracy needed to reflect the impacts of intra-annual changes in land use, land management and climate (Fig. 7), and the inter-annual response of vegetation development and soil erosion to changes in climate forcing, including changes in precipitation sum and intensity (Fig. 9).

The availability of high quality input data and model parameter values are important constraints for many process-based soil erosion models. Although SPHY–MMF requires a wide range of input data, most data were obtained from publicly available global datasets at relatively high spatial and temporal resolution. Whilst climate data may often be an important constraint, long-term, high-resolution, gridded daily climate forcing datasets are becoming increasingly available at national (Berezowski et al., 2016; Kotlarski et al., 2017; Silva et al., 2007; Yin et al., 2015), (sub-)continental (Haylock et al., 2008; Mitchell, 2004; Yatagai et al., 2012; van den Besselaar et al., 2017), and even global scales (Donat et al., 2013; Huffman et al., 2001; Schamm et al., 2016). Although, the model includes a large number of parameters (Table S1), we do not expect this to limit the model's applicability. Apart from the land use specific parameters, all soil erosion model parameters were obtained from the literature, mainly from Morgan and Duzant (2008). Although some parameters may be obtained from calibration, the possibility to use literature values facilitates model application. The model requires a detailed land use map to account for land use specific model parameters that drive soil erosion processes (Table S1), which may be inevitable because soil erosion is inherently land use dependent. For example, we obtained literature values for five land use related model parameters (Table S1), while others were obtained through observations from aerial photographs, expert judgement, and as part of the calibration procedure. The ground cover for cereals, for instance, was obtained by adjusting the stem density, assuming a certain ground cover per stem. Most other input datasets are publicly available and most model parameters may be obtained from literature, which makes the model applicable for any environment.

One of the main limitations of the model is that unrealistic high soil erosion rates are projected where most flow accumulates (Fig. 9). This is mainly due to the high amounts of soil erosion predicted by Eq. (19) for large runoff volumes. Similar behaviour was reported in a daily implementation of the MMF model by Shrestha and Jetten (2018), who subsequently suggested excluding higher order streams from the soil erosion assessment. We correct for the high erosion rates in the river network by including both immediate sediment deposition (Eq. 20) and deposition when the transport capacity is exceeded (Eq. 28). While different erosional processes indeed dominate in channels and streams (i.e. bank erosion and channel incision), which are not captured by Eq. (19), high soil erosion rates may be expected in the river network since many studies stress the large contribution of gully and channel erosion to total sediment yield due to large volumes of accumulated runoff (Poesen, 2018; Poesen et al., 1996; Vanmaercke et al., 2011; de Vente et al., 2008). However, detachment and transport processes in channels are likely different from those on hillslopes and the bed material of rivers differs from the soil texture as used in the model. While Eq. (19) makes a distinction between the three textural classes, it most likely overestimates the erosion in the higher order channels, where in reality the bed material consists of coarser material, such as coarse sand and gravel, which is currently not included in the model. Therefore, in a future update of the model, we suggest the inclusion of a separate channel model to simulate the most relevant channel erosion, transport, and deposition processes more accurately. Recently, several models have been proposed that incorporate fluvial processes into soil erosion models (Arnold et al., 2012; Wei et al., 2017) and landscape evolution models (Baartman et al., 2012; Coulthard et al., 2013), which could serve as a starting point for further model development.

The model is intended to be applied at regional scales and over decadal periods that are most relevant for policy makers. This poses some limitations to the spatial and temporal representation of some essential soil erosion processes. Firstly, the spatial resolution of the model is bounded by the hydrological model. Each cell generates runoff, which is a composite of surface runoff, lateral flow, and base flow (Fig. 1b). Lateral flow and base flow are generally only generated at large spatial scales, imposing a lower limit for the spatial resolution. Therefore, Terink et al. (2015) suggests a spatial resolution in the range of 200 m to 1 km. The lower limit of the spatial resolution may affect the simulation of soil erosion processes. For instance, the detachment of soil particles by runoff is determined by slope, among other variables. The lower limit of the spatial resolution may result in a flattening of the digital elevation model, especially affecting steep slopes. Ultimately, this results in lower soil erosion projections on the most extreme slopes. Secondly, the daily time step limits the simulation of high intensity (sub-)hourly rain storms, which often cause the most soil erosion (Nearing et al., 1990). Therefore, we included a model parameter (α) describing the relation between measured hourly and daily precipitation (see Eq. 8). Nevertheless, given the spatial scale at which the model is applied, running the model with (sub-)hourly time steps may still be infeasible for various reasons. At regional scales, (sub-)hourly precipitation data are often not available, while spatial interpolation of (sub-)hourly precipitation data may introduce large uncertainties, jeopardising the model outcome. Also, climate change assessments would become challenging, since most future climate model output is only available in daily time steps.

5 Conclusions

We have presented a new coupled hydrology, soil erosion, and sediment yield prediction model (SPHY–MMF), and its application to the upper Segura catchment for different climate and land management scenarios. The model is an integration of the MMF soil erosion model into the SPHY hydrological model and simulates most relevant hydrological and soil erosion processes at a daily time step. The model considers soil detachment by raindrop and runoff, uses dynamic vegetation to simulate changes in the canopy cover, simulates saturation excess and infiltration excess surface runoff, simulates soil deposition in the cell of its origin, and routes the sediment through the river network, considering the transport capacity of the flow. The model was successfully applied in a large catchment in southeastern Spain. We have shown that the model is capable of performing scenario assessments of changes in climate and land management. Furthermore, our results show that the model simulates the soil erosion response to intra-annual variability in climate conditions and vegetation development. While multiple challenges to accurately simulating the impacts of environmental change on soil erosion and sediment yield remain, we consider the integrated SPHY–MMF model an important step toward facilitating catchment-scale scenario studies.