Articles | Volume 9, issue 1
Original research article
23 Jan 2023
Original research article |  | 23 Jan 2023

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

Pedro V. G. Batista, Daniel L. Evans, Bernardo M. Cândido, and Peter Fiener

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 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 period of 500 years using measured soil properties from a diversified database of 265 agricultural soil profiles in the UK. For this, we adapted the Modified Morgan–Morgan–Finney 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 the uncertainty in model simulations with a forward error assessment. We found that modelled erosion rates in 99 % 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 displayed a decelerating trend, which depicted an exponential decay in erosion rates over the simulation period. This was largely explained by decreasing silt contents in the soil surface due to selective removal of this more erodible particle size fraction and the presence of clayey or sandy substrata. Moreover, the soil profiles displayed an increased residual stone cover, which armoured the land surface and reduced soil detachment. Contrastingly, the soils with siltier subsurface horizons continuously replenished the plough layer with readily erodible material, which prevented the decline of soil loss rates over time. Although our results are limited by the edaphoclimatic conditions represented in our data, as by our modelling assumptions, we have demonstrated how modelled soil losses can be sensitive to erosion-induced changes in soil properties. These findings are likely to affect how we calculate soil lifespans and make long-term projections of land degradation.

1 Introduction

Rates of soil erosion on agricultural land often exceed the rates at which new soil is formed (Evans et al., 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 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 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 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 a system might develop over time and under assorted conditions is an important step to proactively design and implement effective soil conservation strategies, as different soils are likely to be impacted by erosion in 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 models allow for a representation of multiple mechanisms that influence 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 (LISEM; De Roo et al., 1996), and the Morgan–Morgan–Finney model (MMF; Morgan, 2001; Morgan et al., 1984; Morgan and Duzant, 2008) have been 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 of different soil horizons down a soil profile is constant over the period of a model simulation. As upper soil horizons are removed by erosion, thereby 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 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 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 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 the applications of soil erosion models.

Here, we hypothesise that erosion-induced changes to soil properties affect soil loss rates over long-term periods. To evaluate such a hypothesis, we performed an exploration of the soil erosion feedback system by simulating 500 years of soil losses and surface runoff on 265 agricultural soil profiles in the UK. This allowed us to investigate how soil erosion rates respond to truncation-induced 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 a dynamic soil erosion feedback system will be crucial for improving the calculation of soil lifespans, providing future soil loss projections, and designing long-term soil conservation strategies.

Figure 1Modelling concept: selective soil erosion processes alter topsoil properties, which are then mixed with the underlying substrata as the soil profile thins. Updated soil properties for the plough layer are used as model inputs for the following time step.


2 Materials and methods

2.1 Concept

Our modelling concept is essentially a numerical thought experiment, in which land cover, agricultural management, climate, and topography parameters are held within a constant range, so that any changes in simulated soil losses and surface runoff over a period of 500 years are solely a result of changes in soil properties due to erosion processes (Fig. 1). To perform this experiment, we parameterised a soil erosion model using measured data from 265 agricultural soil profiles spread across the UK (Fig. 2). The abstract spatial scale of the simulations can be perceived as a pedon located on a conventionally 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 0.2 m plough layer (i.e. the average tillage depth in the UK, Townsend et al., 2016), continuously updating soil properties through mass-balance models and pedotransfer functions (PTFs, Fig. 1).

In order to implement our modelling concept, we adapted the Modified Morgan–Morgan–Finney (MMMF) model (Morgan and Duzant, 2008). The model was chosen due to its ability to simulate multiple erosion subprocesses, which is desirable for understanding the specific mechanisms responsible for developing erosion feedback systems. 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 model has a parsimonious parameter set, which facilitates model application using national soil survey datasets. Moreover, MMMF and its derivatives have provided acceptable predictions of annual soil losses for different soils, land covers, and testing sites in the UK (Morgan and Duzant, 2008; Peñuela et al., 2018; Smith et al., 2018). In the following section we provide a brief description of the basic MMMF equations (Sect. 2.2). We subsequently characterise the soil profile database used for the modelling (Sect. 2.3) and describe the model implementation, including the mixing of surface and subsurface horizons (Sect. 2.4).

2.2 MMMF operating equations

The MMMF is a process-oriented conceptual model running on an annual time step, 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 interception by the vegetation cover:

(1) P EF = P ( 1 - PI ) ,

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 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 h−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 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–Finney paper (Morgan, 2001), Eq. (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 Eq. (6) must include leaf drainage depth (see Brandt, 1990; Morgan et al., 1998; 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, generation of saturation-excess runoff is assumed to occur when the mean daily rainfall exceeds the mean daily storage capacity of the soil (SC; mm):

(8) SC = 1000 MS BD HD ET a ET p 0.5 ,

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 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):

(9) Q = P EF e - SC P M L 10 0.1 .

In the sediment phase, the annual detachment of soil particles by raindrop impact (ER; kg m−2) and by surface runoff (EQ; kg m−2) are calculated separately for the clay, silt, and sand texture classes, which are subsequently summed as follows:


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). 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 of stems per unit area (number m−2), vs is the fall velocity for texture class i (0.00002, 0.002 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 unchannelled 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 texture class i:

(15) G = n = 1 i E R i + E Q i 1 - D i 100 .

Figure 2Location of the 265 soil profiles used in this study. Data source: Land Information System (LandIS) (LandIS, 2022).

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:

(16) TC = n = 1 i v a v v v t v b T i 100 Q 2 sin S 10 - 3 .

The average flow velocities for the actual soil conditions (va; m s−1), for the effect of tillage (vt; m s−1), and 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):

(19) n = e ( - 2.11 + 0.03 RFR ) .

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:

(20) If TC i G i ; SL = G i .

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 Eq. (14). The sediment balance becomes

(21) If TC i < G i calculate G1 i = G i 1 - D i 100 , If TC i G1 i ; SL i = TC i ; if TC i < G1 i ; SL i = G1 i .

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):

