Soil erosion is unlikely to drive a future carbon sink in Europe

See allHide authors and affiliations

Science Advances  14 Nov 2018:
Vol. 4, no. 11, eaau3523
DOI: 10.1126/sciadv.aau3523


Understanding of the processes governing soil organic carbon turnover is confounded by the fact that C feedbacks driven by soil erosion have not yet been fully explored at large scale. However, in a changing climate, variation in rainfall erosivity (and hence soil erosion) may change the amount of C displacement, hence inducing feedbacks onto the land C cycle. Using a consistent biogeochemistry-erosion model framework to quantify the impact of future climate on the C cycle, we show that C input increases were offset by higher heterotrophic respiration under climate change. Taking into account all the additional feedbacks and C fluxes due to displacement by erosion, we estimated a net source of 0.92 to 10.1 Tg C year−1 from agricultural soils in the European Union to the atmosphere over the period 2016–2100. These ranges represented a weaker and stronger C source compared to a simulation without erosion (1.8 Tg C year−1), respectively, and were dependent on the erosion-driven C loss parameterization, which is still very uncertain. However, when setting a baseline with current erosion rates, the accelerated erosion scenario resulted in 35% more eroded C, but its feedback on the C cycle was marginal. Our results challenge the idea that higher erosion driven by climate will lead to a C sink in the near future.


The latest projections of soil organic carbon (SOC) changes from different Earth system models (ESMs) show a very wide range of values, adding uncertainty on land-atmospheric carbon (C) feedback under climate change (1). Because the soil is the largest terrestrial C pool, different representations of SOC turnover may amplify the magnitude and even the sign of the C balance (2). Recent studies suggest that the general soil C sink simulated by 2100 is too optimistic, as ESMs underestimate the mean age of SOC (3). Moreover, the level of confidence of current projections and the possibility of improving process descriptions of SOC dynamics may be further limited by the exclusion of relevant processes, such as soil erosion (4). This erosional process not only redistributes SOC across the landscape but also generates different feedbacks (5, 6). In eroding landscape positions, the C that is laterally displaced can be partially replaced with new photosynthates stabilized by the less C-saturated subsoil exposed by erosion, generating a continuous C sink (that is, dynamic replacement) (7). The displaced SOC can be mineralized during transport and preserved in depositional sites with different degrees of efficiency (8). Depending on the strength of these processes, which remain the topic of scientific debate, erosion has been claimed to induce both a net source and net sink of C (911). A recent study, for instance, estimates that erosion has offset one-third of C emissions from land cover change since the onset of agriculture (12). This global-scale analysis highlights how agricultural perturbation of more stable land use (grassland and forest) could generate a C sink over millennia (13). In these approaches, the net effect of soil erosion on soil-atmosphere C exchange is inferred from the measurement of SOC stock change against a reference soil profile, which may integrate the effect of other drivers in addition to the induced erosion effect on the C cycle (the so-called black box model). For instance, primary productivity and soil respiration can be affected by anthropogenic interventions in eroding/depositional areas, leading to a SOC change not entirely driven by the “erosion effect.” Therefore, inventory-based approaches are really useful for tracking past changes but are unlikely to be useful for predicting future change in which the impacts of projected climate on both vertical C fluxes (photosynthesis and respiration) and, indirectly, C redistribution by rainfall erosivity variation are important drivers. As erosion is a widespread phenomenon, especially on current agricultural land use (14), any approach oriented toward projecting future changes requires a framework that can dynamically assess variation of drivers at the same time, integrating current erosion rates as a baseline condition.

The erosion-induced C feedbacks have been studied at the scale of small hillslopes and watersheds (6), but the lack of large-scale data, in conjunction with gaps in process understanding, has prevented the coupling of erosion/distribution processes in large-scale modeling. For explaining past changes, several large-scale data-driven approaches (15) or integrated small-scale modeling (16, 17) studies have been published, but no study to date has attempted to project future changes, which is the subject of this paper.

Biogeochemistry-erosion model framework

