Global soil organic carbon removal by water erosion under climate change and land use change during AD 1850–2005

Erosion is an Earth system process that transports carbon laterally across the land surface and is currently accelerated by anthropogenic activities. Anthropogenic land cover change has accelerated soil erosion rates by rainfall and runoff substantially, mobilizing vast quantities of soil organic carbon (SOC) globally. At timescales of decennia to millennia this mobilized SOC can significantly alter previously estimated carbon emissions from land use change (LUC). However, a full understanding of the impact of erosion on land–atmosphere carbon exchange is still missing. The aim of this study is to better constrain the terrestrial carbon fluxes by developing methods compatible with land surface models (LSMs) in order to explicitly represent the links between soil erosion by rainfall and runoff and carbon dynamics. For this we use an emulator that represents the carbon cycle of a LSM, in combination with the Revised Universal Soil Loss Equation (RUSLE) model. We applied this modeling framework at the global scale to evaluate the effects of potential soil erosion (soil removal only) in the presence of other perturbations of the carbon cycle: elevated atmospheric CO₂, climate variability, and LUC. We find that over the period AD 1850–2005 acceleration of soil erosion leads to a total potential SOC removal flux of 74±18 Pg C, of which 79 %–85 % occurs on agricultural land and grassland. Using our best estimates for soil erosion we find that including soil erosion in the SOC-dynamics scheme results in an increase of 62 % of the cumulative loss of SOC over 1850–2005 due to the combined effects of climate variability, increasing atmospheric CO₂ and LUC. This additional erosional loss decreases the cumulative global carbon sink on land by 2 Pg of carbon for this specific period, with the largest effects found for the tropics, where deforestation and agricultural expansion increased soil erosion rates significantly. We conclude that the potential effect of soil erosion on the global SOC stock is comparable to the effects of climate or LUC. It is thus necessary to include soil erosion in assessments of LUC and evaluations of the terrestrial carbon cycle.

1 Introduction

Carbon emissions from land use change (LUC), recently estimated as 1.0±0.5 Pg C yr⁻¹, form the second largest anthropogenic source of atmospheric CO₂ (Le Quéré et al., 2016). However, their uncertainty range is large, making it difficult to constrain the net land–atmosphere carbon fluxes and reduce the biases in the global carbon budget (Goll et al., 2017; Houghton and Nassikas, 2017; Le Quéré et al., 2016). The absence of soil erosion in assessments of LUC is an important part of this uncertainty, as soil erosion is strongly connected to LUC (Van Oost et al., 2012; Wang et al., 2017).

The expansion of agriculture has accelerated soil erosion by rainfall and runoff significantly, mobilizing around 783±243 Pg of soil organic carbon (SOC) globally over the past 8000 years (Wang et al., 2017). Most of the mobilized SOC is redeposited in alluvial and colluvial soils, where it is stabilized and buried for decades to millennia (Hoffmann et al., 2013a, b; Wang et al., 2017). Together with dynamic replacement of removed SOC by new litter input at the eroding sites, and the progressive exposure of carbon-poor deep soils, this translocated and buried SOC can lead to a net carbon sink at the catchment scale, potentially offsetting a large part of the carbon emissions from LUC (Berhe et al., 2007; Bouchoms et al., 2017; Harden et al., 1999; Hoffmann et al., 2013a; Lal, 2003; Stallard, 1998; Wang et al., 2017).

On eroding sites, soil erosion keeps the SOC stocks below a steady state (Van Oost et al., 2012) and can lead to substantial CO₂ emissions in certain regions (Billings et al., 2010; Worrall et al., 2016; Lal, 2003). CO₂ emission from soil erosion can take place during the breakdown of soil aggregates by erosion and during the transport of the eroded SOC by runoff and later also by streams and rivers.

LUC emissions are usually quantified using bookkeeping models and LSMs that represent the impacts of LUC activities on the terrestrial carbon cycle (Le Quéré et al., 2016). These impacts are represented through processes leading to a local imbalance between net primary productivity (NPP) and heterotrophic respiration, ignoring lateral displacement. Currently, LSMs consider only the carbon fluxes following LUC resulting from changes in vegetation, soil carbon, and sometimes wood products (Van Oost et al., 2012; Stocker et al., 2014). The additional carbon fluxes associated with the human action of LUC from the removal and lateral transport of SOC by erosion are largely ignored.

In addition, the absence of lateral SOC transport by erosion in LSMs complicates the quantification of the human perturbation of the carbon flux from land to inland waters (Regnier et al., 2013). Recent studies have been investigating the dissolved organic carbon (DOC) transfers along the terrestrial–aquatic continuum in order to better quantify CO₂ evasion from inland waters and to constrain the lateral carbon flux from the land to the ocean (Lauerwald et al., 2017; Regnier et al., 2013). They point out that an explicit representation of soil erosion and transport of particulate organic carbon (POC) – in addition to DOC leaching and transport – in future LSMs is essential to be able to better constrain the flux from land to ocean. This is true, since the transfer of POC from eroded SOC forms an important part of the carbon inputs to rivers.

The slow pace of carbon sequestration by soil erosion and deposition (Van Oost et al., 2012; Wang et al., 2017) and the slowly decomposing SOC pools require the simulation of soil erosion at timescales longer than a few decades to fully quantify its impacts on the SOC dynamics. This, and the high spatial resolution that soil erosion models typically require, complicates the introduction of soil erosion and related processes in LSMs that use short time steps (≈30 min) for simulating energy fluxes and require intensive computing resources.

Previous approaches used to explicitly couple soil erosion and SOC turnover have been applying different erosion and carbon dynamics models at different spatial and temporal scales. Some studies coupled process-oriented soil erosion models with carbon turnover models calibrated for specific micro-catchments on timescales of a few decades to a millennium, (Billings et al., 2010; Van Oost et al., 2012; Nadeu et al., 2015; Wang et al., 2015a; Zhao et al., 2016; Bouchoms et al., 2017). Other studies focused on the application of parsimonious erosion–SOC-dynamics models using the Revised Universal Soil Loss Equation (RUSLE) approach together with sediment transport methods at regional or continental spatial scales (Chappell et al., 2015; Lugato et al., 2016; Yue et al., 2016; Zhang et al., 2014). However, the modeling approaches used in these studies apply erosion models that still require many variables and data input that is often not available at the global scale or for the past or the future time period. These models also run on a much higher spatial resolution than LSMs, making it difficult to integrate them with LSMs. The study of Ito (2007) was one of the first studies to couple water erosion to the carbon cycle at the global scale, using a simple modeling approach that combined the RUSLE model with a global ecosystem carbon cycle model. However, there are several unaddressed uncertainties related to his modeling approach, such as the application of the RUSLE at the global scale without adjusting its parameters.

Despite all the differences between the studies that coupled soil erosion to the carbon cycle, they all agree that soil erosion by rainfall and runoff is an essential component of the carbon cycle. Therefore, to better constrain the land–atmosphere and the land–ocean carbon fluxes, it is necessary to develop new LSM-compatible methods that explicitly represent the links between soil erosion and carbon dynamics at regional to global scales and over long timescales. Based on this, our study introduces a 4-D modeling approach that consists of (1) an emulator that simulates the carbon dynamics like in the ORCHIDEE LSM (Krinner et al., 2005), (2) the Adjusted Revised Universal Soil Loss Equation (Adj.RUSLE) model that has been adjusted to simulate global soil removal rates based on coarse-resolution data input from climate models (Naipal et al., 2015), and (3) a spatially explicit representation for LUC. This approach represents explicitly and consistently the links between the perturbation of the terrestrial carbon cycle by elevated atmospheric CO₂ and variability (temperature and precipitation change), the perturbation of the carbon cycle by LUC, and the effect of soil erosion at the global scale.

The main goal of our study is to use this new modeling approach to determine the potential effects of long-term soil erosion by rainfall and runoff without deposition or transport on the global SOC stocks under LUC, climate variability, and increasing atmospheric CO₂ levels. In order to be able to determine if global soil erosion is a net carbon source or sink, it is essential to study first how soil erosion, without deposition or transport, interacts with the terrestrial carbon cycle. Therefore, we also aim to understand the links between the different perturbations to the carbon cycle in the presence of soil erosion and to identify relevant changes in the spatial variability of SOC stocks under erosion.

Figure 1 (a) The structure of the carbon emulator (see variable names in the text, Sect. 2.3). The carbon erosion fluxes are represented by the red arrows and calculated using the soil erosion rates from the Adj.RUSLE model. (b) The land use change module of the emulator.

2 Materials and methods

2.1 Modeling framework concept

