Does soil thinning change soil erodibility? An exploration of long-term erosion feedback systems

. Soil erosion rates on arable land frequently exceed the pace at which new soil is formed. This imbalance leads to soil thinning (i.e., truncation), whereby subsoil horizons and their underlying parent material become progressively closer to the land surface. As soil erosion is a selective process and subsurface horizons often have contrasting properties to the original topsoil, truncation-induced changes to 5 soil properties might affect erosion rates and runoff formation through a soil erosion feedback system. However, the potential interactions between soil erosion and soil truncation are poorly understood due to a lack of empirical data and the neglection of long-term erodibility dynamics in erosion simulation models. Here we present a novel model-based exploration of the soil erosion feedback system over a 500-year period, using measured soil properties from a diversified database of 265 agricultural soil profiles in the 10 United KingdomUK. For such, we adapted the Modified Morgan-Morgan-Finey model (MMMF) to perform a modelling experiment in which topography, climate, land cover, and crop management parameters were held constant throughout the simulation period. As selective soil erosion processes removed topsoil layers, the model gradually mixed subsurface soil horizons into a 0.2 m plough layer and updated soil properties using mass balance mixing models. Further, we estimated Tthe uncertainty in the 15 model simulations was represented with a Monte-Carlo-based forward error assessmentforward error assessment based on Monte Carlo methods. We found that modelled erosion rates in 99 % of the soil profiles in 39 % of the soil profiles were sensitive to truncation-induced changes in soil properties. The soil losses in all except one of the truncation-sensitive profiles

the land surface and reduced soil detachment. Contrastingly, the soils with siltier subsurface horizons 25 continuously replenished the plough layer with readily erodible material, which accelerated prevented the decline of the soil losses loss rates over time.. Ultimately Although our results are limited by the edaphoclimatic conditions represented in our data, as by our modelling assumptions, our resultswe have demonstratede how modelled soil losses can be sensitive to erosion-induced changes in soil properties, which in turn may accelerate or slow down soil thinning. These findings are likely to affect how we 30 calculate soil lifespans and make long-term projections of land degradation.

Introduction
Rates of soil erosion on agricultural land often exceed the rates at which new soil is formed (Evans et al.,35 2019; Montgomery, 2007). This imbalance is one which, left unchecked, can pose a critical threat to the sustainability of global soil resources and their ability to deliver vital ecosystem services across environments and society (Bot et al., 2000;Quinton et al., 2010). Moreover, as soils become thinner (i.e., truncated), the subsoil horizons and their underlying parent material become progressively closer to the land surface. This process might affect physical, chemical and biological topsoil properties (Bouchoms et 40 al., 2019;Papiernik et al., 2009;Vanacker et al., 2019), as well as soil water availability to plants and ultimately crop growth (Herbrich et al., 2018;Öttl et al., 2021;Schneider et al., 2021).
For instance, plot-and catena-based studies report that truncated agricultural soil profiles often display increased clay and/or sand contents in their Ap horizons (Rhoton and Tyler, 1990), compared to the soils from non-eroding positions in the landscape. Moreover, eroded Ap horizons tend to have higher bulk 45 density, lower organic carbon content, and lower water holding capacity (Olson and Nizeyimana, 1988;Stone et al., 1985;Strauss and Klaghofer, 2001). However, such patterns are highly variable and greatly dependent on the properties of the underlying subsoil material being progressively tilled into the Ap horizon (Lowery et al., 1995).
Erosion-induced changes to soil depth and soil properties can therefore influence soil losses and runoff 50 formation through a soil erosion feedback system (Morgan et al., 1984;Vanwalleghem et al., 2017). That is, erosion-induced changes to soil physical properties might affect soil erodibility (i.e., the susceptibility of soil to erosion), which may accelerate or slow down soil losses. Understanding how such system might develop over time and under assorted conditions might beis an important step to proactively design and implement effective soil conservation strategies, as different soils are likely to be impacted by erosion in 55 varied ways (Hoag, 1998). However, the empirical data over decadal to centennial timescales required to explore the feedbacks between soil erosion and soil thinning are currently non-existent. It follows that process-oriented soil erosion models are arguably the only available tool to simulate how erosion processes interact with truncation-induced changes in the soil system.
Process-oriented soil erosion models allow for thea representation of multiple mechanisms that influence 60 soil erosion, from basic processes such as particle detachment by raindrop impact and surface runoff, to more complex interactions between soil properties, hydrological processes, climate, and plant cover (see Merritt et al., 2003 for an overview). This ability to simulate the response of specific soil erosion processes to external stimuli makes process-oriented models useful for exploring what-if scenarios. Hence, models such as the Water Erosion Prediction Project (WEPP; Nearing et al., 1989), the LImburg Soil Erosion Model 65 (LISEM; De Roo et al., 1996), and the Morgan-Morgan-Finey model (MMF; Morgan, 2001;Morgan et al., 1984;Morgan and Duzant, 2008) are frequently used to explore the impacts of land use or climate change on soil erosion (e.g., Anache et al., 2018;Eekhout et al., 2021;Nearing et al., 2005).
To date most soil erosion models and model users assume that the inherent erodibility (i.e., the susceptibility of soil to erosion) of different soil horizons down a soil profile is constant over the period of a model 70 simulation. As upper soil horizons are removed by erosion, exposing the subsurface material, the implicit assumption in soil erosion modelling is that this erodibility is not variable, such that any changes to projected erosion rates are solely a factor of climate, land cover, and topography (e.g., Ciampalini et al., 2020;Eekhout et al., 2021;Panagos et al., 2021). However, since erodibility is a reflection of soil physical, chemical, and biological properties, and given that subsoils typically (although not exclusively) exhibit 75 contrasting soil properties to those observed in upper horizons, it follows that erodibility is not necessarily a constant as a soil profile thins. Furthermore, soil erodibility might change over longer timescales due to the coarsening and armouring of surface soils (Sharmeen and Willgoose, 2007;Willgoose and Sharmeen, 2006) and the depletion of erodible material as a result of extreme soil truncation (Anselmetti et al., 2007).
Although the soil erosion feedback system has been recognised as a key challenge for modelling past and 80 future erosion rates (Vanwalleghem et al., 2017), long-term dynamics of soil erodibility are an underexplored topic in erosion research. Exceptions come from landscape evolution models, which simulate the influence of erosion and deposition on soil development over millennia (Van der Meij et al., 2020;Sommer et al., 2008). However, in areas under severe erosion rates, subsoil horizons can be exposed within a matter of decades (Evans et al., 2020), which might trigger unexpected responses regarding runoff 85 formation and soil losses. Moreover, soil truncation can introduce substantial spatial variability to soil properties, often not accounted for in static soil maps used for a variety of purposes (Świtoniak et al., 2016).
Still, the lack of knowledge about the potential interactions between soil erosion and soil erodibility (i.e., in which timescale are erosion feedback mechanisms developed, how are different soil properties and soil types affected by truncation?) interpose the representation of soil truncation in soil-erosion-model 90 applications.
Here, we hypothesise that erosion-induced changes to soil properties affect soil loss rates over long-term periods. To evaluate such hypothesis, we performed anwe present a novel exploration of the soil erosion feedback system by simulating 500 years of soil losses and surface runoff on 265 agricultural soil profiles 95 in the United KingdomUK. We This allowed us to investigate how soil erosion rates respond to truncationinduced changes in soil properties and unravel the processes potentially driving such responses in different soil types in the UK. To the best of our knowledge, this is the first time soil erosion models have been used to understand the interactions between soil erosion, soil thinning, and soil erodibility, and how these interactions are established in varying soil types. An enhanced understanding of such dynamic soil erosion 100 feedback system will be crucial for improving the calculation of soil lifespans, providing future soil loss projections, and designing long-term soil conservation strategies.

