Evaluating the interaction between sediment fluxes , carbon dynamics and biomass production using an integrated model

During the last centuries, forest clearance has led to an increase of the erosion rates by one to two orders of magnitude. Sustained accelerated soil erosion alters key soil properties such as nutrient, water availability, soil depth and soil 10 texture, which in turns have detrimental effects on crop yields and therefore reduce C input to soils. In this study, we applied a 1D model dynamically linking soil organic carbon (SOC) turnover, soil erosion and crop yield at the profile scale. We extracted a relationship linking crop yield to soil erosion based on available literature and categorized them into three functional forms: high sensitivity to erosion, linear response and low sensitivity to erosion. We tested and validated the model using published observational data from 12 catchments across Europe and the USA. Model evaluation showed that accounting for 15 the erosion-crop yield feedback (i) increased SOC losses by 20% on average and (ii) improved the SOC losses predictions, particularly for higher cumulative soil erosion, compared to the results obtained without the feedback. Cumulative vertical carbon fluxes were reduced by 15 to 71% compared to the no-feedback model, although the large variability highlighted the need to perform site-specific adjustments of the erosion-crop yield relationship. Exploration of parameter sensitivity to SOC parameters and erosion showed that long-term simulations of both SOC loss and vertical C fluxes were primarily influenced 20 by the erosion rate, the yield response to erosion and the depth distribution of the mineralization rate of organic matter. Our simulations further highlighted the increased SOC losses (-3 to -17%) and reduced C uptake from the atmosphere (-30%) in the erosion-crop yield feedback scenario, compared to the no-feedback scenario after 200 years, as well as the importance of the functional form of the erosion-crop yield relationship. Together, this modeling study shows that including the effects of erosion on crop yields has a large potential to reduce uncertainties associated with the estimation of the C budget in landscapes 25 subjected to erosion.


Introduction
The soil system represents one of the most important carbon (C) pools by storing around 1417 PgC in the upper first meter.
As a result, its impact on the global C cycle and climate has been widely recognized and studied (Hiederer and Köchy, 2011;Houghton, 2007;Crowther et al., 2016).The terrestrial carbon cycle is mainly driven by soil-atmosphere exchanges; vegetation takes up carbon from the atmosphere and provides input into the soil in forms of root excretions and plant residues while biologic activity and in-situ mineralization release carbon back to the atmosphere from soils (Houghton, 2007).
Through the removal of natural vegetation for agricultural extension, human activities have an important impact on the soil system, not only by changing the soil C cycle, but also by increasing soil erosion rates by up two orders of magnitude (Vanacker et al., 2013).(Gregorich et al., 1998;Montgomery, 2007).Soil erosion affects vegetation growth and biomass production by changing properties related to soil fertility such as water holding capacity, nutrient status or soil depth (Kirkels et al., 2014;Bakker et al., 2004).Effects of soil erosion on crop productivity have been studied during the past decades for a wide range of pedological and climatic conditions (Kosmas et al., 2001;Bakker et al., 2004;Fenton et al., 2005;Gregorich et al., 1998).Despite highly variable results, these studies have indicated that in absence of fertilizers, yields tend to decrease when soil is subject to erosion (Bakker et al., 2004;den Biggelaar et al., 2003;Larney et al., 2016).Hence, on the long term, reduced biomass production will result in an additional loss of SOC due to decreasing soil C inputs (Gregorich et al., 1998;Doetterl et al., 2016;Kirkels et al., 2014).Although large uncertainties still remain concerning the strength and the mathematical form of the relationship between crop growth and soil erosion general tendencies have now been identified through meta-data analysis (Bakker et al., 2004;Chappell et al., 2012).
In addition to changes in soil C inputs, human-induced erosion also results in the subsequent intensification of lateral SOC losses through erosion and lateral redistribution across the landscape (Van Oost et al., 2005a).Three processes mainly govern SOC dynamics under the impact of erosion: enhanced mineralization during transport, the replacement of eroded C by new photosyntate at eroding sites, and the burial and preservation in depositional areas (Stallard, 1998;Harden et al., 1999;Van Oost et al., 2007b;Lal, 2003;Hoffmann et al., 2009;Wang et al., 2014a).Although each of the described processes is individually relatively well understood, the result of their interaction at the landscape scale is poorly constrained (Kirkels et al., 2014).To our knowledge, the negative feedbacks on soil C inputs in response to erosion, and how it changes over time, have not received much attention when studying erosion-carbon cycling interactions (Harden et al., 1999).Therefore, a dynamic representation of the interactions between soil erosion, crop growth and SOC content is needed in order to better constrain the overall C balance (Chappell et al., 2015;Harden et al., 1999).
Mechanistic SOC models have been developed following two different approaches: the first approach explicitly models vegetation dynamics in relation to soil and climate conditions, but has a limited description of the soil component and its dynamics (Kaplan et al., 2012;Doetterl et al., 2016).This approach is typically useful at larger spatial scales and it represents well the effects of land use and vegetation change on soil properties, including SOC stocks (Kaplan et al., 2010;Kaplan et al., 2012).The second approach couples erosion and SOC models where (i) landscape and vegetation descriptions are limited and the effect of land cover on SOC dynamics is often reduced to a quantification of the total C input to soils (Doetterl et al., 2016;Minasny et al., 2015;Van Oost et al., 2005a;Wang et al., 2014b).Although this approach is considered to enhance our understanding of eroding landscape dynamics, few efforts have been made to link these two approaches to represent the interaction between soil redistribution, SOC dynamics and the feedback of erosion on biomass production (Kirkels et al., 2014;Doetterl et al., 2016).
This study proposes a step in further that direction by explicitly and dynamically linking crop yields, soil properties and SOC dynamics in a process-based soil profile model to explore the effect of soil erosion on SOC stocks and fluxes.The model accounts for SOC displacement and C input into the soil at the profile scale.Specifically, our objectives are (i) to develop an integrated process-based model linking SOC dynamics, crop yields and soil redistribution (ii) to evaluate the performance of the model by confronting model simulations to available observational data, and (iii) to investigate the longer-term (ie centennial) effect of erosion on crop growth and SOC dynamics at the profile scale.

