Modelling long-term, large-scale sediment storage using a simple sediment budget approach

. Currently, the anthropogenic perturbation of the biogeochemical cycles remains unquantiﬁed due to the poor representation of lateral ﬂuxes of carbon and nutrients in Earth System Models (ESMs). This lateral transport of carbon and nutrients between terrestrial ecosystems is strongly affected by accelerated soil erosion rates. However, the quantiﬁcation of global soil erosion by rainfall and runoff, 5 and the resulting redistribution is missing. This study aims at developing new tools and methods to estimate global soil erosion and redistribution by presenting and evaluating a new large-scale coarse resolution sediment budget model that is compatible with ESMs. This model can simulate spatial patterns and long-term trends of soil redistribution in ﬂoodplains and on hillslopes, resulting from external forces such as climate and land use change. We applied the model on the Rhine catchment 10 using climate and land cover data from the Max Planck Institute Earth System Model (MPI-ESM) for the last millennium (850 - 2005 AD). Validation is done using observed Holocene sediment storage data and observed scaling between sediment storage and catchment area. We ﬁnd that the model reproduces the spatial distribution of ﬂoodplain sediment storage and the scaling behavior for ﬂoodplains and hillslopes as found in observations. After analyzing the dependence of the scaling 15 behavior on the main parameters of the model, we argue that the scaling is an emergent feature of the model and mainly dependent


Introduction
Soil erosion by rainfall and the resulting soil redistribution in a landscape play an important role in the cycling of soil carbon and nutrients in ecosystems (Van Oost et al., 2007).On the one hand, vertical fluxes of carbon and nutrients occur due to either mineralization on eroded landscapes and during sediment transport, or due to sequestration in depositional sites (Lal, 2003;Van Oost et al., 2007;Lal, 2005;Quinton et al., 2010).On the other hand, significant lateral fluxes of soil carbon and nutrients can take place when soil redistribution promotes the lateral transport of these elements in and between terrestrial ecosystems (Van Oost et al., 2007;Quinton et al., 2010).
Recent evidence demonstrated that human activities, such as land use change, have accelerated soil erosion rates globally (Van Oost et al., 2012;Wall and Six, 2015).Accelerated soil erosion has the potential to alter not only the vertical fluxes of carbon and nutrients, but also the lateral export of these elements from terrestrial ecosystems to the coastal oceans (Regnier et al., 2013;Stallard;Bauer et al., 2013;Le Quéré et al., 2013).Based on the various effects of soil erosion and redistribution on the carbon and nutrient cycles, these processes can either result in accelerated soil erosion to be a net source (Lal et al., 2004), or in a net uptake or sink of CO 2 (Stallard; Van Oost et al., 2007).
However, data on global soil erosion and redistribution are scarce to non-existing.There exist several modelling approaches to estimate global soil erosion rates (Yang et al., 2003;Ito, 2007;Montgomery, 2007;Doetterl et al., 2012;Naipal et al., 2015).These modelling approaches mainly address the soil detachment process only, and do not simulate soil redistribution by ignoring processes such as sediment deposition and transport.There is, to our knowledge, no spatially explicit model that can simulate soil redistribution on large spatial scales, for the past, present and future.The lack of such kind of large-scale models on soil redistribution substantially limits the understanding of the interaction of soil erosion and redistribution with the global biogeochemical cycles.Therefore, the net global effect of accelerated soil erosion on the vertical and lateral fluxes of soil carbon and nutrients is still unknown.
Consequently, the land components of Earth System Models (ESMs), which are the main tools to investigate the terrestrial carbon cycle and the carbon flux between soil and the atmosphere, ignore the lateral carbon fluxes associated with soil redistribution (Regnier et al., 2013;Van Oost et al., 2012).Therefore, they miss an important aspect of the coupling between land and the ocean.In addition, omitting soil erosion from soil organic carbon cycling schemes results in uncertainties in the soil organic carbon flux with various implications (Chappell et al., 2015).Including soil redistribution processes in ESMs is thus essential to create the possibility to study the full interactions and feedbacks between the soil and atmosphere with respect to the global biogeochemcal cycles, and to better understand the anthropogenic perturbation of these cycles.
The holistic understanding of the interaction and linkages between soil erosion, deposition and transport, can be addressed using the sediment budget approach (Walling and a.L. Collins, 2008).Slaymaker (2003) defined the sediment budget as a mass-balance where the mass of sediments is conserved in the considered system so that the net increase in sediment storage is equal to the excess of inflow over outflow of sediment.However, long-term large-scale sediment budgets are very scarce to non existing.Sediment budgets that have been constructed previously range from small catchments (Verstraeten and Poesen, 2000;Walling et al., 2001) to large river catchments (Milliman and Meade, 1983;Ludwig and Probst, 1998;Syvitski et al., 2003;Slaymaker, 2003).However, these sediment budgets are usually for present day only as they are mostly based on measurements using methods such as sediment tracing or fingerprinting.Also, most of these studies only focus on the sediment delivery from a catchment.They are thus of limited use for assessing the spatial distribution of sediment sources and storage or in predicting long-term sediment yields.Considering explicitly the spatial distribution of these variables within a catchment is not only essential for a proper land management strategy to combat land degradation, but also for a detailed assessment of how soil erosion and redistribution interact with the carbon and nutrient cycles.Furthermore, it is also essential to consider long-term sediment budgets, as they can provide essential information on the forces behind sediment, carbon and nutrient fluxes in a catchment such as human activities and climate change.
There exist different spatial models of suspended sediment flux that also consider the soil redistribution or sediment dynamics in a catchment (Merritt et al., 2003;de Vente and Poesen, 2005;Ward et al., 2009).However, many of them are developed to simulate single events or require input data that is not available for large spatial scales (Wilkinson et al., 2009).There are also partly empirical models which can operate on catchment scale such as the WATEM/SEDEM model, which is used to predict hillslope sediment storage and sediment yields (de Moor and Verstraeten, 2008;Nadeu et al., 2015).Or such as the suspended sediment model from Wilkinson et al. (2009) that also simulates some other processes such as floodplain deposition, gully and riverbank erosion.However, these models are not compatible for a global scale application as they require parameters for which data is not available on a global scale and these type of models also need to be calibrated to measured sediment yields of the studied area (Van Rompaey et al., 2001).Pelletier (2012) proposed a global applicable model for long-term suspended sediment discharge, where he used various environmental controlling parameters to simulate soil detachment and sediment transport.However, in his study he mainly focuses on the sediment discharge and delivery of catchments and his model does not take into account the full dynamics of sediment in a catchment, which would also include the spatial distribution of sediment deposition and storage in the different reservoirs of a catchment.Additionally, he does not consider land use change and thus his approach is limited to natural catchments only.
The overall aim of this study is to make a first step in quantifying large-scale soil erosion and redistribution rates and identifying their drivers, in order to contribute to the representation of sediment dynamics and associated lateral fluxes of carbon and nutrients in ESMs.Therefore, we present and evaluate a new large-scale sediment budget model for the non-Alpine part of the Rhine catchment using the environment of ESMs.We use the model to quantify the spatial variability of floodplain and hillslope sediment storage for the Rhine catchment, and its dependence on climate change and land use change during the last millennium (850-2005 AD).We also investigate the relationship between catchment area and sediment storage on hillslopes or in floodplains to derive a general validation test for our large-scale model.Finally, we discuss the main challenges in modelling large-scale, long-term soil erosion and redistribution, and future perspectives for application in ESMs.