We used the LSM ORCHIDEE-MICT (Guimberteau et al., 2018; Zhu et al., 2016) (in the following simply referred to as ORCHIDEE) to construct a carbon emulator that describes the carbon pools and fluxes exactly as in ORCHIDEE (Fig. 1a). MICT stands for aMeliorated Interactions between Carbon and Temperature, and this version of ORCHIDEE has several major modifications and improvements for especially the high latitudes.

ORCHIDEE has eight biomass pools, four litter pools, of which two are aboveground and two are belowground, and three SOC pools for each land cover type (Fig. 1a). It has been extensively validated using observations on energy, water, and carbon fluxes at various eddy-covariance sites, and with measurements of atmospheric CO₂ concentration (Piao et al., 2009). The land cover types are represented by 12 plant functional types (PFTs) and an additional type for bare soil. A total of 10 PFTs represent natural vegetation and 2 represent agricultural land (C3 and C4 crop). The turnover times for each of the PFT-specific litter and SOC pools depend on their residence time modified by local soil texture, humidity, and temperature conditions (Krinner et al., 2005). The loss of biomass and litter carbon by fire is represented by the parameterization of the Spitfire model from Thonicke et al. (2011) in the full ORCHIDEE model and currently cannot be modified in our version of the emulator. Carbon losses by fire here are considered to contribute directly to the CO₂ emissions from land.

At face value, the emulator merely copies the ORCHIDEE carbon pool dynamics, and for each new atmospheric CO₂ and climate scenario a new run of the original LSM is required to build the emulator. The emulator thus reproduces exactly the carbon pool dynamics of the full LSM. The change in carbon over time for each pool of the original model is represented in the emulator by the following general mass-balance approach:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{dC}}{\mathrm{d}t}}=I\left(t\right)-k\times \mathrm{C}\left(t\right).\end{array}$$

Here, $\frac{\mathrm{dC}}{\mathrm{d}t}$ represents the change in carbon stock of a certain pool over time, calculated by the difference between the incoming flux (I(t)) and the outgoing flux (k×C(t)) to the respective pool, where k is the turnover rate. Although originally calculated by complex equations, the dynamic evolution of each pool can be described using the first-order model of Eq. (1). Complex equations, such as photosynthesis and hydrological processes, are needed to simulate realistically the carbohydrate input to carbon pools and the moisture and temperature conditions controlling litter and soil carbon decomposition over time. All the processes that determine surface and soil temperature and soil moisture are calculated by the ORCHIDEE LSM on a 30 min time step. Such a time step is needed for coupled simulations with a climate model but makes the LSM model CPU intensive. However, there is no need for such high-temporal-resolution calculations of “fast” carbon and energy fluxes to account for erosion-induced effects on SOC stocks. The addition of erosion is here supposed to impact only carbon pools and to have no feedbacks on soil moisture, soil temperature, and photosynthesis. Therefore, we decided to use the emulator concept rather than incorporating erosion processes directly into ORCHIDEE. For each carbon pool the stock and all the incoming and outgoing fluxes are derived at a daily time step from a single simulation performed with the ORCHIDEE LSM. Based on the daily output stock and fluxes, the values of the turnover rates are calculated and archived together with the input fluxes to build the emulator. Then, the emulator can be run to simulate the dynamics of all pools over long timescales without having to recompute carbon fluxes at each time step. In this way the emulator reduces the computation time of the complex ORCHIDEE model significantly and allows us to easily add and study erosion-related processes affecting the carbon dynamics of the soil. Our main objective here is to present a tool able to evaluate erosion-related carbon fluxes at global scale using a state-of-the-art LSM output and to estimate the drivers of carbon erosion at the global scale.

ORCHIDEE also includes crop harvest, defined as the harvest of aboveground biomass of agricultural PFTs and calculated based on the concept of the harvest index (HI) (Krinner et al., 2005). The HI is defined as the yield of crop expressed as a fraction of the total aboveground dry matter production (Hay, 1995). ORCHIDEE uses a fixed HI for crop of 0.45. However, Hay (1995) showed that the HI has increased significantly since 1900 for C3 crops such as wheat. In the emulator, and also in the full ORCHIDEE model, the carbon balance of agricultural lands is sensitive to crop harvest. Based on this we use the findings of Hay (1995) to change the HI of C3 crops to be temporally variable over the period 1850–2005 in the emulator, with values ranging between 0.26 and 0.46. This means that more crop biomass is harvested against what becomes litter. We only changed the HI of C3 plants, because Hay (1995) mentioned that C4 plants, such as maize, had already a high HI at the start of the last century. It should be noted that the harvest index does not vary spatially in our emulator, and harvesting is then done constantly at each time step.

2.2 Net land use change

Land use change is not taken into account in the ORCHIDEE LSM version we are using in this study to build the emulator, but is represented by a net land use change routine in the emulator that includes past agricultural land and grassland expansion over natural PFTs (Fig. 1b). This makes it possible to switch the LUC routine on or off in the emulator or to change LUC scenarios when needed without having to rerun ORCHIDEE. We verified that the LUC routine added to the emulator conserves the mass of all carbon pools for lands in transition to a new land use type. When LUC takes place, the fractions of PFTs in each grid cell are updated every year given prescribed annual maps of agricultural and natural PFTs (Peng et al., 2017). The carbon stocks of the litter and SOC pools of all the shrinking PFTs are then summed and allocated proportionally to the expanding or new PFTs, maintaining the mass balance (Houghton and Nassikas, 2017; Piao et al., 2009). When natural vegetation is reduced by LUC, the heartwood and sapwood biomass pools are harvested and transformed to three wood products with turnover times of 1, 10 and 100 years. The other biomass pools (leafs, roots, sapwood belowground, fruits, heartwood belowground) are transformed to metabolic or structural litter and allocated to the respective litter pools of the expanding PFTs (Piao et al., 2009).

2.3 Soil carbon dynamics

The change in the carbon content of the PFT-specific SOC pools in the emulator without soil erosion can be described with the following differential equations:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{{\mathrm{dSOC}}_{\mathrm{a}}\left(t\right)}{\mathrm{d}t}}={\mathrm{lit}}_{\mathrm{a}}\left(t\right)+k_{\mathrm{pa}}\times {\mathrm{SOC}}_{\mathrm{p}}\left(t\right)+k_{\mathrm{sa}}\times {\mathrm{SOC}}_{\mathrm{s}}\left(t\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }-\left(k_{\mathrm{ap}}+k_{\mathrm{as}}+k{\mathrm{0}}_{\mathrm{a}}\right)\times {\mathrm{SOC}}_{\mathrm{a}}\left(t\right),{\displaystyle }& {\displaystyle \frac{{\mathrm{dSOC}}_{\mathrm{s}}\left(t\right)}{\mathrm{d}t}}={\mathrm{lit}}_{\mathrm{s}}\left(t\right)+k_{\mathrm{as}}\times {\mathrm{SOC}}_{\mathrm{a}}\left(t\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }-\left(k_{\mathrm{sa}}+k_{\mathrm{sp}}+k{\mathrm{0}}_{\mathrm{s}}\right)\times {\mathrm{SOC}}_{\mathrm{s}}\left(t\right),{\displaystyle }& {\displaystyle \frac{{\mathrm{dSOC}}_{\mathrm{p}}\left(t\right)}{\mathrm{d}t}}=k_{\mathrm{ap}}\times {\mathrm{SOC}}_{\mathrm{a}}\left(t\right)+k_{\mathrm{sp}}\times {\mathrm{SOC}}_{\mathrm{s}}\left(t\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }-(k_{\mathrm{pa}}+k{\mathrm{0}}_{\mathrm{p}})\times {\mathrm{SOC}}_{\mathrm{p}}\left(t\right),\end{array}$$

where SOC_a, SOC_s, and SOC_p (g C m⁻²) are the active (unprotected); slow (physically or chemically protected); and passive (biochemically recalcitrant) SOC, respectively. The SOC pools are based on the study of Parton et al. (1987) and are defined by their residence times. The active SOC pool has the lowest residence time (∼ 1.5 years) and the passive the highest (∼ 1000 years). lit_a and lit_s (g C m⁻² day⁻¹) are the litter input rates to the active and slow SOC pools, respectively; k0_a, k0_s, and k0_p (day⁻¹) are the respiration rates of the active, slow, and passive pools, respectively; k_as, k_ap, k_pa, k_sa, k_sp are the coefficients determining the flux from the active to the slow pool, from the active to the passive pool, from the passive to the active pool, from the slow to the active pool, and from the slow to the passive pool, respectively (Fig. 1a).

