Modelling a century of soil redistribution processes and carbon delivery from small watersheds using a multi-class sediment transport model

Over the last few decades, soil erosion and carbon redistribution modelling has received a lot of attention due to large uncertainties and conflicting results. For a physically based representation of event dynamics, coupled soil and carbon erosion models have been developed. However, there is a lack of research utilizing models which physically represent preferential erosion and transport of different carbon fractions (i.e. mineral bound carbon, carbon encapsulated by aggregates and particulate organic carbon). Furthermore, most of the models that have a high temporal resolution are applied to relatively short time series (< 10 yr−1), which might not cover the episodic nature of soil erosion. We applied the event-based multi-class sediment transport (MCST) model to a 100-year time series of rainfall observation. The study area was a small agricultural catchment (3 ha) located in the Belgium loess belt about 15 km southwest of Leuven, with a rolling topography of slopes up to 14 %. Our modelling analysis indicates (i) that interrill erosion is a selective process which entrains primary particles, while (ii) rill erosion is non-selective and entrains aggregates, (iii) that particulate organic matter is predominantly encapsulated in aggregates, and (iv) that the export enrichment in carbon is highest during events dominated by interrill erosion and decreases with event size.


Introduction
Numerical models of soil detachment, transport and deposition are important tools for improving our understanding of soil systems and the linkages between the terrestrial and aquatic ecosystems. At present, a wide range of simulation models are available. Conceptual models, such as the RUSLE (Römkens et al., 1997), focus largely on the prediction of long-term sediment production under various environmental and management conditions. In parallel, physically-based models have been 5 developed to simulate the routing of soil over complex topographies, taking hydrological and sediment-sorting processes into consideration (e.g., WEPP: Nearing et al., 1989;EROSION-3D: Schmidt et al., 1991;LISEM: De Roo et al., 1996). These models operate over relatively short time scales, typically one to several events, and are concerned with modelling the detachment and movement of mineral particles. Over the last few decades, they have been instrumental in improving our understanding of erosion processes and currently serve as tools for landscape management. 10 Erosion-induced changes in biogeochemical cycles, in particular carbon fluxes between soils, the aquatic environment and the atmosphere, have received considerable attention over the past two decades (e.g., Stallard, 1998;Renwick et al., 2004;Quinton et al., 2010). However, large uncertainties and conflicting results remain (e.g., Lal, 2003;Van Oost et al., 2007;Kuhn et al., 2009), and this has spurred renewed interest in the application of soil erosion models. To date, few soil and carbon erosion models integrate transport processes. There have been attempts to address this issue using single-point models with varying 15 degrees of complexity (Harden et al., 1999;Manies et al., 2001;Liu et al., 2003;Billings et al., 2010). These models apply prescribed carbon erosion and/or deposition rates and simulate the resulting effects on the SOC profile using CENTURY (Parton et al., 1988) parameterizations. Recently, spatially explicit models that combine erosion models with models of carbon dynamics have been developed (e.g., Changing Relief and Evolving Ecosystems Project (CREEP): Rosenbloom et al., 2001;Yoo et al., 2005;SPEROS-C: Van Oost et al., 2005a;Fiener et al., 2015). Both CREEP and the model presented by Yoo et al. 20 (2005) focus on long-term landscape development (i.e. millennial scale) and diffusive geomorphic processes that occur on undisturbed grasslands. The CREEP model also simulates textural differentiation and preferential transport of the finer fractions by surface wash. Compared to CREEP, SPEROS-C focuses on shorter timescales (i.e. years to decennia), puts a greater emphasis on erosion process description and emphasizes different processes, such as interrill erosion, concentrated flow erosion (e.g. rills and gully erosion), and diffusive processes due to tillage erosion (Van Oost et al., 2005a;Dlugoß et al., 25 2012).
Although these model concepts have facilitated an improved qualitative understanding of carbon erosion and erosion-induced changes in soil carbon storage, they are largely based on unverified assumptions and simplified process descriptions: First, carbon erosion is mostly approximated as being proportional to the bulk carbon:sediment ratio of topsoils. However, both experimental and modelling studies have clearly shown that erosion preferentially removes and exports soil organic carbon 30 (SOC; e.g., Polyakov and Lal, 2004;Schiettecatte et al., 2008a b;Kuhn et al., 2010). This preferential transport results from the fact that SOC is not distributed uniformly throughout the soil, but instead consists of several fractions, characterized by different densities and particle sizes. For example, some soil organic carbon is bound to the fine mineral fraction, some is Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License. encapsulated in soil aggregates, while another SOC fraction exists as mineral-free particulate organic carbon (POC) and has a much lower density (e.g., John et al., 2005;Von Lützow et al., 2007). This differentiation is particularly relevant for the C cycle, since for example the C fraction with the highest potential mobilization and transport capacity (i.e. POC due to its low density) is also a very labile fraction (Haynes, 2005). Thus, carbon erosion models should always consider the differential behaviour of sediment particles and SOC fractions when simulating erosion and transport processes. Second, carbon erosion 5 simulation models need to consider relatively long time scales, i.e. several years to decades, as carbon erosion fluxes are relatively small when compared to rates of soil C turnover (Fiener et al., 2015). Current models addressing erosion (e.g., CENTURY: Parton et al., 1988;EPIC: Williams, 1995;WaTEM: Van Oost et al., 2000;EDCM: Liu et al., 2003) use a constant average annual soil erosion rate by assuming uniformity. However, empirical observations indicate that soil erosion and sediment delivery are to a large extent controlled by extreme events (e.g. Fiener and Auerswald, 2007). This calls into question 10 whether the effects of erosion on biogeochemical cycles can reasonably be derived from continuous average long-term erosion rates. Event size also influences the extent to which selective transport takes place in erosion processes. For example, interrill erosion, which is a selective process (e.g., Kuhn et al., 2010), is more pronounced during smaller erosion events. As a result, there is more enrichment of fine soil fractions, including carbon associated with clay particles, during small events compared to large ones. It is therefore important that carbon erosion models correctly represent the different processes that control 15 selectivity. Furthermore, analysis on the relative contribution of low intensity (but high frequency) and more extreme (but low frequency) erosion events is required to understand the long-term effect on soil carbon dynamics. An important limitation of current approaches is therefore the frequent use of the USLE (Wischmeier and Smith, 1978) as a basis for erosion prediction.
The USLE was not designed to estimate frequency distributions of soil erosion but is in fact designed to average out variability (Wischmeier and Smith, 1978). 20 The Multi-Class Sediment Transport (MCST) model (Van Oost et al., 2004;Fiener et al., 2008) enables the simulation of carbon erosion, deposition and export. It also has the potential to be modified in order to address some of the issues discussed above. For example, assessing the effect of event magnitude on soil carbon erosion can be achieved by utilizing a range of rainfall input data. In addition, the MCST model uses selective transport to route sediment over areas of erosion and net deposition. Therefore, preferential export of finer particles, as well as the deposition of coarser particles, is represented by the 25 model. Overall, by combining a model for the dynamic simulation of rainfall-runoff processes with a multi-class sediment transport model that considers different carbon fractions, the MCST allows for the examination of sediment and carbon mobilization and export in a temporal framework.
The main objective of this paper is to improve our mechanistic understanding of sediment redistribution and carbon delivery using numerical experiments. To this end, the MCST model is modified to incorporate the natural long-term variability of soil 30 and soil organic carbon erosion. Existing empirical observations, covering a century of carbon delivery, will be used to assess the model behaviour and to identify potential deficiencies in model process descriptions. Finally, the long-term role of event size on both rates and patterns of soil and carbon erosion will be evaluated and discussed.