(22) SL = n = 1 i SL i .

2.3 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 only selected the profiles under agricultural land cover and those that had complete information on the key soil properties used for modelling.

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 to 1.96 m (median depth is 0.60 m) and are typically composed by four characteristic horizons: an A, E, B, and C horizon. More information about how each horizon was surveyed and differentiated in the field can be found in Hodgson (1997).

Table 1Descriptive statistics of the soil properties from the 265 agricultural soil profiles (depths between 0.22 and 1.96 m) from the LandIS database used in the simulations.

Download Print Version | Download XLSX

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 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 soil moisture at field 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 3 times, while both clay and silt decrease by more than 5 times.

Figure 3Boxplots of the key soil properties for each horizon of the soil profiles used in this study. Horizons which were not classified, or which occurred less than 5 times in the dataset are not shown in the figure. Organic carbon and rock fragment values underwent a square root transformation to improve the visualisation of the data.


2.4 Model implementation

Our modelling framework consists of an application of the MMMF model for each of the 265 soil profiles over a period of 500 years (Fig. 4). Whilst the variability of the soil properties across profiles and horizons was incorporated into the model, all modelling units were parameterised the same for their climatic, land cover, and topographic variables. This was performed to test the sensitivity of modelled soil losses to erosion-induced changes to soil properties in different soil types.

Figure 4Flowchart of the modelling framework.


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 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 (GC). This is because we found that the suggested value of 1.5 m was excessive for winter cereals, especially considering that PH in Eq. (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 height of 0.61 m for wheat varieties in the UK (Berry et al., 2015), we adopted a baseline PH value of 0.4 m. This was also performed to constrain disproportionate estimates of leaf drainage kinetic energy, as trial model runs indicated a high sensitivity of the model outputs to parameter PH. Moreover, GC was calculated using the number of plants per unit area (NV) and the average diameter of plant elements at ground surface (ø; m), in order to avoid inconsistencies between sampled GC values and the remaining plant parameters.

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). Texture-dependent parameter values were taken from the model guide, considering the soils' particle size distribution. A Monte Carlo simulation with 100 iterations per year was included to provide 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.

Table 2Parameter values which were applied to all soil profiles and sampled in the Monte Carlo simulation.

 The parameter value was held constant during the Monte Carlo simulation.

Download Print Version | Download XLSX

In order to simulate soil thinning, the soil losses (SL; kg m2 yr−1) estimated with the MMMF model were converted into SLm (m yr−1) using soil bulk density (BD; Mg m−3):

(23) SL m = SL BD 10 - 3 .

Next, the model reduced the depth of the upmost soil horizon based on the amount of eroded soil in the previous time step.