Basic model concept
The main purpose of the sediment budget model presented here, is to estimate large-scale longterm floodplain and hillslope sediment storage and lateral fluxes of sediment.The model should, therefore, be spatially explicit and capable of estimating erosion, deposition and sediment transport processes.For this purpose we will use a grid cell based approach.Compatibility of this new model with ESMs is important for a future extension of the model to include the carbon and nutrient cycling.Furthermore, it is essential to distinguish between floodplain and hillslope systems due to the distinct differences in sediment dynamics between these systems (Hoffmann et al., 2013).Human activities usually lead to a stronger increase in sediment deposits on hillslopes compared to floodplains, and an overall decreased export of sediment out of a catchment, despite increased soil erosion (de Moor and Verstraeten, 2008).In this way, sediments stored in floodplains and on hillslopes over long timescales can significantly delay or alter the human induced changes to the carbon and nutrient cycles (Hoffmann et al., 2013).
Before we can define a model that satisfies the above mentioned conditions we have to make some basic assumptions.First, as it is difficult to disentangle the floodplains and hillslopes in available soil data sets, we assume that each grid cell contains both a hillslope and a floodplain reservoir.When estimating large-scale sediment storage with the aim of predicting the effects of soil redistribution on the biogeochemical cycles, the focus is to get the large-scale spatial patterns right, rather than accurate numbers for the sediment storage and fluxes.Second, we assume that the sediment deposition and transport behave differently between the floodplain and hillslope reservoirs on the timescale of the last millennium.Third, erosion is considered to mainly take place on hillslopes, where part of the eroded sediment is directly transported from hillslopes and deposited in the floodplains.
The underlying model framework (Fig. 1a) that consists out of the erosion, deposition and sediment transport modules, is based on the sediment mass-balance method.The change in sediment storage (M ) within a certain unit of time and space is given by the difference between sediment input and sediment output (Slaymaker, 2003).For sediment stored in floodplains (M a ), this leads to: Here, D a is the sediment deposition rate in floodplains, and L is the sediment loss.Equation 1 can be seen as a representation of the net soil redistribution flux, and is approximated by the following as function of time: where D a (t) is the time-dependent input rate in the model, which is independent from M a (t).
k * M a (t) is the loss term of the floodplain reservoir, and k is the specific rate for floodplains.
The specific rate is the inverse of the residence time (1/τ ) for floodplain sediment.τ is defined as the time (in years) a soil particle stays in the floodplain reservoir of a certain grid cell.Here, we assume that τ is independent of time for timescales in the order of a few thousands of years.We also assume that τ is increasing exponentially with catchment area or weighted flow-accumulation (F lowAcc): a τ and b τ are the adjustment parameters of the model that relate the residence time to F lowAcc.
F lowAcc is defined as the number of grid cells located upstream that flow into a certain grid cell.As each grid cell represents a certain catchment area, the value of τ will be dependent on the location of the grid cell in the catchment.The presented relationship between τ and catchment area in equation 3 is based on the fact that large catchment areas are usually characterized by low slopes, which mainly result in a low connectivity making the system capable of storing sediment for long time periods.
The opposite is true for small catchment areas, where the connectivity is usually high, resulting in short residence times for sediment (Hoffmann, 2015).
The deposition rate (D a ) in the floodplain reservoir can be defined as a certain fraction of the erosion rate (E).In this way equation 2 can be rewritten as: Where f is the dimensionless floodplain deposition fraction ranging between 0 and 1.
) is computed according to the adjusted Revised Universal Soil Loss Equation (RUSLE) model (Naipal et al., 2015), which computes annual averaged rill and interril erosion rates and is formulated as a product of a rainfall erosivity factor (R, MJ mm ha −1 h −1 yr −1 ), a slope steepness factor (S,dimensionless), a soil erodibility factor (K, t ha h ha −1 MJ −1 mm −1 ), and a land cover factor (C,dimensionless): The underlying RUSLE model stems from the original Universal Soil Loss Equation (USLE) model developed by USDA (USA Department of Agriculture), which is based on a large set of experiments on soil loss due to water erosion from agricultural plots in the United States (Renard et al., 1997).These experiments covered a large variety of agricultural practices, soil types and climate conditions, making it a potentially suitable tool on a regional to global scale.Although, RUSLE was originally developed for agricultural land, model parameters for other land cover types such as forest and grassland have also been estimated using observational data (Dissmeyer and Foster, 1981;Millward and Mersey, 1999;Lu et al., 2004).
In the adjusted RUSLE model, as presented above, the effects of the slope-length (L factor) and support practice (P factor) are excluded.In the original RUSLE model (Renard et al., 1997), these factors are part of the model, however, on large to global scale there is too little data available on these factors.Including them in the model would only result in additional uncertainties, while we try to keep the model simple, to be able to capture and quantify the main processes and drivers behind large-scale sediment mobilization.We do, however, agree that leaving these two factors out could introduce some biases in erosion rates, especially in agricultural areas.
f is calculated by a simple growth function where deposition is a function of the average percent topographical slope (θ) and the main land cover type in a grid cell: where a f and b f are adjustment parameters that relate f to the average slope depending on the land cover type.θ max is the maximum percent slope.According to equation 6, an increase in the overall average slope of a grid cell leads to a larger transport of eroded soil from the hillslopes to the floodplains, leading to an increased deposition rate to the floodplain reservoir of that specific grid cell.Hereby, we assume the increase of f to be exponential.
The effect of the land cover type on f in our model represents mainly the interaction of the landscape connectivity with sediment transport.The connectivity of a natural landscape, consisting out of mainly forest, is largely based on the vegetation density and morphological structures (Gumiere et al., 2011;Bracken and Croke, 2007).In crop-and grassland, however, the landscape connectivity is strongly affected by anthropogenic structures.Several recent studies (Hoffmann et al., 2013;de Moor and Verstraeten, 2008;Gumiere et al., 2011) show that these anthropogenic structures and activities reduce the sediment transport from hillslopes to the floodplains.In this way, the stored hillslope sediment is disconnected from the fluvial system on timescales of 100 to a few 1000 years.
Based on this, we assume in our model that for crop-and grassland the sediment connectivity is disturbed.A disturbed sediment connectivity will result in a larger fraction of eroded soil that remains on the hillslopes compared to the fraction that flows along the hillslopes and is deposited in the floodplains.For natural landscapes we assume a better sediment connectivity, meaning that an equal or larger fraction of the eroded soil will be deposited in the floodplains compared to the fraction that remains on the hillslope.Here we ignore morphological conditions that can cause deconnectivity in the landscape.
After calculating erosion and deposition, the sediment is transported between grid cells based on the multiple flow sediment routing scheme such as presented by Quinn et al. (1991) (Fig. 1b).In the multiple flow routing scheme the weight (W , dimensionless), which specifies the part of the flow that comes in from a neighboring grid cell, is calculated as: where c is the contour length and is 0.5 in the cardinal direction and 0.354 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 cells where k and l can be either -1, 0 or 1. θ is calculated here as: where, h is the elevation in meters derived from a digital elevation model, d is the grid size in meters.
The floodplain sediment storage rate (t ha −1 year −1 ) of a grid cell (i, j) is then a function of the deposition rate in that grid cell, the loss from that grid cell and the incoming sediment from the neighboring grid cells, and is calculated at each time step t as: For hillslopes the change in sediment storage is assumed to be equal to the input rate (Eq.10), because we assume that the stored hillslope sediment has an infinite residence time on the timescale of the last millennium in accordance with the study of Hoffmann (2015).This means that the hillslope sediment storage will increase linearly with time (Eq.11).The hillslope sediment deposition rate (D c ) is here defined as the remaining part of the eroded soil that has not be been transferred to the floodplain directly (1-f ).The equations for the hillslope sediment storage rate (M c , t ha −1 year −1 ) are represented by: And as a function of time: The modelling approach as presented by the equations above focuses on the net soil redistribution by separately modelling the main processes of soil redistribution, which are erosion, deposition and transport.In the following paragraphs we will show how this dynamical modelling approach performs when applied on the Rhine catchment.