Concept
Our modelling concept is essentially a numerical thought experiment, in which land cover, agricultural 105 management, climate, and topography parameters are held within a constant range, so that any changes in simulated soil losses and surface runoff over a 500-year period are solely a result of changes in soil properties due to erosion processes (Fig. 1). To perform this experiment, we parametrised a soil erosion model using measured data from 265 agricultural soil profiles spread across the UK (Fig. 2)nited Kingdom (Fig. 2). The abstract spatial scale of the simulations can be perceived as a pedon located on conventionally 110 tilled hillslope with winter cereals. For simplicity, we assume this spatial unit does not receive runoff and sediment input from upslope. As the original topsoil of each profile/pedon is successively removed by erosion, our model gradually mixes the subsurface horizons into a 20 cm0.2 m plough layer (i.e., the average tillage depth in the UK (Townsend et al., 2016)), continuously updating soil properties through massbalance models and pedotransfer functions (Fig. 1). 115  In order to implement our modelling concept, we adapted the Modified Morgan-Morgan-Finey model (MMMF) (Morgan and Duzant, 2008). The model was chosen due to its ability to simulate multiple erosion subprocesses (e.g., runoff formation, detachment by raindrop impact and runoff, particle size selectivity), 125 which is desirable for understanding the specific mechanisms responsible for developing erosion feedback systems with a parsimonious parameter set. That is, MMMF represents particle-size selectivity during erosion, transport, and deposition, incorporates the effects of stone cover on soil detachability, and simulates particle detachment by both raindrop impact and surface runoff. In addition, the MMMF has a parsimonious parameter set, which facilitates model application using national soil survey datasets. 130 ImportantlyMoreover, the MMMMF and its derivatives haves provided acceptable predictions of annual soil losses for different soils, land covers, and testing sites in England the UK (Morgan and Duzant, 2008;Peñuela et al., 2018;Smith et al., 2018). (Morgan and Duzant, 2008). In the following we provide a brief description of the basic MMMF equations (Section 2.2). We subsequently characterise the soil profile database used for the modelling (Section 2.3) and describe the model implementation, including the mixing 135 of surface and subsurface horizons (Section 2.4).