The soil texture of the 0.2 m plough layer was updated after each time step k using a mass-balance model considering the 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 p in time step k (Mp; kg m−2) and the mass of the subsoil layer that will be incorporated by tillage in time step k+1 (Ms; kg m−2):


where BDp is the bulk density (Mg m−3) of the original plough layer in time step k and BDs is the bulk density (Mg m−3) of the subsoil layer s being incorporated by tillage in time step k+1. Please note that the estimation of BD for each time step is described in detail below.

The masses of each particle size fraction i in the plough layer for time step k, after selective removal (Mpi; kg m−2), and in the subsoil layer being incorporated (Msi; kg m−2) are then calculated as


where Tip is the percentage of each particle size fraction i in the original plough layer for time step k, Tis is the percentage of each particle size fraction i in the subsoil layer being incorporated by tillage for time step k+1, and SLik (kg m−2 yr−1) is the soil loss for particle size fraction i in time step k.

The percentage of each textural class Ti in the new plough layer for time step k+1 is calculated as

(28) T i k + 1 = M p i + M s i M p - SL k + M s 100 .

Rock fragments were assumed not to be removed from the soil matrix; therefore, the stone cover (if present) (ST; %) undergoes a residual increment. As such, we used the volumetric percentage of rock fragments as a proxy for the stone cover model parameter:

(29) ST k + 1 = ST p 200 + ST s M s BD s M p - SL k BD p + M s BD s ,

where 200 is the volume (L) of the 0.2 m plough layer in 1 m2 of soil.

If the upmost horizon depth was greater than the 0.2 m plough depth, we mixed 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 stable as the selective removal associated with finer soil fractions was not simulated. However, if the upper horizon was thinner than 0.2 m for any given time step, 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:

(30) OC k + 1 = OC p M p - SL k M p - SL k + M s + OC s M s M p - SL k + M s

For every time step, soil bulk density and soil moisture at field capacity were estimated using pedotransfer 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 (Figs. A1, A2 in the Appendix). 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 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 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. 0.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 to integrate soil formation rates into the model calculations.

Figure 5Soil erosion trends over 500 years of model simulations for the profiles from the UK SOILPITS dataset. The dark solid line represents the median of all profiles and simulations, whereas the dark- and light-grey shaded areas are 50 % and 95 % prediction intervals, respectively.


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, whereas the ones with a negative correlation (r<0, p<0.00001) were assumed to display a decelerating feedback trend. The remaining profiles were not considered sensitive to truncation and were assumed to present a stable erosion progression. It is of note that we imposed a more restrictive significance level to screen out the profiles with very slight responses to soil truncation, considering this is a fully controlled modelling experiment.

Figure 6Soil erosion trends over 500 years of model simulations for two representative profiles from the UK SOILPITS dataset. Coloured symbols are the median of the simulations per year, and the solid lines are local regression functions adjusted from the data.


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 time steps of the simulations. In this case, the differences between soil parameter values for the initial and final time steps were used as explanatory variables. All model 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 that such factors would likely have an impact on 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.

3 Results

From the 265 soil profiles in the UK SOILPITS database, 262 (99 %) displayed a significant correlation (p<0.00001) between soil truncation and annual soil losses, of which 261 presented a decelerating feedback trend (Pearson's r<0). Considering the median of the simulations from all soil profiles, the temporal evolution of erosion rates 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 there is only a decrease of 3 % in the last 100 years.

The changes in simulated soil losses between the initial and final model time steps 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 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. Figure 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 simulated changes in erosion rates were largely explained by alterations in soil texture. This is demonstrated in Fig. 7, which displays a random forest importance ranking for predicting the differences in erosion rates between the initial and final time steps 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 (Fig. 7).

Figure 7Random forest importance ranking for predicting the differences in erosion rates between the initial and final time steps of the simulations, for each soil profile. 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 prediction error).


The sensitivity of the simulated changes in soil erosion rates to the erosion-induced changes to soil properties can be further visualised by comparing the difference in model parameter values with the variation in soil losses over 10-year rolling means (Fig. 8). Positive and negative changes in single parameter values did not yield consistent responses regarding the simulated soil losses (e.g. a 1 % decrease in sand content can lead to both accelerating and decelerating erosion rates over the rolling means). However, the direction of the changes in model parameters explains the general decelerating trend for the profiles (Fig. 8). For instance, the profile trajectories were characterised 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 a 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 centred pattern along the y axis (Fig. 8).