Model implementation and parameter estimation
The resolution of the sediment budget model is 5 arcmin.The main reason for choosing this particular model resolution is based on the assumption that this resolution is optimal when considering that each grid cell contains a floodplain and hillslope fraction.Here, a higher resolution could lead to cases where this assumption is not met.Also, the 5 arcmin resolution fits well with the resolution of the adjusted RUSLE model.
The sediment budget model uses climate and land cover data from simulations of MPI-ESM that have been performed under the Coupled Model Intercomparison Project Phase 5 (CMIP5) framework (Hurrell and Visbeck, 2011;Taylor et al., 2009).As this data was given at a resolution of approximately 1.875 degrees, we had to downscale the data to the resolution of the sediment budget model.For the period 1850-2005 AD three ensemble members from MPI-ESM (r1i1p1, r2i1p1, r3i1p1) were available, while for the period 850-1850 AD only one ensemble member (r1i1p1) was available.This data existed on a 6 hourly, monthly or yearly time step for the last millennium.
Calculation of soil erosion according to the adjusted RUSLE model is mostly based on the methods presented in the study of Naipal et al. (2015).However, the calculation of the R and C factors had to be adapted due to the very coarse resolution of the data from MPI-ESM or the lack of data on certain parameters of the model.A detailed description of erosion estimation with the adjusted RUSLE model in combination with data from the MPI-ESM model is presented in the supply material.
Additionally, due to the overestimation of erosion rates by the adjusted RUSLE model in the Alps, we defined a mean soil erosion rate of 20 t ha −1 year −1 for this region based on high resolution erosion data from Bosco et al. (2008).
We chose f to range between 0.5 and 0.8 for forest, and between 0.2 and 0.5 for crop-and grassland.
These numbers are based on findings from the study of de Moor and Verstraeten (2008), where they show an approximately equal deposition rate in floodplains as on hillslopes before agricultural activities started in the Geul river catchment in The Netherlands.However, for present-day they show that much more sediment is trapped on hillslopes than is transferred to the floodplains.Based on the chosen ranges for f and equation 6 we calculated a f and b f for forest to be 0.5 and 0.47, respectively, and for crop-and grassland to be 0.2 and 0.917, respectively.This means that for low slopes (<±0.2 %) in a forest an equal amount of sediment is deposited in floodplains as on hillslopes, while for crop-and grassland only 20 % of the eroded soil from the hillslopes will reach the floodplains.
The floodplain residence time is made to range between the median and maximum residence time of floodplain sediment in the Rhine catchment of 260 and 1500 years, respectively.This is in ac-cordance with the residence times derived from observed sediment storage in the Rhine catchment.a long land use history, where agricultural activities started about 7500 years ago (Houben et al., 2006;Hoffmann et al., 2007;Dix et al., 2016;Kalis and Merkt, 2003), the Rhine catchment serves as a good case study to investigate the impact of human activities on erosion and sediment yields through history.The Rhine catchment (Fig. 2) has a size of ∼ 185000 km 2 with a main river channel length of ∼ 1320 km and drains large parts of the area between the European Alps and the north sea.It has a complex topography where the elevation ranges between -180 and 1967 m with a mean topographical percent slope of 0.07, where percent slopes can go up to 1.2.It consists out of two large sedimentary catchments (ie, upper Rhine Graben and the lower Rhine Embayment-Southern North Sea Basin) that serve as large floodplain sinks for sediment, and some upland areas, such as the Black Forest and the European Alps that serve as major sediment production areas.
From observed Holocene sediment storage Hoffmann et al. (2013) derived scaling relationships between sediment storage S (10 9 kg = 1 Mt) and catchment area A (km 2 ) for floodplains and hillslopes.They found that for floodplains the sediment storage increases in a non-linear way with catchment area, while for hillslopes this increase is linear.With these scaling relationships, for the first time, a direct comparison is made between the behavior of soil redistribution on hillslopes and in floodplains at large spatial scales.This is an essential difference between hillslopes and floodplains that large-scale sediment budget models like ours need to capture in order to reliably simulate the spatial distribution of sediment on such a scale.The scaling relationships, given by equation 12for hillslopes and equation 13 for floodplains, will be used as a simple validation test for our coarse resolution sediment budget model.
Here, A ref is an arbitrary chosen reference area, in this case 10 3 km 2 .The observation data contains 41 hillslope and 36 floodplain sediment storage values, derived from a large number of auger and bore holes that are used to measure sediment thickness related to human induced soil erosion.
With the estimated scaling exponents (Eq.12 and Eq.13) Hoffmann et al. (2013) showed that even for large catchments (in the order of 10 5 km 2 ) hillslopes store an equal amount of sediment as floodplains.They pointed out that this is a substantial sink that needs to be considered in sediment budgets of large catchments.tons was found for the whole Rhine catchment, of which 50 % is stored in the Rhine Graben and the delta.The spatial variability of this observed sediment storage in floodplains will be a second validation test for our model.
Finally, Hoffmann et al. (2007) also found an average erosion rate of 0.55 ± 0.16 t ha −1 year −1 for the last 10000 years, with extreme minimum and maximum values of 0.3 and 2.9 t ha −1 year −1 .
However, Hoffmann et al. (2013) included also hillslope sediment storage and calculated a total sediment storage of 126 ± 41 Gt for the Rhine catchment, which requires a minimum Holocene erosion rate of approximately 1.2 ± 0.32 t ha −1 year −1 .This shows that hillslopes are not only the main sources of eroded sediment but can be major millennial-scale sinks for eroded sediment that comes from agriculture.We will use the average erosion erosion rates from the above mentioned studies as a comparison to the rates derived from our sediment budget model.