Erosion effect on crop yields: Metadata analysis
This study used the data compiled by Bakker et al. (2004) to constrain the relationship between erosion intensity and crop yields as this is one of the most comprehensive meta-analysis.Only paired-plot experiment were selected as Bakker et al. (2004) pointed out that this method is the most realistic at estimating the effect of erosion on crop yields.
Following the analysis of Bakker et al. (2004), three functional forms of erosioncrop yield relationships are possible (Fig. 1) : a rapid and nonlinear decrease of crop yield as a function of soil truncation, a linear decrease corresponding to the influence of soil depth, and a slow and nonlinear decrease due to reduction of water availability (Bakker et al., 2004).
We explored the full range of constraints of soil truncation on crop yields using the following equation: where Ydr is the relative yield (compared to a reference yield of 1 for no-erosion), Tr is the cumulative soil truncation since the start of cultivation (m), α is the maximum yield reduction and B is the power law exponent linking the relative crop yield to soil erosion.
In our simulations, only the power-law exponent B varies, allowing to fully consider the wide range of relationships that are reported in literature.The concave form (with B < 0.9) can be related to the nutrient depletion case whereby the removal of nutrient-rich topsoil layers by erosion quickly affects crop growth (Bakker et al., 2004).As the depth distribution of nutrients typically follow an exponential evolution with depth, the constraints on crop growth become less important for more intense soil truncation (Bakker et al., 2004).The linear form (with 0.9 ≤ B ≤ 1.1) is related to soil depth limiting crop yields where soil depth limits the space available for root growth.In this case, physical hindrance decreases crop yields as soon as the root growth is limited by a compacted soil layer.The convex form is related to water availability, where crops that are typically not very sensitive to water limitation are not affected by soil truncation until a certain threshold beyond which crop yields are reduced quickly (Bakker et al., 2004).Based on an analysis of the Bakker et al. (2004) data, alpha was set to 0.7 (Fig. 1).
Finally, soil depth is indirectly taken into account as the input clay content profile is given in absolute percentage of volume (ie % of clay in a given volume of soil including coarse fragment).(e.g.Van Oost et al, 2005).ICBM is a two-pool carbon model simulating SOC transfer from the roots, residue and manure to a 'young' C pool, transfer from the 'young' pool to an 'old' C pool and C mineralization in both pools (Andren and Katterer, 1997).
SOC fluxes are described by the following equations: Where Y (Mg C ha -1 ) and O (Mg C ha -1 ) are respectively the young and old SOC pools and ky (yr -1 ) and ko (yr -1 ) their turnover rates (Andren and Katterer, 1997).i stands for the total carbon input which is the sum of the input from the crops (ic) and manure (im).The transfer from the young pool to the old pool is proportional to the humification factor (h) and the climatic and edaphic conditions which are condensed in the r coefficient (Andren and Katterer, 1997).
The humification factor is estimated as follows: with ic(z) and im(z) the C input from crop and manure at the depth z, hc and hm the humification coefficient for respectively crops and manure, and cl(z) the clay content at depth z (%).
The climate factor r, is corrected for the local climate using a Q10 relationship based on temperature (Andren and Katterer, 1997).
where T is the mean annual temperature (°C).
The model is depth-explicit as it considers a depth-dependent C input and mineralization rate (Nadeu et al., 2015;Van Oost et al., 2005b;Wang et al., 2014a).While manure and residue input only affect the topsoil layers, the carbon input from plant roots is distributed throughout the soil profile using the following relationship: With () the relative root-derived C input at depth z, z the soil depth (m),   is the depth of the top soil where ploughing is assumed to homogenize the SOC content and δ is the root density coefficient.
The turnover rates of the SOC pools at each depth are computed as an exponential function: With ur a dimensionless coefficient of depth attenuation,  0 (yr -1 ) the turnover rate at the soil surface.and  (yr -1 ) represents the SOC turnover rate at depth z.
The model starts with a prescribed SOC and soil profile having an initially defined clay content distribution.
Where i(t) is the C input at the time t, i(0) the initial C input and Ydr the relative yield at time t compared to initial yield.The implications of this assumption are discussed further.