The SOC pools are not vertically discretized in the version of the ORCHIDEE LSM used to build the emulator, so we implemented a simple vertical discretization scheme for the SOC pools in the emulator based on the concept of Wang et al. (2015a, b). In this scheme the carbon dynamics of each soil layer are calculated separately, based on layer-dependent litter input and respiration rates (Fig. 1a). The vertical discretization scheme of the emulator does not change the total input and respiration as simulated by ORCHIDEE in the case where erosion and land use change processes are switched off. We apply the same scheme for all three SOC pools, assuming that each SOC pool is equally distributed across all layers of the soil profile, while the ratios between the pools per soil layer are equal to those from ORCHIDEE. We base this assumption on the fact that there is very little information or data to constrain the pool ratios globally, mainly because the three SOC pools cannot be directly related to measurements (Elliott et al., 1996). Furthermore, neither the emulator nor the ORCHIDEE LSM model we used include soil processes that may affect these pool ratios with depth, such as vertical mixing by soil organisms, diffusion, leaching, changes in soil texture (SOC protection and stabilization by clay particles), and limitations by oxygen and by access to deep organic matter by microbes. It is also uncertain how sensitive SOC is to these processes. For example, the study of Huang et al. (2018), who implemented a matrix-based approach to assess the sensitivity of SOC, showed that equilibrium SOC stocks are more sensitive to input than to mixing for soils in the temperate and high-latitude regions.

In the vertical discretization scheme of the emulator, the soil profile is divided into thin layers of 1 cm thickness down to a depth of 2 m, which is the soil depth used by ORCHIDEE to calculate SOC. The first 10 cm of the soil profile are referred to as the “topsoil”, where we assume that the SOC content is homogeneously distributed. The rest of the soil profile is referred to as the subsoil. The topsoil receives carbon from above- and belowground litter, which is homogeneously distributed across the soil layers of the topsoil.

The belowground litter input for the active SOC pool is the sum of a fraction of the belowground structural and metabolic litter pools of ORCHIDEE being recalculated by the emulator to change with depth, while the belowground litter input for the slow SOC pool is equal to a fraction of the belowground structural litter pool only. This setting is consistent with the structure of the SOC module of the ORCHIDEE LSM to ensure that the emulator reproduces the same C pool dynamics as the LSM. The litter fractions are based on the Century model as introduced by Parton et al. (1987) and later implemented inside ORCHIDEE (Krinner et al., 2005). We assume that the subsoil receives carbon only from belowground litter and that this input decreases exponentially with depth following the root profile of ORCHIDEE. This discretization of the total belowground litter input (lit_be) is the same for both SOC pools and can then be represented as

$$\begin{array}{}\text{(5)}& {\mathrm{lit}}_{\mathrm{be}}={\int }_{z=\mathrm{0}}^{z=z_{max}}{I}_{\mathrm{0}\mathrm{be}}\times e^{-r\times z}\mathrm{d}z,\end{array}$$

where I_0be is the belowground litter input to the surface layer and is equal to

$$\begin{array}{}\text{(6)}& I_{\mathrm{0}\mathrm{be}}={\mathrm{lit}}_{\mathrm{be}}\times {\displaystyle \frac{r}{\mathrm{1}-e^{-r\times z_{max}}}}.\end{array}$$

The homogeneously distributed belowground litter input (I_be) to the layers of the topsoil is equal to

$$\begin{array}{}\text{(7a)}& \sum _{z=\mathrm{0}}^{z=\mathrm{10}}{I}_{\mathrm{0}\mathrm{be}}\times e^{-r\times z}\times \mathrm{0.1}.\end{array}$$

The belowground litter input to the layers of the subsoil is equal to

$$\begin{array}{}\text{(7b)}& I_{\mathrm{be}}\left(z\right)=I_{\mathrm{0}\mathrm{be}}\times e^{-r\times z},\end{array}$$

where z_max is the maximum soil depth equal to 2 m and dz is the soil layer discretization (1 cm); r is the PFT-specific vertical root-density attenuation coefficient as used in ORCHIDEE.

The SOC respiration rates of ORCHIDEE are determined by soil temperature, moisture, and texture. In the emulator we conserve the total respiration of ORCHIDEE (when LUC or erosion is absent) and assume that the respiration is homogeneously distributed across the layers of the topsoil. For the rest of the soil profile the respiration rates of all three SOC pools decrease exponentially with depth:

$$\begin{array}{}\text{(8)}& k_i\left(z\right)=k_{\mathrm{0}i}\times e^{-r_{\mathrm{e}}z}.\end{array}$$

Here k_0 i is the SOC respiration rate at the surface layer for each SOC pool (i=a, s, p), and r_e (m⁻¹) is a coefficient representing the impact of external factors, such as oxygen availability, which reduce SOC mineralization rate with depth (z). To ensure that the total soil respiration of the emulator is similar to that of the ORCHIDEE LSM model for each grid cell, each PFT, and each SOC pool, we have calibrated the exponent “r_e” and variable “k_0i” of Eq. (8) for each grid cell and PFT under equilibrium conditions. First we selected a default value for r_e between 0 and 5 and calculated the respiration rate of the surface soil layer (k₀) when all SOC pools are in an equilibrium state, with the following equation:

$$\begin{array}{}\text{(9)}& {\mathrm{SOC}}_{\mathrm{orchidee}}={\sum }_{z=\mathrm{0}}^{z=n}{\displaystyle \frac{L\left(z\right)}{k_{\mathrm{0}}\times e^{-r_{\mathrm{e}}\times z}}},\end{array}$$

where SOC_orchidee is the total equilibrium SOC stock derived from ORCHIDEE for a certain grid cell and PFT. L(z) is the total litter input to the soil for a certain soil layer discretized according to the root profile. Then we derived the equilibrium SOC stocks per soil layer as

$$\begin{array}{}\text{(10)}& \mathrm{SOC}\left(z\right)={\displaystyle \frac{L\left(z\right)}{k_{\mathrm{0}}\times e^{-r_{\mathrm{e}}\times z}}}.\end{array}$$

Assuming that the ratios between the active, slow, and passive SOC pools do not change with depth and are equal to the ratios derived from ORCHIDEE, we calculated the SOC stocks of each pool with the following equation:

$$\begin{array}{}\text{(11)}& \mathrm{1}+{\displaystyle \frac{{\mathrm{soil}}_{\mathrm{s}}\left(z\right)}{{\mathrm{soil}}_{\mathrm{a}}\left(z\right)}}+{\displaystyle \frac{{\mathrm{soil}}_{\mathrm{p}}\left(z\right)}{{\mathrm{soil}}_{\mathrm{a}}\left(z\right)}}={\displaystyle \frac{\mathrm{SOC}\left(z\right)}{{\mathrm{soil}}_{\mathrm{a}}\left(z\right)}},\end{array}$$

where soil_a(z), soil_s(z), and soil_p(z) are the emulator-derived active, slow, and passive SOC stocks per soil layer, grid cell, and PFT. Now, for the equilibrium state the input is equal to the output, allowing us to derive k_0a, k_0s, and k_0p with the following equations:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\sum }_{z=\mathrm{0}}^{z=n}\left({\displaystyle \frac{L_{\mathrm{a}}\left(z\right)+k_{\mathrm{sa}}\times {\mathrm{soil}}_{\mathrm{s}}\left(z\right)+k_{\mathrm{pa}}\times {\mathrm{soil}}_{\mathrm{p}}\left(z\right)}{k_{\mathrm{0}\mathrm{a}}\times e^{-r_{\mathrm{e}}\times z}+k_{\mathrm{as}}+k_{\mathrm{ap}}}}\right)\\ \text{(12a)}& {\displaystyle }& {\displaystyle }={\mathrm{SOC}}_{\mathrm{a}},\text{(12b)}& {\displaystyle }& {\displaystyle }{\sum }_{z=\mathrm{0}}^{z=n}\left({\displaystyle \frac{L_{\mathrm{s}}\left(z\right)+k_{\mathrm{as}}\times {\mathrm{soil}}_{\mathrm{a}}\left(z\right)}{k_{\mathrm{0}s}\times e^{-r_{\mathrm{e}}\times z}+k_{\mathrm{sa}}+k_{\mathrm{sp}}}}\right)={\mathrm{SOC}}_{\mathrm{s}},\text{(12c)}& {\displaystyle }& {\displaystyle }{\sum }_{z=\mathrm{0}}^{z=n}\left({\displaystyle \frac{k_{\mathrm{sp}}\times soi{l}_{\mathrm{s}}\left(z\right)+k_{\mathrm{ap}}\times soi{l}_{\mathrm{a}}\left(z\right)}{k_{\mathrm{0}\mathrm{p}}\times e^{-r_{\mathrm{e}}\times z}+k_{\mathrm{pa}}}}\right)={\mathrm{SOC}}_{\mathrm{p}},\end{array}$$