Simulation setup
In order to simulate sediment storage for a certain catchment, an initial state of that catchment has to be assumed.Here we assume the initial state to be the equilibrium state of a catchment, defined as the state of a catchment where the sediment input is equal to the sediment output, and thus the sediment yield at the outlet of the river is constant in time.External forces working on a catchment such as, land use activities or deglaciation, can bring the catchment out of equilibrium into a transient state.
In the case of the Rhine catchment the period directly after the Last Glaciation Maximum (LGM) could be of major importance due to strong erosion that was triggered by the retreating ice sheets.
From today's observations on sediment yields or erosion rates we cannot determine when the Rhine catchment was in an equilibrium state.Additionally, there are no observations of sediment storage before the start of agricultural activities in the Rhine catchment.This poses a problem in simulating and interpreting present-day absolute values of sediment storage and yields with our sediment budget model.
In order to still being able to interpret the simulated results for the Rhine catchment, we will only focus on the change in sediment storage due to land use and climate change since 850 AD.Consid-ering mainly the changes induced by external forcing, it is not necessary to know if the system was in an equilibrium or transient state at 850 AD.Based on this reasoning, we use the environmental conditions of the period between 850 and 950 AD to determine the equilibrium state of the model.
In the rest of this study, we will refer to the period between 850 and 950 AD as the 'default equilibrium state' that we define based on the mean environmental conditions between 850 and 950 AD, while one should keep in mind that this is not the 'real' equilibrium state of the catchment.The period 850 to 950 AD is used here as the equilibrium state due to reasons related to data availability, and because human impact in this time period is still small compared to present day.
Hence, our simulation setup structure is generally defined by an equilibrium simulation based on the mean climate and land cover conditions of the period between 850 and 950 AD, followed by a transient simulation for the last millennium.
We performed three equilibrium simulations, one based on the mean climate and land cover conditions of the period 850 to 950 AD, and the two others based on the mean climate and land cover conditions of the mid-Holocene period (6000 years ago) from the mid-Holocene experiment of the MPI-ESM (Table 1).The reason for performing an equilibrium simulation for the mid-Holocene period is to investigate how different initial conditions for climate and land cover influence the overall sediment storage change during the last millennium.
In the equilibrium simulations the erosion and deposition rates are kept constant and the model is the equilibrium state of 850 -950 AD (Table 1).In the climate change simulation the land cover was fixed to the mean conditions of 850 -950 AD during the whole period of the last millennium, while the climate was variable.In the land use change simulation the climate was fixed to the mean conditions of 850 -950 AD, while the land cover was variable (Table 1).
3 Application of the sediment budget model