Here, we use a consistent biogeochemistry-erosion model framework to quantify the impact of future climate on the C cycle, taking into account accelerated water erosion through variation in rainfall erosivity (18). The model couples the process-based biogeochemistry model CENTURY (19, 20) and the RUSLE2015 erosion model (21), the latter based on the Revised Universal Soil Loss Equation (RUSLE). The strength of this approach is that (i) any SOC change related to land use or management activity, including the effect of climate on plant productivity and soil respiration, propagates as a different amount (and type) of displaced SOC, potentially inducing a dynamic feedback, and (ii) land use and management directly affect the soil erosion rates, giving a bidirectional feedback on C balance.

This integrated model has already been used to quantify the current impact of erosion in the agricultural soil of the European Union (EU) at a resolution of 1 km2 (fig. S1) (22), showing a high level of accuracy in simulating dynamic replacement in comparison with previous studies and the overall SOC budget (section S2: figs. S2 to S5). Some of the previous assumptions on sediment distribution that are required for running the model at continental scale were further tested using a delivery/sedimentary model, as detailed in section S3 (figs. S6 and S7 and table S1).

The model was run for the Representative Concentration Pathway (RCP) 4.5 climate scenario as simulated by the HadGEM2 climate model, with current erosion (CE) and accelerated erosion (AE), the latter induced by the rainfall erosivity change as calculated recently by Panagos et al. (18). When the model incorporates erosion, the C balance accounts for C displacement components (see Materials and Methods) and does not generate CO2 fluxes until the C is mineralized. As we were interested in the net C exchange between soil and atmosphere, the net soil C flux was calculatedEmbedded Image(1)where CI and CH are the so-called vertical fluxes representing the atmospheric CO2 fixed in litter inputs and released by soil heterotrophic respiration, respectively, while CTm and CBm are the C mineralized during transport and in the burial deposits. There is limited consensus in the literature about C mineralization during transport and its preservation at depositional sites (23). To provide a range of model uncertainty, we ran the model with two configurations according to previous model sensitivity analysis (22): the first enhancing the erosion-induced C sink effect (e) and the second minimizing it (r) (see Materials and Methods). Last, we calculated the contribution of accelerated erosion (ΔCAE) on the net C exchange as the difference between the scenario with accelerated and current erosion (both under e and r uncertainty configurations)Embedded Image(2)


The effect of climate change on the terrestrial C cycle is predominantly explained by productivity-driven changes in ESMs (24), with potential significant biases when omitting the land management effect (25). Our model framework, which includes a detailed description of land use and agricultural management, also showed an increase in C input by 2100 (Fig. 1A), driven by the higher productivity under the HadGEM2 projection observable at mid-high latitude. This increase was almost offset by enhanced soil respiration, which showed a similar continental response to C input fluxes (Fig. 1B). By the end of the century, the model also predicted an overall increase of C displaced by erosion, with an opposite trend in western countries (Fig. 1C) and some large areas such as Andalusia and central France. This pattern partially overlapped with the change in rainfall erosivity (fig. S8), as eroded C is dependent not only on soil erosion but also on the evolution of SOC stock due to trends in productivity and soil respiration. These complex feedbacks emerged in the net soil flux figure (Fig. 1D), where no clear continental patterns of land C sink and source were predicted over the period 2016–2100.

Fig. 1 Response of C fluxes (Mg C km−2 year−1) under the RCP 4.5 (implemented in HadGEM2) climate scenario for the accelerated soil erosion simulation (AE_e).

(A to C) Average difference in C input, soil heterotrophic respiration, and eroded C, respectively, between the 2090–2100 and 2000–2010 time frames. (D) Cumulative change in net soil C flux over the period 2016–2100 (Mg C km−2) in EU agricultural soils. In (D), negative values represent a C source to the atmosphere, while positive values are a C sink.

Looking at the continental C budget in the AE_e scenario (Fig. 2), we estimated a small C sink from agricultural land equal to 0.09 Tg C year−1 (CI − CH). Despite the fact that eroded C was more than one order of magnitude lower than the vertical fluxes, its contribution was important; the additional losses generated by C displacement (CTm and CBm) offset the vertical C sink, leading to a C source from land equal to 0.92 Tg C year−1 over the period 2016–2100. When we ran the model without erosion, the C source from land was even higher and equal to 1.8 Tg C year−1 (table S2). This result suggests that climate change was the overarching driver on net CO2 losses, which were partially offset by soil erosion through dynamic replacement (fig. S5), under this conservative model configuration.