Model implementation
The initial soil profile has a thickness of 1 meter and is represented by 100 layers of 1 cm, each layer being characterized by its own clay content, SOC input and turnover rate.As we assume a constant bulk density, the amount of clay and SOC vertically transferred between layers is proportional to the amount of erosion (upward transfer) (Van Oost et al., 2005a).The SOC content in the profile is then updated each year as a response to the vertical advection of matter, the C input at the soil surface and the clay content following erosion or deposition.The model keeps track of the SOC and clay content per layer, and tracks the evolution of the crop yield over time.
After performing a model spin-up without erosion for the C pools to reach equilibrium, we performed transient simulations where the soil profile is modified by erosion.During the simulation period, erosion rates are assumed to be constant through time.We presented the results in terms of the total SOC content evolution for the 1m profile and the net vertical C balance with the atmosphere.The annual vertical balance, i.e. the net vertical exchange of C between the soil and the atmosphere, of a soil layer at depth z and time t is calculated as follows: Where Cv(z,t) is the amount of carbon exchanged between the soil layer at a depth z and the atmosphere at time t.Positive vertical carbon fluxes denote C fluxes from the atmosphere to the soil while negative vertical carbon fluxes represent a C emission to the atmosphere.We evaluated the cumulative balance by integrating the vertical carbon fluxes over the entire duration of the simulation.

Model parametrization and calibration
In this study, we use published estimates of lateral SOC losses and net vertical C exchange that were derived from 137 Cs as a tracer for erosion while the C balance was derived from space for time substitutions (Van Oost et al, 2007).We used data on erosion rates and duration of the erosional disturbance for ten study sites in Europe and the US to estimate SOC loss and vertical C exchange.Eight sites were located in Europe and 2 sites in the US.The European sites represent a diversity of soil and/or climate conditions.Belgian, English and Danish sites are located in temperate climates and vary from each other by their climate, erosion rate and soil properties: from loamy soils with relatively high erosion rates in Belgium toward more clayloam soils and slightly lower erosion rates in for the English and Danish sites (Table 1b) (Quine and Zhang, 2002;Heckrath et al., 2005)..The two Belgian sites are sampled in nearby catchments and differ from each other by their topographic position: Belgium 1 being on steeper slope with higher erosion rate than Belgium 2. American sites were sampled in catchments in Iowa and are characterized with fine-textured loamy to silty soils and temperate continental climate (Ritchie et al., 2007).Due to their deep soils and humid climate, all of the aforementioned sites are therefore more prone to experience nutrient depletion constraint (B < 0.9) than soil thinning or water availability shortage.Mediterranean sites were characterized by warm and drier climate, clayey soils, high erosion rate (except for the Greeke case) and similar cultivation periods (Table 1).The Spanish sites differ by their topographic characteristics with steeper slope, shallower soils and higher erosion rate in the first site (Van Oost et al, 2007).Greek, Spanish and Portuguese sites experienced intense soil thinning (Greece and Spain, 0.9 < B < 1.1) or a mix of soil thinning and water availability constraints (Portugal) (Kosmas et al, 2001;Van Oost et al, 2007).For the Spanish and Portuguese sites, reference SOC data (ie representing stable profiles) were taken from the national surveys while for the other sites local data was used (Van Oost et al, 2007).
We estimated the initial SOC profile using site-specific observations and the reported mean annual temperature and clay content for each site.Based on this information, the following parameters were optimized for each site: C inputs at the surface (ic, Eq 4), the root-derived C input at depth z ((), Eq 6) and the depth attenuation of C mineralization (ur, Eq 7).We optimized the values by minimizing the relative root mean square error (RRMSE) between observed and simulated SOC profile (Eq 10).
Where N is the layer number, Ci,s is the simulated carbon content of the layer i (%) and Ci,o is the carbon content observed at the depth of the mid-point layer i. N varies for each site due to data sampling (Fig 3).reference erosion rate and +/-2% around the published clay content (Table 1, Van Oost et al, 2007).We produced two sets of 1000 simulations: with and without erosion-yield feedbacks (respectively designated by the FB and CTL abbreviations) to evaluate the effect of the erosion-crop yield relationship.SOC losses were calculated for the same periods as reported in the empirical study (Table 1).