Scaling test
In order to validate the sediment budget model we tested if the model can reproduce the scaling relationships found by Hoffmann et al. (2013) for the non-Alpine part of the Rhine catchment (Eq.12 and 13).For this we chose the grid cells in the Rhine catchment that correspond to the observation points from Hoffmann et al. (2013).Observation points that fell outside the Rhine catchment, were not considered.When considering only the selected grid cells and applying the same scaling approach as in the study of Hoffmann et al. (2013), we find average scaling exponents of 1.2 ± 0.04 and 1.05 ± 0.07 for floodplains and hillslopes, respectively (Table 2).These values fall in the range of floodplain and hillslope scaling exponents of 1.23 ± 0.06 and 1.08 ± 0.07, respectively, found by Hoffmann et al. (2013).The uncertainty in the scaling exponents is mainly due to the regression, while the uncertainty due to different ensemble simulations is very small (Table 2).These results indicate that our model reproduces the characteristic differences in scaling between floodplains and hillslopes as found by Hoffmann et al. ( 2013) (Fig. 3a and b).One should note that the grid resolution of the model limits the prediction of sediment storage to grid points with a catchment area When considering all the grid cells of the Rhine catchment we find a scaling exponent for floodplain storage of 1.33±0.02(Table 3).This is somewhat higher than the value found when only the selected grid cells are used, which can be explained by the inclusion of grid cells located in the Alpine region of the Rhine catchment.Including the Alpine region leads then to a stronger gradient in sediment storage and catchment area between the Alps and the Rhine delta.In the Alpine region the model predicts much less sediment storage due to the low residence time and high sediment connectivity, while for the Rhine delta the sediment storage is large due to high residence times and low sediment connectivity.For hillslope storage the scaling exponent is also slightly higher when including all grid cells in the scaling approach (Table 3).This can also be explained by including the Alpine region, where the model predicts more sediment storage on hillslopes in contrast to the rest of the Rhine catchment as a result of high erosion rates.
Furthermore, when including all grid cells in the scaling approach the spread in the data increases, which is clear from the lower r-value of the regression.The small difference between the scaling exponents when considering all grid cells and the scaling exponents when considering only selected grid cells indicates that the selected observation points from Hoffmann et al. ( 2013) are robust and representative for the catchment.The relatively small difference can be partly attributed to biases in simulated erosion and deposition rates and the floodplain residence times.
Finally, we find that keeping either the climate or land cover constant throughout the last millennium has very little impact on the scaling exponent for floodplain storage.Here, the climate change simulation results in a slightly higher and the land use change simulation in a slightly lower scaling exponent.The different forcings have a stronger impact on the scaling for hillslope storage, as hillslope storage is only dependent on erosion and deposition rates.In the climate change simulation the scaling exponent for hillslope storage increases by 3.8 %, while in the land use change simulation a small decrease of 0.1 % is found.This decrease can result from the fact that most land use change took place in the lower parts of the Rhine catchment resulting in an increased sediment storage there.
In contrast, the land use conditions in the Alpine region did not change that rapidly, resulting in a decreased difference in sediment storage on hillslopes between the upper and lower areas of the catchment.
With the above results we show that the scaling relationships are a general feature for the entire Rhine catchment and are independent of the selected observation points.As the Rhine catchment is a large catchment with a complex topography, this indicates that the scaling relationships might also be applicable for other large river catchments.

Origin of scaling between sediment storage and catchment area
We also performed a sensitivity study to test the robustness of the scaling relationships derived with the model.For this we investigated the dependence of the scaling on the three main variables of the model, namely, the residence time, erosion and deposition.First, we investigated the dependence of the scaling exponent for floodplain storage on the residence time.We chose different median residence times for floodplain sediment in the Rhine catchment, while keeping the maximum residence time fixed.Changing the median residence time by a factor of 10, from 50 to 500 years, results in a decrease of 21.8 % in the scaling exponent for floodplains in the transient simulation (Table 4).
When the median floodplain residence time is increased, the range in the residence time decreases.
This leads to a decreased difference in sediment loss between grid cells with small and large catchment areas, which then leads to a decrease in the scaling exponent.We find that when the residence time is increased by 5.2 % (from 50 to 260 years) the scaling exponent decreases by 18.2 %, while an increase in the residence time of 1.9 % (from 260 to 500 years) results only in a decrease of the scaling exponent of 4.4 %.This indicates that the scaling exponent for floodplain storage does not change linearly with the residence time, and points out that the model behaves in a non-linear way.Applying the same approach for the equilibrium simulation results in a similar behavior for the scaling exponent.However, here the 10 fold change in the residence time leads to a slightly larger change in the scaling exponent.
Next, we investigated the dependence of the scaling exponents for floodplain and hillslope storage on erosion.We changed the spatial variability of erosion in the Rhine catchment by changing the spatial variability of the R factor.We increased the R values in the Alpine region and decreased the R values in the rest of the catchment.This results in a larger difference between the sediment storage in small catchment areas and sediment storage in large catchment areas.Although the resulting scaling exponent for floodplain storage is still much higher than the scaling exponent for hillslope storage, both scaling exponents increase significantly.
For the deposition we find a minor to neglecting effect on the scaling parameters.
Overall we find that changing erosion and residence time does not change the basic property of the scaling, which is that floodplain storage increases in a non-linear way with catchment area while hillslope storage increases linearly with catchment area.As the residence time is determined by flow-accumulation and flow-accumulation determines the spatial variability of floodplain sediment storage, we expect that the scaling parameters for floodplain sediment storage are also mainly determined by flow-accumulation.Erosion is mainly determined by the slope, and slope determines the spatial variability of hillslope sediment storage.Therefore, we expect that the slope determines the scaling parameters for hillslope sediment storage.Based on this we argue that the scaling for both floodplain and hillslope storage is an emergent property of the model and that the scaling parameters are controlled by the underlying topography.