Fig. 2 Cumulative C budget (Tg C year−1) over the period 2016–2100 at the EU level in the accelerated (AE) and current (CE) soil erosion scenarios.

The numbers in brackets [e; r] are the outcomes of the two model configurations: enhanced erosion-induced C sink (e), with the mineralization during transport set to 2%, the burial efficiency to 95%, and the enrichment factor to 2, and reduced erosion-induced C sink (r), where the same parameters were set to 10%, 20%, and 1, respectively. Dark arrows represent C displacements, while blue arrows represent C fluxes as C exchanges with the atmosphere (CO2-C). For the net soil flux, negative values represent a C source to the atmosphere, while positive values represent a C sink (see Materials and Methods for the full C balance component details). The agricultural area simulated (arable crops, grassland, and permanent crops) covers 1.88 Mkm2.

In the AE_r scenario, the halving of the enrichment factor caused a direct reduction of eroded C, which amounted to 63% than the AE_e (Fig. 2). Despite the lower SOC loss from the soil profile by erosion, the concomitant lower dynamic replacement led to a positive feedback on the vertical fluxes (CI < CH with a value of −0.56 Tg C year−1). In this model configuration, we estimated a net soil C flux of −10.1 Tg C year−1, largely driven by the lower burial efficiency of buried SOC in the model setup (that is, 0.2). Although this value represents a quite extreme configuration of our model framework, a recent study indicated burial efficiencies of 0.22 and 0.14 in colluvium and alluvium stores of temperate zones, respectively (12). In this scenario, the C feedbacks generated by erosion resulted in a net soil C flux five times higher than that driven by climate change alone (−10.1 versus −1.8 Tg C year−1).

The C displacement between eroding and depositional areas also induced different C feedbacks mediated by the associated nitrogen (N) (that is, in the organic matter), in turn affecting crop productivity at different landscape positions (see fig. S9 and related discussion) (26). At the EU level, we calculated a CI loss of 0.05 Tg C year−1 per Tg of C eroded, which suggests relatively small reductions in productivity over the period 2016–2100, as supported by other studies (27). Conversely, the effect may be much higher locally and the magnitude of our estimate may be conservative, if the increase in rainfall erosivity will be associated with a higher frequency of extreme events (28). In CENTURY, the unique feedback on plant productivity due to erosion is related to the displacement of N (associated with SOC) and its turnover, while the model does not simulate mechanistically physicochemical soil changes (for example, in structure, texture, and loss of all nutrients) due to rain-driven intensive erosion events or the direct impact on the vegetation (29). Direct and indirect effects of erosion on plant productivity may then be very important (30). Despite the fact that the direct plant response to erosion is only modulated by N in our model, this is the first attempt to dynamically include lateral C fluxes in a process-based model framework at continental scale.

Soil erosion is a widespread process that currently affects EU agricultural soils (21); therefore, the evaluation of future C feedbacks driven by its change requires a baseline reflecting the current conditions. To quantify the potential C feedback of accelerated erosion in the future, we ran the model, keeping erosion rates constant at current levels (current erosion: CE scenarios). With this run, we estimated that the future variation in rainfall erosivity led to 35% more C displaced by erosion in both configurations (e and r; Fig. 2). Overall, the climate-induced erosion feedback (ΔCnetAE) was marginal in the e scenario but was more pronounced in r, with a loss of 2.34 Tg C year−1 (Fig. 2 and table S2). Considering the contrasting ranges of parameterization in the two configurations, it is likely that we are catching the magnitude of this feedback, although there are no data against which to compare the model.