Figure 8Hexagonal heatmaps relating the changes in model parameters to soil loss responses over 10-year rolling means for the soil profiles. Colours represent the number (n) of cases in each hexagon.


The main processes driving the soil erosion feedback systems were particle detachment by raindrop impact and silt sediment supply (R2=0.85 and 0.66, respectively; Fig. 9), while changes in runoff amounts, detachment by runoff, and runoff transport capacity had a narrow effect on the simulated soil losses (Fig. 9). 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. 9).

Figure 9Hexagonal heatmaps and linear regression lines (p<0.00001) of the changes in soil loss over 10-year rolling means for the soil profiles. Colours represent the number (n) of cases in each hexagon.


Annual runoff depths and soil losses were correlated in 187 (71 %) of the 265 soil profiles (p<0.00001); however, these correlations did not amount to causation. The positive correlations (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. The negative correlations (143 profiles; 54 %) occurred due to selective erosion processes, soil organic carbon depletion, and the presence of sandy soil substrata. That is, as the silty material was removed, enriching the topsoil with sand and carbon-depleted subsoil, particle 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. 10).

Figure 10Annual changes in model parameter values and simulated runoff depths and soil losses over 500 years for the soil profiles from the UK SOILPITS database. The solid dark line represents the median of all profiles and simulations, whereas the shaded dark- and light-grey areas are 50 % and 95 % prediction intervals, respectively.


4 Discussion

Our model simulations underline the strong interaction between current soil erosion dynamics and the erosion history of soil profiles. In particular, 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 soil profiles presented a decelerating feedback trend, 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 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 depletion of readily erodible material 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 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 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 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 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 a 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 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 only 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).

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, saturation-excess 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 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. For instance, if upslope run-on is considered, increases in overland flow due 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 the formation of infiltration-excess runoff, 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., 2011), the truncation-induced changes we have simulated would likely alter the susceptibility of surface soils to crusting and, consequently, to the formation of infiltration-excess runoff. Another caveat in the model structure worth highlighting is that, different to model outputs, 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 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., 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 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 mechanisms (Doetterl et al., 2016), this would lead to even greater truncation-induced depletions in soil organic carbon and therefore 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.

Although not all the complex interactions between soil erosion and soil thinning could be described in our model, 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 (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 period of 50 years 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-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 eroding hillslope segment with 6 or 10.5 % slope gradient). It is worth highlighting that (i) the plot data described in Benaud et al. (2020) are mostly derived from slopes below 10 % and (ii) measured erosion rates at field or catchment scale in Britain are much lower than our model simulations.

Moreover, our model outputs further help to identify where more empirical evidence would help constraining 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 between erosion and soil thinning. Interactions between soil truncation, water availability to plants, and crop growth – and how these could in turn reduce soil cover and organic carbon input – might 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 described 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 unravelling 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 rates for eroding hillslope 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 M 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 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 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.

5 Conclusions

Here, we explored the soil erosion feedback system in 265 agricultural soil profiles in the UK. 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 period of 500 years. We found that (i) soil erosion rates in 99 % of the soil profiles were sensitive to soil truncation, (ii) the truncation-sensitive profiles essentially displayed a decelerating trend, and (iii) 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 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 the so-called intrinsic properties of surface soils might change in a matter of decades in areas under accelerated erosion rates. Moreover, 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. 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.

Appendix A

Figure A1Pedotransfer 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 SOILPITS dataset used in this study


Figure A2Pedotransfer 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 SOILPITS dataset used in this study.


Code availability

The model code is available online at: (Batista et al., 2022).

Data availability

The soil profile data used in this research were restricted under licence. For further information on data accessibility, please contact

Author contributions

All authors were part of the conceptualisation of this research. PVGB wrote the model code and the paper with contributions from all authors.

Competing interests

At least one of the (co-)authors is a member of the editorial board of SOIL.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Pedro V. G. Batista would like to thank Franz Conen and Diego Tassinari for the discussions about soil truncation and Leonardo Ferreira for the coding advice. The authors kindly thank Isabella Teles for preparing Fig. 1 and Claudia Mignani for her comments on an earlier draft of this paper. We are also highly thankful to Enrico Balugani, Andres Peñuela Fernandez, and Joris Eekhout for their referee comments, which largely improved the quality of our work.

Financial support