Model validation 5
We performed a model validation using empirical observations on SOC losses and cumulative vertical C fluxes (Table 1) (Van Oost et al., 2007).In a first step, we calculated the observed SOC losses and cumulative vertical C fluxes for each sites based on the data of Van Oost et al (2007).The carbon balances were derived from soil erosion measurements using 137 Cs as a tracer.
The carbon budgets therefore integrate over the period 1954-1996 (the date of sampling).We computed the observed vertical C fluxes by summing the annual rates provided by Van Oost et al (2007) over the period of interest while the relative SOC losses were derived from the annual lateral C fluxes and the initial C profile.In a second step, we ran site-specific simulations using the 1000-parameters sets obtained from the calibration (see section 2.2.5 and Table 1).Finally, we confronted the results to the observations and evaluated performances of our model for both CTL and FB scenarios.We further quantified the effect of adding the feedback in the model by comparing the performances (i.e.RRMSE) of each scenario.

Long-term experiment
After the model evaluation, we further explored the behavior of the model by running additional simulations where we focused on the potential impact of including the yield reduction feedback in landscape modelling on longer timescales (i.e.c. 200 years).We built another 1000 model parameters sets in which the distribution of root depth and carbon mineralization rate, the initial clay content (cl), erosion rate (E) and yield response to erosion (exponent B) vary in a large but realistic range (Table 2).In order to explore the effect of clay content, we linearly scale the initial clay content profile (cl) creating an effective range of clay content from 5% to 45% in the top soil (Table 2).The values for all these parameters were randomly generated assuming a uniform distribution.We used the Fourier Amplitude Sensitivity Test (FAST) to assess the contribution of each individual parameter to the global variance of the results (Cukier et al., 1973;Cukier et al., 1975).We performed the FAST analysis with all the parameters varying between scenarios with the observed erosion rates (derived from Table 1a) and a uniform erosion rate of 1 mm.yr -1 in Table 2b.
Finally, using the same set of 1000 randomly-generated parameters, we evaluated the impact of the erosion-yield feedback on the SOC content and vertical balance by comparing the results produced by the model in FB configuration (erosion-yield feedback) and in CTL configuration (no erosion-yield feedback).
Table 2: Selected parameters, range of tested values and results of the FAST analysis.The FAST analysis results can be interpreted as the relative contribution of each parameter variability to the total variance of the selected output, i.e.ie the relative SOC loss compared to the initial SOC content and the cumulative vertical C fluxes at the end of the 200 years transient simulations.The asterisk designate the parameter accounting of more than 15% of the total variance.

Model evaluation 5
In this section, we perform a validation of the model by comparing our calibrated-simulation results with available data from the literature.Observations in Table 1 are taken from Van Oost et al. (2007a).The range of parameters used for calibration is displayed in Table 1 and Fig 3 shows the comparison between the optimized SOC profiles and the observations.All estimated initial SOC profile are close to the observations for each of the sites, with a RRMSE ranging between 0.01 %C to 0.09%C (Fig. 3).The observed SOC losses relative to the initial SOC content varied between 0.09 +/-0.02(UK) and 0.41 +/-0.08 10 (USA), with typical values around 0.15 +/-0.05 of the initial content.Figure 4a  Inclusion of the erosion-yield feedback (FB) increased SOC losses by 20% on average but showed a strong spatial variability related to the feedback nature.Simulations for sites characterized by soil thinning or water availability limitation did not result in substantial SOC losses (Table 1): for example, the Greek and Spanish sites exhibited only around 1% increase in SOC losses relative to the control scenario (Fig. 4a).On the contrary, sites where nutrient depletion (B < 0.9) was the dominant factor for yield reduction clearly showed an increase of both mean SOC loss and of the associated standard deviation compared to the CTL scenario.
Adding the feedback while using the same parameters improved the overall accuracy with an RRMSE of 0.63 for FB compared to 0.93 for CTL when all sites are taken into account (Table 1).The model performances are highly site-dependent: the addition of the feedback increased the prediction accuracy for the cases where cumulative soil truncation was substantial (e.g.. Belgium 1, Portugal, USA) while it decreased the accuracy for cases with small soil truncation (Fig 4a , Table 1).The FB scenario was able to reproduce the observed trend in which higher cumulative soil truncation leads to higher SOC losses.This shows that the addition of the feedback is important, particularly when soil truncation has been intense.However, it should be noted that the margins of error of both CTL and feedback modeled SOC losses are overlapping, except for the American sites.
The observed values for the vertical C fluxes ranged from 0.030 kgC.m -2 to 0.279 kgC.m -2 (Fig. 4b).The simulations without feedback produced results, which are of the same order of magnitude as the results from Van Oost et al ( 2007) with an RRMSE of 0.826 kgC.m-².Despite an accurate prediction of the observed trend, simulations tended to underestimate vertical carbon fluxes, particularly for the higher cumulative soil truncation cases (Fig 4b , Table 1).Including the feedback reduced vertical carbon fluxes by 34% in average relative to the CTL scenario and reduced the accuracy of the model with a RRMSE of 1.026 kgC.m-² (Fig 4b , Table 1)