Despite the fact that accelerated erosion induced a negligible C source continentally under the AE_e scenario (0.01 Tg C year−1), we observed a marked regional variability (Fig. 3A). In 35% of the cells (that is, 66 Mha), we estimated an induced C sink by the accelerated erosion, completely offset by sources for the remaining agricultural areas. In the AE_r scenario, the erosion-induced C sink was present in only 13% of the agricultural soils (Fig. 3B), while lower values than those in AE_e were estimated all over the continent.

Fig. 3 Climate-induced erosion feedback on net soil C flux (Mg C km−2 year−1).

The values represent the difference in net soil C flux (ΔCnetAE) between the accelerated erosion (AE) and current erosion (CE) scenarios (Eq. 2) over the period 2016–2100. (A and B) The e and r model configurations. Negative values indicate an erosion-induced source of C to the atmosphere driven by a variation in rainfall erosivity, while positive values show an erosion-induced C sink.

Projections of land C fluxes in current Earth system simulations under climate change are particularly affected by C turnover in soils (31) and are still missing the erosion feedbacks over decadal scales (32). This study may be a preliminary step toward coupling landscape dynamics in the biogeochemistry component of land schemes. Our results highlight how the interaction among soil variability, local management, climate, and geophysical processes leads to a variable C response, both in space and in sign of emissions. Upscaling of inventory-based approaches to large areas should therefore be done with care so that processes that have inherent variability are not oversimplified. Furthermore, our results challenge the suggestion that erosion will lead to a C sink in the near future, but they also claim for new researches investigating poorly known landscape C processes.


Our modeling framework has some limitations. Soil erosion by water is a selective process that redistributes different size particles across the landscape, in turn potentially changing the environmental conditions affecting plant growth and SOC decomposition (33). In addition, the displaced C is associated to different mineral/particle fractions (for example, C bound to mineral particles, C in aggregates, and particulate C) (34), and the turnover of these different fractions is not yet implemented into the large-scale models, which are still currently based on the compartment approach of “conceptual” C pools (35). We estimate that the development of full sediment delivery models coupled with new generations of biogeochemical models based on measurable C fractions will take another decade. To do this, we will need more information and increased process understanding to transfer microscale processes, such as aggregate breakdown, to millions of hectares. Meanwhile, we could account for lateral C fluxes with feasible approaches, as we have described here. Future land C projections rely on complex Earth system and biogeochemical models that do not incorporate the erosion process, lowering the level of confidence of our predictions. We have attempted to include the main lateral and vertical C fluxes as driven by a climate change projection to increase the confidence of future projections of soil C erosion on the C cycle.


Model framework and configuration

The model framework was based on coupling the process-based biogeochemistry model CENTURY to the RUSLE2015 erosion model. This model integration has been presented in previous studies, and a summary flowchart is given as fig. S1. We summarized the model setup and main assumptions as follows.

1) The CENTURY model was ran at a resolution of 1 km2 for the agricultural soil of the EU, using the soil erosion from RUSLE2015 model as input for CENTURY. Starting from 1900, the erosion process was implemented, keeping the climate, soil, and topographic factors (R, K, and LS, respectively) constant. While we considered K and LS factors quite invariable on a centennial scale, the C factor associated with the crop type was dynamically varied with crop rotations and land use changes. The simulated land use was based on the CORINE Land Cover 1990, 2000, and 2006, supplemented with Eurostat (Statistical Office of the European Communities)/FAO (Food and Agriculture Organization of the United Nations) statistics to build up crop rotations and implement consistent agronomic inputs (fertilization, irrigation, etc.). Before 1990, we assumed the same land use but with different agricultural techniques characterized by lower productivity crops, lower rates of mineral N, and different rotation schemes.

2) Originally, we assumed that each 1-km2 grid cell is composed of an eroding area and a depositional area, the latter retaining a proportion of eroded C. The partition was based on the study of Van Oost et al. (11), who found that 53 to 95% of eroded SOC were retained in the catchment and redeposited in a limited area (14 to 35%) within the same catchment. Taking the central values, we assumed that 25% of our grid cells were depositional areas, which received 70% of eroded soil. The remaining 30% was accounted for as leaving the grid cell as potential sediments and C discharged to riverine systems. These assumptions on sediment distributions, necessary to work at continental scale, were further tested using a delivery/sedimentary model in regional simulations, as detailed in section S3. After the intercomparison between the original and the sedimentary model–driven configuration, we set a new sediment delivery ratio of 0.11, as most of the sediments were predicted to be retained in land under the delivery/sedimentary model runs.