where L_a is the total litter input to the active SOC pool and L_s is the total litter input to the slow SOC pool. SOC_a, SOC_s, and SOC_p are the total active, slow, and passive SOC stocks per grid cell and PFT, respectively, derived from ORCHIDEE. k_as, k_ap, k_sa, k_sp, k_pa are the coefficients determining the fluxes between the SOC pools. After the derivation of k_0i we tested if the difference in the SOC stocks between the emulator and the original ORCHIDEE LSM is less than 1 g m⁻² per grid cell and PFT. If this was not the case we increased or decreased the value of “r_e” and repeated the calibration cycle. If we did not find an optimized value for both r_e and k_0i that meet this criteria, we used values that minimized the difference in SOC stocks between the emulator and the original ORCHIDEE LSM.

For the transient period (without LUC or erosion) we assumed a time-constant r_e, where the values are equal to those at equilibrium. Using the mass-balance approach we calculated the daily values for k_0a, k_0s, k_0p per grid cell and PFT with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{{\mathrm{dSOC}}_{\mathrm{a}}}{\mathrm{d}t}}={\sum }_{z=\mathrm{0}}^{z=n}\left(L_{\mathrm{a}}\left(z,t\right)\right.\\ {\displaystyle }& {\displaystyle }+k_{\mathrm{sa}}\times {\mathrm{soil}}_{\mathrm{s}}\left(z,t-\mathrm{1}\right)+k_{\mathrm{pa}}\times {\mathrm{soil}}_{\mathrm{p}}\left(z,t-\mathrm{1}\right)\\ \text{(13a)}& {\displaystyle }& {\displaystyle }\left.-\left(k_{\mathrm{0}\mathrm{a}}\left(t\right)\times e^{-r_{\mathrm{e}}\times z}+k_{\mathrm{as}}+k_{\mathrm{ap}}\right)\times {\mathrm{soil}}_{\mathrm{a}}\left(z,t-\mathrm{1}\right)\right),{\displaystyle }& {\displaystyle \frac{{\mathrm{dSOC}}_{\mathrm{s}}}{\mathrm{d}t}}={\sum }_{z=\mathrm{0}}^{z=n}\left(L_{\mathrm{s}}\left(z,t\right)+k_{\mathrm{as}}\times {\mathrm{soil}}_{\mathrm{a}}\left(z,t-\mathrm{1}\right)\right.\\ \text{(13b)}& {\displaystyle }& {\displaystyle }\left.-\left(k_{\mathrm{0}\mathrm{s}}\left(t\right)\times e^{-r_{\mathrm{e}}\times z}+k_{\mathrm{sa}}+k_{\mathrm{sp}}\right)\times {\mathrm{soil}}_{\mathrm{s}}\left(z,t-\mathrm{1}\right)\right),{\displaystyle }& {\displaystyle \frac{{\mathrm{dSOC}}_{\mathrm{p}}}{\mathrm{d}t}}={\sum }_{z=\mathrm{0}}^{z=n}\left(k_{\mathrm{sp}}\times {\mathrm{soil}}_{\mathrm{s}}\left(z,t-\mathrm{1}\right)\right.\\ {\displaystyle }& {\displaystyle }+k_{\mathrm{ap}}\times {\mathrm{soil}}_{\mathrm{a}}\left(z,t-\mathrm{1}\right)-\left(k_{\mathrm{0}\mathrm{p}}\left(t\right)\times e^{-r_{\mathrm{e}}\times z}+k_{\mathrm{pa}}\right)\\ \text{(13c)}& {\displaystyle }& {\displaystyle }\left.\times {\mathrm{soil}}_{\mathrm{p}}\left(z,t-\mathrm{1}\right)\right).\end{array}$$

In case there was no solution for the k_0i at a certain time step we took the values from the previous time step.

The annual average soil erosion rate (E, t ha⁻¹ year⁻¹) is calculated by the Adj.RUSLE (Naipal et al., 2015, 2016) according to

$$\begin{array}{}\text{(14)}& E=S\times R\times K\times C,\end{array}$$

where R is the rainfall erosivity factor (MJ mm ha⁻¹ h⁻¹ year⁻¹), K is the soil erodibility factor (t ha h ha⁻¹ MJ⁻¹ mm⁻¹), C is the land cover factor (dimensionless), and S is the slope steepness factor (dimensionless). The S factor is calculated using the slope from a 1 km resolution digital elevation model (DEM) that has been scaled using the fractal method to a resolution of 150 m (Naipal et al., 2015). In this way the spatial variability of a high-resolution slope dataset can be captured. The computation of the R factor has been adjusted to use coarse-resolution input data on precipitation and to provide reasonable global erosivity values. For this Naipal et al. (2015) derived regression equations for different climate zones of the Köppen–Geiger climate classification (Peel et al., 2007). The results from the Adj.RUSLE model have been tested against empirical large-scale assessments of soil erosion and rainfall erosivity (Naipal et al., 2015, 2016). The original RUSLE model as described by Renard et al. (1997) also includes the slope-length (L) and support-practice (P) factors. Although these factors can strongly affect soil erosion in certain regions, the Adj.RUSLE does not include these factors due to several reasons. Firstly, Doetterl et al. (2012) showed that these factors do not significantly contribute to the variation in soil erosion at the continental to global scales, in comparison to the other RUSLE factors. Secondly, data on the L and P factors and methods to estimate them at the global scale are very limited. Thus, including them in global soil erosion estimations would result in large uncertainties. Finally, the focus of this study is to show the effects of potential soil erosion on the terrestrial carbon cycle, without the explicit effect of management practices such as covered by the P factor. For more information on the validation of our erosivity values and a more detailed description of the calculation of each of the RUSLE factors see Supplement Sect. S1.

The Adj.RUSLE model is not imbedded in the C emulator but is run separately on a 5 arcmin spatial resolution and at a yearly timestep. The resulting soil erosion rates are then read by the C emulator at each time step and used to calculate the daily SOC erosion rate of a certain SOC pool i (Ce_i in g C m⁻² day⁻¹) at the surface layer by

$$\begin{array}{}\text{(15)}& {\mathrm{Ce}}_i={\mathrm{SOC}}_i\times {\displaystyle \frac{\frac{E}{\mathrm{365}}\times \mathrm{100}}{{\mathrm{BD}}_{\mathrm{top}}\times \mathrm{d}z\times {\mathrm{10}}^{\mathrm{6}}}},\end{array}$$

where BD_top is the bulk density of the surface layer (g cm⁻³). We assume that the enrichment ratio, i.e., the volume ratio of the carbon content in the eroded soil to that of the source soil material, is equal to 1 here, which implies that our estimates of SOC mobilization are likely conservative (Chappell et al., 2015; Nadeu et al., 2015).

When erosion takes place, the surface layer is truncated by the erosion height, and at the same time an amount of SOC corresponding to this erosion height is removed. As we assume that the soil layer thickness does not change, part of the SOC of the next soil layer is allocated to the surface layer proportional to the erosion height and the SOC concentration (per volume) of the next layer. In this way, SOC from all the following soil layers moves upward and becomes exposed to erosion in the surface layer at some point in time (Fig. 1a). To preserve mass balance, we assume that there is no SOC below the 2 m soil profile represented in the emulator and new substrate replacing the material of the last soil layer is SOC free, so that SOC in the bottom layer will decrease towards zero after erosion has started.

Figure 2 Spatial patterns of the difference in forest, crop, and grassland area between 1851 and 2005 represented as a fraction of a grid cell. And spatial patterns of the change in average annual precipitation between 1851 and 2005 in mm yr⁻¹, calculated as the total change in precipitation over the period 1850–2005 and divided by the number of years in this period.

2.4 Input datasets

2.4.1 For ORCHIDEE