Methodology
The MCST model (Van Oost et al., 2004;Fiener et al., 2008) combines a soil infiltration component with a kinematic wave routine to produce continuous series of runoff events. The event-based soil erosion component describes detachment as a function of rainfall characteristics, slope and discharge, while transport and deposition are simulated using the Hairsine and Rose (1992a b) equations. The two-dimensional implementation in a regular grid (1x1 m 2 to 5x5 m 2 ) uses a digital elevation 5 model to route overland flow and sediment redistribution. A detailed model description can be found in Van Oost et al. (2004) and Fiener et al. (2008); here we focus on its main features and modifications made in order to continuously simulate longterm (up to centuries) soil and carbon erosion.

Modelling surface runoff
The model calculates rainfall excess at a fine temporal resolution (minutes to hours) using a modified Curve Number approach. 10 The original version of the MCST model simulates single rainfall events but is converted into a continuous simulation model as follows: The input of the model is a continuous rainfall series with a time resolution of 10 minutes. A rainfall-runoff event is identified as a period (i) in which rainfall depth exceeds 2 mm in 24 h (<1% of total runoff excluded) and (ii) which is separated by at least 72 h without rainfall. Accordingly, a rainfall-runoff event is not necessarily defined by a single hydrograph, but might contain multiple runoff peaks. A moving window of 24 h is used to estimate cumulative rainfall (Pi,cum) 15 and cumulative abstractions for each time-step (i.e. initial abstraction (Ia,cum) and continuing abstraction (Fa,cum) of the curve number method). The excess rainfall hyetograph (Ri) at time step is calculated as: where , (mm) is the cumulative excess rainfall during the last 24 hr. 20 is a scaling factor for rainfall intensity which is calculated as: 10 is the maximum 10-minute rainfall intensity (mm/h). 25 Flow discharges for each grid cell and time-step are calculated by numerically solving the kinematic wave equations (Van Oost et al., 2004). For sheet flow, cross-sectional flow area is calculated assuming a homogeneous flow depth in each raster cell, while for concentrated flow, a relationship between discharge and cross-sectional flow area is used (Govers, 1992). To distinguish between sheet and concentrated flow, a critical shear velocity of 3.5 cm s -1 for rill initiation, based on flume experiments conducted by Govers (1985), is used. The model keeps track of changes in the pattern of concentrated flow and 30 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License. rill network development. Finally, sediment movement is described by utilizing an event-based steady-state sediment continuity equation proposed by Yu et al. (1997).