3) The replacement of eroded soil in the fixed profile comes by “recruitment” from the subsoil layers (SSLs), characterized by a SOC composition defined quantitatively and qualitatively as a partition among the three CENTURY SOC pools (active, slow, and passive). These pools are functionally defined on the basis of mean residence time, as opposed to measured fractions, and thus cannot be constrained by measurements applicable at the EU scale. The most rational approach was to adopt the calibration of Harden et al. (7), which implicitly related the subsoil composition to the topsoil composition; accordingly, the subsoil SOC pools at time t = 0 (that is, 1900) were assumed to be composed as a fraction of the top soil pools, as followsEmbedded Image

The lower horizon pools were decremented at each time step using the following relationships for each of the three SOC pools (i)Embedded Imagewhere t is the time and FLOST is the fraction of topsoil lost to erosion [for the sake of clarity, SSLi(t − 1) × FLOST(t) is the C moving up from subsoil annually].

Although CENTURY does not explicitly simulate depositional processes, they can be mimicked while ensuring conservation of mass (soil and C) within the system boundaries. First, we calculated the amount of soil transported from the eroding to the depositional fraction of each grid cell. Second, we simulated the burial effect of the deposition event as a “negative” erosion event, which moves a SOC fraction below the fixed soil profile (C burial flux), proportional to the amount of soil deposited on the surface. Last, the associated amount of deposited C was controlled by setting the incoming SOC pool in the model. This procedure required an iterative process, running the model year by year, first in the eroding area and then in the depositional area of each pixel. A check of the C balance closure was done at the end of simulations.

Future rainfall erosivity and erosion

The methodology for the estimation of the future rainfall erosivity was presented in a very recent study (18). Summarizing, this made use of the REDES (Rainfall Erosivity Database at European Scale) and a statistical approach (Gaussian process regression), used to spatially interpolate rainfall erosivity data with climatic scenarios. Using the HadGEM2 general circulation model climate projections (for RCP 4.5), available at the WorldClim repository, the rainfall erosivity was calculated for two time periods, namely, 2050 and 2070. Therefore, the soil erosion was recalculated, keeping the current factors of RUSLE, except the rainfall erosivity. The ratio between the future and the current soil erosion by water is reported in fig. S8.

Climate change scenarios with and without accelerated erosion

