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:

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:

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

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

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

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

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:

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:

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

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:

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:

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

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

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

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 2Spatial 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.

Download

2.4 Input datasets

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 1Description 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.

Download Print Version | Download XLSX

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 3Conceptual 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.

Download

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.

Download

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 2Area-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 %.

Download Print Version | Download XLSX

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 3Model 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.

Download Print Version | Download XLSX

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.

Download

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 4Model estimates per continent of changes in SOC stocks since 1851 from simulations S1, S2, S1–S2, S3, S4, and S3–S4.

Download Print Version | Download XLSX

Table 5Statistics 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.

Download Print Version | Download XLSX

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

Download

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

Download

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 6Comparison 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.

Download Print Version | Download XLSX

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 7Comparison 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 %.

Download Print Version | Download XLSX

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.

Download

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

Download

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