Modelling erosion and deposition
Experimental research has shown that the Hairsine-Rose model provides an accurate physically based description of sediment transport and deposition for multiple sediment classes that differ in terms of settling velocities (Beuselinck et al., 2002a b). 5 Transport of soil by overland flow is characterized by simultaneous re-entrainment and deposition (i.e. temporary settlement) of sediments: where is the mass rate of deposition per unit area of size class (kg s -1 m -2 ), is the rate of sediment re-entrainment for settling velocity class (kg s -1 m -2 ), is the mean sediment concentration (settling velocity class ) (kg m -3 ), is the ratio of 10 the sediment concentration adjacent to the bed to the mean sediment concentration, is the settling velocity of sediment size class (m s -1 ), is the fractional shielding of the soil by the deposited layer, is the fraction of stream power used for reentrainment, is gravity (m s -2 ), is the sediment density of settling velocity class (kg m -3 ), is the water density (kg m -3 ), Ω is the stream power (W m -2 ), Ω is the critical stream power (W m -2 ), is the depth of the water flow (m), is the mass of sediment class in the deposited layer (kg m -2 ), is the total mass of the deposited layer per unit area (kg m -2 ). 15 If the local stream power (Ω) is less than a critical threshold (Ω ), re-entrainment does not occur and deposition of size class is a function of its specific settling velocity (Beuselinck et al., 1999;Hairsine et al., 2002). If the local stream power exceeds this threshold value, a shielding factor is calculated to decide whether net erosion or deposition occurs (Hairsine and Rose, 1992a): 20 If H ≥ 1, then net deposition, characterized by steady state flow and re-entrainment of previously deposited sediment, occurs.
If < 1, net erosion occurs, and soil detachment is modelled as: where and are the rill detachment rate and the interrill sediment transport to the rill (kg m -2 s -1 ), respectively, is the rill erodibility factor, is the interrill erodibility factor, is the rill discharge (m 3 s -1 ), is the local slope gradient, is the maximum 10 min rainfall intensity, is a slope factor and , and are calibration exponents.
Rill erosion is considered to be unselective, i.e. the sediment particle size distribution of the eroded material equals the 30 distribution of the source material at the source location. In contrast, interrill erosion is simulated as a selective process: Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License. assuming steady state flow conditions, Eq. (4) is used to estimate the particle size distribution of the sediment detached by interrill erosion and Eq. (7) is used to estimate the transport for sediment delivered to the rill network (or that leaves a grid cell when there is no incised rill). This approach is consistent with empirical observations showing that the enrichment of finer sediment particles and SOC in suspended sediment is mainly controlled by the transport capacity of the flow (Schiettecatte et al., 2008a). 5 The MCST model keeps track of spatio-temporal changes in particle size distribution of the eroded and deposited topsoil sediment. The model considers 10 different size fractions.