The model was run with the current climate up to 2015 and with the same HadGEM2-ES_rcp45 scenario for the 2016–2100 time frame for consistency. Monthly precipitation and maximum and minimum air temperature were downloaded from the World Climate Research Programme–Coordinated Regional Downscaling Experiment (WCR-CORDEX) portal ( The present land use and agricultural management remained the same to isolate the climate effect on both productivity/turnover and erosion C feedbacks.

In the accelerated erosion (AE) scenario, we linearly interpolated the soil erosion rates since 2016, considering the two soil erosion projections at 2050 and 2070. The name “accelerated erosion” refers to the average effect of rainfall erosivity change on increasing soil erosion. Higher and lower rainfall erosivity patterns are present at the EU level, as depicted in fig. S8. In the current erosion (CE) scenario, we kept the actual soil erosion rates.

Given the fact that there is limited consensus on C preservation/decomposition upon its displacement, we ran the model with two configurations, enhancing (e model setup) and minimizing (r model setup) the erosion as an induced C sink (or source). In particular, in the e configuration, the enrichment factor (that is, the C concentration in eroded sediments with respect to the bulk soils) was set to 2, the C mineralization during sediment transport to 2% yearly, and the burial efficiency (that is, the amount of C preserved when moved downward into the soil profile) to 95%. In the r configuration, the three parameters were set to 1, 10, and 20%, respectively.

According to previous sensitivity analyses, a higher enrichment factor promotes a higher dynamic replacement of C (inducing a C sink in eroding position), while lower mineralization and higher burial efficiency preserve the displaced C. Considering the ranges of parameters from different studies, our configuration should reflect the outer boundary of the erosion effect in inducing a sink/source of C.

C budget

The CENTURY model was parameterized to simulate the soil C and N turnover in a fixed profile of 0 to 30 cm. When the model incorporates soil erosion (22), an upward C flux from a subsoil pool follows the erosion event (as described in point 3 of the “Model framework and configuration” section). Moreover, the C balance accounts for C displacement components, some of them representing C moving outside the land domain (such as C to rivers, lakes, and ocean). For the purpose of the study, we set our boundary conditions to the land and, despite estimating a lump term accounting for C delivered to the riverine system, we did not calculate its decomposition (to CO2) related to the turnover outside the land.

For the eroding area, the C balance isEmbedded Image(3)where ΔSOC is the SOC stock change in the fixed profile, CI is the C input through the remaining net primary productivity (NPP; after C exportation by harvest but including roots and manure), CH is the heterotrophic respiration, CS is the incoming SOC from a deeper layer, CE is the lateral C flux by sediment transport, and DOC is the C exported as dissolved organic C.

For the depositional area, the balance isEmbedded Image(4)where CD is the C deposited coming from eroding areas and CB is the C that is moved out (that is, buried) of the simulated profile as a consequence of soil deposition.

However, the net changes in SOC storage do not directly equate to the net CO2 exchange from soil to the atmosphere due to the combination of vertical C fluxes (CI − CH) of eroding and depositional areas and the mineralized fraction of SOC upon displacement. Therefore, the net soil C flux (as CO2) in the grid cell is given byEmbedded Image(5)where [CI − CH] is the net vertical C flux, CTm represents the fraction of eroded SOC mineralized during the transport, and CBm is the fraction of buried C that is mineralized in depositional area. Because the uncertainty in simulating C turnover in deeper layers increases with depth (due to cumulative sediment deposition), the cumulative CB fluxes, over the period 2016–2100, were assumed to be preserved from decomposition with two contrasting burial efficiency rates over a 100-year horizon (0.95 and 0.2). The corresponding CO2 flux to the atmosphere was then calculated as CB × (1 − burial efficiency).


Supplementary material for this article is available at

Section S1. Model framework

Section S2. Model validation

Section S3. Simplified versus delivery sediment soil budget

Fig. S1. Modeling framework of the integration between the CENTURY biogeochemistry model application and the RUSLE2015 erosion model.

Fig. S2. Distribution of LUCAS sampling location for all land uses.

Fig. S3. Comparison of cumulative SOC stock distribution between simulated and LUCAS values grouped at the national level.

Fig. S4. Comparison of cumulative SOC stock distribution between simulated and LUCAS values grouped for environmental factors.

Fig. S5. Comparison of modeled erosion-induced vertical C exchange against inventory-based data.

Fig. S6. Selected regions where the model configuration intercomparison was run.

Fig. S7. Simulated C balance components of CENTURY driven by the WATEM model in comparison with the standard original and corrected sedimentary parameterization.

Fig. S8. Change in soil erosion by 2050 and 2070 due to projected changes in rainfall erosivity.

Fig. S9. Distribution of change in NPP between depositional and eroding areas for the different scenarios simulated.

Table S1. Comparison of C budget between different couplings of CENTURY with erosion/delivery sediment parameterizations (Mg C km−1 year−1) in the selected regions.

Table S2. Cumulative C budget (Tg C year−1) over the period 2016–2100 at the EU level for scenarios with accelerated (AE) and current (CE) soil erosion.

References (3644)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank I. Biavetti for graphical support in designing Fig. 2. Funding: The work was carried out as part of the JRC’s Institutional Work Programme under the Natural Capital Soil Project (Project 702), Work Package 5037 “Soil for Climate Change.” Author contributions: E.L. developed the research concepts and conducted the modeling and analyses. E.L., P.S., and P.B. performed the data interpretation. E.L. took a lead on writing the paper with contributions from P.S. All authors assisted in the revision of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are available at the European Soil Data Centre (ESDAC) of the European Commission Joint Research Centre ( Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article