Last millennium sediment storage
We estimate an average soil erosion rate of 2.8 ± 0.002 t ha −1 year −1 for the last millennium for the entire Rhine catchment.We find that this value is twice as high as the 1.2 ± 0.32 t ha −1 year −1 , which was estimated as the minimum average soil erosion rate for the Holocene by Hoffmann et al. (2013).
The average soil erosion rate for the last millennium results in a mean floodplain and hillslope sediment storage change for the last millennium of 11.95 ± 0.01 and 29.68 ± 0.03 Gt, respectively (Table 5).Altogether, floodplain and hillslope storage result in 41.63 ± 0.02 Gt of sediment, which can be considered as the contribution of climate and land use change to sediment storage in the last millennium.It is, however, hard to say what the range in the change of sediment storage should be for this period, as there are no related studies for this specific time period.The total sediment storage we find is lower than the total Holocene sediment storage of 126 ± 41 Gt found by Hoffmann et al. (2007) for the Rhine catchment.This is logical as we consider only the last millennium and not the past 7500 years as in the study of Hoffmann et al. (2007).Our results show that the sediment storage of the last millennium form 25 to 50 % of the total sediment storage of the last 7500 years.This indicates that the average sediment storage rate during the last millennium is higher than the average rate during the last 7500 years.This also supports the findings from previous studies (Bork, 1989;Notebaert et al., 2011), which show that land use change has a significant and long-term impact on erosion and sediment mobilization.
Furthermore, Hoffmann et al. (2013) found a floodplain to hillslope ratio of about 0.88, indicating that during the Holocene more sediment was stored on hillslopes than in floodplains.We find with our model a floodplain to hillslope ratio of about 0.46, confirming that more sediment is stored on hillslopes.However, the floodplain to hillslope ratio from our model indicates a much larger difference in sediment storage between floodplains and hillslopes than in the study of Hoffmann et al. (2013).This can be attributed to the lack of an explicit representation of the size and location of floodplains in the model, and the simple representation of the sediment deposition processes for floodplains and hillslopes.
We also analyzed the spatial variability of the simulated sediment storage in floodplains, and find that the model reproduces the spatial variability well when compared to the observed values from Hoffmann et al. (2007) for the Holocene (Fig. 4).We find a correlation coefficient of 0.77, where sediment storage in floodplains increased with catchment area.Furthermore, we find that most floodplain sediment is stored in the Mosel sub-catchment, in contrast to the observations that show the largest storage in the Upper-Rhine sub-catchment (Table 6).This can be related to the fact that different dynamical processes, which are not captured with our model, play a role in the Upper-Rhine catchment.Melting ice sheets for example can produce a lot of erosion that is not captured by our model and in this way the total stored sediment in the catchment could be underestimated.Furthermore, the Mosel sub-catchment has a highly complex topography, which may indicate that our sediment budget model is too coarse for an accurate representation of floodplain storage for this catchment.
For hillslope sediment storage we find a similar spatial trend as for the floodplain sediment storage, with some more variation between the minimum and maximum values (Table 6).Also here, the Mosel sub-catchment stores the most sediment.Furthermore, when comparing floodplain to hillslope sediment storage we find that the floodplain to hillslope ratio varies significantly between the various sub-catchments.The highest ratio of 0.48 is found for the Lower Rhine sub-catchment, while the lowest ratio of 0.14 is found for the Emscher sub-catchment.The ratios seem not to be correlated with slope or catchment area and can be assumed as independent features of the model.
The sediment budget model presented here, has been developed to simulate long-term trends and to determine the main drivers behind these trends.Figure 5 shows the land use change and the 10 year-mean precipitation averaged over the Rhine catchment for the last millennium.There are two interesting periods, which are 1350-1400 AD and 1750-1950 AD, that show increased precipitation amounts correlating with a sudden increase in land use change (increase in crop and pasture).Both periods lead to maxima in the erosion timeseries of 2.8 t ha −1 year −1 and 4.3 t ha −1 year −1 , respectively (Fig. 6a and 6b).These rates correspond to increased erosion rates during the 14th and 18th century found by (Bork, 1989;Lang et al., 2003) for Germany.
We find the strongest increase in the sediment storage rate for floodplains during the period 1750-1850 AD, and for hillslopes during the period 1850-1950 AD.For hillslopes this maximum sediment storage rate corresponds to a maximum increase in the deposition rate, which is a result of a maximum increase in land use change.Land use change leads to a sediment disconnectivity in the landscape, which prevents the sediment stored on hillslopes of reaching the fluvial system on the timescale of the last millennium.In contrast to hillslopes, the maximum sediment storage rate for floodplains is a result of the interplay between deposition and sediment loss from the catchment.
In the period 1750-1850 AD land use change started to increase in the Alpine region.This region did not experience such a strong change in land-use before 1750 AD compared to the downstream regions of the catchment.During the period 1750-1850 AD, the deposition to floodplains increased significantly due to the increased erosion rates as a result of land use change.As land use change started to impact the Alpine region, steep slopes and short residence times lead to a strong sediment flux downstream.However, due to the long residence time of the areas located downstream, the sediment loss from the total catchment did not increase as much, leading to an increased sediment storage in the floodplains.This is in accordance with the findings of Asselman et al. (2003), who found that due to an inefficient sediment delivery, an increase in soil erosion in the Alps will have a little effect on sediment load downstream the Rhine catchment.
Furthermore, if we disentangle the effects of land use and climate on the sediment storage in floodplains and on hillslopes, we find that land use change explains most of the change in sediment storage.For floodplains climate change also has a non-negligible impact on the temporal variability of sediment storage.For example, in the periods 1350-1400 AD and 1750-1950 AD, the sediment storage rate is increased due to increased precipitation that lead to a strong sediment flux downstream.
If the land use conditions of the period 850-950 AD are kept constant, the total change in sediment storage in floodplains and on hillslopes during the last millennium is 2.9 and 15.4 Gt, respectively.This is four and two times, respectively, less than the change in floodplain and hillslope sediment storage when land use change is variable (Fig. 7a and 7b).Here, the overall sediment storage still increases due to the overall increased trend in precipitation during the last millennium.If only the climate conditions are kept constant, the resulting change in sediment storage in floodplains and on hillslopes is 10 and 27.4 Gt, respectively.