This research has been supported by the German Federal Ministry of Food and Agriculture (Grant No. 28DK118B20) and by QR GCRF funding.

Review statement

This paper was edited by Nikolaus J. Kuhn and reviewed by Enrico Balugani, Joris Eekhout, and Andres Peñuela Fernandez.


Anache, J. A. A., Flanagan, D. C., Srivastava, A., and Wendland, E. C.: Land use and climate change impacts on runoff and soil erosion at the hillslope scale in the Brazilian Cerrado, Sci. Total Environ., 622, 140–151,, 2018. 

Anselmetti, F. S., Hodell, D. A., Ariztequi, D., Brenner, M., and Rosenmeier, M. F.: Quantification of soil erosion rates related to ancient Maya deforestation, Geology, 35, 915–918,, 2007. 

Auerswald, K., Fiener, P., Martin, W., and Elhaus, D.: Use and misuse of the K factor equation in soil erosion modeling: An alternative equation for determining USLE nomograph soil erodibility values, Catena, 118, 220–225,, 2014. 

Batista, P. V. G., Evans, D. L., Cândido, B. M., and Fiener, P.: Erosion Feedback System – Soil Thinning MMMF Model (2.0), Zenodo [code],, 2022. 

Benaud, P., Anderson, K., Evans, M., Farrow, L., Glendell, M., James, M., Quine, T., Quinton, J., Rawlins, B., Rickson, J., and Brazier, R.: National-scale geodata describe widespread accelerated soil erosion., Geoderma, 371, 114378,, 2020. 

Berry, P. M., Kendall, S., Rutterford, Z., Orford, S., and Griffiths, S.: Historical analysis of the effects of breeding on the height of winter wheat (Triticum aestivum) and consequences for lodging, Euphytica, 203, 375–383,, 2015. 

Beven, K. J.: Environmental Modelling: An Uncertain Future, Routledge, Oxon, ISBN 10: 0-415-46302-5, 2009. 

Beven, K. J.: Rainfall-Runoff Modelling, 2nd ed., John Wiley & Sons, Chichester, ISBN 13: 9780470714591, 2012. 

Bonetti, J. A., Anghinoni, I., Moraes, M. T., and Fink, J. R.: Resilience of soils with different texture, mineralogy and organic matter under long-term conservation systems, Soil Tillage Res., 174, 104–112,, 2017. 

Bot, A. J., Nachtergaele, F. O., and Young, A.: Land resource potential and constraints at regional and country levels, in: World Soil Resources Reports, edited by: FAO, 90, 1–114, Rome, 2000. 

Bouchoms, S., Wang, Z., Vanacker, V., and Van Oost, K.: Evaluating the effects of soil erosion and productivity decline on soil carbon dynamics using a model-based approach, SOIL, 5, 367–382,, 2019. 

Brandt, C. J.: Simulation of the size distribution and erosivity of raindrops and throughfall drops, Earth Surf. Proc. Land., 15, 687–698,, 1990. 

Ciampalini, R., Constantine, J. A., Walker-Springett, K. J., Hales, T. C., Ormerod, S. J., and Hall, I. R.: Modelling soil erosion responses to climate change in three catchments of Great Britain, Sci. Total Environ., 749, 141657,, 2020. 

Cousin, I., Nicoullaud, B., and Coutadeur, C.: Influence of rock fragments on the water retention and water percolation in a calcareous soil, Catena, 53, 97–114,, 2003. 

De Roo, A. P. J., Wesseling, C. G., and Ritsema, C. J.: Lisem: a single-event physically based hydrological and soil erosion model for drainage basins, I: theory, input and output, Hydrol. Process., 10, 1107–1117, 1996. 

Doetterl, S., Berhe, A. A., Nadeu, E., Wang, Z., Sommer, M., and Fiener, P.: Erosion, deposition and soil carbon: A review of process-level controls, experimental tools and models to address C cycling in dynamic landscapes, Earth-Sci. Rev., 154, 102–122,, 2016. 

Dunne, T. and Black, R. D.: An experimental investigation of runoff production in permeable soils, Water Resour. Res., 6, 478–490, 1970. 

Eekhout, J. P. C., Millares-Valenzuela, A., Martínez-Salvador, A., García-Lorenzo, R., Pérez-Cutillas, P., Conesa-García, C., and de Vente, J.: A process-based soil erosion model ensemble to assess model uncertainty in climate-change impact assessments, L. Degrad. Dev., 32, 2409–2422,, 2021. 

