CE-DYNAM (v1), a spatially explicit, process-based carbon erosion scheme for the use in Earth system models

. Soil erosion by rainfall and runoff is an important process behind the redistribution of soil organic carbon (SOC) over land, hereby impacting the exchange of carbon (C) between land, atmosphere and rivers. However, the net role of soil erosion in the global C cycle is still unclear as it involves small-scale SOC removal, transport and re-deposition processes that can only be addressed over selected small regions with complex models and measurements. This leads to uncertainties in future projections of SOC stocks and complicates the evaluation of strategies to mitigate climate change through increased SOC sequestration. In this study we present the parsimonious process-based Carbon Erosion DYNAMics model (CE-DYNAM) that links sediment dynamics resulting from water erosion with the C cycle along a cascade of hillslopes, floodplains and rivers. The model simulates horizontal soil and C transfers triggered by erosion across landscapes and the resulting changes in land-atmosphere CO 2 fluxes at a resolution of about 8 km at the catchment scale. CE-DYNAM is the result of the coupling of a previously developed coarse-resolution sediment budget model and the ecosystem C cycle and erosion removal model derived from the ORCHIDEE land surface model. CE-DYNAM is driven by spatially explicit historical land use change, climate forcing, and global atmospheric CO 2 concentrations affecting ecosystem productivity, erosion rates and residence times of sediment and C in deposition sites. The main features of CE-DYNAM are (1) the spatially explicit simulation of sediment and C fluxes linking hillslopes and floodplains, (2) the relative low number of parameters that allow running the model at large spatial scales and over long-time scales, and (3) its compatibility with global land surface models, hereby, providing opportunities to study the effect of soil erosion under global changes. We present the model structure, concepts, limitations and evaluation at the scale of the Rhine catchment for the period 1850-2005 AD. Model results are validated against independent estimates of gross and net soil and C erosion rates, and the spatial variability of SOC stocks from high-resolution modeling studies and observational datasets . We show that despite local differences, the resulting soil and C erosion rates, and SOC stocks from CE-DYNAM are comparable to high-resolution estimates and observations at sub-basin level. 39 41 42 We find that soil erosion mobilized around 66 ± 28 Tg (10 12 g) of C under changing climate and land use over the non-Alpine region of the Rhine catchment, assuming that the erosion loop of the C cycle was in near steady-state by 1850. This caused a net C sink equal to 2.1 - 2.7% of the Net Primary Productivity of the non-Alpine region over 1850-2005 AD. This sink is a result of the dynamic replacement of C on eroding sites that increases in this period due to rising atmospheric CO 2 concentrations enhancing the litter C input to the soil from primary production.