Uncertainty assessment and limitations of the modelling approach
As shown in the previous sections, the average modelled erosion rate for the last millennium of the Rhine catchment is overestimated when compared to the average erosion rate for the Holocene from the study of Hoffmann et al. (2013).As we consider in this study only the last millennium, where human impacts through land use change are strongest pronounced, it is logical that our estimated average soil erosion rate is higher.For present day, we estimate an average soil erosion rate of 3.3 t ha −1 year −1 for the non-Alpine part of the Rhine catchment, which is also overestimated when compared to other studies.Cerdan et al. (2010) found for the non-Alpine part of the Rhine catchment a value of 1.5 t ha −1 year −1 , while Auerswald et al. (2009) found for Germany a value of 2.7 t ha −1 year −1 .
Comparing our simulated erosion rates for present day with high resolution estimates from Cerdan et al. ( 2010), we find that our rates are overestimated for the entire Rhine catchment.We expect that the overestimation is mainly due to uncertainties related to the coarse input datasets on climate and land cover, and biases in the adjusted RUSLE model.For example, we find that precipitation is generally overestimated by MPI-ESM for the Rhine catchment.Even after introducing a correction factor, which partly adjusted the R value estimation to values from present-day observational datasets, biases related to the R factor remain.It is, therefore, also important to test the sensitivity of the sediment budget model with input data on precipitation and land cover from other ESMs.
Additionally, coarse resolution land cover fractions and leaf area index (LAI) from MPI-ESM also affect the total erosion rates.Using coarse resolution data to calculate the C factor of the adjusted RUSLE model results in discrepancies between the C and S factors.For example, consider a coarse resolution grid cell with a complex topography where cropland is located in flat areas and forest in the steeper areas.Even though the C factor is calculated correctly as combination of cropland and forest fractions, it is applied to the whole grid cell.This leads to an overestimation of erosion rates for flat areas, as erosion is in the first order controlled by the slope through the S factor.We attempted to correct for this by introducing slope classes for each coarse grid cell with resolution of MPI-ESM (1.875 degrees).The cropland was then allocated to the flatter areas, while other land cover types were allocated to the steeper areas.However, this only had a minor effect on the overall erosion rates, indicating that this is not the major source for the overestimation in erosion rates.
Additionally, the absence of the seasonality in the C factor results in discrepancies between the C and R factors.
Neglecting the support practice (P ) and slope-length (L) factors in agricultural regions, where they may play an important role, results in an overestimation of the increases of soil erosion, especially during the 1950's.However, we expect that this does not affect the overall trends.This assumption is also supported by Doetterl et al. (2012), who shows that the L and P factors explain only up to 22% of the variability in water erosion rates on cropland in the USA.
Also, biases in the adjusted RUSLE model, such as the unadjusted C and K factors and the low performance of the model in mountainous areas, have an equally important effect on the total erosion rates.
Another large uncertainty in our sediment budget model, besides the biases in erosion rates, is the choice of the equilibrium state.We find a decreasing trend in the floodplain sediment storage in the transient simulation when using the equilibrium state based on the mean conditions of 6000 BP.
This can be attributed to the different spatial distribution of erosion and the average high erosion rate for the mid-Holocene of 7.8 t ha −1 year −1 .When switching from the equilibrium state to the transient state, the erosion rates drop and its spatial distribution changes significantly.This leads to a decreased sediment flux from upstream areas and overall decreased sediment production rates that result in a decrease in sediment storage in the floodplains.For hillslopes we find that the equilibrium state has minimal to no influence on the total sediment storage for the last millennium.
The initial conditions determine the amount and spatial distribution of erosion in the catchment during the time that the model runs to equilibrium.Therefore, the equilibrium state that is then reached, largely determines the spatial distribution, trend, and amount of the sediment storage during the transient period.
Furthermore, the different ensemble simulations for the period 1850-2005 AD do not differ strongly in precipitation and land cover/land use, and, therefore, do not contribute much to the uncertainty in the overall erosion rates and sediment storage.This period is also too short to find significant effects on the sediment storage change using different ensemble simulations.
There are also some limitations to the model.The sediment yield cannot be accurately simulated for catchments where the initial state of the catchment is uncertain.However, with accurate data input on climate and land cover, the model can be made applicable for tropical catchments on the timescale of the last millennium, after adjusting the model parameters for these catchments.This is because we expect the effect of the last glaciation to be minimal on tropical catchments.In combination with little human activities during 850-950 AD, assuming an equilibrium state for these catchments in this time period seems reasonable.This can be tested in a future application of the model on other large catchments.We show that the model reproduces the scaling behavior between catchment area and sediment storage found in observed data from Hoffmann et al. (2013).The scaling behavior shows that the floodplain storage increases non-linearly with catchment area in contrast to hillslope storage.The scaling exponents can be modified by changing the spatial distribution of erosion or by changing the residence time for floodplains.However, the main feature of the scaling behavior is not changed.
Based on this we conclude that the scaling behavior is an emergent feature of the model and mainly dependent on the underlying topography.
We find a mean soil erosion rate of 2.8±0.002t ha −1 year −1 for the last millennium (850 -2005AD).
This is an overestimation when compared to the minimum Holocene erosion rate of 1.2±0.32t ha −1 year −1 from Hoffmann et al. (2013).Also for present day the erosion rates from our model are overestimated.We argue that this is mainly due to the coarse resolution input data on climate and land cover, and the fact that the land cover factor of the erosion model is not adjusted for a coarse resolution application.Additionally, the absence of the seasonality in the C factor plays a role, and other biases of the adjusted RUSLE model, such as the neglection of the land management and slope-length factors.However, we aim with the sediment budget model to distinguish between the floodplain and hillslope sediment storage, simulate their long-term behavior, and more specifically estimate the spatial distributions of sediment rather than the total amounts.For this objective a coarse estimation of erosion is sufficient.Furthermore, the model reproduces the overall spatial distribution of sediment storage in floodplains during the last millennium.However, there are some outliers, such as the Mosel sub-catchment for which the model simulates a too high sediment storage.This could be a result of biases in the erosion rates and the fact that our model is limited to the last millennium.We also found that the hillslope storage of the sub-catchments show a similar spatial pattern as the floodplain storage.
When analyzing the timeseries of erosion and storage during the last millennium we find that the model reproduces the timing of the maxima in erosion rates as found in the study of Bork (1989).
We also find that land use change is the main driver behind the trends in erosion and sediment storage for both floodplains and hillslopes.For floodplains, however, climate change has a non-negligible impact on the temporal variability of sediment storage.When keeping the land cover constant to the conditions in the period 850 to 950AD, we find that the sediment storage still increases due to an increased trend in precipitation during the last millennium.
We conclude that our sediment budget model is a promising tool for estimating large-scale long-term sediment redistribution.An advantage of this model is its capability to use the framework of ESMs to predict trends in sediment storage and yields for the past, present and future.
The next steps in quantifying soil redistribution on the global scale are the application of the sediment budget model on other large catchments, and validation of the model with existing data on net soil redistribution, sediment storage or yields.Furthermore, in order to make the soil redistribution model better applicable on a global scale and to prevent conflict with the underlying assumption of the simultaneous presence of floodplains and hillslopes in each grid box, the model needs to be made independent of grid resolution.
Finally to have a complete picture of the net soil redistribution and the feedbacks on the carbon and nutrient cycles, it is essential to model also other types of soil erosion, such as wind erosion (Chappell et al., 2015), tillage erosion (Van Oost et al., 2009) and gully erosion (Poesen et al., 2003).Table 1.Simulation specifications for the application of the sediment budget model for the Rhine catchment.
For each experiment with the sediment budget model the type of simulation (equilibrium or transient), the time period, and the initial conditions on which the simulation is based on, are given.Furthermore, we also provide the number of simulations we made with the model for a certain type of simulation, and the experiment from MPI-ESM that we used to derive the input data to force the sediment budget model.Table 2. Summary of regression results of the scaling of at the end of the equilibrium and transient simulations.
Here we consider only the grid cells that correspond to the observation points from Hoffmann et al. (2013) and fall into the borders of the Rhine catchment.The r-value is the Pearson correlation coefficient, and the slope and intercept are the scaling parameters.