Evans, D. L., Quinton, J. N., Tye, A. M., Rodés, Á., Davies, J. A. C., Mudd, S. M., and Quine, T. A.: Arable soil formation and erosion: a hillslope-based cosmogenic nuclide study in the United Kingdom, SOIL, 5, 253–263,, 2019. 

Evans, D. L., Quinton, J. N., Davies, J. A. C., Zhao, J., and Govers, G.: Soil lifespans and how they can be extended by land use and management change, Environ. Res. Lett., 15, 0940b2,, 2020. 

Fernández, C. and Vega, J. A.: Evaluation of the rusle and disturbed wepp erosion models for predicting soil loss in the first year after wildfire in NW Spain, Environ. Res., 165, 279–285,, 2018. 

Fiener, P., Govers, G., and Van Oost, K.: Evaluation of a dynamic multi-class sediment transport model in a catchment, Earth Surf. Proc. Land., 33, 1639–1660, 2008. 

Fiener, P., Auerswald, K., and Van Oost, K.: Spatio-temporal patterns in land use and management affecting surface runoff response of agricultural catchments-A review, Earth-Sci. Rev., 106, 92–104,, 2011. 

Finke, P. A.: Modeling the genesis of luvisols as a function of topographic position in loess parent material, Quat. Int., 265, 3–17,, 2012. 

Govers, G., Van Oost, K., and Poesen, J.: Responses of a semi-arid landscape to human disturbance: A simulation study of the interaction between rock fragment cover, soil erosion and land use change, Geoderma, 133, 19–31,, 2006. 

Hamza, M. A. and Anderson, W. K.: Soil compaction in cropping systems A review of the nature, causes and possible solutions, Soil Tillage Res., 82, 121–145,, 2005. 

Hao, H., Wei, Y., Cao, D., Guo, Z., and Shi, Z.: Vegetation restoration and fine roots promote soil infiltrability in heavy-textured soils, Soil Tillage Res., 198, 104542,, 2020. 

Herbrich, M., Gerke, H. H., and Sommer, M.: Root development of winter wheat in erosion-affected soils depending on the position in a hummocky ground moraine soil landscape, J. Plant Nutr. Soil Sci., 181, 147–157,, 2018. 

Hodgson, J. M.: Soil Survey field handbook: describing and sampling soil profiles, Cranfield, Cranfield University, ISBN 0901128821, 1997. 

Hoag, D. L.: The intertemporal impact of soil erosion on non-uniform soil profiles: A new direction in analyzing erosion impacts, Agr. Syst., 56, 415–429,, 1998. 

Kendon, M., McCarthy, M., Jevrejeva, S., Matthews, A., Sparks, T., and Garforth, J.: State of the UK Climate 2020, Int. J. Climatol., 41, 1–76,, 2021. 

Koiter, A. J., Owens, P. N., Petticrew, E. L., and Lobb, D. A.: The role of soil surface properties on the particle size and carbon selectivity of interrill erosion in agricultural landscapes, Catena, 153, 194–206,, 2017. 

LandIS: The Land Information System,, last access: 18 March 2022. 

Le Bissonnais, Y.: Aggregate stability and assessment of soil crustability and erodibility: I. Theory and methodology, Eur. J. Soil Sci., 67, 11–21, 12311, 2016. 

Lewis, D. T. and Witte, D. A.: Properties and Classification of an Eroded Soil in Southeastern Nebraska, Soil Sci. Soc. Am. J., 44, 583–586,, 1980. 

Lowery, B., Swan, J., Schumacher, T., and Jones, A.: Physical properties of selected soils by erosion class, J. Soil Water Conserv., 50, 306–311, 1995. 

Merritt, W. S., Letcher, R. A., and Jakeman, A. J.: A review of erosion and sediment transport models, Environ. Model. Softw., 18, 761–799,, 2003. 

Montgomery, D. R.: Soil erosion and agricultural sustainability, P. Natl. Acad. Sci. USA, 104, 13268–13272,, 2007. 

Moraes, J. M., de Schuler, A. E., Dunne, T., Figueiredo, R. O., and Victoria, R. L.: Water storage and runoff processes in plinthic soils under forest and pasture in Eastern Amazonia, Hydrol. Process., 20, 2509–2526,, 2010. 