Long-term simulations: Sensitivity analysis
Moving from the model validation to the long term modelling, the following sections present the results of the 200 years simulations using the wider range of randomly generated parameters (see model implementation).
We performed a model sensitivity analysis to explore the model behavior.The results of the FAST procedure are reported in Table 2a and 2b.Changes in SOC stock were almost entirely controlled by (i) the soil erosion rate (70% of the total variance) and (ii) the functional form of the erosion feedback on yields (22% of the total variance) (Table 2a).Similar observations are valid for the cumulative vertical balance, although it is more sensitive to yield reduction than erosion rate.Mineralization distribution with depth was the third major factor influencing SOC losses and the cumulative C balance, accounting for ~ 10% of the variability.It should be noted that clay content and root depth distribution played only a minor role When variability in erosion rate was excluded from the analysis, both SOC loss and the vertical carbon fluxes were mainly sensitive to the functional form of the feedback between erosion and yield (70% of the variance) and the mineralization rate distribution with depth (21%) (Table 2b).The root depth distribution had an effect, although much weaker, on SOC loss only (13% of the variance).The model simulations showed a strong positive correlation between SOC loss and erosion rate as well as with the functional form of the yield reduction (color scale, concave if b<1, linear when b ~1 and convex when b > 1) (Fig. 5a).The variability of the relative SOC loss to yield reduction was substantially larger for concave than for convex relationships.Although the same observations can be made for the vertical carbon fluxes, as the FAST analysis indicated, the influence of erosion was less important than the shape of the yield decrease.

Long-term SOC stock loss
Simulated SOC losses after 200 years ranged from 0.02 to 0.77 of the initial stock, depending on the erosion rate and the feedback type.For the FB scenario, the average SOC loss for the 1000 simulations equaled 0.38 with a standard deviation of 0.18 (Fig 5 and Table 3).In cases of erosion rates lower than 0.5 mm.yr -1 , the simulated SOC loss is limited to 0.20 of the initial content but increases almost linearly up to 0.2 to 0.5 at an erosion rate of 1.5 mm.yr -1 .Higher soil truncation rates exhibit a smaller variability in SOC losses with respective range of relative SOC loss of 0.32 to 0.60 and 0.40 to 0.70, for an erosion rate of 2 mm.yr -1 and 3 mm.yr - .
For the CTL scenario (Fig 5b), SOC losses ranged from 0.02 up to 0.57 of the initial carbon content depending on the feedback strength and the soil truncation, with a mean loss of 0.31 and a standard deviation of 0.14 (Fig 5b and Table 3).When including feedback in our simulations (FB scenario), SOC losses increased, particularly for the high sensitivity cases (i.e.nutrient depletion and soil thinning).Comparing the CTL with the FB scenario, FB further enhanced SOC loss by an additional 3% to  3).Furthermore, adding the feedbacks created more variability as indicated by the increase of the standard deviations shows.
When considering a period of 200 years, nutrient depletion (B<0.9)resulted in the strongest SOC losses with an average loss of 0.43+/-0.18(range of 0.05 to 0.74) compared to the results obtained in the CTL scenario (Fig 5 and Table 3).In contrast, physical hindrance (B~1,) had a weaker effect (0.36 +/-0.16 of average relative SOC loss, ranging from 0.04 to 0.68) while water availability (B >1.1,) had virtually no effect on the mean relative SOC loss with 0.34 +/-0.15,ranging from 0.02 to 0.65 to compare to the (Fig 5b, Table 3).

Vertical carbon fluxes
We present the cumulative amount of carbon uptake from the atmosphere to the soil (positive values) or emitted by the soil (negative values) due to erosion after 200 years of transient simulations (Fig. 7).Provided that C input remains unaffected by erosion in our simulations (CTL scenario), a higher erosion rate increased the total C uptake due to the enhanced dynamic replacement, i.e. positive vertical carbon fluxes (Fig 6).For the CTL scenario, vertical carbon fluxes increased almost linearly by 0.28 kgC.m -2 for each additional 1 mm.yr -1 of soil erosion.It should be noted that the variability increases as well when increasing the cumulative soil truncation, because of depth-dependency of the SOC mineralization rate.
As expected, the FB scenario resulted in substantially lower values for the vertical carbon fluxes (Fig. 6).However, most of the simulations still resulted in net C uptake with an average value of 0.41 +/-0.21kgC.m -2 (-30% compared to the CTL scenario) (Table 3, Fig 6, Fig 7).While higher erosion rates resulted in higher values of vertical carbon fluxes, the FB scenario induced a much larger variability, relative to the CTL scenario, particularly in the case of nutrient depletion and physical hindrance (Fig 6).Furthermore, the simulations showed the large effect of nutrient depletion and physical hindrance on the vertical carbon fluxes; relative to CTL, FB reduced the vertical carbon fluxes by 71% and 45%, respectively (Table 3, Fig 6).
Conversely, convex relationships resulted in estimates which were very similar to the CTL scenario in terms of mean response, range and variability (Table 3, Fig 7).