Furthermore,
Wittmann and von Blanckenburg (2009) found a residence time of 600 years for floodplain sediments at Rees in the Rhine catchment, which falls in the range of the floodplain residence times of our study.We used the median and maximum residence times and the maximum flowaccumulation of the Rhine catchment to determine the a τ and b τ in equation 3. The exact values for a τ and b τ are -922442.54and 165886.77,respectively.2.3 Criteria for model evaluationHoffmann et al. (2013) compiled published data on sediment storage for regions in Central Europe, mainly for the Rhine catchment, where human induced soil erosion took place.Combined with

Furthermore,
Hoffmann et al. (2007) established a Holocene sediment budget for sediments in the floodplains and the delta of the non-Alpine part of the Rhine catchment.They derived sediment thickness of Holocene deposits from borehole data that consists out of 563 drillings and available geological maps.This was then multiplied with floodplain areas to calculate floodplain volumes.Sediments on hillslopes were not addressed in this study.A total floodplain sediment mass of 53.5 ± 12.4 × 10 9 run with a yearly time step till the total floodplain sediment storage of a catchment does not change more than 1 ton per year.The floodplain and hillslope sediment storage at equilibrium were then used as a starting point for the transient simulation that covers the period 850 -2005 AD.In the transient simulation erosion and deposition rates are averaged over time steps of 100 and 50 years, based on the time resolution of the rainfall erosivity factor (R) that is part of the erosion module.We performed five 'default' transient simulations, two based on the mid-Holocene equilibrium states, and three based on the equilibrium state of 850 -950 AD.The different ensemble simulations were used to investigate the uncertainty in the resulting sediment storage due to the input data of MPI-ESM.Additionally, we also performed a climate change and land use change simulation based on

Furthermore, a more
concrete parameterization for the residence time and deposition of floodplain sediment, and a possible new parameterization for the residence time of hillslope sediment could lead to an improvement of the model.Finally, more validation with long-term sediment storage from other catchments, especially tropical catchments, would be an important contribution in making the model applicable on the global scale.4ConclusionsIn this study we introduced a new model to simulate long-term, large-scale soil erosion and redistribution based on the sediment mass-balance approach.The main objective here was to develop a sediment budget model that is compatible with Earth System Models (ESMs), to simulate largescale spatial patterns of soil erosion and redistribution for floodplains and hillslopes following climate change and land use change.We applied this sediment budget model on the Rhine catchment as a first attempt to investigate its behavior and validate the model with observed data on sediment storage and erosion rates.
The simulated erosion rates result in a change in floodplain and hillslope sediment storage during the last millennium of 11.95 ± 0.03 and 29.68 ± 0.01 Gt, respectively.Based on this and the observed data we estimate that the climate and land use changes during the last millennium contribute between 25 -50% to the total sediment storage for the past 7500 years.Disentangling the contribution from climate change and land use change to the change in sediment storage during the last millennium for the Rhine catchment, we find that land use change contributes the most to the total change in sediment storage.

Figure 3 .
Figure 3. Scaling of floodplain (a) and hillslope (b) sediment storage from the transient simulation in the non-Alpine part of the Rhine catchment.The black dots and black trend line correspond to the observed sediment storage values from Hoffmann et al. (2013).The colored dots and colored trend line correspond to modelled sediment storage values that correspond to the observation points from Hoffmann et al. (2013) and fall into the borders of the Rhine catchment.

Figure 4 .
Figure 4. Observed versus modelled floodplain sediment storage for Rhine sub-catchments.The values are in percentage (actual storage divided by the sum times 100).Data on the observed sediment storage is taken from Hoffmann et al. (2007).RMSE is the root mean square error.

Figure 5 .
Figure 5. Land cover and precipitation variability averaged over the Rhine catchment for the last millennium.The red line is the 10 year mean total precipitation for the Rhine catchment.The background colors are land cover types, starting from the darkest grey to the lightest: forest, bare soil, grass, crop and pasture.Land cover and precipitation data is from MPI-ESM.

Figure 6 .
Figure 6.(a) Timeseries of simulated average erosion (black line), average deposition (green line) and the total change in sediment storage (blue line) with respect to 850-950 AD for floodplains in the last millennium in the Rhine catchment.(b) Timeseries of simulated average erosion (black line), average deposition (green line) and the total change in sediment storage (blue line) with respect to 850-950 AD for hillslopes in the last millennium in the Rhine catchment.

Figure 7 .
Figure 7. Simulated change in (a) floodplain and (b) hillslope sediment storage for the Rhine catchment during the last millennium.Shown are the sediment storage for the climate change simulation, where land cover is set to the conditions of the period 850-950 AD (CC -blue line), the sediment storage for the land use change simulation, where the climate is set to the conditions of the period 850-950 AD (LUC -red line), and the sediment storage where both climate and land cover change during the last millennium (CC and LUC -black line).

Table 3 .
Summary of regression results of the scaling of sediment storage after the equilibrium and transient simulations.Here we consider all grid cells in the Rhine catchment area.The r-value is the Pearson correlation coefficient, and the slope and intercept are the scaling parameters.

Table 4 .
Summary of regression results of the sensitivity analysis on floodplain sediment storage scaling.Here we consider only the previously mentioned selected grid cells in the Rhine catchment area.τ is the residence time of floodplain sediment.The r-value is the Pearson correlation coefficient, and the slope and intercept are the scaling parameters.

Table 5 .
Summary of sediment storage M (Gt), erosion (E) and deposition (D) rates in t ha −1 year −1 , and the related uncertainty ranges for the Rhine catchment for the period 850-2005 AD.The uncertainty values represent the range in the mean values due to different ensemble simulations.
Mean M Ensemble uncertainty M Mean E Ensemble uncertainty E Mean D Ensemble uncertainty D

Table 6 .
Hoffmann et al. (2007)ed sediment storage (Gt) for Rhine sub-catchments.The catchment area is given in km 2 .Data on the observed sediment storage is taken fromHoffmann et al. (2007).