MMMF operating equations
The MMMF is a process-oriented conceptual model running on an annual timestep, in which soil erosion processes are separated into a water phase and a sediment phase (Morgan et al., 1984;Morgan and Duzant, 2008). In the water phase, effective annual rainfall (PEF; mm) is calculated considering the effect of 140 interception by the vegetation cover: Where P is the mean annual rainfall (mm) and PI is the average rainfall interception (proportion 0-1) afforded by the vegetation cover. For annual crops, PI and all other land cover parameters are taken as an approximate average over the growing season.
Annual leaf drainage (LD; mm) and direct throughfall (DT; mm) are separated as a function of the average 145 canopy cover of the vegetation (CC; proportion 0-1): The kinetic energy of direct throughfall KEDT is calculated with the typical value of erosive rainfall intensity for a given location (I; mm hr -1 ) and the amount of annual direct throughfall, whereas the kinetic energy of leaf drainage KELD is a function of the average plant height for the growing season (PH; m) and the amount of the annual leaf drainage. Total kinetic energy of the effective annual rainfall (KE; J m -2 ) is then calculated 150 as the sum of the throughfall and leaf drainage components.
Of note is that in the original MMMF publication (Morgan and Duzant, 2008), as well as in the Revised Morgan-Morgan-Finey paper (Morgan, 2001), equation 6 does not include the LD parameter. However, this would lead to the KELD being expressed in J mm -1 m -2 , instead of J m -2 . Hence, we highlight that a correct application of equation 6 must include leaf drainage depth (see Brandt, 1990;Morgan et al., 1998;155 Peñuela et al., 2018).
Saturation excess overland flow typically occurs in climates with low intensity precipitation and without a pronounced seasonal rainfall regime (Morgan and Duzant, 2008), specifically in areas with shallow soils and impermeable bedrocks (Beven, 2012). In the MMMF model, saturation excess runoff generation is assumed to occur when the mean daily rainfall exceeds the mean daily storage capacity of the soil (SC; 160 mm): Where MS is soil moisture at field capacity (% w w -1 ), BD is bulk density (Mg m -3 ), HD is effective hydrological depth (m) (i.e., the land cover dependent soil depth in which storage capacity controls the generation of runoff), and ETa/ETp is the ratio of actual to potential evapotranspiration. These parameters 165 represent an approximate average for the cropping season.
The annual runoff generation (Q; mm) is then estimated as a function of annual effective rainfall, the ratio between storage capacity and mean daily rainfall (PM, mm), and the slope length of the spatial modelling element (L; m): In the sediment phase, the annual detachment of soil particles by raindrop impact (ER kg m -2 ) and by surface 170 runoff (EQ; kg m -2 ) are calculated separately for the clay, silt, and sand texture classes, which are subsequently summed: Where KR is the detachability of the soil by raindrop impact (J m -2 ), T is the percentage of texture class i, ST is stone cover (proportion 0-1), KQ is the detachability of the soil by runoff (J m -2 ), GC is the average proportion of the soil covered by vegetation during the growing season (0-1), and S is slope angle (degrees). 175 Soil detachability values for each texture class i (clay, silt, and sand) are taken from Quansah (1982).
The immediate deposition of detached sediments (D; %) (i.e., the percentage of sediments not delivered to the runoff for transport) is estimated as a function of the average annual flow velocity, in our case for vegetated conditions (vv; m s -1 ), and the particle fall number (FN): Where g is the gravitational acceleration (9.81 m s -2 ), ø is the diameter of plant stems (m), NV is the number 180 of stems per unit area (number m -2 ), vs is the fall velocity for texture class i (0.00002 m s -1 , 0.002 m s -1 , and 0.02 m s -1 for clay, silt, and sand, respectively), and d is the hydraulic radius of the flow (0.005 m for unchanneled flow, 0.01 m for shallow rills, and 0.25 m for deeper rills). Again, in this case, land cover parameter values describe an average over the cropping season.
The total detached material delivered annually to transport (G; kg m -2 ) is modelled separately for each soil 185 texture class i: The annual transport capacity of the surface runoff (TC; kg m -2 ) is calculated as a function of annual runoff volume (Q; mm), slope, and the effect of plant cover/tillage on flow velocities, for each particle size class i: The average flow velocities for the actual soil conditions (va; m s -1 ), for the effect of tillage (vt; m s -1 ), and 190 for the standard bare soil condition (vb; m s -1 ) are calculated using the Manning equation: Where n is Manning's roughness coefficient. For the tilled conditions, Manning's n is estimated as a function of an implement-dependent surface roughness parameter (RFR) taken from Morgan (2005): The annual soil loss (SL; kg m -2 ) is calculated by comparing the annual transport capacity (TC; kg m -2 ) and the annual sediment delivered to the runoff (G; kg m -2 ), for each texture class i: 195 If the amount of sediment delivered to the runoff is greater than the transport capacity, the excess sediment will be deposited until G = TC. Such deposition is modelled using the settling velocities and fall numbers described in equation 14. The sediment balance becomes: Finally, the soil losses for the clay, silt, and sand texture classes are summed to produce total estimates of annual soil losses (SL; kg m -2 ): 200