Morgan, R. P. C.: A simple approach to soil loss prediction: A revised Morgan-Morgan-Finney model, Catena, 44, 305–322,, 2001. 

Morgan, R. P. C.: Soil Erosion & Conservation, 3rd ed., Blackwell, Oxford, ISBN 1-4051-1781-8, 2005. 

Morgan, R. P. C. and Duzant, J. H.: Modified MMF (Morgan–Morgan–Finney) model for evaluating effects of crops and vegetation cover on soil erosion, Earth Surf. Proc. Land., 34, 613–628,, 2008. 

Morgan, R. P. C., Morgan, D. D. V., and Finney, H. J.: A predictive model for the assessment of soil erosion risk, J. Agric. Eng. Res., 30, 245–253,, 1984. 

Morgan, R. P. C., Quinton, J. N., Smith, R. E., Govers, G., Poesen, J. W. A., Auerswald, K., Chisci, G., Torri, D., and Styczen, M. E.: The European soil erosion model (EUROSEM): a dynamic approach for predicting sediment transport from fields and small catchments, Earth Surf. Proc. Land., 23, 527–544,< 527::AID-ESP868>3.0.CO;2-5, 1998. 

Nearing, M. A., Jetten, V., Baffaut, C., Cerdan, O., Couturier, a., Hernandez, M., Le Bissonnais, Y., Nichols, M. H., Nunes, J. P., Renschler, C. S., Souchère, V., and van Oost, K.: Modeling response of soil erosion and runoff to changes in precipitation and cover, Catena, 61, 131–154,, 2005. 

Nearing, M. A., Foster, G. R., Lane, L. J., and Finkner, S. C.: Process-based soil erosion model for USDA-water erosion prediction project technology, Trans. Am. Soc. Agric. Eng., 32, 1587–1593,, 1989. 

Olson, K. R. and Nizeyimana, E.: Effects of Soil Erosion on Corn Yields of Seven Illinois Soils, J. Prod. Agric., 1, 13–19,, 1988. 

Öttl, L. K., Wilken, F., Auerswald, K., Sommer, M., Wehrhan, M., and Fiener, P.: Tillage erosion as an important driver of in-field biomass patterns in an intensively used hummocky landscape, L. Degrad. Dev., 32, 3077–3091,, 2021. 

Panagos, P., Ballabio, C., Himics, M., Scarpa, S., Matthews, F., Bogonos, M., Poesen, J., and Borrelli, P.: Projections of soil loss by water erosion in Europe by 2050, Environ. Sci. Policy, 124, 380–392,, 2021. 

Papiernik, S. K., Schumacher, T. E., Lobb, D. A., Lindstrom, M. J., Lieser, M. L., Eynard, A., and Schumacher, J. A.: Soil properties and productivity as affected by topsoil movement within an eroded landform, Soil Tillage Res., 102, 67–77,, 2009. 

Parsons, A. J., Abrahams, A. D., and Luk, S.-H: Size characteristics of sediment in interrill overland flow on a semiarid hillslope, Southern Arizona, Earth Surf. Proc. Land., 16, 143–152,, 1991. 

Peñuela, A., Sellami, H., and Smith, H. G.: A model for catchment soil erosion management in humid agricultural environments, Earth Surf. Proc. Land., 622, 608–622,, 2018. 

Quansah, C.: Laboratory experimentation for the statistical derivation of equations for soil erosion modelling and soil conservation design, PhD Thesis, 1982. 

Quinton, J. N., Govers, G., Van Oost, K., and Bardgett, R. D.: The impact of agricultural soil erosion on biogeochemical cycling, Nat. Geosci., 3, 311–314,, 2010. 

Radziuk, H. and Switoniak, M.: Soil erodibility factor (K) in soils under varying stages of truncation, Soil Sci. Annu., 72, 134621,, 2021. 

Rhoton, F. E. and Tyler, D. D.: Erosion-Induced Changes in the Properties of a Fragipan Soil, Soil Sci. Soc. Am. J., 54, 223–228,, 1990. 

Rieke-Zapp, D., Poesen, J., and Nearing, M. A.: Effects of rock fragments incorporated in the soil matrix on concentrated, Earth Surf. Proc. Land., 32, 1063–41076,, 2007. 