We used 6-hourly climate data supplied by the CRU-NCEP (version 5.3.2) global database (https://crudata.uea.ac.uk/cru/data/ncep/; last access: 16 July 2018) available at 0.5^∘ resolution to perform simulations with the full ORCHIDEE model for constructing the emulator. CRU-NCEP climate data were only available for the period 1901–2012. To be able to run ORCHIDEE for the period 1850–1900, we randomly projected the climate forcing after 1900 to the years before 1900. The random projection of the climate data is necessary to avoid the risk of including the effects of extreme climate conditions multiple times when only a certain decade is used repeatedly.

The historical changes in PFT fractions were derived from the historical annual PFT maps of Peng et al. (2017). These PFT maps were available at a resolution of 0.25^∘ (Fig. 2) and were regridded to the resolution of the ORCHIDEE emulator, which is 2.5×3.75, using the nearest-neighbor approach.

2.4.2 For the Adj.RUSLE

Due to the resolution of the Adj.RUSLE, which is 5 arcmin (∼0.0833^∘), all the RUSLE factors had to be regridded or calculated at this specific resolution before calculating the soil erosion rates.

The land cover fractions from the historical 0.25^∘ PFT maps were used in combination with the LAI values from ORCHIDEE at the resolution of $\mathrm{2.5}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$ to derive the values for the C factor of the RUSLE model. We first regridded the yearly average LAI to the resolution of the PFT maps before calculating the land cover factor of RUSLE (C factor) at the resolution of 0.25^∘. The C values were then regridded using the nearest-neighbor method to the resolution of the Adj.RUSLE model. We used the nearest-neighbor approach here, because the C factor is strongly dependent on the land cover class.

Daily precipitation data for the period 1850–2005 to calculate soil erosion rates is derived from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP), product ISIMIP2b (Frieler et al., 2017). These data are based on model output of the Coupled Model Intercomparison Project Phase 5 (CMIP5 output of IPSL-CM5A-LR (Taylor et al., 2012), which are bias corrected using observational datasets and the method of Hempel et al. (2013) and made available at a resolution of 0.5^∘ (Fig. 2). We chose these data as input to the Adj.RUSLE model, because the dataset extended to 1850, in contrast to the CRU-NCEP data. Also, this dataset being bias corrected provides a better distribution of extreme events and frequencies of dry and wet days (Frieler et al., 2017), which is important for the calculation of rainfall erosivity (R factor). The ISIMIP precipitation data were regridded using the bilinear interpolation method to the resolution of the Adj.RUSLE model, before being used to calculate the R factor. This was necessary because the erosivity equations from the Adj.RUSLE model are calibrated at this specific resolution (Naipal et al., 2015).

Data on soil bulk density and other soil parameters to calculate the soil erodibility factor (K), available at the resolution of 1 km, have been taken from the Global Soil Dataset for use in Earth System Models (GSDE) (Shangguan et al., 2014). The K factor has been calculated at the resolution of 1 km before being regridded to 5 arcmin using the bilinear interpolation method. We also used the SOC concentration in the soil from GSDE, which was derived using the “aggregating first” approach, to compare to our SOC stocks from simulations with the emulator. Finally, the slope steepness factor (S), which was originally estimated at the resolution of 1 km, was also regridded to the resolution of 5 arcmin using the bilinear interpolation method.

Using the above-mentioned data, soil erosion rates were first calculated at the resolution of 5 arcmin and afterwards aggregated to the coarse resolution of the emulator ($\mathrm{2.5}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$) to calculate daily SOC erosion rates.

Table 1 Description of the simulations used in this study. The S1 simulation is also the control simulation (CTR). S1 and S2 minimum and maximum are the sensitivity simulations using minimum or maximum soil erosion rates. “Best” stands for the best estimated soil erosion rates using optimal values for the R and C factors of the Adj.RUSLE model.

2.5 Model simulations

To be able to understand and estimate the different direct and indirect effects of soil erosion on the SOC dynamics, we propose a factorial simulation framework (Fig. 3 and Table 1). This framework allows isolating or combining the main processes that link soil erosion to the SOC pool, namely the influence from climate variability, LUC, and atmospheric CO₂ increase. The different model simulations described in this section will be based on this framework.

We performed two different simulations with the full ORCHIDEE model to produce the required data input for the emulator for the period 1850–2005. For this we first performed a spinup with ORCHIDEE to get steady-state carbon pools for the year 1850. We chose the period 1850–2005 based on the ISIMIP2b precipitation data availability and the fact that this period underwent a significant intensification in agriculture globally and a substantial rise in atmospheric CO₂ concentrations. In the first simulation of ORCHIDEE the global atmospheric CO₂ concentration was fixed to the year 1850 to calculate time-varying NPP not impacted by CO₂ fertilization and subsequent carbon pools, while in the second simulation the atmospheric CO₂ concentration was made variable. In both simulations, climate is variable and derived from CRU-NCEP (Fig. 3).

Figure 3 Conceptual diagram of SOC affected by erosion in the presence of other perturbations of the carbon cycle, namely climate variability, increasing atmospheric CO₂ concentrations, and land use change. A separation of these components and of the role of erosion is obtained with the factorial simulations (S1–S8), presented in Table 1.

Furthermore, we performed seven simulations with the Adj.RUSLE model to precalculate the soil erosion rates that will be used as input to the ORCHIDEE emulator. Three of the seven erosion simulations used best estimates for each model parameter, and the rest used either the minimum or maximum values for the R and C factors to derive an uncertainty range for our soil erosion rates and to analyze the sensitivity of the emulator. In the first simulation with the best estimated model parameters we kept the climate and land cover variable through time (the “CC + LUC” simulation). In the second simulation we only varied the climate through time and kept land cover fractions fixed to 1850 (the “CC” simulation, Fig. 3). In the third simulation we only varied the land cover through time and kept the climate constant to the average cyclic variability of the period 1850–1859 (the “LUC” simulation, Fig. 3).The erosion simulations with either minimum or maximum model parameters were either a CC + LUC or a CC type of simulation.

Figure 4 (a) Global annual soil erosion rates, (b) global annual SOC erosion rates, (c) agricultural annual soil erosion rates, (d) agricultural annual SOC erosion rates, (e) grassland annual soil erosion rates, (f) grassland annual SOC erosion rates, (g) forest annual soil erosion rates, and (h) forest annual SOC erosion rates over the period 1850–2005 for scenarios with only LUC (green lines), the scenario with only climate and CO₂ change (blue line), and the scenario with LUC, climate, and CO₂ change (red line). In (a), (c), (e), and (g) the dashed green line is the difference between the CTR and S2 simulations, while the straight green line is the LUC-only simulation with the Adj.RUSLE model.

From the two simulations of ORCHIDEE (with variable and constant CO₂) and the seven soil erosion simulations of the Adj.RUSLE, we constructed four versions of the emulator to perform eight main simulations and four sensitivity simulations. The different simulations and their description are given in Table 1 and Fig. 3. In the simulations without LUC (S2, S4, S6, and S8), the PFT fractions and the harvest index are constant and equal to those in the year 1850. In the simulations with LUC (S1, S3, S5, and S7) the harvest index increases and the PFT fraction changes with time during 1850–2005. In each emulator simulation we first calculated the equilibrium carbon stocks analytically before calculating the change in the carbon stocks in time depending on the perturbations during the transient period (1851–2005). In simulations with erosion, the equilibrium state of the SOC pools has been calculated using the average erosion rates of the period 1850–1859, assuming erosion to be constant before 1850 and a steady-state condition where erosion fluxes are equal to input from litter.

Table 2 Area-weighted average and standard deviation of soil and SOC erosion rates per land cover type for the year AD 2005; the uncertainty range for soil erosion rates is 25 %–53  % and for SOC erosion rates 16 %–50 %.

3 Results

3.1 Erosion versus no erosion

After including soil erosion in the ORCHIDEE emulator we obtain a total global soil loss flux of 47.6±10 Pg C year⁻¹ for the year 2005, of which 20 % to 29 % is attributed to agricultural land and 51 % to 55 % to grassland. This global soil loss flux (here “loss” meaning horizontal removal of SOC by erosion) leads to a total SOC loss flux of 0.52±0.14 Pg C year⁻¹, of which 26 % to 33 % is attributed to agricultural land and 54 % to 64 % to grassland (CTR, Fig. 4). Grassland and agricultural land thus have much larger annual average soil and SOC erosion rates compared to forest (Table 2).

Table 3 Model estimates per continent of area-weighted average annual soil erosion and SOC erosion rates for the year 2005, their spatial standard deviations, and the changes in average soil and SOC erosion rates since 1851; the uncertainty range for soil and SOC erosion rates is 2 %–36 % and 3 %–52 %, respectively. The uncertainty range for the changes in soil and SOC erosion rates since 1851 is 3 %–83 % and 11 %–166 %, respectively.

The total soil and SOC losses in the year 2005 show an increase of 11 %–19 % and 23 %–35 %, respectively, compared to 1850 (CTR, Fig. 4), with the largest increases found in the tropics (Fig. 5b, d). The largest increase in soil and SOC erosion during 1850–2005 is found in South America (Table 3) despite the significant decreases in simulated precipitation leading to less intense erosion rates in this region. One should keep in mind that, due to uncertainties in the simulated LUC and climate variability for certain regions and the assumptions made in our modeling framework, these trends in soil and SOC erosion rates are linked to some uncertainty. However, it is difficult to assess this uncertainty, mainly due to the lack of observations for the past in regions such as the tropics and the lack of model testing in these regions.

Figure 5 (a) Average annual soil erosion rates at a 5 arcmin resolution in the year 2005, (b) change in average annual soil erosion rates over the period 2005–1850, (c) average annual SOC erosion rates at a resolution of 2.5×3.75^∘ in 2005, (d) change in average annual SOC erosion rates over the period 2005–1850, and (e) difference in SOC stocks at a resolution of 2.5×3.75^∘ between the year 2005 and 1850 (CTR simulation). For the SOC stocks positive values (green color) indicate a gain, while negative values (red color) indicate a loss. For the erosion rates positive values (red color) indicate an increase over 1850–2005, while negative values (green color) indicate a decrease over 1850–2005.

We found that the total soil erosion flux on agricultural land increased by 55 %–58 % in the year 2005 compared to 1850, while the SOC erosion flux increased by 11 %–70 % (Fig. 4) and led to a cumulative SOC removal of 22±5 Pg C (CTR). On grassland the soil erosion flux increased only by 8 %–20 %, while the SOC erosion flux increased by 44 %–54 % (Fig. 4) and led to a cumulative SOC mobilization of 38±7 Pg C since 1850. It is evident that on agricultural land the uncertainty range of soil erosion leads to a large uncertainty range in SOC erosion compared to grassland. The increase in SOC erosion is much larger than the increase in soil erosion for grasslands because in our model LUC (without erosion) leads to a significant increase in SOC on grassland amplifying the increasing trend in SOC erosion for grassland. This simulated increase in SOC stocks on grasslands after LUC is not unrealistic, as it is observed from paired chronosequences worldwide where grasslands have higher SOC densities than forests, for instance (Li et al., 2017).

In total 7183±1662 Pg of soil and 74±18 Pg of SOC is mobilized across all PFTs by erosion during the period 1850–2005, which is equal to approximately 46 %–74 % of the total net flux of carbon lost as CO₂ to the atmosphere due to LUC (net LUC flux) over the same period estimated by our study (S1–S2). In this study, we do not address the fate of this large amount of eroded SOC, be it partly sequestered (Wang et al., 2017) or released to the atmosphere as CO₂.

3.2 Validation of model results

We calculated a total global SOC stock for the year 2005 in the absence of soil erosion (S3) of 1284 Pg C, which is a factor of 0.73 lower than the total SOC stock from GSDE (Shangguan et al., 2014) for a soil depth of 2 m (Table 5). Including soil erosion (S1, minimum, maximum) leads to a total SOC stock of 1001±58 Pg C for the year 2005 (Table 5). We also find that including soil erosion in the SOC-dynamics scheme slightly improves the root mean square error (RMSE) between the simulated SOC stocks and those from GSDE for the top 30 cm of the soil profile. This improvement in the RMSE occurs especially in highly erosive regions. Furthermore, the total SOC stock of agricultural land is significantly lower than that of the GSDE, because we assume a steady-state landscape at 1850, where soil erosion losses are equal to the carbon input to the soil. We did not perform a more in-depth comparison with SOC global observations as our emulator and the original ORCHIDEE LSM do not include various soil processes that have been proven to affect SOC substantially such as vertical mixing, diffusion, priming, changes in soil texture, and carbon-rich organic soils formation. The ORCHIDEE LSM model we use to build the emulator also lacks processes such as nitrogen and phosphorus limitations, which affect the productivity and SOC decomposition (Goll et al., 2017; Guenet et al., 2016). The emulator also simulates only the removal of SOC but not the subsequent SOC transport and deposition after erosion, and there is a general uncertainty in the simulation of underlying processes that govern the SOC dynamics (Todd-Brown et al., 2014). Finally, large uncertainties in the global soil databases (Hengl et al., 2014; Scharlemann et al., 2014; Tifafi et al., 2018) complicate the exact quantification of the uncertainties of the resulting SOC dynamics simulated by our emulator.

Table 4 Model estimates per continent of changes in SOC stocks since 1851 from simulations S1, S2, S1–S2, S3, S4, and S3–S4.

Table 5 Statistics of a grid cell by grid cell comparison of global SOC stocks between GSDE soil database and simulations S1 (with erosion) and S3 (without erosion). RMSE is the root mean square error and r value is the correlation coefficient of the linear regression between GSDE and S1 or S3.

Using the Adj.RUSLE model to estimate agricultural soil loss by water erosion for the year 2005 resulted in a global soil loss flux of 12.28±4.62 Pg year⁻¹ (Fig. 4). This flux is paralleled by a SOC loss flux of 0.16±0.06 Pg C year⁻¹ after including soil erosion in the CTR simulation (Fig. 4). This soil loss flux is in the same order of magnitude as earlier high-resolution assessments of this flux, while the SOC removal flux is slightly lower compared to previously published high-resolution estimates, but within the uncertainty (Table 6). We also find a fair agreement between our model estimates of recent agricultural soil and SOC erosion fluxes per continent and the high-resolution estimates (excluding tillage erosion) from the study of Doetterl et al. (2012) (Table 7). However, the continental SOC erosion fluxes from our study are generally lower, because of the lower SOC stocks on agricultural land. Only South America shows a higher SOC flux for the present day compared to the high-resolution estimates of Doetterl et al. (2012), which is the result of the simulated high productivity of crops in the tropics.

Furthermore, we find a cumulative soil loss of 1888±753 Pg and cumulative SOC removal flux of 22±5 Pg C from agricultural land over the entire time period (CTR simulation). This soil loss flux lies in the range of 2480±720 Pg found by Wang et al. (2017) for the same time period, while the SOC removal flux is significantly lower than the 63±19 Pg C found by Wang et al. (2017). Wang et al. (2017) used only recent climate data in his study while we explicitly include the effects of changes in precipitation and temperature on global soil erosion rates and the SOC stocks in our study, which may partly explain this difference.

Figure 6 Cumulative SOC stock changes during 1850–2005 for (a) simulations with variable atmospheric CO₂ concentration and (b) for simulations with a constant CO₂ concentration, implied by variable land cover alone (dash–dotted lines), by variable climate (dashed lines), and variable land cover and climate (straight lines), without erosion (black lines) and with erosion (red lines).

4 Discussions

4.1 Significance of including soil erosion in the ORCHIDEE emulator

To estimate the net effect of soil erosion on the global SOC stocks under all perturbations we compare the cumulative SOC stock change from simulation S3 (no erosion; Table 1) with that of the CTR simulation with all factors included, that is, land use, CO₂, and climate. When considering our best estimated soil erosion rates and assuming that the SOC mobilized by soil erosion in the CTR simulation is all respired, we find an overall global SOC stock decrease that is 62 % larger compared to a world without soil erosion during the period 1850–2005 (Fig. 6a). This assumption is certainly an extreme and unrealistic assumption, as in reality a fraction of the mobilized SOC will remain stored on land, but we take this assumption as an extreme scenario. Including soil erosion in the SOC cycling scheme under the previously mentioned assumption thus reduces the global land C sink, with the largest impact observed for Asia, where the decrease in the total SOC stock is 156 % larger when the effects of soil erosion are taken into account (Table 4, Fig. 8a). Some regions, such as Western Europe show instead a smaller SOC loss when erosion is taken into account. This is because we assumed a steady state in 1850, where carbon losses by erosion are equal to the carbon input by litter. And as soil erosion decreased during 1850–2005 in Western Europe – mainly due to a decreasing trend in precipitation since 1965, less intense expansion of agricultural lands and grasslands (Fig. 5b), and agricultural abandonment – it partly offsets the decrease of SOC by LUC (Fig. 8).

Figure 7 Cumulative SOC stock changes per PFT during 1850–2005 implied by variable land cover, climate, and CO₂, without erosion (grey colors) and with erosion (red colors).

The significantly smaller increase in SOC stocks on agricultural land when the best estimated soil erosion rates are taken into account (Fig. 7) explains the larger decrease in the global SOC stock during 1850–2005 (S1) compared to a world without soil erosion (Fig. 6a). Due to the slow response of the global SOC stocks to perturbations, this impact of soil erosion can be even larger at longer timescales. The effect of soil erosion on the SOC stocks is also influenced by the mechanism where removal of SOC causes a sink in soils that tend to return to equilibrium.

Table 6 Comparison of our model estimates of agricultural soil and SOC loss fluxes for the year 2005 with high-resolution model/observation estimates.

^* Quinton et al. (2010) included also pasture land in their study.

Furthermore, we find that the variability in the temporal trend of global SOC erosion is mainly determined by the variability in soil erosion rates and less by climate and rising atmospheric CO₂ that are affecting SOC stocks (Fig. 4). Also, the spatial variability in SOC erosion rates for the year 2005 and the spatial variability in the change in SOC erosion during 1850–2005 follow closely the spatial variability of soil erosion rates (Fig. 5b, d). This can be explained by the slow response of the SOC pools to changes in NPP and decomposition caused by CO₂ and climate in contrast to the fast response of soil erosion to changes in land cover and climate.

Table 7 Comparison of our model best estimates of agricultural soil erosion and SOC erosion rates for the year 2005 with best model/observation estimates from Doetterl et al. (2012) per continent. The uncertainty range for the present-day sediment and SOC fluxes is 18 %–62 % and 5 %–51 %, respectively. The uncertainty range in the values of Doetterl et al. (2012) is 30 %–70 %.

4.2 LUC versus precipitation and temperature change

Although the variability in the temporal trends of soil and SOC erosion is dominated by the variability in precipitation changes, the overall trend follows the increase in agricultural land and grassland. The global decrease in precipitation in many regions worldwide, especially in the Amazon, as simulated by ISIMIP2b, leads to a slight decrease in soil and SOC erosion rates (Fig. 4). At the same time precipitation is very variable and might not lead to a significant global net change in soil erosion rates over the total period 1850–2005. This result might be contradictory to the fact that major soil erosion events are caused by storms. But in this study we only simulate rill and interrill erosion, which are usually slow processes. In addition, previous studies (Lal, 2003; Montgomery, 2007; Van Oost et al., 2007) have shown that land use change is usually the main driver behind accelerated rates of these types of soil erosion. Our study confirms this observation.

Figure 8 (a) Difference between the changes in SOC stocks over the period 1850–2005 under all perturbations including soil erosion and the changes in SOC stocks excluding soil erosion, S1–S3; (b) difference between the changes in SOC stocks under LUC including soil erosion and the changes in SOC stocks excluding soil erosion, S1–S2 to S3–S4; (c) difference between the changes in SOC stocks under a variable climate and CO₂ increase including soil erosion and the changes in SOC stocks excluding soil erosion, S2–S4.

If we separate the effects of LUC and climate variability co-varying with soil erosion, we find that in the LUC erosion scenario with constant climate (see Sect. 2.5) the total global soil loss from erosion increases by a factor of 1.27 since 1850, while in the CC erosion scenario with constant LUC at the level of 1850 the soil loss flux from erosion decreases by a factor of 1.12 (Fig. 4). Analyzing the effects of LUC and climate variability separately on SOC erosion we find that in the LUC-only scenario (S2–S1) the total global SOC loss increases by a factor of 1.35 since 1850, while in the climate-change-only scenario (S2) SOC loss decreases by a factor of 1.12 (Fig. 4). This shows that LUC slightly dominates the trend in both soil and SOC erosion fluxes on the global scale during 1850–2005.

For soil erosion, however, LUC dominates the temporal trend less than for SOC erosion. This effect is especially clear for grasslands, where we find that climate variability offsets a large part of the increase in soil erosion rates by LUC, but not in the case of SOC erosion. This is the due to the fact that LUC has a much stronger effect on the carbon content in the soil than the effect of climate and CO₂ change on the timescale of the last 200 years. Also, intense soil erosion is typically found in mountainous areas where climate variability has significant impacts, while at the same time these regions are usually poor in SOC due to unfavorable environmental conditions for plant productivity.

Regionally, there are significant differences in the relative contributions of LUC versus climate variability to the total soil erosion flux (Figs. 2 and 5). In the tropics in South America, Africa, and Asia, where intense LUC (deforestation and expansion of agricultural areas) took place during 1850–2005, a clear increase in soil erosion rates is found even in areas with a significant decrease in precipitation due to a higher agricultural area being exposed to erosion. However, in regions where agriculture is already established and has a long history, precipitation changes seem to have more impact than LUC on soil erosion rates. A combination of our assumption that erosion rates are in steady state with carbon input to the soil at 1850 and minimal agricultural expansion during the last 200 years may be the reasons for this observation.

We also find that summing up the changes in soil erosion rates due to LUC alone and the changes in soil erosion due to climate variability alone do not exactly match the results in the changes in soil erosion obtained when LUC and climate variability are combined (Fig. 4). The nonlinear differences between soil erosion rates calculated with changing land cover fractions in combination with a constant climate (LUC) and soil erosion rates calculated by subtracting the erosion simulation CC from CC + LUC are significant for agricultural land but much smaller for other PFTs and at the global scale. It implies that the LUC effect on erosion depends on the background climate. This is important to keep in mind when evaluating the LUC effect on SOC stocks in the presence of soil erosion.

The decrease in global SOC stocks in simulation S3 is due to the various effects of LUC (without erosion) (Fig. 6a). During 1850–2005 LUC has led to a decrease in natural vegetation and an increase in agricultural land. At the global scale, the replacement of natural PFTs by crops results in increased SOC decomposition and decreased carbon input to the soil by litter fall due to harvest and a lower productivity. Regionally this effect of LUC may be different, depending for example on the natural PFTs that are replaced. Furthermore, the increase in carbon input into the soil after LUC due to increased litter fall when natural vegetation is removed may play a role, but this effect is only temporary. In addition, wood harvest after deforestation and crop harvest contribute to the decreased carbon input to the soils.

Figure 9 Contribution to the cumulative global SOC stock change over 1850–2005 by CO₂ fertilization (red), by the effect of precipitation and temperature change on the carbon cycle (dark blue), by the effect of precipitation change on soil erosion ( aqua), by the LUC effect on the carbon cycle (dark green), and by the LUC effect on soil erosion (light green).

We find that the global SOC stock decreases by 17 Pg C due to LUC only during 1850–2005 (Figs. 6a, S3–S4). The overall change in carbon over this period summed up over all biomass, litter, SOC, and wood-product pools due to LUC without erosion is a loss of 102 Pg C, which lies in the range of cumulative carbon emissions by LUC from estimates of previous studies (Houghton and Nassikas, 2017; Li et al., 2017; Piao et al., 2009). When we use our best estimated soil erosion rates in the SOC-dynamics scheme of the emulator we find that the LUC effect on the global SOC stock is amplified by 4 Pg C or a factor of 1.2 (S1–S2, Fig. 6a). The main reason behind this is the increase in soil erosion rates by expanding agricultural lands and grasslands that limit the increase in the global agricultural land and grassland SOC stock due to LUC (Fig. 7). This leads to a total change in the overall carbon stock on land due to an LUC of −106 Pg. Regionally the amplification of the LUC effect on SOC stocks by the increase in soil erosion ranges between a factor of 0.9 and 1.6 (Table 4).

Regionally, changes in precipitation can amplify or offset a large part of the increase in soil erosion due to LUC (Figs. 2, 5e). Globally we find that the decrease in global total precipitation, especially in the Amazon after AD 1960, partly offsets the increase in soil loss due to land use change (Fig. 4). It should be noted that the uncertainty in precipitation from global climate models for the Amazon is significant, making this result uncertain (Mehran et al., 2014). Furthermore, we find that precipitation and temperature changes lead to a small net decrease in SOC stocks at the global scale since 1950 (Fig. 9, S8). This is likely related to the decreased productivity under drought stress (Piao et al., 2009). However, soil erosion offsets this decrease by a small net increase of 2 Pg in SOC stocks, mostly due to the decreasing trend in precipitation globally after AD 1950.

4.3 Effects of atmospheric CO₂ increase

In the ORCHIDEE model, increasing CO₂ leads to a fertilization effect as it increases the NPP and results in an increase in biomass production on land for most PFTs, depending on the temperature and moisture conditions (Arneth et al., 2017; Piao et al., 2009). Figure 9 shows the contribution of this fertilization effect to the cumulative SOC stock change during 1850–2005 (S4–S8), which is in the same order of magnitude as the effect of LUC excluding soil erosion. Together with climate variability the atmospheric CO₂ increase offsets all the carbon losses by LUC in our model and leads even to a net cumulative sink of carbon on land over this period of about 30 Pg C (S3). This value is calculated by summing up the changes in all the biomass, litter, and SOC pools and is in line with other assessments that found a net carbon balance that is close to neutral over 1850–2005 (Arora et al., 2011; Ciais et al., 2013; Khatiwala et al., 2009).

In the presence of soil erosion, climate variability and the atmospheric CO₂ increase lead to a net cumulative sink of carbon over land of about 28 Pg C (S1), which is still within the uncertainty of assessed estimates (Arora et al., 2011; Ciais et al., 2013; Khatiwala et al., 2009). Soil erosion can thus slightly change the sink strength by influencing the net effect of LUC on the terrestrial carbon balance.

When the CO₂ fertilization effect is absent (S5, S7), we find that the temporal trend in the cumulative change in global SOC stocks is largely determined by the effect of LUC (Fig. 6b) and leads to a cumulative source of carbon on land of 76 Pg C. LUC alone leads to a cumulative decrease in SOC stocks of −14 Pg C (S7–S8), which is 3 Pg C less than the decrease in SOC stocks due to LUC in the presence of increasing atmospheric CO₂ concentrations (S3–S4). The overall change in carbon over 1850–2005 summed up over all biomass, litter, and SOC pools due to LUC alone is −74 Pg C in absence of increasing CO₂ (S7–S8), which is 27 Pg C less than the LUC effect on carbon stocks under variable atmospheric CO₂ (S3–S4). LUC has indeed a smaller effect on carbon stocks in the absence of increasing CO₂ concentrations as expected, because the productivity of the vegetation is lower (lower NPP), resulting in less biomass that can be removed by deforestation.

The previously calculated global total soil erosion flux of 47.6 Pg year⁻¹ leads to an annual SOC erosion flux of 0.48±0.13 Pg C year⁻¹ in the year 2005 in the absence of increasing atmospheric CO₂ (S5), which is about 0.04 Pg C year⁻¹less than the SOC erosion flux under increasing CO₂ (S1). The global cumulative SOC erosion over the entire time period in the absence of increasing atmospheric CO₂ is about 4.78 Pg C less (S5). Although these changes in SOC are small, the effect of LUC on the SOC stocks is amplified by erosion with a factor of 1.26 in absence of increasing CO₂ (S5–S6), which is slightly larger than the effect of LUC with increasing CO₂ (S1–S2). This means that the LUC effect in combination with soil erosion has a stronger effect on SOC stocks losses under constant atmospheric CO₂ conditions, because the CO₂ fertilization effect does not replenish SOC in agricultural lands everywhere.

Finally, it should be mentioned here that the absence of nutrients in the current version of the ORCHIDEE model may result in an overestimation of the CO₂ fertilization effect on NPP and may introduce biases in the effect of erosion on SOC stocks under increasing atmospheric CO₂ concentrations. Soil erosion may also lead to significant losses of nutrients in the real world, especially in agricultural areas. For a more complete quantification of the effects of soil erosion on the carbon cycle, nutrients have to be included in future studies.

4.4 Model limitations, uncertainties, and next steps

One of the uncertainties in our modeling approach is related to the application of the Adj.RUSLE model at the global scale and the estimation of the model parameters for various different environmental conditions and biomes. Although the Adj.RUSLE model was extensively validated using large high-resolution datasets, we calculated an uncertainty range for the R and C factors of the model to investigate the sensitivity of the emulator to the uncertainty in soil erosion rates. In Sect. 4.1 we show that including soil erosion in the emulator decreases the land carbon sink due to the large SOC losses on agricultural land triggered by erosion that reduce the SOC stocks significantly. Without soil erosion (S3) the global agricultural SOC stock increases by 60 Pg C due to agricultural expansion, while soil erosion reduces this increase by 11 Pg C in the minimum soil erosion scenario (S1 minimum) and by 18 Pg C in the maximum soil erosion scenario (S1 maximum). Thus, LUC results in a smaller increase in the global agricultural SOC stock under all soil erosion scenarios, while the magnitude of this effect is region dependent. The larger the soil erosion rates, the less carbon can be stored on agricultural land.

Furthermore, the aggregation of the high-resolution soil erosion rates from the Adj.RUSLE model to the resolution of the emulator can induce some uncertainties, as we might not capture correctly the hotspots of carbon erosion and their effects on the local SOC dynamics in these regions. However, the aggregation was needed to be consistent with the coarse resolution of ORCHIDEE and to limit the computational power of the emulator.

In addition, our soil erosion model is limited to water erosion only. This might result in biases for regions where other types of soil erosion are dominant such as tillage erosion (Van Oost et al., 2007), gully erosion, and landslides (Hilton et al., 2008, 2011; Valentin et al., 2005).

Although our erosion model runs on a daily time step, the soil erosion rates are calculated on a yearly time step, and thus we might miss extreme climate events triggering large soil losses. In addition, the Adj.RUSLE is not trained for extreme events. The effect of precipitation and temperature change on the SOC stocks under soil erosion might thus be larger than in our model simulations.

Concerning the reconstructed PFT maps, only expansion and abandonment of agriculture is taken into account, but not soil conservation measures as implemented in Australia and the US to prevent erosion (Chappell et al., 2012; Houghton et al., 1999). Regarding the land use change method that we applied, we only account for net land use change and do not account for shifting cultivation or distinguish between areas that have already seen LUC. Forest regrowth and forest age are also not considered, which could bring uncertainties in our estimates of LUC emissions (Yue et al., 2018). To show the potential uncertainty in our results due to uncertainties in underlying land use data we performed four additional simulations (S1 to S4) with a new PFT map using the same methods and data as by Peng et al. (2017), however, where the forest area is not constrained with historical data from Houghton (2003, 2008) and present-day data from satellite land cover products. In the following we will refer to the new PFT map as the unconstrained PFT map. In the unconstrained PFT map there is a stronger decrease in forest area over the period 1850–2005. Also, the grassland shows an increasing trend, while in the PFT map with constrained forest the grassland shows globally a slight decreasing trend (Fig. S2 of the Supplement).

After calculating soil erosion with the unconstrained PFT map we find that the differences in global average soil erosion rates between the different PFT maps are small (Fig. S3a). This can be related to the fact that the C factor of the Adj.RUSLE model is similar for forest and dense natural grass. As the change in global agricultural area is not significantly different between the two PFT maps, the overall soil erosion rates are similar. We expect, however, that the differences in soil erosion rates between the PFT maps can be larger in areas where the change in forest area is substantial over the historical period.

In contrast to the soil erosion rates, the two PFT maps result in significant differences in the SOC erosion rates and cumulative changes in SOC stocks during the transient period (Fig. S4). The global SOC stock in the equilibrium state without soil erosion (S3) is 8 % higher when the unconstrained PFT map is used, due to a larger global forest area in this map at 1850. The higher global SOC stock of the unconstrained PFT map leads to higher SOC erosion rates (Fig. S3b). According to the unconstrained PFT map, soil erosion leads to a total SOC removal of 79 Pg C (S1) over the period 1850–2005, which is 5 Pg C larger than the total SOC removal by soil erosion for the constrained PFT map.

Interestingly, due to the unconstrained PFT map, the global cumulative SOC stock change over 1850–2005 under soil erosion, climate change, and LUC (S1) is 60 % smaller than the stock derived using the constrained PFT map. This is most likely due to the higher forest area at the start of the period 1850–2005, leading to a larger increase in SOC stocks by increasing atmospheric CO₂ concentrations. We find the global LUC effect on the SOC stocks of both PFT maps to be similar (Fig. S3c, d).

5 Conclusions

In this study we introduced a 4-D modeling approach where we coupled soil erosion to the carbon cycle of ORCHIDEE and analyzed the potential effects of soil erosion, without sediment deposition or transport, on the global SOC stocks over the period 1850–2005. To calculate global potential soil erosion rates we used the Adj.RUSLE model that includes scaling approaches to calculate erosion at a coarse spatial and temporal resolution. The SOC dynamics are represented by an emulator that imitates the behavior of the carbon cycle of the ORCHIDEE LSM and enables us to easily couple our soil erosion model to the carbon cycle and calculate the effects of soil erosion under different climatic and land use conditions. Although our modeling approach is rather coarse and fairly simple, we found a fair agreement of our soil loss and SOC loss fluxes for the year 2005 with high-resolution estimates from other studies.

When applying the model on the time period AD 1850–2005 we found a total soil loss flux of 7183±1662, where soil erosion rates increased the most on agricultural land. This potential soil loss flux mobilized 74±18 Pg of SOC across all PFTs, which compares to 46 %–74 % of the total net flux of carbon lost as CO₂ to the atmosphere due to LUC estimated by our study for the same time period. When assuming that all this mobilized SOC is respired we find that the overall SOC change over the period 1850–2005 would increase by 62 % and reduce the land carbon sink by 2 Pg of carbon. The effect of soil erosion on the cumulative SOC change between AD 1850 and 2005 differs significantly between regions, where the largest decrease in SOC due to soil erosion is found in Asia. The expansion of agricultural land and grassland is the main driver behind the decreasing SOC stocks by soil erosion. Including soil erosion in the SOC dynamics amplifies the decrease in SOC stocks due to LUC by a factor of 1.2. Overall, the potential effects of soil erosion on the global SOC stocks show that soil erosion needs to be included in future assessments of the terrestrial carbon cycle and LUC.