Soil database
The soil profile data used in the model were retrieved from the UK SOILPITS dataset, which is one of many datasets held within the Land Information System (LandIS) operated by the Soil and Agrifood Institute at Cranfield University, UK (LandIS, 2022). The UK SOILPITS dataset represents a compilation of a series of soil profile surveys conducted across the UK since 1984. We selected only the profiles under 205 agricultural land cover, and which had complete information on the key soil properties used for modelling for all described horizons. Table 1 presents a descriptive summary of the data representing each whole soil profile; that is, data for each horizon from each profile has been bulked together. The profiles range in thickness from 0.22 m to 1.96 m (median depth is 0.60 m) and are typically composed by four characteristic horizons: an A, E, B, 210 and C horizon. More information about how each horizon was surveyed and differentiated in the field can be found in Hodgson (1997).  Figure 3 demonstrates the variability of key soil properties within the four characteristic soil horizons, compiling all soil profiles used in the dataset (by 'key', we mean those properties which are employed directly or indirectly as input variables in our model). There is considerable overlap between each horizon, largely due to the heterogeneity of soil types represented in the dataset. However, some distinctive patterns 220 can be discerned. For example, between the A and B horizon, bulk density tends to increase, while organic carbon tends to decrease. Some 79 profiles were also observed to have an E horizon directly below the A horizon. This was distinguished by the presence of a mineral layer with less organic carbon and clay content than the underlying B horizon indicating downward and/or lateral translocation into the subsoil. Another notable boundary lies between the B and C horizon, where the median volumetric water content at field 225 capacity reduces by more than 2.5 times. This may be reflective of some distinctive textural changes between these two horizons: the median sand content increases thrice, while both clay and silt decrease by more than five times. Organic carbon and rock fragment values undergone a square root transformation to improve the visualisation of the data.