Schneider, S. K., Cavers, C. G., Duke, S. E., Schumacher, J. A., Schumacher, T. E., and Lobb, D. A.: Crop responses to topsoil replacement within eroded landscapes, Agron. J., 113, 2938–2949,, 2021. 

Sharmeen, S. and Willgoose, G. R.: A one-dimensional model for simulating armouring and erosion on hillslopes: 2. Long term erosion and armouring predictions for two contrasting mine spoils, Earth Surf. Proc. Land., 32, 1437–1453,, 2007. 

Smith, H. G., Peñuela, A., Sangster, H., Sellami, H., Boyle, J., Chiverrell, R., Schillereff, D., and Riley, M.: Simulating a century of soil erosion for agricultural catchment management, Earth Surf. Proc. Land., 43, 2089–2105,, 2018. 

Smith, R. E. and Goodrich, D. C.: Rainfall Excess Overland Flow, Encyclopedia of Hydrological Sciences,, 2005. 

Sommer, M., Gerke, H. H., and Deumlich, D.: Modelling soil landscape genesis – A “time split” approach for hummocky agricultural landscapes, Geoderma, 145, 480–493,, 2008. 

Stone, J. R., Gilliam, J. W., Cassel, D. K., Daniels, R. B., Nelson, L. A., and Kleiss, H. J.: Effect of Erosion and Landscape Position on the Productivity of Piedmont Soils, Soil Sci. Soc. Am. J., 49, 987–991,, 1985. 

Strauss, P. and Klaghofer, E.: Effects of soil erosion on soil characteristics and productivity, Bodenkultur, 52, 147–153, 2001. 

Świtoniak, M.: Use of soil profile truncation to estimate influence of accelerated erosion on soil cover transformation in young morainic landscapes, North-Eastern Poland, Catena, 116, 173–184,, 2014. 

Świtoniak, M., Mroczek, P., and Bednarek, R.: Luvisols or Cambisols? Micromorphological study of soil truncation in young morainic landscapes – Case study: Brodnica and Chełmno Lake Districts (North Poland), Catena, 137, 583–595,, 2016. 

Tanner, S., Katra, I., Argaman, E., and Ben-Hur, M.: Erodibility of waste (Loess) soils from construction sites under water and wind erosional forces, Sci. Total Environ., 616, 1524–1532,, 2018. 

Townsend, T. J., Ramsden, S. J., and Wilson, P.: How do we cultivate in England? Tillage practices in crop production systems, Soil Use Manag., 32, 106–117,, 2016. 

Vanacker, V., Ameijeiras-Mariño, Y., Schoonejans, J., Cornélis, J. T., Minella, J. P. G., Lamouline, F., Vermeire, M. L., Campforts, B., Robinet, J., Van de Broek, M., Delmelle, P., and Opfergelt, S.: Land use impacts on soil erosion and rejuvenation in Southern Brazil, Catena, 178, 256–266,, 2019. 

Vanwalleghem, T., Gómez, J. A., Infante Amate, J., González de Molina, M., Vanderlinden, K., Guzmán, G., Laguna, A., and Giráldez, J. V.: Impact of historical land use and soil management change on soil erosion and agricultural sustainability during the Anthropocene, Anthropocene, 17, 13–29,, 2017. 

van der Meij, W. M., Temme, A. J. A. M., Wallinga, J., and Sommer, M.: Modeling soil and landscape evolution – the effect of rainfall and land-use change on soil and landscape patterns, SOIL, 6, 337–358,, 2020.  

Veihe, A., Rey, J., Quinton, J. N., Strauss, P., Sancho, F. M., and Somarriba, M.: Modelling of event-based soil erosion in Costa Rica, Nicaragua and Mexico: Evaluation of the EUROSEM model, Catena, 44, 187–203,, 2001. 

Willgoose, G. R. and Sharmeen, S.: A One-dimensional model for simulating armouring and erosion on hillslopes: 1. Model development and event-scale dynamics, Earth Surf. Proc. Land., 31, 970–991,, 2006. 

Short summary
Most agricultural soils erode faster than new soil is formed, which leads to soil thinning. Here, we used a model simulation to investigate how soil erosion and soil thinning can alter topsoil properties and change its susceptibility to erosion. We found that soil profiles are sensitive to erosion-induced changes in the soil system, which mostly slow down soil thinning. These findings are likely to impact how we estimate soil lifespans and simulate long-term erosion dynamics.