Inclusion of the feedback on the SOC losses
As shown by our simulations, accounting for erosion-crop yield feedbacks increases SOC losses by 3% to 17% relative to the CTL scenario, depending on the cumulative amount of soil truncation and the functional form of the erosion-crop yield feedback.Our results are supported by the findings reported by previous studies (Berhe et al (2005), Gregorich et al. (1998), Lal (2004), Quinton et al. (2010) or Starr et al. (2000)) that suggested that the modfication of key soil properties (such as nutrients, water retention capacity) by erosion degrades the agronomic quality of the soil and decreases the SOC stock triggered by a C input decline.
The model evaluation showed that, for both the CTL and FB scenarios, model predictions were close to the observations for sites having experienced relatively small soil truncation (i.e.short cultivation period or low erosion rates) (Fig. 4).The FB scenario resulted in an overall better prediction because it was able to predict the large relative SOC losses for the cases where intensive erosion took place.However, the addition of the feedback in the model led to contrasting results.On the one hand, and in line with the sensitivity analysis and the feedback nature, SOC losses were higher for sites where yields were more sensitive to erosion due do nutrient depletion or soil thinning when compared to the CTL scenario.The FB scenario showed an increase in the model performance when SOC losses were important (e.g.USA, Denmark, Belgium 1) (Table 2,Fig 4).In contrast, the predictive power of the model was smaller for sites with low SOC losses (e.g.Belgium 2, Spain 2, UK).On the other hand, adding yield feedbacks related to water availability and physical hindrance had little effect on the model results (Fig 4).Model validation, however, indicated that results with and without the erosion feedback on yields were similar.Despite the increased SOC losses, both aggregated statistics and ranges for similar erosion rates did not exhibits clear differences between the scenarios expect for (i) nutrient-limited environments (B < 0.9) under moderate to high erosion rate or (ii) a high erosion rates cases.Based on the results of the FAST analysis (Table 3), where the strong impact of cumulative soil truncation and the feedback nature on SOC losses and fluxes was highlighted, we argue that the small differences between the scenarios in the model validation are mainly due to the short timescale during which the selected sites have been eroding.
Longer timescale simulations (200 years) improved the understanding of the link between erosion-yield and SOC losses and vertical fluxes.The sensitivity analysis highlighted a strong influence of soil erosion rate and yield reduction rate while exhibiting a weak link between SOC response (SOC stock loss and vertical carbon fluxes) and the initial SOC content or the profile shape, as determined by the clay content, the mineralization rate and the root depth distribution (Table 2).Furthermore, differences between the CTL and FB scenarios as well as between the feedbacks nature were larger than those observed in the validation.The addition of the feedback had the most impact for higher erosion rate ,the high sensitivity settings and, to a lesser extent, the physical hindrance.Nutrient depletion had the largest impact owing to their immediate response to erosion and lead to potentially high losses.In contrast, water availability constraints are less distinguishable from the CTL scenario, except in the cases where B is close to 1 and total soil truncation is high (Table 3, Fig  The observed timescale dependency of predictions performance of respective model scenarios can be explained by the nonlinear evolution of the SOC content of a profile exposed to erosion.Without feedback, the SOC stock follows a non-linear evolution divided into two phases.Given the exponential form of the SOC depth profile, a quick initial decrease of the SOC content is followed by a stabilization of SOC content to a steady state level due to an equilibrium between the C input, C uptake from the atmosphere, the lateral export and the C mineralization (Bouchoms et al, 2017, Kirkels et al., 2014;Kuhn et al., 2009;Liu et al., 2003).Under continuous erosion, the rate of C export from a profile is decreasing over time owing to the differential SOC distribution between subsoil and topsoil (Kirkels et al., 2014;Kuhn et al., 2009;Liu et al., 2003).Hence, the fast initial decrease of the SOC stock is linked to the erosion of a SOC-rich topsoil, whereby a small sediment flux may carry a relatively large amount of SOC (Kirkels et al., 2014).In the later stages of the transient simulation, i.e.where the SOC-poor subsoil is exposed to erosion, the SOC loss is smaller for a similar amount of soil truncation (Kirkels et al., 2014).At this stage, the impact of the erosion-yield feedback becomes more important and drives the SOC stock decline.Depending on the erosion rate, the first phase can last for several decades before a steady state could be reach.High sensitivity feedbacks on SOC input such as nutrient depletion and physical hindrance tends to increase the SOC losses in the first phase while the effect of low sensitivity effect may be partially masked until the later stages of the above-described SOC loss processes.

SOC dynamics in eroding landscape: discussion of the addition of erosion-yield feedback on the vertical C balance
In eroding landscapes, several studies have highlighted that a fraction of the erosional SOC loss is replaced by new photosynthate, thereby creating a local atmospheric carbon sink (Harden et al., 1999;Berhe et al., 2007;Van Oost et al., 2007a).
Although much weaker than the C release rate from land cover conversion or SOC lateral export, this so-called dynamic replacement operates on long time scales and can be sustained as long as (i) new C-depleted subsoil material is exposed to the surface and (ii) new C inputs, mainly from plants are available (Doetterl et al., 2016;Kirkels et al., 2014).Both conditions can be questioned here, particularly for landscapes having experienced intense cultivation, and hence erosion, for several centuries.
We do not address the first one here as it not the main focus of this study.As for the second one, when not including the feedback between erosion and yield, one assumes a constant C input to the soil and the absence of soil resource limitations (e.g.water, nutrients, rooting depth).
In their meta-data analysis, Bakker et al. (2004) highlighted that deeply truncated soils exhibit a large reduction in crop productivity.Our results showed that reducing C input in response to long-term erosion actually decreased the SOC stocks by 5% to 74% in the intense eroding cases ( Simulations pointed out that intense sustained erosion coupled with sufficient C input reduction can turn the soil profiles into a net C source for the atmosphere when the C input to the soil profile become smaller than the mineralization rate, due to the decreasing yields.FAST analysis (Table 2) indicates that most of the uncertainty in our simulations is related to the combination of yield reduction strength, erosion rate and mineralization rate, which stressed the importance of (i) accurate SOC profile calibration and linking erosion and vegetation.
Although the model provided estimations of vertical carbon fluxes, which were of the correct order of magnitude and represented the relative differences between the sites well, there is an overall underestimation of the net vertical balance.
Although we derived the functional form of the feedbacks from Bakker et al (2004), yields reduction effects on SOC and vertical carbon fluxes as presented in Fig. 5 should be carefully interpreted.SOC content and cumulative vertical fluxes seem to be much more sensitive to high sensitivity relationships (B < 0.9) than threshold ones (convex, B > 1.1), although this observation is a direct consequence of the equation nature.With only the exponent varying, the yield reduction functions intersect only when the soil truncation is zero or equals 1 meter; it can be questioned to what extent this assumption is representative of real environmental situations.Furthermore, the investigated soil truncation range (60cm) was not sufficient to pass the threshold of yield degradation in convex cases (B > 1.1).The observations on which the feedback relationships were based, were produced by the comparative plots methods, which although better at representing the erosion effect on crop yields than desurfacing or transect methods, may not be as accurate as local on-site evaluations (Bakker et al, 2004).
Furthermore, the assumption of a linear decrease of C input in response to yield may not be accurate in some cases.Particularly, in shallow soils environment, plants tend to adapt their root morphology and increase their root density in response to limited rooting-depth, leading to a slower decline of both C inputs and C stocks over time (Bardgett et al., 2014;Jin et al., 2017), Kosmas et al, 2001).It can be argued that a better fit of the equation to specific conditions would improve the prediction power of the model particularly in the context of nutrient depletion situations where our model overestimated the losses and underestimated the C uptake for low rates of soil truncation.Furthermore, it should be noted that agricultural practices have changed drastically during the last 60 years.Mechanization and intense usage of amendments and fertilizers tend to counterbalance the yield loss consecutive to long term erosion (Gregorich et al., 1998;Doetterl et al., 2016;Kirkels et al., 2014).
Therefore, SOC content and yield evolution may have been partially decoupled whereby, without soil depth constraints, soil erosion has not altered the yields while still driving the SOC losses in eroding landscapes (Meersmans et al., 2009;Bakker et al., 2007;Fenton et al., 2005).

Conclusion
Based on a meta-data analysis, we extracted three main functional relationships linking soil truncation and crop yields discretized by the value of the exponent: two non-linear relationships corresponding to a high-sensitivity to erosion (i.e.nutrient depletion), a low-sensitivity threshold behavior (i.e.water availability) and a linear relationship representing soil thinning.We implemented the erosion-yield feedback in a simple but depth-explicit model of SOC dynamics.The integration of the erosion- SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1 :
Figure 1: Relative crop yields decrease as a function of soil truncation, as shown from yields desurfacing experiments.Datapoints are taken from the metadata analyses by Bakker et al. (2004).The three shaded areas represent the space of the relationships investigated in our study.Dark blue, cyan and orange shades denote respectively the nutrient depletion (B =< 0.9, concave relationships), physical hindrance (0.9 < B < 1.1, linear relationship) and water availability cases (B > 1.1, convex relationships).
Carbon turnover is then coupled to the clay content depth profile through a depth-dependent humification factor.Yields are updated each year following Eq.1, in relation to the cumulative soil truncation.Crop yields affect the SOC content by modifying the amount of SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.soil C inputs.Under the absence of site specific data on the relation between yield and soil C inputs, we here assume a linear relationship between the crop yield and the C input:

Figure 2 :
Figure 2: conceptual framework of processes represented in the model.Black arrows designate processes included in earlier versions of depth-discretized ICBM and red arrows highlight the added explicit processes in this study.
SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.We explicitly accounted for the uncertainties associated with the estimation of model parameters during the transient simulations by building a set of 1000 model parameters sets for each site.The parameter sets combine fixed values (for temperature) and randomly generated parameters in a prescribed range assuming a uniform distribution: ur and  were allowed to vary by +/-2% around the optimum value.Erosion rate, clay content and yield reduction exponent (when available in observations) were constrained around the reported values per sites, with respective tolerances of +/-0.05 mm.yr -1 around the

Figure 3 :
Figure 3: Measured (blue) and simulated (red) SOC profiles for the specific sites (Van Oost et al., 2007).Data of the Spanish sites were not available and the same SOC profile has been used for both US sites.
clearly show a cluster of SOC losses in this range (Table 1b, Fig 4a).Site-specific simulations produce SOC losses, which are in line with the ones estimated by Van Oost SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.et al (2007).Without feedback, the model produced relative SOC losses ranging from only 0.05+/-0.01(Greece) to 0.19 +/-0.01 (USA) of the initial carbon content.

Figure 4 :
Figure 4: Modelled against observed (A) relative SOC losses and (B) vertical C fluxes.Colors denote the different scenarios: Control scenario (CTL, red) and feedback scenario (FB, blue).Error bars represent one standard deviation from the mean for both observed and modeled values.
Fig. 6 also shows that, the vertical carbon fluxes can be close to zero or even 5 negative for those cases where erosion resulted in a large reduction in C input (and this represents a net emission of C to the atmosphere).

Figure 5 :
Figure 5: (a) Relative SOC loss after 200 years for the feedback scenario as a function of erosion rate (mm.yr-1).(b) Relative SOC loss for the FB scenario against the relative SOC loss of the CTL scenario, for all the erosion rates.The 10 colors scale represents the exponent B value, standing for the yield constraint form: B<0.9 exhibits a high sensitivity to small truncation (concave relationship) and B>1.1 shows low sensitivity to small soil truncation (threshold relationship, convex).If 0.9<B<1.1,yields decrease linearly with soil truncation.

Figure 6 :
Figure 6: (a) Cumulative vertical C fluxes (kgC.m -²) after 200 years for the feedback scenario as a function of erosion rate (mm.yr-1).(b) Cumulative vertical C balance (kgC.m -²) for the FB scenario against the cumulative vertical C balance (kgC.m -²) loss of the CTL scenario, for all the erosion rates.The colors scale represents the exponent B value, standing for the yield constraint form: B<0.9 exhibits a high sensitivity to small truncation (concave relationship) and B>1.1 shows low sensitivity to small soil truncation (threshold relationship, convex).If 0.9<B<1.1,yields decrease linearly with soil truncation.

Figure 7 :
Figure 7: Comparison the cumulative vertical C balance (kgC.m -²) between the CTL dataset (grey) and the FB dataset (colors) after 200 years of transient simulations.Positive values indicate a net flux from the atmosphere to the soil.Green, blue, light blue and orange represent respectively the FB scenario results, nutrient depletion, physical hindrance 5 and water availability.Nutrient depletion, physical hindrance and water availability statistics are subsets of the FB scenario results dataset.

Fig 5 )
while producing results in line with observed SOC losses (seeabove and Fig   4).AsHarden et al. (1999) andDoetterl et al. (2016) reported, taking into account the erosion feedback on productivity leads to a better estimation of the C budget and particularly the dynamic replacement, which could be overestimated when ignoring the erosion-yield feedback, particularly when considering longer timescales.Our study supports these assumptions, finding that, compared to the CTL scenario, cumulative C balance decreased on average by 15% to 71% after 200 years depending on the feedback nature (Fig 6 and Fig 7).SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.

Table 1 :
(a) Characteristics of the study sites used for the model validation.Site selection, observed range of relative SOC loss and cumulative vertical balance are from Van Oost et al (2007).Modeled range of SOC losses and cumulative vertical C balance for the no-feedback (CTL) and the erosion-crop yield feedback-scenario (FB).RRSME is calculated over the whole 1m profile between observed and optimized SOC profile, (b) RRMSE of the CTL and FB scenario for each location and as well as the RRMSE of each scenario, including all observations (all).SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.
Table 2a include the effect of the erosion intensity with a large range of erosion tested while Table 2b assumes a constant erosion rate of 1 mm.yr 1 .SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.

Table 3 :
Relative SOC loss and cumulative C fluxes (kgC.m -2 ) after 200 years of transient simulations for the CTL dataset and the FB dataset.Results are given for the whole FB dataset and for the sub-sets of the FB dataset corresponding to nutrient depletion, physical hindrance and water availability.SOIL Discuss., https://doi.org/10.5194/soil-2018-23Manuscript under review for journal SOIL Discussion started: 26 July 2018 c Author(s) 2018.CC BY 4.0 License.