Model implementation
Our modelling framework consists of an application of the MMMF model for each of the 265 soil profiles over a 500-year period (Fig. 4). Whilst the variability of the soil properties across profiles and horizons was incorporated into the model, all modelling units were parametrised the same for their climatic, land cover, and topographic variables. This was performed to test the sensitivity of modelled soil losses to erosion-240 induced changes to soil properties in different soil types. We selected rainfall-associated parameter values based on UK average climatic variables for the 1991-2020 period (Kendon et al., 2021). For the land cover parameters, which represent an approximate average over 245 the crop growing season, we took the guide values recommended by Morgan and Duzant (2008) for conventionally tilled winter cereals, except for plant height (PH) and ground cover (CG). This is because we found the suggested value of 1.5 m excessive for winter cereals, especially considering that PH in equation 6 does not represent the height of the top of the canopy, but rather the height of the fall of drops from the plant, which depends on the thickness of the canopy (Brandt, 1990). Assuming an average absolute These values represent an approximate average for the crop growing season.
In addition, we assumed a 10 m slope length and 6° slope gradient for the spatial element of the simulation unit. For the soil parameters, we used the measured properties from the soil profile database (Table 1). In addition, Ttexture-dependent parameters values were taken for the model guide, considering the soils' particle-size distribution. A Monte Carlo simulation with 100 iterations per year was included to provide 260 a forward error assessment of the model outputs (Beven, 2009). Model parameters were sampled from a normal distribution with a 10 % standard deviation to partially account for measurement errors and the uncertainty in parameter estimation. The constant parameter distributions used in all simulations are displayed in Table 2. In order to simulate soil thinning, As the MMMF model is applied with an annual resolution, we used the soil bulk density to convert the modelled soil losses from their native units to m yr -1 , after each timestepthe soil losses (SL; kg m 2 yr -1 ) estimated with the MMMF were converted into SLm (m yr -1 ) using soil bulk 270 density (BD; Mg m -3 ): Next, the model reduced the depth of the upmost soil horizon considering based on the amount of eroded soil in the previous timestepyear.
Since the model calculates soil losses separately for each particle size fraction, theThe soil texture of the 20 cm0.2 m plough layer was updated after each timestep k using a mass-balance model considering the 275 amount of fresh subsoil being incorporated into the plough layer and the selective removal of different particle size fractions. This requires estimating the mass of the original plough layer in timestep k (Mk; kg m -2 ) and the mass of the subsoil layer that will be incorporated by tillage in timestep k +1 (Ms; kg -2 ): Where BDk is the bulk density (Mg m -3 ) of the original plough layer in timestep k and BDs is the bulk density (Mg m -3 ) of the subsoil layer s being incorporated by tillage in timestep k+1. Please note that the estimation 280 of BD for each timestep is described in detail below. The masses of each particle size fraction i in the plough layer for timestep k, after selective removal (Mki; kg m -2 ), and in the subsoil layer being incorporated (Msi; kg m -2 ) are then calculated as: Where SLik (kg m -2 yr -1 ) is the soil loss for particle size fraction i in timestep k. 285 The percentage of each textural class Ti in the new plough layer for timestep k +1 is calculated as: Rock fragments were assumed not to be removed from the soil matrix, and therefore the stone cover (if present) (ST; %) undergoes a residual increment. , considering the soil losses from the previous timestepFor such, we used the volumetric percentage of rock fragments as a proxy for the stone cover model parameter: Where 200 is the volume (L) of the 0.2 m plough layer in 1 m 2 of soil. 290 . Of note, the selective removal of soil organic carbon associated to finer soil fractions was not simulated.
If the upper upmost horizon depth was greater than the 20 cm0.2 m plough depth, we assumed that no mixing would occur with the underlying horizonsmixed the eroded plough layer with fresh material from this same upmost soil horizon using the mass balance model described above (Eqs. 24-29) to recalculate soil texture and the percentage of rock fragments. Accordingly, soil organic carbon was assumed to remain 295 stable as the selective removal associated to finer soil fractions was not simulated. However, if the upper horizon was thinner than 20 cm0.2 m for any given timestep, the mass balance model (Eqs. 24-29) mixed the material in the plough layer with the underlying soil horizon. In this case, soil organic carbon (OC; %) values were also updated with a mass balance model: For every timestep, soil bulk density and soil moisture at field capacity were estimated using pedotransfer 300 functions (PTFs). We established the PTFs by fitting a linear regression of bulk density and soil moisture at field capacity as a function of sand (%) and organic carbon content (%) for the A horizons in our soil profile dataset (Fig. A1, A2). This was performed because i) we assumed that bulk density and soil moisture at field capacity would be affected by the changes in soil texture due to selective particle size removal; and ii) we presupposed that, as the subsoil horizons get incorporated into the plough layer and become closer 305 to the surface, their bulk density and soil moisture at field capacity would become more characteristic of an A horizon, due to tillage and organic matter input from plant biomass. Similarly, we established a pragmatic lower limit for soil organic carbon content for different soil texture classes, based on the lowest values observed in the A horizons from our dataset. That is, we assumed the organic carbon would not decrease indefinitely with soil truncation due to the continuous input from plant material and potentially other 310 farming practices. This assumption is based on observations that even heavily eroded arable soils typically contain some type of Ap horizon (Świtoniak, 2014). The successive soil thinning and mixing processes continued for 500 years or until the plough layer reached the end of the lowermost soil horizon (i.e., 20 cm0.2 m above the bedrock). We assumed that soil losses would outpace soil formation within the simulated system (see Evans et al., 2019), and since we did not focus on calculating soil lifespans, it was not necessary 315 to integrate soil formation rates into the model calculations.
The sensitivity of the simulated erosion rates to soil thinning was assessed with the correlation (Pearson's r) between soil truncation (i.e., the cumulative annual reduction in soil depth) and annual soil losses (here taken as the median of the Monte Carlo simulations per year). Soil profiles exhibiting a positive correlation (r > 0, p<0.00001) were assumed to display an accelerating erosion feedback trend for the simulation period, 320 whereas the ones with a negative correlation (r < 0, p<0.00001) were assumed to display a deaccelerating decelerating feedback trend. The remaining profiles were not considered sensitive to truncation and were assumed to present a stable erosion progression. Of note, we imposed a more restrictive significance level in order to screen out the profiles with very slight responses to soil truncation, considering this is a fully controlled modelling experiment. 325 In order to understand the processes driving the soil erosion feedback system, we used a random forest analysis to rank the importance of model parameters for predicting the changes in soil erosion rates between the first and final timesteps of the simulations. the above-mentioned erosion trends (i.e., deaccelerating, stable, accelerating). In this case, the differences between soil parameter values for the initial and final timesteps were used to predict the trends for each of soil the as explanatory variablesprofiles. All model 330 simulations and statistical analyses were performed in R (R Core Team, 2022), and the model code is available as supplementary material (Batista et al., 2022). Importantly, we did not consider all potential changes to the modelled systems. That is, we did not consider any feedbacks between soil thinning and crop development, nor the effects of climate change on rainfall, temperature, and farming practices. Although we are aware such factors would likely have an impact on 335 the model simulations, our aim here is to analyse the sensitivity of modelled soil losses to erosion-induced changes to soil properties. This involves making fixed assumptions about other system components. In addition, we would like to highlight that our model simulations should not be mistaken as projections of future erosion rates in Britain, due to all the above-mentioned reasons.

Considering the median of the simulations from all soil profiles, the temporal evolution of erosion rates 345
was characterised by an exponential decay function (Fig. 5), which meant erosion rates initially decelerated at a higher pace until reaching a horizontal asymptote at approximately 250 years into the model runs.
Based on the adjusted decay function, the initial soil loss rates (intercept = 0.81 kg m -2 yr -1 ) decrease in 6 % and 10 % over the first 50 and 100 years of the simulations, respectively; but only 3 % in the last 100 years. 350 Figure 4 illustrates the typical behaviour of these trends using data from two representative profiles with similar topsoil but different subsoil properties at the beginning of the simulation. The changes in simulated soil losses between the initial and final model timesteps per profile followed a normal distribution, which signified an average 20 ± 8 % decrease in erosion rates. The steepest declines in soil losses (maximum 46 %) typically occurred in soil profiles with a negative gradient in silt contents in 360 their subsoil horizons. Contrarily, the soil profiles with a stable erosion progression or with a gentler decay steepness were associated with the presence of silty substrata. Fig. 6 illustrates the typical behaviour of these trends using data from two representative soil profiles with similar topsoil but with different subsoil properties at the beginning of the simulations. It follows that simulatedThe changes in erosion rates were mostly largely explained by alterations in soil texture. This is demonstrated in Fig. 75, which displays a random forest importance ranking for predicting 370 the feedback trends in model outputs (i.e., deaccelerating, stable, accelerating differences in erosion rates between the initial and final timesteps of the simulations, for each soil profile). The random forest analysis described how soil erosion responses were highly influenced by variations in silt content and stone cover.
Changes in soil moisture at field capacity and bulk density had a lower impact on the modelled soil losses, whereas the simulated variations in bulk density had an almost null effect on model outputs (Fig. 75). 375 Figure 75. Random forest importance ranking for predicting the differences in erosion rates between the initial and final timesteps of the simulations, for each soil profile. the erosion trends for each soil profile (i.e., deaccelerating, stable, accelerating). Feature importance is represented by the relative increase of the mean squared error (MSE) of the random forest (i.e., how much removing the feature increased the 380 prediction error).
The soil erosion feedback system sensitivity of the simulated changes in soil erosion rates to the erosioninduced changes to soil properties can be further visualised comparing the difference in model parameter values with the variation in soil losses over 10-year rolling means (Fig. 86). Positive and negative changes in single parameter values did not yield consistent responses regarding the simulated soil losses (e.g., a 1 385 % decrease in sand content can lead to both accelerating and deaccelerating erosion rates over the rolling means). However, the direction of the changes in model parameters explains the feedback trends for the individual soil profilesthe general decelerating trend for the profiles (Fig. 86). For instance, the profiles with an accelerating erosion trend are characterised almost exclusively by increases in silt content within the 10-year rolling means. ContrarilyFor instance, for the profile trajectoriess were characterisedin which 390 soil losses show a deaccelerating trend over the whole simulation period have almost exclusively by decreasing contents of silt and increasing stone cover within the rolling means, which led to net decreases in soil losses over the whole simulation period. Clay and sand contents mostly increased within the 10-year rolling means, although such pattern was less pronounced compared to the silt and stone cover progressions.
Less important parameters, such as bulk density and soil moisture at field capacity, displayed a more 395 scattered centred pattern, along the y-axis regardless of the erosion trends (Fig. 86).  The main processes driving the soil erosion feedback systems were particle detachment by raindrop impact and silt sediment supply (R 2 = 0.8590 and 0.6684, respectively; Fig. 97), while changes in runoff amounts, detachment by runoff, and runoff transport capacity had a narrow effect on the simulated soil losses (R 2 = 0.01; Fig. 97). That is, only acute changes in discharge seem to have produced a sufficient response in detachment by runoff and in transport capacity to influence the net erosion rates (Fig. 97). 405 Figure 97. Hexagonal heatmaps and linear regression lines (p < 0.00001) of the changes in soil loss over 10-year rolling means for the soil profiles. Fill colours represent the number (n) of cases in each hexagon. Importantly, aAnnual runoff depths and soil losses were correlated in only 18731 (7112 %) of the 265 soil profiles (p < 0.00001); however, these correlations did not amount to causation. The positive correlations 410 (44 profiles; 17 %) occurred for instance when the uprise of clayey subsurface horizons led to a reduction in both runoff amounts (due to an increase in soil moisture storage capacity) and soil detachment (Fig. 8).
The negative correlations (143 profiles; 54 %) typically occurred due to selective erosion processes, soil organic carbon depletion, and the presence of sandy soil substrata in profiles with silty/loamy topsoils and sandy substrata. That is, aAs the silty material was removed, enriching the topsoil with sand and carbon-415 depleted subsoil, the erosion ratesparticle detachment by raindrop impact declined, whereas runoff amounts increased due to a reduction in soil moisture at field capacity. However, this increase in runoff was not sufficient to accelerate soil losses (Fig. 108).

Discussion
Our model simulations underline the strong interaction between current soil erosion dynamics and the erosion history of soil profiles. That is, the simulations demonstrate how different soils are likely to have contrasting responses to soil thinning, depending on the properties of the surface material, as well as those of the underlying soil horizons. In our database, most of the truncation-sensitive soil profiles presented a 430 deaccelerating decelerating feedback trends, which reflects both the characteristics of these soils and, importantly, our basic modelling assumptions.
For instance, as silt detachability in the MMMF model is assumed to be much higher than for other particle size fractions, silt was preferentially removed from the soil matrix. In addition, as silt contents typically remain stable or decrease in the subsurface horizons of the soil profiles in our database (Fig. 2), silt was 435 often not replenished by the underlying substrata being mixed into the plough layer. Such behaviour is overall consistent with empirical observations of selective particle size removal by interrill erosion processes (Koiter et al., 2017) and the progressive coarsening depletion of readily erodible material of from eroding surface soils (Parsons et al., 1991). However, it is worth highlighting that the soil detachability values used in MMMF have a limited empirical basis (Quansah, 1982), and the erodibility of the clay 440 particles might be poorly described due to the neglection of other variables, such as aggregate stability (Morgan and Duzant, 2008).
The residual accumulation of rock fragments further contributed to the reduction in erosion rates for the profiles with a deaccelerating trend in the model simulations, as stone cover is assumed to armour the land surface and to reduce soil detachment. Specifically, even a small number of rock fragments in the soil 445 matrix can disperse the overland flow and dissipate its energy, reducing rill incision and soil losses (Rieke-Zapp et al., 2007). Moreover, decreases in water erosion rates due to a residual increment of stone cover have previously been simulated by Govers et al. (2006), who warned, however, that the accumulation of rock fragments in arable soils depends on tillage practices. Notwithstanding, increases in rock fragment contents might also affect soil hydraulic conductivity and water holding capacity (Cousin et al., 2003), and 450 therefore influence runoff formation. None of these potential interactions were represented in our model simulations, and might therefore warrant further scrutiny.
The few soil profiles displaying an acceleratinga stable or slightly accelerating erosion trend were characterised by the presence of sandy or loamy surface horizons over a siltier substratum, which successively supplied the plough layer with readily erodible material as the original topsoil was removed 455 by erosion. Moreover, these profiles were defined by the absence of a surface stone cover and by subsoil horizons with very limited amounts of rock fragments. Although accelerating erosion trends were relatively infrequent in our datasetonly simulated for one profile in our modelling experiment, such behaviour might be expected in soils with highly incrementing silt contents in their C horizons, which are common, for instance, where loess is the parent material (Finke, 2012;Świtoniak et al., 2016). 460 Moreover, potential accelerating erosion trends might have been under-detected by the model simulations, as we did not consider how truncation can decrease the soil moisture storage capacity due to a reduction in soil depth, and how this might affect runoff generation (Dunne and Black, 1970;Morgan et al., 1984). That is, as soils become shallower, excess saturation overland flow might increase, depending on the permeability of the bedrock or the presence of an impeding horizon (Beven, 2012;Moraes et al., 2010). We 465 did not consider these processes in our simulation due to the absence of an explicit soil thickness parameter in the MMMF equations and the lack of data regarding the permeability of the soil profiles parent materials.
Although simulated increases in runoff formation were generally not sufficient to increase particle detachment by overland flow and soil losses in the model outputs (Figs. 9 and 10), different results would be conceivable at a landscape scale. That is, if upslope run-on is considered, increases in overland flow due 470 to soil truncation might lead to rill initiation in flow accumulation zones, which would largely increase the simulated erosion rates.
Furthermore, very different runoff responses to soil truncation can be expected in areas where infiltration excess is the dominant overland flow mechanism. While saturation excess is common under British edaphoclimatic conditions, most subhumid and semiarid zones are prone to infiltration excess runoff 475 formation, which is a process primarily controlled at the soil surface (Smith and Goodrich, 2005). Under such circumstances, erosion-induced changes to topsoil properties might have an even greater interaction with runoff generation. In particular, soil crusting slows down infiltrability and rapidly increases the overland flow, leading to greater soil losses (Le Bissonnais, 2016;Fiener et al., 2008;Veihe et al., 2001).
As the development of soil crusts is influenced by soil texture and organic carbon content (Fiener et al., 480 2011), the truncation-induced changes we have simulated would likely alter the susceptibility of surface soils to crusting and, consequently, to infiltration excess runoff formation. Another caveat in the model structure worth highlighting is that, differently to what is exhibited in Fig. model outputs8, an uprise of subsurface clayey material might in fact increase infiltration excess overland flow, due to the lower infiltrability and saturated hydraulic conductivity of heavy-textured soils (Hao et al., 2020), which are also 485 more susceptible to compaction, depending on their mineralogy (Bonetti et al., 2017;Hamza and Anderson, 2005).
In addition, accelerating erosion responses to soil truncation might have been more frequent if we assumed an increase in topsoil erodibility due to the lower aggregate stability and looser structure of the carbon depleted subsurface material being incorporated into the plough layer (Le Bissonnais, 2016;Doetterl et al., 490 2016;Tanner et al., 2018). That is, as the soil detachability coefficients from MMMF do not take soil organic carbon into account, the sensitivity of topsoil erodibility to soil truncation was likely downplayed.
For instance, (Radziuk and Switoniak, (2021) used an equation from the Erosion Productivity Impact Calculator (EPIC) model to estimate the erodibility of Luvisols at different truncation stages in Poland.
They found that more truncated soils had higher erodibility, due to decreases in soil organic carbon and 495 sand content in the eroded Ap horizons.
Our model simulations may have particularly underestimated erosion responses to soil thinning for the profiles which displayed an accumulation of sand and a depletion in clay and silt contents. This progressive coarsening should lead to lower soil water availability, and therefore, lower soil cover, lower crop biomass production, and less organic carbon input from plants. As sandier soils already have less carbon stabilisation 500 mechanisms (Doetterl et al., 2016), this would lead to even greater truncation-induced depletions in soil organic carbon and therefore to increases in erodibility (Auerswald et al., 2014;Fernández and Vega, 2018).
In general, as organic carbon was only an indirect model input via the PTFs for estimating bulk density and soil moisture at field capacity, the interplays between soil thinning, soil organic matter, and soil erodibility were likely underrepresented here. 505 Although not all the complex interactions between soil erosion and soil thinning could be described in our simulationsmodel, the simulated trajectories of erosion-induced changes to soil properties (Fig. 10) are consistent with field-observations. That is, measured data from eroded Ap horizons typically indicate (i) incrementing clay, sand, and rock fragment contents, (ii) increasing soil bulk density, (iii), depleting soil organic carbon contents, and (iv) decreasing soil water holding capacity with increasing erosion severity 510 (Lowery et al., 1995;Rhoton and Tyler, 1990;Stone et al., 1985;Strauss and Klaghofer, 2001). Moreover, the decelerating erosion trend in our model outputs corroborate the results from (Govers et al., (2006), who simulated an exponential decay in sediment production from arable hillslopes over a 50-year period following land use intensification, due to a rapid emergence of a soil with high stone cover. Importantly, the soil losses estimated from our model outputs (median = 0.67 kg m -2 yr -1 ; interquartile range = 0.58 -515 0.80 kg m -2 yr -1 ) are encompassed by the median and the upper quartile of measured erosion rates from arable land at plot scale in the UK (~ 0.1 -1 kg m -2 yr -1 ) (see Benaud et al., 2020). We compare model outputs to plot-scale measurements due to their similarity with our modelling spatial unit (i.e., 10 m erodinghillslope segment with 6 ° or 10.5 % slope gradient). It is worth highlighting that (i) the plot data described in Benaud et al. (2020) mostly derives from slopes below 10 %, (ii) measured erosion rates at field or 520 catchment scale in Britain are much lower than our model simulations. -well encompassed by, , Moreover, our model outputs further help to we can identify where more empirical evidence would help constraining our modelling assumptions and improve process representation. For instance, investigating how truncation affects the aggregate stability of surface soils (due to changes in soil texture, mineralogy, and organic carbon dynamics) might be an important step in order to further understand the feedbacks 525 between erosion and soil thinning. Interactions between soil truncation, water availability to plants, and crop growthand how these could in turn reduce soil cover and organic carbon inputmight also warrant further investigations. Similarly, exploring the responses of different runoff generation mechanisms to soil thinning and erosion-induced changes to soil properties should be beneficial to increase our understanding of the erosion feedback system. Modelling-wise, an important next step would be to adapt the framework 530 describe here into spatially distributed soil erosion models, in order to evaluate the effects of soil erosion feedback systems at landscape or watershed scale. This would allow for an unrevealing of different feedback processes at different positions in the landscape, which might be affected by erosion-induced changes to soil properties in different ways.
Finally, although our model simulations indicate that soil thinning has a decelerating effect on soil loss 535 rates for eroding hillslopes segments in the UK, the results also demonstrate that some of these soils' physical properties might become restrictive for agriculture before the profiles are completely truncated.
That is, some simulated erosion progressions create dense Ap horizons (max. bulk density of 1.78 Mg m -3 ) with excessive rock fragment contents (max. 36 %), very high clay (max. 89 %) or sand contents (max. 99 %), and low water holding capacity (min. soil moisture at field capacity of 0.14 % w w -1 ). For such 540 scenarios, rehabilitation techniques based on topsoil replacement might be necessary to sustain crop production (Schneider et al., 2021). These significant changes in soil properties, simulated in relatively short time periods (Fig. 10), further indicate how accelerated erosion in agricultural landscapes can potentially affect soil classification (Lewis and Witte, 1980;Olson and Nizeyimana, 1988;Świtoniak et al., 2016), which is important to consider when interpreting decades-old soil maps. Since soil erosion rates in 545 Britain are much lower compared to other regions of the world (Benaud et al., 2020), our model simulations represent a somewhat conservative scenario of erosion feedback systems. Erosion-induced changes to soil properties, and their implications for modelling, should be particularly relevant in areas under severe erosion rates.

Conclusions 550
Here we explored the soil erosion feedback system in 265 agricultural soil profiles in the United Kingdom.
In particular, we simulated how selective erosion processes and the incorporation of different subsoil horizons into the plough layer affected the erodibility of surface soils. We further analysed how these processes could change erosion rates during a 500-year period. We found that i) soil erosion rates in 39 99 % of the soil profiles were sensitive to soil truncation, ii) that the most of these truncation-sensitive profiles 555 essentially (75 %) displayed a deaccelerating decelerating trend, and iii) that changes in soil texture and stone cover were the main drivers of the modelled soil erosion feedbacks loops. Importantly, we found that different soils had drastically different simulated responses to soil truncation, depending on the properties of the surface material as well as those of the underlying soil horizons.
Ultimately, our findings highlight the dynamic nature of the soil as a three-dimensional body. That is, even 560 the so-called intrinsic properties of surface soils might change in a matter of decades in areas under accelerated erosion rates. MoreoverIn specific, erosion-induced changes to soil properties can have a significant impact on the rates with which soils are eroded, which in turn affects the calculation of soil lifespans and model-based erosion projections. To date, this soil erosion feedback system has been largely overlooked in soil erosion models, which might have led to spurious estimates of long-term erosion rates. 565 Therefore, understanding how erosion-induced changes to soil properties reverberate with erosion itself will be crucial for improving long-term model predictions, investigating the resilience of different soils to erosion disturbances, and for developing appropriate soil conservation strategies for a changing world. Figure A1. Pedotransfer function for estimating bulk density (BD) as a function of sand and organic carbon (OC) content. The regression was fit using only the data for the A horizons in the 265 agricultural soil profiles from the UK SOIL PITS dataset used in this study. Figure A2. Pedotransfer function for estimating soil moisture at field capacity (MS) as a function of sand and organic carbon (OC) content. The regression was fit using only the data for the A horizons in the 265 agricultural soil profiles from the UK SOIL PITS dataset used in this study.

Data availability
The soil profile data used in this research is restricted under licence. For further information on data accessibility, please contact nsridata@cranfield.ac.uk.

Author contribution
All authors were part of the conceptualisation of this research. PVGB wrote the model code and the 585 manuscript with contributions from all authors.

Competing interests
PF is a topical editor for SOIL. The authors have no other competing interests to declare.

Acknowledgments
PVGB was supported the project entitled "Dynamic Agricultural Weather Indicators for Extreme Weather 590 PVGB would like to thank Franz Conen and Diego Tassinari for the discussions about soil truncation and 595 soil science, and Leonardo Ferreira for the coding advice. The authors kindly thank Claudia Mignani for her comments on an earlier draft of this manuscript. We are also highly thankful to Isabella Teles for preparing Figure 1 of this manuscript.