In the ORCHIDEE LSM, terrestrial C is represented by eight biomass pools, four litter pools and three SOC pools. Each of the pools varies in space, time and over the twelve Plant Functional Types (PFTs). An extra PFT is used to represent bare soil. Anthropogenic and natural disturbances (as a result of climatic changes) to the C pools include fire, crop harvest, changes to the Gross Primary Productivity (GPP), litterfall, autotrophic and heterotrophic respiration (Krinner et al., 2005;Guimberteau et al., 2018). The C-cycle processes are represented by a C emulator that reproduces for each PFT all C pools and fluxes between the pools exactly as in ORCHIDEE in absence of erosion. A net land use change scheme is included in the emulator with mass-conservative bookkeeping of SOC and C input when a PFT is changed into another from anthropogenic land use change (Naipal et al., 2018). The sediment budget model has been added in the emulator to simulate large-scale long-term soil and SOC redistribution by water erosion using coarse-resolution precipitation, land-cover and LAI data from Earth System Models (Naipal et al. 2015(Naipal et al. , 2016. The C emulator including erosion removal was developed by Naipal et al. (2018) to reproduce the SOC vertical profile, removal of soil and SOC starting from the topsoil, and compensatory SOC storage from litter input. As soil erosion is assumed not to change soil and hydraulic parameters but only the SOC dynamics, the emulator allows substituting for the ORCHIDEE model and performing simulations on time scales of millennia with a daily time step, which would be a very computationally expensive or nearly impossible with the full LSM. The concept and all equations of the emulator are described in Naipal et al. (2018). The following subsections describe the different components of the CE-DYNAM that couples the C and soil removal scheme (Naipal et al., 2018) with the horizontal transport and burial of eroded soil and C (Naipal et al., 2016).

The soil erosion scheme
The potential gross soil erosion rates are calculated by the Adjusted Revised Universal Soil Loss Equation (Adj. RUSLE) model (Naipal et al., 2015), which is based on the Revised Universal Soil Loss Equation (RUSLE) (Renard et al., 1997) and part of the sediment budget model (Naipal et al., 2016) (Fig 1). In the Adj. RUSLE the yearly average soil erosion rate is a product of rainfall erosivity ( R ), slope steepness ( S ), land cover and management ( Cm ) and soil erodibility ( K ): Note that the original RUSLE model further includes a slope-length factor ( L ), which gives the length of a field in the direction of steepest descent, and a support practice factor ( P ), which accounts for management practices to mitigate soil erosion. These two factors have been excluded here, because their quantification still includes many uncertainties and is not practical for applications at regional to global scales. These factors are largely affected by local man-made structures (such as field size) and management practices, which are difficult to assess for present day and whose changes over the past are even more uncertain. In addition, we focus in this study on the potential effect of soil erosion on the C budget without erosion-control (EC) practices. Naipal et al. (2015) have developed a methodology to derive the S and R factors from 5 arcmin resolution (5 x 5 arcminute raster) data on elevation and precipitation, hereby preserving the high-resolution spatial variability in slope and temporal variability in erosivity. In the rest of the manuscript we will refer to X by X km/arcminute raster cells alway with X km/arcmin resolution. Despite the comparatively coarse resolution of the erosion model, the so derived R factor was shown to compare well with the corresponding high-resolution product published by Panagos et al. (2017). In the study of Naipal et al. (2016), where the soil erosion model was applied for the last millennium, the change in climate was taken into account in the calculation of the R factor. For this study, we assume that the climate zones as defined by the Koeppen-Geiger climate classification have not changed drastically since 1850 AD.

The sediment deposition and transport scheme
The sediment deposition and transport scheme is adapted from the sediment budget model described by Naipal et al. (2016), which was calibrated and validated for the Rhine catchment (Fig 1). In the sediment budget model rivers and streams are not explicitly simulated. Instead each grid cell contains a floodplain fraction to ensure sediment transport between the grid cells (transport from one grid cell to another can only follow the connectivity of floodplains). It should be noted that global soil databases do not identify floodplain soil as a separate soil class, although national soil databases might. Because we aim to present a carbon erosion model that should be also applicable for other similar catchments, we followed a two-step methodology to derive floodplains in the Rhine catchment. For this purpose we used hydrological parameters and existing data on hillslopes and valleys. First, grid cells were identified that consisted entirely out of floodplains. For this we used the gridded global data set of soil at 5 arcminute resolution, with intact regolith, and sedimentary deposit thicknesses of Pelletier et al. (2016) (Table 1), and identified lowlands and hillslopes based on soil thickness and depth to bedrock. The lowlands were classified as grid cells that contain only floodplains and no hillslopes.
Second, we calculated the floodplain fraction ( A fl ) of a grid cell i that has both hillslopes and floodplains as a function of stream length and width based on the methodology developed by Hoffmann et al. (2007): Where ,L stream is the stream length derived from the HydroSHEDS database (Lehner and Grill, 2013) (Table 1). Where, A upstream is the upstream catchment area, and a is equal to 60.8, and b is equal to 0.3.
The parameters a and b have been derived using the scaling behavior of floodplain width as estimated from measurements on the Rhine (Hoffmann et al., 2007). The sediment deposition on hillslopes ( D hs ) and floodplains ( D fl ) is calculated as a function of the gross soil removal rates ( E ) according to Naipal et al. (2016) with the following equations: Where, f is the floodplain deposition factor at 8 km resolution that determines the fraction of eroded material transported and deposited in the floodplain fraction of a grid cell. a f and b f are constants that relate f to the average topographical slope ( θ ) of a grid cell depending on the type of land cover . θ max us the maximum topographical slope of the entire Rhine catchment.
The parameters a f and b f are chosen in such a way that f varies between 0.2 and 0.5 for cropland, reflecting the decreased sediment connectivity between hillslopes and floodplains created by man made structures such as ditches and hedges. For natural vegetation such as forests and natural grassland, a f and b f are chosen in a way that f varies between 0.5 and 0.8 assuming that in these landscapes hillslopes and floodplains are well-connected. This assumption on the reduced sediment connectivity for agricultural landscapes is supported by several previous studies on the effect of erosion on sediment yield (Hoffmann et al., 2013a;de Moor and Verstraeten, 2008;Gumiere et al., 2011;Wang et al., 2015). These studies showed that man-made activities on agricultural landscapes result in a trapping of eroded soil in colluvial deposition sites, reducing the sediment transport from hillslopes to floodplains. The model parameter f has been calibrated for the Rhine catchment by Naipal et al. (2016), where the ranges mentioned above are found to produce a ratio between hillslope and floodplain sediment storage that was comparable to observations. The studies of Wang et al. (2010; identified a range for the hillslope sediment delivery to be between 50 and 80 %, which is similar to the range in the (1-f) factor in our model. In each case and within the defined boundaries, the slope gradient determines the final value of f . Eroded material that has not been deposited in the floodplains stays on the hillslopes and is assumed to be deposited at the foot of the hillslopes as colluvial sediment.
The floodplain fractions of the grid cells are connected through an 8 km resolution flow routing network (Naipal et al., 2016), where the rivers and streams are indirectly included in the floodplain area but not explicitly simulated. By routing the sediment and C through the floodplain fractions of grid cells we lump together the slow process of riverbank erosion by river dynamics (time scale ≈ a few years to thousands of years), and the rather fast process of transport of eroded material by the rivers (time scale ≈ days). The rate by which sediment and SOC leave the floodplain of a grid cell to go to the floodplain of an adjacent grid cell is determined by the sediment residence time. The sediment residence time ( τ ) is a function of the upstream contributing area ( Flowacc ): The study of Hoffmann et al. (2008) showed that the majority of floodplain sediments have a residence time that ranges between 0 and 2000 years, with a median of 50 years. The constants a τ and b τ are chosen in such a way that basin τ varies between the 5th and 95th percentile of those observations, with a median for the whole catchment of 50 years. These constants are uniform for the whole basin. These constants need to be calibrated based on local data of sediment ages before CE-DYNAM can be applied to other catchments.
Floodplain C storage follows the same residence time as sediment on top of the actual decomposition rate of C in a grid cell of ORCHIDEE. The routing of sediment and C between the grid cells follows a multiple-flow routing scheme. In this scheme the flow coming from a certain grid cell is distributed across all lower-lying neighbors based on a weight ( W , dimensionless) that is calculated as a function of the contour length ( c ): Where c is 0.5xgrid size (m) in the cardinal direction and 0.354xgrid size (m) in the diagonal direction. ( i, j ) is the grid cell in consideration where i counts grid cells in the latitude direction and j in the longitude direction. i + k and j + l specify the neighboring grid cell where k and l can be either -1, 0 or 1. θ is calculated as the division between the difference in elevation ( h ) give in meters difference and the grid cell size ( d ), also in meters: The sediment and C routing is done continuously at a daily time-step to preserve the numerical stability of the model. More detailed explanation of the methods presented in this section can be found in the study of Naipal et al. (2016).

Litter dynamics
The four litter pools in the emulator are an below-and an above-ground litter pool, each split into a metabolic and structural pool with different turnover rates as implemented in ORCHIDEE (Krinner et al., 2005). The belowground litter pools consist mostly out of root residues. Both the biomass and litter pools have a loss flux due to fire as incorporated in ORCHIDEE by the Spitfire model of Thonicke et al. (2010). The litter that is not respired or burnt is transferred to the SOC pools based on the Century model (Parton et al., 1987), which was modified by Naipal et al. (2018) to include a vertical discretization scheme for SOC.
The vertical discretization scheme was introduced in the emulator to account for a declining C input and SOC respiration with depth, and consists of 20 soil layers with each 10 cm thickness. The litter to soil fluxes from aboveground litter pools are all attributed to the top 10 cm of the soil profile. The litter to soil fluxes from belowground litter pools are distributed exponentially over the whole soil profile according to: Where I 0be is the below-ground litter input to the surface soil layer and r is the PFT-specific vertical root-density attenuation coefficient as used in ORCHIDEE. The sum of all layer-dependent litter to soil fractions is equal to the total litter to soil flux as calculated by ORCHIDEE. The vertical SOC profile is modified by erosion and the resulting deposition fluxes, which is discussed in detail in the following sections.

Crop harvest and yield
We adjusted the representation of crop harvest from ORCHIDEE by assuming a variable harvest index for C3 plants that increases during the historical period as shown in the study of Hay (1995) for wheat and barley, which are also the main C3 crops in the Rhine catchment. The harvest index is defined by the ratio of harvested grain biomass to above-ground dry matter production (Krinner et al., 2005). In this study the harvest index increases linearly between 0.26 and 0.46 (Naipal et al. 2018) consistent with the average values of Hay (1995). We also found that in certain cases the cropland Net Primary Productivity (NPP) was too high during the entire period of 1850-2005, especially in the early part of the 20 th Century. This is because the cropland photosynthetic rates were adjusted in ORCHIDEE to give a cropland NPP representative of present day values that are higher than for the low input agriculture of the early 20 th Century. To derive a more realistic NPP for crop and barley in the Rhine catchment we used the long-term crop yield data obtained from a dataset on 120000 yield observations over the 20 th century in Northeast French Départements (NUTS3 administrative division) (Schauberger et al., 2018). According to the yield data assembled by Schauberger et al. (2018), yields in Northeast France (covers part of the Rhine catchment) for these crops increased fourfold during the last century. Note that crop residues like straw constituted a larger fraction of the total biomass in 1850 than in 2005, but those residues were likely collected and used for animal feed, housing fuel. We did not account for this harvest of residue in the simulation of SOC.

SOC dynamics without erosion
The change in the C content of the PFT-specific SOC pools in the emulator without soil erosion was described by Naipal et al. (2018) (Fig 1) as follows: Where, SOC a , SOC s , and SOC p (g C m -2 ) are the active, slow and passive SOC, respectively. The distinction of these SOC pools, defined by their residence times, are based on the study of Parton et al. (1987). The active SOC pool has the lowest residence time (1 -5 years) and the passive the highest (200-1500 years). lit a and lit s (g C m -2 day -1 ) are the daily litter input rates to the active and slow SOC pools, respectively; k0 a , k0 s and k0 p (day -1 ) 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.
The vertical C discretization scheme in the emulator assumes that the SOC respiration rates decrease exponentially with depth: Where k i is the respiration rate at a soil depth z and re (m -1 ) is a coefficient representing the impact of external factors, such as oxygen availability that decreases with depth. k 0 is the respiration rate of the surface soil layer for a certain SOC pool i .
The variable re is determined in such a way that the total soil respiration of a certain pool over the entire soil profile without erosion is similar to the output of the full ORCHIDEE model. Detailed description of how this is done can be found in the study of Naipal et al. (2018).

Net C erosion on hillslopes
In the model we assume that soil erosion takes place on hillslopes, and not in the floodplains due to the usually low topographical slope of floodplains. The factor (1-f) determines the fraction of the eroded soil that is deposited in the colluvial reservoirs (Fig 1). Soil erosion always removes a fraction of the SOC stock in the upper soil layer depending on the erosion rate and bulk density of the soil. The next soil layer contains less C and therefore at the following time-step less C will be eroded under the same erosion rate. To account for this effect, the SOC profile evolution is dynamically tracked in the model and updated at a daily time step, conform with the method of Wang et al. (2015). First, a fraction of the C from each soil pool in proportion to the erosion height is removed from the surface layer. Then, at the same erosion rate, SOC from the subsoil layer becomes the surface layer, maintaining the soil layer thickness in the vertical discretization scheme. Similarly, the SOC from the subsoil later also moves upward one layer. The removal of C by erosion also triggers a compensatory C sink due to the reduction in SOC respiration on eroding land. This compensatory C sink and reduced C erosion over time will ultimately lead to an equilibrium state. The change in C content due to net erosion (the eroded sediment/C that leaves the hillslopes after deposition) of the PFT-specific pools for hillslopes can be represented by the following equations: Where dSOC HSi (z,t) is the change in hillslope SOC of a component pool i at a depth z and at time step t . The daily net erosion fraction k E (dimensionless) is calculated as following: Where, E is the gross soil erosion rate (t ha 2 year -1 ), f is the floodplain deposition factor, BD is the average bulk density of the soil profile (g cm -3 ), dz is the soil thickness (=0.1 m), and EF is the C enrichment factor that is set to 1 by default. A model sensitivity analysis will be performed (see section 4.3) with EF>1 to represent a higher C concentration in eroded soil compared to the original soil as a result of the selectivity of erosion.
Hillslope erosion without the deposition term has already been tested and applied at the global scale as part of the C removal model presented by Naipal et al. (2018).

C deposition and transport in floodplains
The SOC-profile dynamics of floodplains are controlled by: (1) C input from the hillslopes, (2) C import by lateral transport from the floodplain fractions of upstream grid cells, and (3) C export to the floodplain fractions of downstream g grid cells (Fig 1). First, the net erosion flux from the surface layer of the hillslope fraction of the grid cell ( k E × SOC HS (z=0)) is incorporated in the surface layer of the floodplain. At the same deposition rate, the SOC of the surface layer of the floodplain is incorporated in the subsoil layer. Similarly, a fraction of the SOC of the subsoil layer is moved downward one layer. We will refer to this process as the 'downward' moving of C in the soil layer profile. It should be noted that C selectivity during transport and deposition is not taken into account here, meaning that the C pools of the deposited material are the same as the eroded material from the topsoil of eroding areas. At the same time a fraction of the C of the surface layer proportional to the sediment residence time ( τ ) is exported out of the catchment following the sediment routing scheme, resulting in the 'upward' moving of the C from the subsoil layers. This process represents the river bank erosion and resulting POC export by the water network, although rivers and streams are not explicitly represented in the model. As we do not have information on the sub-grid spatial distribution of land cover fractions we first sum the exported C flux over all PFTs before assigning the flux proportionally to the land cover fractions of the receiving downstream-located grid cells. The C that is imported from the neighboring grid cells follows the same procedure as the deposition of eroded material, and results in a 'downward' moving of the C in the soil profile. The change in C content due to deposition and routing of the PFT-specific SOC pools for floodplains can be represented by the following equations: , for z=0 (17) Where n is the neighboring grid cell that flows into the current grid cell, dSOC FLi (z,t) is the change in floodplain SOC of a component pool i at a depth z and at time step t, and SOC HS is the hillslope SOC stock. k D is the deposition rate and equal to: Where AREA HS is the hillslope area and AREA FL is the floodplain area (m 2 ). is the import rate per C pool i from k i out neighboring grid cells (dimensionless) and can be calculated as: Where, W is the weight index of equation 7.
The first term of equation 16  ), and part is moved 'downwards' in the soil profile as a result of C deposition ( k D ) and the incoming later C from upstream grid cells ( ). k i out

The land use change bookkeeping model
The land use change bookkeeping scheme includes the yearly changes in forest, grassland and cropland areas in each grid cell as reconstructed by Peng et al. (2017) (Table 1). Peng et al. (2017) derived historical changes in PFT fractions based on LUHv2 land use dataset (Hurtt et al., 2011), historical forest area data from Houghton, and present day forest area from ESA CCI satellite land cover (European Space Agency, ESA, 2014). By using different transition rules and independent forest data to constrain the changes in crop and urban PFTs they derived the most suitable historical PFT maps.
When land use change takes place, the litter and SOC pools of all shrinking PFTs are summed and allocated proportionally to the expanding PFTs, maintaining the mass-balance. In this way the litter pools and SOC stocks get impacted by different input and respiration rates for each soil layer. When forest is reduced, three wood products with decay rates of 1, 10 and 100 years are formed and harvested. The biomass pools of other shrinking land cover types are transformed to litter and allocated to the expanding PFTs. More details on the land use scheme are described in the study of Naipal et al. (2018).

Study-Area
The model is tested for the Rhine catchment (Fig 2), which has a total basin area of about 185,000 km 2 covering five different countries in Central Europe. Its large size is beneficial for the application of a coarse-resolution model such as CE-DYNAM to study large-scale regional dynamics in the C cycle due to soil erosion. The Rhine catchment has a contrasting topography, with steep slopes larger than 20% upstream in the Alps, and large, wide and flat floodplains at the foot of the Alps, the upper Rhine and the lower Rhine. The floodplains store large amounts of sediment and C that originate from eroding hillslopes upstream. These sediment storages provide the possibility to study the long-term effect of erosion on hillslope and floodplain dynamics. Furthermore, the Rhine catchment has been experiencing different stages of land use change over the Holocene, with land degradation dating back to more than 5500 years ago (Dotterweich, 2013). In contrast, during the last two decades there has been a general afforestation and soil erosion has been decreasing. These land use changes and changes in erosion make an interesting and important case to study the effect of anthropogenic activities on the C cycle in Europe. In addition, the Rhine catchment has been the focus of many erosion studies providing observations on erosion and sediment dynamics that can be used for model validation (Asselman, 1999;Asselman et al., 2003;Erkens, 2009;Hoffmann et al., 2007Hoffmann et al., , 2008Hoffmann et al., , 2013aHoffmann et al., , 2013bNaipal et al., 2016). The global sediment budget model that forms the basis for the sediment dynamics scheme of CE-DYNAM has been validated and calibrated for the Rhine catchment with observations on sediment storage from Hoffmann et al. (2013a) and the derived scaling relationships between sediment storage and basin area (Naipal et al., 2016). Hoffmann et al. (2008Hoffmann et al. ( , 2013a) did an inventory of 41 hillslope and 36 floodplain sediment and SOC deposits related to soil erosion over the last 7500 years. The floodplain sediment observations consist mostly out of organic material (gyttja, peat) and fine sediments (fine sand, loam, silt) in overbank deposits (Hoffmann et al., 2008).
These fine sediments are a result of long-term soil erosion on the hillslopes. Hoffmann et al. (2013a) found that the sediment and SOC deposits were quantitatively related to the basin size according to certain scaling functions, where floodplain deposits increased in a non-linear way with basin size while the hillslope deposits showed a linear increase with basin size. We use these relationships to validate the spatial variability in SOC storage of floodplains and hillslopes simulated by CE-DYNAM. The scaling relationships have the form of a simple power law: Where M is the sediment storage or the SOC storage, a is the storage (Mt) related to an arbitrary chosen area A ref , while b is the scaling exponent.

Input data and model simulations
To create the C emulator that forms the underlying C cycle part of CE-DYNAM, we first ran the full ORCHIDEE model for the period 1850-2005 at a coarse resolution of 2.5°degrees latitude and 3.75° degrees longitude, and output all C pools and fluxes. The pools and fluxes were then archived together and used to derive the turnover rates to build the emulator.
The SOC scheme of the emulator that has been modified to account for soil erosion processes has been made to run at a spatial resolution of 5 arcmin, similar to the original global sediment budget model. Then, we performed three main simulations with CE-DYNAM for the Rhine catchment. Simulation S0: The baseline simulation or no-erosion simulation, where SOC dynamics are similar to the full ORCHIDEE model. Simulation S1: The erosion -only simulation, where the hillslopes erode and all eroded C is respired to the atmosphere without reaching the colluvial and alluvial deposition sites.
Simulation S2: The simulation with full sediment dynamics where hillslopes and floodplains are connected and can store or lose C. We ran the emulator for 3000 years at a daily time step with the initial climate and land cover of the period We also performed seven additional sensitivity simulations and four additional uncertainty simulations. Simulation S1_EF and S2_EF are performed to test the model assumption of C enrichment during erosion. Here, we changed the enrichment factor EF to two, based on the study of Lugato et al. (2018). Simulations S2_Tmin and S2_Tmax are performed to test the rate of C transport between floodplains. Here we modified the mean sediment residence time for the Rhine catchment to a minimum of 60 years (50 % lower than the current value), and to a maximum of 128 years (50 % higher than the current value), respectively. However, we kept the maximum sediment residence time at 1500 years. Simulations S0_RM, S1_RM and S2_RM are performed to test the model assumption on crop residue management, where we assumed that all above-ground crop litter is harvested.
For the uncertainty analysis we performed simulations S1_min and S2_min based on a minimum soil erosion scenario, and S1_max and S2_max based a maximum soil erosion scenario. These soil erosion scenarios are derived from the uncertainty ranges in the rainfall erosivity and land cover factors of the erosion model. All the model simulations are summarized in table 2.

Validation methods and data
We performed a detailed model validation of the sediment and the C parts of the model according to the following steps: (1) validation of soil erosion rates using observational and high-resolution model estimates for Germany and Europe, (2) validation of C erosion rates using high-resolution model estimates for Europe from Lugato et al. (2018), (3) validation of the spatial variability of hillslope and floodplain C storage using observational results from Hoffmann et al. (2013a), (4) validation of SOC stocks using observational data from a global soil database and a European land use survey.
The validation of the soil erosion module has been done before in the studies of Naipal et al. (2015Naipal et al. ( , 2016. However, we do it again in this study due to different input datasets. In addition, the validation includes soil erosion data from new global soil erosion studies such as Borrelli et al. (2018) and Panagos et al. (2015). For the validation of gross soil erosion rates we used the high-resolution model estimates of Panagos et al. (2015), who applied the RUSLE2015 model at a 100 m resolution at European scale for the year 2010. Similarly to the Adj.RUSLE, RUSLE2015 is also derived from the original RUSLE model. However, in contrast to our model, RUSLE2015 does include the erosion factors L and P . Furthermore, our model uses more coarsely resolved input datasets (Table 1), for which the equations for the R and S factors have been modified. The extensive validation of the Adj. RUSLE model in this study and previous studies (Naipal et al., 2015(Naipal et al., , 2016(Naipal et al., , 2018, shows that despite its coarse resolution, it is applicable at large spatial scales. Thus, even though both Adj.RUSLE and RUSLE2015 are derived from the same erosion model, the differences between the models are large, which justifies our model comparison. Furthermore, we used independent high-resolution erosion estimates from the study of Cerdan et al. (2010), available at a 1 km resolution at European scale, which were based on an extensive database of measured erosion rates under natural rainfall in Europe. For the comparison we aggregated the high-resolution model results of both datasets to the resolution of CE-DYNAM. We also used the potential soil erosion map of the Federal Institute for Geosciences and Natural Resources of Germany (Bug and Stolz, 2014). This map presents the yearly average soil erosion rates at a 250 m resolution on agricultural land derived from a USLE-based approach, with some modifications to the erosion factors and input data. Before validating our model results we aggregated these high-resolution erosion rates to the coarser resolution of our model.
Validation of our net soil erosion rates is done based on the 100 m resolution net soil erosion rates derived with the WATEM-SEDEM model (Borrelli et al., 2018). WATEM-SEDEM simulates soil removal by water erosion based on the USLE approach, sediment transport and deposition based on the transport capacity. The model has been extensively employed to estimate net fluxes of sediments across hillslopes at catchment-and regional-scales.
For the validation of C erosion rates, we used the high-resolution model results from Lugato et al. (2018), where they coupled the RUSLE2015 erosion model to the Century biogeochemistry model. These model results were available at a resolution of 1 km, where each grid cell was composed of an erosion and deposition fraction. The C erosion rates provided by Lugato et al. (2018) were multiplied with the erosion fraction of a 1 km grid cell. Then, the C erosion rates were aggregated to the resolution of CE-DYNAM. Lugato et al. (2018) provided an enhanced and a reduced erosion-induced C sink uncertainty scenario, based on different assumptions for C enrichment, burial and C mineralization during transport. In CE-DYNAM the C erosion rates from simulation S1 are multiplied with the hillslope area to get the total C erosion flux of a grid cell. As the study of Lugato et al. (2018) considers only agricultural areas, we considered only the crop fraction of a grid cell. It should be noted that the SOC dynamics scheme of CE-DYNAM, which is derived from ORCHIDEE LSM, is also based on the Century model. However, there are large differences between the Century model used by Lugato et al. (2018) and the C dynamics scheme of ORCHIDEE used in this study. For example, in the Century model the crop productivity is mediated by nitrogen availability, which is not the case in the ORCHIDEE version used for this study. The Century model also includes some management practices such as crop rotations, which are not represented in ORCHIDEE.
The Century model runs at a much higher resolution and is calibrated for agricultural land, while ORCHIDEE also simulates forest, grasslands and bare soil. In this way, the final SOC stocks derived with CE-DYNAM are also a result of erosion from other land cover types and land use changes. This is an important feature for land use change, which is not included in the Century model. Furthermore, the ORCHIDEE LSM has been used in many global intercomparisons and extensively evaluated for C budgets (Mueller et al., 2019;Todd-Brown et al., 2013). Finally, ORCHIDEE also includes the last century change in crop production calibrated against data . For the validation of the spatial variability of the SOC stocks of hillslopes and floodplains we used the scaling relationships between basin area and SOC storage derived by Hoffmann et al. (2013a). The study by Naipal et al. (2016) found that the global sediment budget model is able to reproduce the scaling behaviour of sediment storage. After analyzing the dependence of this scaling behavior, they argue it is an emergent feature of the model and mainly dependent on the underlying topography. This indicates that the scaling features of floodplain and hillslope sediment and C storage should also be applicable to a more recent time period. In order to evaluate the ability of CE-DYNAM to reproduce this scaling behavior, we selected the grid cells that contained the points of observation of the study of Hoffmann et al. (2013a) and . The LUCAS topsoil SOC stocks, available at a high spatial resolution of 500 m, were calculated using the LUCAS SOC content for Europe (de Brogniez et al., 2015) and soil bulk density derived from soil texture datasets (Ballabio et al., 2016).

Results
Due to large uncertainties in the model and validation data for the Alpine region we only present and discuss the model and validation results for the non-Alpine part of the Rhine catchment.

Model validation
In this section we present the model validation results using the methods and data described in detail in the previous section.
We find that the quantile distribution of the simulated gross soil erosion rates compares well to the distributions of other observational and high-resolution modelling studies (Cerdan et al., 2010, Panagos et al., 2015, Bug et al., 2014, although CE-DYNAM usually underestimates the very large soil erosion rates such as is found by Cerdan et al. (2010) (Fig 3A, B, C). This is due to the coarse spatial and temporal resolution of CE-DYNAM, and the lack of the slope-length factor (L) (Cerdan et al. (2010) assumed a constant slope length of a 100m). It should be noted that our study, Cerdan et al. (2010) and Bug et al. (2014) simulated potential soil erosion rates, not accounting for EC practices represented by the P-factor. We also find that the quantile distribution of the simulated net soil erosion from hillslopes compares well with the distribution from the high-resolution modelling study of Borrelli et al. (2018) (Fig 3D). In addition we performed a spatial comparison of our simulated gross and net erosion rates to those of the studies mentioned above. For this purpose we delineated 13 sub-basins in the Rhine catchment ( Fig S3). Table 3 summarizes the resulting goodness-of-fit statistics of this comparison and shows that our erosion model is generally in good agreement with the other studies at sub-basin level, except for net soil erosion. The estimation of the net soil erosion between our study and the study of Borrelli et al. (2018) is done in different ways, which may explain the difference in the results. In our study the deposition of sediment in hillslopes is explicitly calculated as a function of the slope, and vegetation type/cover. Borrelli et al. (2018) used the transport capacity concept (Van Rompaey et al., 2001). Both methods have their uncertainties when applied at large spatial scales.
The method in our study has been designed and calibrated to be used at a large spatial scale, and at coarse resolution, while the method of Borrelli et al. (2018) was originally designed to be applied at spatial scales <100m.
We find that the quantile distributions of our simulated agricultural C erosion and deposition rates are similar to those of the high-resolution modelling study of Lugato et al. (2018) (Figs 4A-D). Also the spatial variability of the C erosion rates at sub-basin level is in good comparison to the validation data (Table 4). However, the linear regression between soil erosion and C erosion rates of our study lies at the lower end of the relationships derived from the enhanced and reduced erosion scenarios of Lugato et al. (2018) (Fig 5). On the one hand, our study does not include EC practices, leading to substantially larger simulated soil erosion rates in regions with EC. Figure 5 shows that our simulated erosion rates are in general larger than the erosion rates from Lugato et al. (2018), which may be explained by this mechanism. On the other hand, the C erosion rates of our study are lower than those of Lugato et al. (2018), due to the coarse spatial resolution of our underlying C-scheme derived from the ORCHIDEE LSM. The decreased spread in our simulated values is also a result of the coarse resolution of our model. Accounting for erosion, deposition and transport of SOC leads to a better representation of the simulated topsoil C stocks per land cover type when compared to SOC stocks of the LUCAS database (Fig 6). The simulated SOC stocks of the top 20 cm of the soil profile fall within the quantile range of the LUCAS SOC stocks for cropland and forest (Fig 6). Although the topsoil SOC stocks for grassland improved, a large uncertainty range remains. Furthermore, we find that in both the erosion and no-erosion simulation the SOC stocks for grassland are higher than for forest. This is also observed in the study of Wiesmeier et al. (2012), where they found considerable higher SOC stocks for grassland with a median of 11.8 kg C m -2 compared to forest based on the analysis of 1460 soil profiles in South-Germany. Furthermore, the comparison of the simulated total SOC stocks to those of the LUCAS and GSDE databases at sub-basin level shows a good model performance with respect to the spatial variability in topsoil SOC stocks (Table 5) (Table 6). However, when deriving the scaling relationships at sub-basin level instead of using individual grid cells we do not find a significant difference in scaling between floodplains and hillslopes (Table 6).

Model application
We find an average annual soil erosion rate of 1.44 ± 0.82 t ha -1 year -1 over the period 1850-2005, which is about half of the average erosion rate simulated for the last millennium (Naipal et al., 2016) and falls within the range of the average erosion rates of the Holocene (Hoffmann et al., 2013). This soil erosion flux mobilized around 66 ± 28 Tg of C over the same time period, of which on average 57 % is deposited in colluvial reservoirs, 43% is deposited in alluvial reservoirs, and 0.2% is exported out of the catchment.
The lower average annual soil erosion rate over the study period compared to the last millennium is a result of the general afforestation in the non-Alpine part of the Rhine catchment that started around 1910 AD according to the data on land cover and land use (Peng et al., 2017;Fig 7B). This land cover data also shows that forest increases by 24 % over the period 1910-2005, mostly as a result of grassland to forest conversion. Cropland decreases by 6 % over the period 1920 and 1970, and is relatively stable afterwards. This afforestation leads to a long-term decreasing trend in gross soil and SOC erosion rates on hillslopes (Fig 7C). The temporal variability in the soil and C erosion rates is a result of direct changes in precipitation, such as the temporary increase in erosion rates over the period 1940-1960 ( Fig 7A). Furthermore, we find that the temporal variability in C erosion rates follows the soil erosion rates closely, indicating that soil erosion dominates the variations in C erosion over this time-period, while increased SOC stocks due to CO 2 fertilization and afforestation play a secondary role as a slowly varying trend. It should be noted that the correlation between soil and C erosion might be affected by processes not properly captured by the model such as the selectivity of erosion including the enrichment of C in eroded material.
The cumulative C erosion removal flux of 66 ± 28 Tg of C leads to a cumulative net C sink for the whole Rhine region of 216 ± 23 Tg C (Fig 7D). This is about 2.1 -2.7 % of the cumulative NPP and of the same magnitude as the cumulative land C sink of the Rhine without erosion. It should be noted that these are potential fluxes, assuming that the photosynthetic replacement of C is not affected by the degradation of soil due to the removal of nutrients, declining water-holding capacity and other negative changes to the soil structure and texture (processes not covered by our model).
The breaking point in figure 7D around 1910 AD is a result of the climate data used as input.
To better understand the erosion-induced net C flux, we analyze the erosion-induced C exchange with the atmosphere by explicitly track the movement of eroded C through all reservoirs (for example between eroding hillslopes and colluvial reservoirs), we make use of the changes in SOC stocks and Net Ecosystem Productivity (NEP), which is the difference between NPP and heterotrophic respiration, of the three main simulations (S0, S1, S2) to derive the erosion-induced vertical C fluxes. By subtracting the NEP of hillslopes (NEP HS ) of the no-erosion simulation (S0) from the erosion-only simulation (S1), we derive the additional photosynthetic replacement of SOC on eroding sites (Eq. 21): Where, E rep is the potential dynamic Photosynthetic replacement of C on eroding sites (assuming no feedback of erosion on NPP). Part of the eroded C that is transported to and deposited in colluvial reservoirs can be respired or buried (Eq. 22).
The difference between NEP of simulation S2 and S1 is the NEP caused by the deposition of eroded C in colluvial areas and equal to the difference between the burial and respiration of C in colluvial sites. As we do not explicitly track the respiration of deposited material in the model, we can only derive the net respiration or net burial of colluvial deposits (Rc net ) with the following equation: The same concept can be applied for the net respiration of floodplains: Where, NEP FL is the floodplain NEP, and Ra net is the net respiration or net burial of alluvial deposits. Positive values for Ra net or Rc net indicate a net burial (respiration S2 < respiration S0/S1) of the deposited material.
We find that the dynamic replacement of C on eroding sites increased by 17-33% at the end of the period despite decreasing soil erosion rates (Figs 8A& B). This increase in the photosynthetic replacement of C is due to the globally increasing CO 2 concentrations that lead to the CO 2 fertilization effect, amplified by the afforestation trend in the Rhine over this period. Without this fertilization effect, soil erosion and deposition would be likely a weaker C sink or even a C source over the period 1850-2005 (Figs S4A & B). This CO 2 fertilization effect promotes a 100% replacement of the eroded C on hillslopes and even leads to a C sink on hillslopes at the end of the study period ( Fig 8B). Furthermore, we find that the yearly average gross C erosion flux from eroding sites decreases by 10 -34 %, while the yearly deposition fluxes in colluvial and alluvial sites decreases by 20 % and 19 -47 %, respectively. The decrease in the deposition flux to floodplains is compensated by a better sediment connectivity between hillslopes and floodplains due to afforestation.
Forests have less man-made structures that can prevent the erosion fluxes from reaching the floodplains, which is represented by a higher floodplain deposition ' f ' factor in the model. The decrease in the erosion flux also leads to a decreased POC export of the catchment at the end of the study period.
We also find that both the colluvial and alluvial reservoirs show a net respiration flux throughout the time period (Figs 8A & B). This is consistent with previous studies who found that deposition sites can be areas of increased CO 2 emissions (Billings et al., 2019;Van Oost et al., 2012). However, there is a slight difference in the respiration of deposited C between the start and end of the transient period. The respiration of deposited SOC in colluvial sites increases with time while the respiration of deposited SOC in alluvial sites shows rather a decreasing trend. These changes in SOC respiration of deposited material depends on (1) the amount of deposited material, (2) increasing temperatures over 1850-2005 for the entire catchment, and (3) the constant removal of C-rich topsoil and its deposition in alluvial and colluvial reservoirs, which makes the deposited sediments generally richer in C than soils on erosion-neutral sites, providing more substrate for respiration. The largest increase in total respiration of alluvial and colluvial deposits takes place in hilly regions due to the initial increase in erosion rates resulting in large deposits of C. Overall, we find that the increased respiration of deposited material slightly offsets the increased dynamic C replacement, however, the dynamic C replacement on eroding sites still dominates the erosion-induced C sink.

Discussion
In this section we discuss some of the most important model limitations, uncertainties and assumptions.

Initial conditions and past global changes
Initial climate and land cover/use conditions, and the length of the transient period are essential parameters that determine the resulting spatial distribution of soil and C. Landscapes are in a constant transient state due to global changes, such as climate change, land use change, accelerated soil erosion. However, we assumed an equilibrium state so that we can quantify the changes during the transient period. The longer the transient period that covers the essential historical environmental changes, the more accurate are the present-day distribution of SOC stocks, sediment storages, and related fluxes. This is especially true when analyzing the redistribution of soil and C as a result of erosion, deposition and transport, as these soil processes can be very slow. For example, the study of Naipal et al. (2016) showed that by simulating the soil erosion processes for the last millennium a spatial distribution of sediment storages that is similar to observations can be found. In this study we simulated the steady state based on the initial conditions of the period 1850-1860 due to constraints in data availability on precipitation and temperature. By focusing only on the period 1850-2005 we miss the effects of significant land use changes in the past that coincided with times of strong precipitation such as in the 14 th and 18 th century (Bork et al., 2003). These major anthropogenic changes in the last Holocene substantially affected the present-day spatial distribution and size of sediment storage and SOC stocks. The absolute value of the SOC storage from the S2 simulations of the non-Alpine region of the Rhine catchment for the year 2005 is in the range of 2.74 -2.99 Pg of C, which is larger than the 1.7 ± 0.6 Pg of C that Hoffmann et al. (2013a) measured. It should be noted that the ORCHIDEE model (S0 simulation) already overestimates the total SOC stock of the non-Alpine region of the Rhine (2.43 Pg of C), when the initial conditions of the period 1850-1860 are used. Due to the fact that we miss the climate and land use changes before the year 1850, we find that floodplains store less SOC than hillslopes. Although this is in contrast to the findings of Hoffmann et al. (2013a), the difference in SOC stocks between floodplains and hillslopes from the S2 simulations is better than the difference derived from the S0 simulation. We find that floodplains store 1.28 -1.72 and hillslopes 1.7 -2 Pg of C when erosion and deposition processes are taken into account, compared to 0.69 Pg of C for floodplains and 2.29 Pg of C for hillslopes when these processes are lacking.
We also find that floodplains have an overall higher C concentration (12 kg m -2 ) compared to hillslopes (9 kg m -2 ) at the end of the transient period (Fig 9A), which is in line with the findings of Hoffmann et al. (2013a) and what can be derived from global soil databases. This is a result of higher SOC concentrations in deeper soil layers of floodplains compared to hillslopes (Figs 9 A&B), as is also shown in the study of Hoffman et al. (2013). To be closer to the observational difference between floodplains and hillslopes we would need to consider the period before 1850, extreme climate events, and a higher plant productivity in floodplains resulting from favorable soil nutrient and hydrological conditions.

Model advantages and limitations
Although we parameterized and applied CE-DYNAM for the Rhine catchment, it is intended to be made applicable to other large catchments. CE-DYNAM combines soil erosion processes, for which small scale differences in topography are of utter importance, with a state-of-the-art representation of large-scale SOC dynamics driven by land use and environmental factors (climate, atmospheric CO 2 ) as simulated by the ORCHIDEE LSM. The flexible structure of CE-DYNAM makes the model adaptable to the SOC dynamics of other LSMs. In this way it is possible to study the main processes behind the linkages between soil erosion and the global C cycle.
CE-DYNAM explicitly accounts for hillslope and floodplains re-deposition, which is to our knowledge unique for a large-scale C erosion model and highly novel. However, it still lacks important processes affecting the C dynamics in floodplains. The model does not account for a slower respiration rate due to low-oxygen conditions, physical and chemical stabilization (Berhe et al., 2008;Martínez-mena et al., 2019) or a higher NPP for nutrient-rich floodplains (Van Oost et al., 2012;Hoffmann et al., 2013). The oxidation and preservation of C in deposition environments, especially in alluvial reservoirs remain highly uncertain (Billings et al., 2019). Due to its simplistic nature and coarse-resolution, CE-DYNAM does not resolve rivers and streams explicitly but assumes that they are included in the floodplain part of the grid cells. As a result, CE-DYNAM does not differentiate between eroded hillslope soil that reaches the water network directly (where the residence time of suspended sediment is in the order of days), and the sediment that is first retained in the floodplains before it reaches the water network due to fluvial erosion (sediment residence time is in the order of a few to thousands of years). CE-DYNAM has been developed and calibrated to simulate long-term changes in sediment and C storage on land and not the short-term variations in sediment and POC fluxes carried by rivers. This limits the application of CE-DYNAM in its current form to accurately quantify sediment and POC fluxes of rivers and streams, and to compare them to observations. This sediment export rate leads to a yearly sediment bound POC export of about 2x10 8 g C year -1 2005. This POC flux is also two orders of magnitude lower than the 2.6x10 10 g C year -1 given by the GlobalNEWS2 model (Mayorga et al., 2010) or the 1.5x10 11 g C year -1 found by Beusen et al. (2005), which is mainly a result of the underestimated simulated sediment export rate.
Furthermore, CE-DYNAM does not simulate fluvial erosion as a complex function of the channel geometry, riverbank erodibility and sheer stress (Dröge et al., 1992), due to the lack of data on these parameters at the regional scale, and to keep a balance between model complexity and its computational ability. Also, our model does not resolve erosion of the deposited river sediment by flooding events. This simplified model concept for fluvial erosion contributes to the underestimation of sediment and C export in floodplains. Finally, with the current model setup we do not account for large soil erosion events before 1850 or extreme precipitation events that may have a long-term effect on the sediment export rate of the Rhine.
Although we underestimate the riverine sediment and POC fluxes, we find that the spatial variability in sediment storage and SOC stocks of the sub-basins are within or close to observational uncertainty ranges (Table 5, 6;Naipal et al., 2016).
We also find that the C density in the topsoil layers of floodplain soils located downstream of the Rhine and the C concentration of the POC flux are realistic. We find a C concentration of ~3.3 % in the exported fine sediments. Abril et al. significantly to the higher SOC stocks of floodplains compared to hillslopes, and to the export of DOC and POC to rivers (Van Oost et al., 2012;Hoffmann et al., 2013a).
In a future study we aim to improve the sediment and POC export, and account for a higher floodplain plant productivity by using a nutrient-enabled version of the ORCHIDEE LSM (Goll et al., 2017). Furthermore, the model does not take into account the full effects of the selectivity of erosion, often expressed as the enrichment ratio, where the C content of eroding soil or the deposited sediment can be different from that of the original soil. The enrichment varies substantially across landscapes, while the importance of erosion selectivity for C is still under debate (Nadeu et al., 2015;Wang et al., 2010). However, we did a simple sensitivity test to study the effect of C enrichment by erosion (section 4.3).
CE-DYNAM does not account for different ratios between the SOC pools (active, slow, passive) with depth due to the limitation in information to constrain these fractions for floodplains and hillslopes. However, this can be potentially important for respiration of C in depositional sites and during transport. Studies show that the labile C is decomposed first during sediment transport and directly after deposition, leaving behind the more recalcitrant C in deposition sites (Berhe et al., 2007;Billings et al., 2019). Due to the simplistic nature of our coarse-resolution model and the lack of data on oxidation of eroded C during transport we did not include C respiration during transport in the model.
The current SOC scheme of CE-DYNAM does also not account for different residence times of SOC as a function of landscape position along a hillslope. The SOC decomposition rates can vary significantly along a hillslope due to changes in soil moisture, temperature, aggregation, and the transport of minerals and nutrients (Doetterl et al., 2016). Currently, these processes are not resolved in coarse resolution LSMs, contributing to the uncertainty in the large-scale linkage between soil erosion and SOC dynamics.
Furthermore, there is no feedback between soil erosion and plant productivity in the model. To account for this feedback, soil erosion processes would need to be explicitly included in a LSM, such as ORCHIDEE, which would increase the computational complexity of the simulations substantially. The lack of this feedback results in an unlimited dynamic replacement of C on eroding sites.
Currently, the erosion scheme of CE-DYNAM does not include the L (slope-length) and P (support-practice) factors. This might induce some bias in the results, especially for agricultural land. In a future study we aim to make CE-DYNAM better applicable for agricultural land, where these factors play an important role. For this purpose we will focus on the development of new methods that can quantify the L and P factors reliably at the global scale, and will need to re-calibrate the Adj.RUSLE model. Our decision of leaving out the L and P factors from the erosion equation in this study is based on the global study of Doetterl et al. (2012), which showed that the S, R, C and K factors explain approximately 78 % of the total erosion rates on cropland in the USA. This indicates that on cropland the L and P factors, which are related to agriculture and land management, contribute only for 22 % to the overall erosion rates. This percentage is comparable to the uncertainty range in the estimation of the S, R, C and K factors at the regional scale from coarse resolution data. Renard and Ferreira (1993) also mention that the soil loss estimates are less sensitive to slope-length than to most other factors.
Furthermore, various studies argue that the estimation of the L factor for large areas is complicated and thus can induce significant uncertainty in soil erosion rates calculated based on coarse resolution data (Foster et al., 2002;Kinnell, 2007).
Especially, for natural landscapes, such as forest, the estimation of the L factor is not straightforward as these natural landscapes usually include steep slopes (Elliot, 2004). In order to stay consistent with the estimation of potential soil erosion for all land cover types, we removed the L factor from the equation. The Adj.RUSLE has been already successfully validated at the regional scale, without the L and P factors where the spatial variability of soil erosion rates compares well to other high resolution modeling studies and observational data and the absolute values fall within the uncertainty ranges of those validation data (Naipal et al., 2015;Naipal et al., 2016;Naipal et al., 2018;and this study). Finally, the aim of this study was to develop and validate a C erosion scheme for applications at the global scale, where the estimation of the L and P factors is limited. By showing that the erosion rates from the Adj.RUSLE and CE-DYNAM are within the uncertainty of other data and modelling studies, we assume that it will be applicable for other large catchments in the temperate region.
Finally, CE-DYNAM considers only the rather 'slow' rill and interrill soil erosion processes, and does not take into account severe erosion processes such as gully erosion and landslides, which are bound to extreme precipitation events.
The daily timestep of CE-DYNAM and the current setup of the sediment budget module allows only for long-term yearly average changes in erosion and deposition rates and cannot be applied to estimate episodic erosion and deposition events.

Sensitivity analysis
We analyzed the effects of the following model assumptions: (1) C enrichment during erosion, (2) the floodplain sediment residence time, and (3) crop residue management.
To test the C enrichment we increased the EF (Eq. 15) from 1 to 2, assuming a strong enrichment of C during erosion (section 2.11). We find that this enrichment results in a gross C erosion flux that is 1.61 times larger than the flux without enrichment (Table 7). This leads also to a larger dynamic replacement of C on eroding sites in combination with a larger burial in depositional sites, which is in accordance with the study of Lugato et al. (2018). The resulting C sink from the enrichment simulation is 1.25 times larger than the sink under default conditions (Table 7).
To test the potential effects of a different sediment residence time on the SOC dynamics, we performed a sensitivity study where we changed the basin average sediment residence time to be 50% higher or 50% lower but keeping the maximum sediment residence time at 1500 years (section 2.11). By changing the average sediment residence time and keeping the maximum fixed, the grid cells with the lowest residence times will undergo the largest changes in the residence time and consequently in the floodplain SOC storage and export. The higher the residence time, the longer the deposited soil C will reside in the floodplains, where it can either be respired or buried in deeper soil layers. Therefore, we find that the effects of the sediment residence time on the SOC dynamics are non-linear. Under default conditions we find the highest SOC storage. A 50 % higher average sediment residence time leads to the lowest total SOC storage, with a decrease of 30 % compared to default conditions, while the erosional C sink is reduced by 20 % (Table 7). This could be explained by a higher C decomposition flux for floodplains due to the long residence time of C in deposition areas. Especially in mountainous regions where the soil erosion flux is large and removes a large part of the labile C, a higher sediment residence time will lead to higher C emissions due to decomposition in floodplains. The turnover seems to dominate over the C burial in deeper layers and export. A 50 % lower average sediment residence time also leads to a decrease (of 8 %) in the total SOC storage and a decrease of 6 % in the erosional C sink compared to default conditions (Table 7). Also here, the largest changes are found in mountainous regions where a low sediment residence time leads to a large export of C, which is then deposited in lower lying, more extensive floodplains. Thus, increasing or decreasing the residence time leads to a smaller total SOC storage, resulting from different spatial distributions of this SOC storage. The POC flux under the high sediment residence time scenario is substantially higher than under default conditions (Table 7).
To test the effects of crop residue management we harvested all above-ground crop residues (section 2.11). We find that the total litter C stock is about 15 % smaller than the default case by the end of the year 2005. This leads to a total change in the transient SOC stocks that is 20 % smaller under no erosion (S0), and 26 % smaller under erosion (S2) ( Table 7). Our findings confirm that soil management practices such as residue management have a substantial effect on the SOC dynamics.

Conclusions
We presented a novel spatially-explicit and process-based C erosion dynamics model, CE-DYNAM, which simulates the redistribution of soil and C over land as a result of water erosion and estimates the implications for C budgets at catchment scale. We demonstrated that CE-DYNAM captures the spatial variability in soil erosion, C erosion and SOC stocks of the non-Alpine region of the Rhine catchment when compared to high-resolution estimates and observations. We also showed that the quantile ranges of erosion and deposition rates and C stocks fall within the uncertainty ranges of previous estimates at basin or sub-basin level. Furthermore, we demonstrated the model ability to disentangle vertical C fluxes resulting from the redistribution of C over land and develop C budgets that shed light on the role of erosion in the C cycle. The simple structure of CE-DYNAM and the relative low amount of parameters make it possible to run several simulations to investigate the role of individual processes on the C cycle such as the removal by erosion only, or the role of deposition and transport. Its compatibility with land surface models makes it possible to investigate the long-term and large-scale effect of erosion processes under various global changes such as increasing atmospheric CO 2 concentrations, changes to precipitation and temperature, and land use change. The application of CE-DYNAM for the Rhine catchment for the period 1850-2005 AD reveals three key findings:

25
• Soil erosion leads to a cumulative net C sink of 216±23 Tg by the end of the period, which is in the same order of magnitude as the cumulative land C sink of the Rhine without erosion. This C sink is a result of an increasing dynamic replacement of C on eroding sites due to the CO 2 fertilization effect, despite decreasing soil and C erosion rates over the largest part of the catchment. We conclude that it is important to take into account global changes such as climate change in order to better quantify the net effect of erosion on the C cycle.
• After performing a sensitivity analysis on key model parameters we find that the C enrichment by erosion, crop residue management and the residence time of floodplain sediment can substantially change the overall values of C fluxes and SOC storages. However, the main findings, such as soil erosion being a net C sink for the Rhine catchment, remain.
• Initial climate and land cover conditions and the transient period over which erosion under global changes takes place are essential for determining if soil erosion is a net C sink or source and to what extent.
Altogether, these results indicate that despite model uncertainties related to the relative coarse spatial resolution, missing or simplified processes, CE-DYNAM represents an important step forwards into integrating soil erosion processes and sediment dynamics in Earth system models. The next step would be to improve CE-DYNAM with respect to riverine sediment and POC export fluxes and management practices.

Code and data availability
The source code of CE-DYNAM is included as a supplement to this paper. Model data can be accessed from the Zenodo repository under the doi:10.5281/zenodo.2642452 (not published yet). For the other data sets that are listed in Table 1, it is encouraged to contact the first authors of the original references.

Author contributions
VN built and implemented the mode. YW provided the basic structure for the model. All authors contributed in the interpretation of the results and wrote the paper.

Competing interests
The authors declare that they have no conflict of interest.    Table 4: Goodness-of-fit results of the comparison of the simulated gross and net C erosion rates to those of the study of Lugato et al. (2018) in the enhanced and reduced scenario, taking into account 13 sub-basins of the Rhine. RMSE is the root mean square error in tons year -1 . Ce stands for gross C erosion, while Cd stands for net C erosion.  and those of the Global Soil dataset for Earth System Modeling (GSDE) and from the LUCAS database. The regression is done after aggregating the data at sub-basin level for the 13 sub-basins that were delineated in the Rhine catchment.
RMSE is the root mean square error given in Tg of C per year, while the r-value is the spatial correlation coefficient.