Model implementation
For our modelling based analysis, we combined data from different sources into a virtual catchment data-set: (i) All basic data (i.e. digital elevation model, soils) were taken from a small first-order catchment in central Belgium, located about 15 km 10 southwest of Leuven. The site has a size of 3 ha with a mean and maximum slope of 7% and 14%, respectively. The catchment consists of diverging convex hillslopes and a central concavity where ephemeral gullying and sediment deposition are frequently observed ( Fig. 1; Desmet and Govers, 1997). Soils in the catchment are loess-derived, silty-loamy Luvisols, with a clay, silt and sand content of 14%, 82% and 4%, respectively (Desmet and Govers, 1997). (ii) For the 100 yrs. modelling period, high resolution rainfall data (1989-1997; 10- This range resulted in runoff volumes that are consistent with field observations (Gillijns et al., 2005). (ii) Two annual tillage operations are assumed to erase the network of rills and ephemeral gullies which may have evolved during preceding erosion 20 events. (iii) The rill and interrill erodibility parameter values, as well as the slope and discharge exponents (Eq. 6 and 7), were assumed to be constant over time and space. Therefore, spatio-temporal variability of soil moisture is not accounted for. The parameter values are taken from flume and plot-scale experiments, conducted using soils from the Belgium loess belt (Table   1; Van Oost et al., 2004). With these parameters, MCST has already shown to be able to predict the spatial patterns and rates of sediment detachment and transport in the test catchment (Van Oost et al., 2004). 25 In simulation studies, the particle size distribution is typically derived from dispersed sediment samples and therefore reflects the settling velocities of the primary particles of the sediment. However, sediment transport and deposition can also occur in the form of aggregates, particularly for fine textured soils, as is the case in our study area . Therefore, we considered the particle size distributions of both aggregated soil and primary particles in our simulations. We consider two erosion scenarios: (i) erosion scenario 1, in which both rill and interrill erosion detach primary particles and (ii) erosion scenario 30 2, in which interrill erosion detaches primary particles and rill erosion detaches aggregates. Following detachment, the model simulates the transport of primary particles or aggregates based on the type of detachment that they underwent. The particle size distributions of primary particles and aggregated soil were taken from direct measurements (n=81) in the Belgian loess Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License.
belt conducted by Beuselinck et al. (1999). The corresponding settling velocities were calculated according to the model of Dietrich (1982), using a density of 2.6 and 1.3 kg m -3 for primary particles and aggregates, respectively. The density of primary particles is assumed to be close to quartz, whereas a pore space of 50% is assumed for aggregates. The settling velocity distributions (Fig. 2) show that the aggregated sediments are dominated by fractions with settling velocities between 10 -4 and 10 -3 m s -1 , i.e. silt sized particles. In contrast, erosion scenario 1, which solely considers primary particles, shows very low 5 settling velocities relative to aggregated sediments (Fig. 2). This results from differences in particle size between the two fractions: aggregated soils contain fewer clay and silt sized particles, because particles of this size tend to be occluded in micro and macro aggregates. As a result, aggregates have larger particle sizes and faster settling velocities.
The implementation of SOC characteristics in the model is based on a SOC fractionation study that was carried out with similar soils (Luvisols) from the Belgian loess belt (Doetterl et al., 2012). In the study of Doetterl et al. (2012), soil samples were 10 taken at 11 locations along a topographic gradient, from non-eroded to eroding and depositional sites. The results showed that 85% (±10%) of the total SOC was associated with the mineral fraction (clay and silt size), while the remaining 15% (±3%) was particulate organic matter (POM; Fig. 3). The POM fraction was enriched in SOC relative to the bulk soil (enrichment factor of 1.87; Doetterl et al., 2012). To our knowledge, no detailed information is available on the allocation of SOC in particle size fractions from 2 to 63 µm. However, others have found that 55 to 82% of mineral-bound SOC is generally associated with 15 clay-sized particles in temperate arable soils, while the remaining 18 to 45% is associated with finer silt fractions (2-20 µm; Von Lützow et al., 2007). Based on these observations, and given the constraints imposed by the model structure, we consider two types of SOC: (i) mineral-bound SOC, which represents 90% of the total and is associated with the finest sediment class (< 2 µm) and (ii) a POM fraction, which represents 10% of total SOC and is considered a separate class in the model, with a particle size of 250 µm and a density of 1 kg m -3 (Doetterl et al., 2012). 20

Model evaluation
We evaluate the performance of the model by comparing the predicted characteristics with those that were continuously observed in the Kinderveld (250 ha) and Ganspoel (117 ha) agricultural catchments for two observation periods of 3 years each (6 years total observation; Van Oost et al., 2005b). The two catchments are situated approx. 15 km from our study site and are very similar to our site in terms of soil properties and geomorphology. We were unable to directly apply our model to 25 these 2 agricultural catchments, as our model has high data requirements, which could not be met due to large uncertainties in input data, or in some cases because the data was simply not available. Rather than providing an evaluation on an event-basis, we evaluate the model performance by looking at the characteristics of sediment and carbon delivery, in response to a range of erosion event-sizes. This provides a first, but stringent, test of model structure and assumptions.

Frequency analysis 30
Because the time duration of a simulated runoff event may differ significantly and can contain multiple hydrographs, it is not possible to calculate event-based frequency distributions. Therefore, the event-based model results are aggregated on a monthly The recurrence interval is expressed in years when is multiplied by the number of modelled years.
To calculate the frequency of exceedance, monthly soil erosion values were ranked in increasing order, and a rank is given 5 to each modelled soil erosion event. The exceedance probability for event is given by: where is the total number of events during the period.

Rainfall/runoff 10
Application of the rainfall/runoff model over a period of 100 years resulted in 792 individual rainfall/runoff events. The temporal variability of rainfall events is relatively low, as more than 70% of total rainfall is associated with events with a recurrence interval of less than 1 year. Extreme rainfall events do occur, but their relative contribution to total rainfall is limited (i.e. events with a recurrence interval ≥ 2 years contribute less than 18%). The model simulates that, integrated over the period of simulation, about 10% of the total rainfall does result in surface runoff. This is consistent with field observations in the 15 study area, where an average of 8% was reported by Steegen et al. (2000Steegen et al. ( , 2001. In contrast, the simulated temporal variability in runoff is high, and events with a larger recurrence interval, i.e. ≥ 2 years, make up more than 36% of total runoff. The variability in runoff is higher than that of rainfall because it is controlled by multiple factors, including rainfall amount and intensity, vegetation characteristics, soil surface conditions and the presence and/or absence of a rill/ephemeral gully network at the beginning of an event. 20

Interrill and rill/ephemeral gully erosion
Interrill erosion, modelled here as a function of slope and rainfall intensity, accounts for 14% of total sediment mobilization over a 100 year period. Rainfall intensity is the main factor controlling interrill erosion and explains about 70% (p < 0.001) of the variability. In contrast, incised (i.e. rill/ephemeral gully) erosion was modelled as a function of slope and discharge and is therefore mainly controlled by surface runoff (R² = 88%; p < 0.001). The simulated relative contribution of interrill erosion 25 depends on the suspended sediment concentration (SSC) and the sediment delivery ratio (SDR). For events with low values for SSC and SDR, the contribution of interrill erosion can account for up to 100% of total sediment mobilization, but this contribution declines rapidly with increasing SSC and SDR (Fig. 4). This reflects the role of event size, whereby only larger events produce significant amounts of concentrated erosion once the hydraulic threshold for rill initiation is exceeded.

Sediment and carbon mobilization and export under different scenarios
We evaluated two erosion scenarios where different assumptions about particle size distribution are made. In erosion scenario 1, sediment is detached as primary particles, whereas in erosion scenario 2, soil is detached and transported in aggregated form during rill erosion. The simulated long-term enrichment ratio of the deposits for the fine fraction (< 2 μm), which results from selective erosion and transport processes, was found to be 0.02 and 0.6 for erosion scenario 1 and 2, respectively. For erosion 5 scenario 1, this implies that the deposition of clay particles is virtually non-existent and also suggests a very efficient export of clay materials from first-order catchments. However, these results are not consistent with field observations. Data derived from the Belgian Soil Map (Baeyens, 1959) shows only small differences between the primary particle size distributions of colluvial and non-eroded agricultural soils in our study area. The reported enrichment ratio for the clay fraction of colluvial soils is 0.8 (Baeyens, 1959). The colluvial sediment is thus only slightly depleted in clay when compared to the source material. 10 Based on this analysis, we consider the results of the simulations for erosion scenario 1 to not be physically valid. In contrast, the results of erosion scenario 2 are qualitatively similar to the observations (Baeyens, 1959), which suggests that erosion scenario 2 more accurately depicts erosion processes in our study area. This implies that the assumptions made for erosion scenario 2 are appropriate, i.e. interrill erosion detaches and transports primary particles, whereas rill erosion is unselective and detaches and transports soil aggregates. The model tends to slightly over-estimate the export of the finest fractions 15 (enrichment of colluvial sediments: 0.6 model versus 0.8 field observations). Beuselinck et al. (2002b)  The clay enrichment ratios at the outlet of the catchment (i.e. clay exported sediment / clay source material) for the simulated events range between 1 and 4.8 (Fig. 5). These ratios are strongly related to the SSC: high enrichment ratios occur when SSC is low, i.e. during small events with a low recurrence interval. In contrast, low enrichment ratios (i.e. close to 1) are associated with events characterized by a high SSC. These findings are in line with other studies (e.g., Schiettecatte et al., 2008ab, Wang et al., 2010 and emphasize the importance of event size. The contribution of interrill erosion is higher for small events and, 25 since interrill erosion is modelled as the detachment and export of individual sediment particles, this results in a higher clay enrichment ratio. Vice versa, the contribution of interrill erosion is small for large events, resulting in enrichment ratios close to 1, since concentrated erosion is assumed to be unselective. The simulated range of enrichment ratios and the relationship between those ratios and SSC are both very similar to that which was observed at the Kinderveld and Ganspoel catchment (Fig. 5;Wang et al., 2010). 30 However, over a simulation period of 100 years, the flux-weighted predicted clay enrichment ratio in exported sediments was found to be 1.4, which is lower than the field observed ratio of 1.5 -2.6 for a 6 year period (Wang et al., 2010). Possible explanations for this discrepancy are (i) that the model underestimates the detachment and/or breakdown of aggregates into Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License. finer fractions during transport, or (ii) that concentrated erosion also results in selective detachment of sediments. Nonetheless, it is not unreasonable that the much longer simulation period resulted in lower clay enrichment ratios, since the contribution of extreme events, with high export rates but low enrichment ratios, is more important over longer time scales.
The simulated enrichment of carbon is directly related to the preferential export of the clay fraction through by interrill erosion.
The simulated carbon enrichment ratios are higher than the clay enrichment ratios previously discussed and range between 1 5 and 9 (Fig. 5). Exported sediment is more enriched in carbon than it is in clay due to the fact that the clay fraction is itself enriched in carbon relative to the bulk soil. The simulated relationship between SSC and carbon enrichment is similar to what was found for clay enrichment, i.e. enrichment is higher when SSC is low. This is again consistent with the Kinderveld and Ganspoel field observations (Wang et al., 2010).
It should be noted that the enrichment of exported clay and carbon was simulated assuming that interrill erosion resulted in the 10 detachment of primary particles while concentrated erosion resulted in the detachment of aggregates. Alternatively, clay and POM fractions could be considered as individual classes in the model. However, due to very low settling velocities, nearly the entire mobilized clay and POM fractions are exported from the catchment when this is simulated. This is not in line with field observations or with experiments that show that the transport of fine-textured sediments mainly occurs in the form of aggregates (e.g., Beuselinck et al., 2000, Wang et al., 2013. 15

Frequency and magnitude of erosion and delivery of soil constituents
Based on the 100 year modelling period, we analysed the effect of event frequency and magnitude of erosion on mobilisation and delivery of bulk sediment, clay and SOC on a monthly basis (Fig. 6). We found that for within catchment erosion, a large number of relatively small events (recurrence interval < 1.4. yrs.) accounts for about half of the cumulative erosion, while larger events (> 10 yrs. recurrence) account for only about 15%. 20 The sediment delivery ratio (SDR), which is the fraction of eroded soil that is transported to the catchment outlet, was 0.18 over the 100 year simulation period, while the mean erosion rate was 12.5t ha -1 yr -1 . Figure 6 clearly shows that larger events play a more important role in determining SDR than they do in determining soil erosion. Approximately 59% of the total sediment delivery comes from events with a recurrence interval less than or equal to 10 years (Fig. 6). This is explained by the fact that sediment delivery is not linearly related to runoff amount: once the hydraulic threshold is exceeded (i.e. an extensive 25 network for concentrated flow is established, and the rill/ephemeral gully network is directly connected with the outlet of the catchment) sediment delivery rates can be very large. The simulation of this key mechanism requires the introduction of (i) differentiated hydrological behaviours for sheet and concentrated flow and (ii) rill/ephemeral gully network development tracking. Our simulations show that the highest export rates occurred when the rill/ephemeral gully network was already well established at the beginning of an event. 30 The importance of event size for simulating clay and SOC delivery is also shown in Figure 6. Compared to bulk sediment, the delivery of clay and SOC is less driven by rare, large events, since small events with more interrill erosion already deliver relatively large amounts of clay and SOC. In general, the model results underline the importance of a more process based Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License.
analysis of SOC redistribution, as the effects of small erosion events, e.g. upon aquatic ecosystems, are underestimated when modelling only mean bulk erosion rates.
In order to qualitatively evaluate our predicted temporal patterns for sediment delivery, we compared our results to studies that continuously measured export from small catchments. In one such study, which was carried out in small agricultural catchments in the Belgian loess belt, Steegen et al. (2000) measured sediment delivery over a 3 year period in two first-order 5 streams. The authors found that a single event contributed to more than 40% of the total sediment delivery during the observation period and that the sum of rare and extreme events accounts for 46%. Even more extreme results were reported from a small loess catchment (3.7 ha) in Southern Germany, where a short-term series of single runoff events accounted for up to 67% of total sediment export (> 0.5 mm runoff) over an eight-year period (Fiener et al., 2008). Although a quantitative comparison of the model results with these empirical observations is not possible, as empirical observations typically cover 10 far fewer than 100 years, this analysis strongly indicates that the mechanisms incorporated into the MCST model (i) allow for a quantitative representation of the relative importance of both small and large events and (ii) account for event size related spatio-temporal variability of both soil erosion and sediment delivery distributions.

Conclusions
In this study, we incorporated preferential erosion and transport of sediment and soil organic carbon (SOC) fractions into a 15 numerical model of surface runoff and sediment transport. In doing so, we were able to predict the export of these different classes of sediment and SOC from small hilly watersheds, located in a temperate region with fine-textured soils. The model predictions were only consistent with field observations when (i) interrill erosion was simulated as a process that entrains primary particles, (ii) rill erosion is unselective and (iii) low-density POM is encapsulated within soil aggregates and cannot be entrained by interrill erosion. These results suggest that carbon enrichment at the outlet of small watersheds occurs as a 20 result of the selective interrill transport of clay and fine-silt associated carbon. Based on the application of the model over a period of 100 years, we conclude that sediment export is a highly episodic process. Our results show that 63% of the total sediment export was caused by 20 single events with a rainfall recurrence > 5 years. This highlights the need to consider sufficiently long timescales when addressing the impact of lateral fluxes of sediment and nutrients on soil processes. However, the dominance of large events is less pronounced in the case of soil organic carbon delivery, where only 44% of total delivery 25 is caused by extreme events (recurrence interval > 5 years). This reduced importance is associated with the selective process of interrill erosion and transport.
This study highlights the need for an event-based analysis of carbon erosion and delivery in order to assess the overall effect of soil organic carbon redistribution on the terrestrial carbon balance. Moreover, the episodic nature of soil organic carbon redistribution is particularly important when considering the effects of SOC input to surface water bodies . 30 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-32, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 1 July 2016 c Author(s) 2016. CC-BY 3.0 License.