A millennium of arable land use – the long-term impact of tillage and water erosion on landscape-scale carbon dynamics

In the last decades, soils and their agricultural management have received great scientific and political attention due to their potential to act as a sink of atmospheric carbon dioxide (CO₂). Agricultural management has strong potential to accelerate soil redistribution, and, therefore, it is questioned if soil redistribution processes affect this potential CO₂ sink function. Most studies analysing the effect of soil redistribution upon soil organic carbon (SOC) dynamics focus on water erosion and analyse only relatively small catchments and relatively short time spans of several years to decades. The aim of this study is to widen this perspective by including tillage erosion as another important driver of soil redistribution and by performing a model-based analysis in a 200 km² sized arable region of northeastern Germany for the period since the conversion from forest to arable land (approx. 1000 years ago). The spatially explicit soil redistribution and carbon (C) turnover model SPEROS-C was applied to simulate lateral soil and SOC redistribution and SOC turnover. The model parameterisation uncertainty was estimated by simulating different realisations of the development of agricultural management over the past millennium. The results indicate that, in young moraine areas, which are relatively dry but have been intensively used for agriculture for centuries, SOC patterns and dynamics are substantially affected by tillage-induced soil redistribution processes. To understand the landscape-scale effect of these redistribution processes on SOC dynamics, it is essential to account for long-term changes following land conversion as typical soil-erosion-induced processes, e.g. dynamic replacement, only take place after former forest soils reach a new equilibrium following conversion. Overall, it was estimated that, after 1000 years of arable land use, SOC redistribution by tillage and water results in a current-day landscape-scale C sink of up to 0.66 ‰ yr⁻¹ of the current SOC stocks.

1 Introduction

Soils play an important role in the global carbon (C) cycle (Bellamy et al., 2005; Berhe et al., 2008; Lal, 2004) and have received great scientific (e.g. Amelung et al., 2020; Bellassen et al., 2022; Van Oost et al., 2007) and political attention as one of the cornerstones to tackle climate change, e.g. the 4 ‰ Initiative (Minasny et al., 2017), Article 3.4 of the Kyoto Protocol (United Nations Framework Convention on Climate Change, 1998), and the IPCC Special Report (Intergovernmental Panel on Climate Change, 2019).

A substantial loss of soil organic carbon (SOC) to the atmosphere before industrialisation is generally associated with the conversion of (natural) forest sites to cropland (Lal, 2019; Le Quéré et al., 2016; Sanderman et al., 2017). However, tillage operations and water erosion lead to an accelerated lateral redistribution of SOC within agricultural landscapes (Montgomery, 2007). In consequence, the spatial variability of SOC within the soils of arable landscapes increase, but this also creates complex interactions between changing SOC profiles, site-specific C mineralisation and sequestration, and potential losses to aquatic ecosystems (Doetterl et al., 2016). In a nutshell, (i) the removal of SOC-rich topsoil at erosional areas stimulates dynamic replacement of C via fresh photosynthates and the uplift of more reactive subsoil minerals (Harden et al., 1999; Stallard, 1998). (ii) During transport by different erosion agents, some SOC might be mineralised due to erosion-induced aggregate breakdown (Doetterl et al., 2016); however, this has a relatively short-lived effect due to the episodic nature of erosion processes (Van Oost and Six, 2023). (iii) At depositional sites, SOC is buried in deeper soil layers and hence is protected from fast mineralisation (Berhe et al., 2008; Rumpel and Kögel-Knabner, 2011; Stallard, 1998). (iv) In the case of water erosion, SOC will also partly enter aquatic ecosystems, where it is either buried in sedimentary deposits or mineralised during fluvial transport (Aufdenkampe et al., 2011; Battin et al., 2009).

The impact of soil redistribution on C dynamics has been assessed in various studies as reviewed in, e.g. Doetterl et al. (2016), Kirkels et al. (2014), and Van Oost and Six (2023). Such studies have often benefited from a strong modelling component, which has been explored by both process-oriented models and more conceptual approaches. Most process-oriented studies focused on water-erosion-prone micro-catchments, where field surveys regarding spatial patterns of SOC and erosion, or general erosion monitoring, can be used for model development and testing (e.g. Doetterl et al., 2012; Van Oost et al., 2005; Wilken et al., 2017a). The focus on small erosion-prone catchments has several implications: (i) results can only be partially generalised as these small-scale water erosion studies tend to be located in steeper areas; (ii) water erosion studies are often associated with loess-burden soils (e.g. Dlugoß et al., 2012; Li et al., 2007; Wilken et al., 2017a), which, although highly erodible, are also deep and display a low sensitivity to soil truncation regarding crop productivity; and (iii) the focus on water erosion makes it difficult to close the C balance as the fate of SOC after leaving the micro-catchment is open to debate (Aufdenkampe et al., 2011; Battin et al., 2009; Van Oost and Six, 2023). Apart from these process-oriented studies, there are also regional (Lugato et al., 2018; Nadeu et al., 2015) and even global (Naipal et al., 2018; Van Oost et al., 2007) model-based estimates of the effect of soil redistribution on SOC stocks, which are based on coupled conceptual soil erosion and C turnover models. These (water erosion) modelling studies give valuable insights for large areas but are mostly focused on current erosion and C turnover (e.g. Nadeu et al., 2015; Van Oost et al., 2007), while long-term effects of erosion-induced C dynamics after centuries or even millennia of land management are ignored. In consequence, model results might overestimate the effect of intensive modern agriculture as they typically only take the last 50 to 100 years into account (e.g. Dlugoß et al., 2012; Nadeu et al., 2015; Wilken et al., 2017b). Only a few of these regional studies addressed longer timescales (e.g. Bouchoms et al., 2017; Wang et al., 2017), which is a prerequisite to compare today's SOC profiles with model outputs in regions with a long agricultural land use history.

However, such long-term regional erosion and C turnover modelling is obviously challenged by the rapid decline in data accessibility and quality when moving back for centuries or even millennia. Apart from natural factors (e.g. climate, topography, soil cover, soil development), it is most challenging to reconstruct factors governed by agricultural practices (e.g. crop rotations, productivity, modification of soil cover, tillage methods). Moreover, estimates of initial (undisturbed) soil conditions (especially SOC stock profiles) are required to initiate long-term modelling. The existing long-term modelling studies (Bouchoms et al., 2017; Wang et al., 2017) used undisturbed soil profiles from long-term arable land, while, to our knowledge, studies accounting for the decline in SOC following conversion from forest to arable land in combination with erosion-induced C fluxes have not been carried out. Moreover, tillage erosion has been shown to be the main soil redistribution process in different parts of the world (e.g. Gerontidis et al., 2001; Lobb et al., 1995; Van Oost et al., 2003), and ignoring its effects on long-term C dynamics might lead to spurious conclusions.

Any large-scale and long-term study faces the challenge of assumption-based input data. Hence, the aim of large-scale and long-term modelling must be to simulate plausible patterns instead of process-based reconstruction. The aims of this study are (i) to simulate long-term changes (1000 years) in soil profiles in an agricultural landscape heavily affected by tillage erosion and less affected by water erosion, (ii) to perform a model-based soil redistribution and SOC turnover analysis for a larger area (about 200 km²) in order to avoid a systematic bias typically found in small-scale studies focussing on erosion processes in steep areas, and (iii) to model the long-term effect of soil redistribution when moving from a SOC-rich forest soil to a heavily eroded arable soil after 1000 years of cultivation.

2 Materials and methods

2.1 Study area

The study area covers an area of 196 km² and is located in the Quillow river catchment about 100 km north of Berlin in northeastern Germany (Fig. 1). It represents a typical ground moraine landscape formed after the retreat of the Weichselian glaciers ca. 20 000–15 000 years ago (shaded area in Fig. 1; Lüthgens et al., 2011). The area is characterised by a hilly topography with short summit–footslope distances (on average 35 m) and a mean slope (± standard deviation) of ca. 4.4 % ± 3.7 %. A large number of kettle holes that were formed by the delayed melting of bigger ice blocks (Anderson, 1998) are typical landscape elements. Drainage is only possible via sub-surface flow from the kettle holes (Lischeid et al., 2017). The kettle holes can be filled with water or (degraded) peat or are covered by colluvial material resulting from arable land use over centuries (Van der Meij et al., 2019).

Figure 1 The study area is located north of Berlin in the young moraine landscape of northeastern Germany, indicated by the grey area of the inset map (a). Location of the two test sites A and B (black letters), as well as of the four non-eroded soil profiles used for calibration (yellow circles) within the study area. Thereby, the yellow circle close to test site B represents two profiles. Topographic position index (TPI) in metres and kettle holes of test site A (b) and B (c) with 2 m contour lines (black lines). Positive TPI values indicate hilltops and knolls, while negative TPI values represent depressions.

The land cover of the study area is dominated by arable land and pasture (ca. 70 %), followed by wetlands and lakes (ca. 16 %), while only a small part is made up of forest (ca. 11 %) and settlements (ca. 3 %; Heinrich et al., 2018). Some parts of the study area have been used for agriculture since Neolithic times (ca. 5500 BCE; Behre, 2008), while it is assumed that agricultural land use became widespread approximately in 1000 CE (Behre, 2008; Herrmann, 1985). Intensive mechanisation of agriculture started in the second half of the 20th century. This was accompanied by a substantial increase in field sizes during the socialistic era of the German Democratic Republic (Bayerl, 2006), resulting in recent average field sizes of 21 ha (± 20 ha). The region is characterised by a relatively dry subcontinental climate, with an average annual air temperature of 9.4 °C and a mean annual precipitation of 466 mm (20-year average, 2001–2020, DWD meteorological station at Grünow; DWD Climate Data Center (CDC), 2018, 2021).

The soil pattern of the region follows the heterogeneity of Pleistocene deposits and has been strongly modified by soil redistribution over the past centuries (Deumlich et al., 2010; Sommer et al., 2008; Koszinski et al., 2013). Nowadays, non-eroded soils can only be found in ca. 20 % of the arable land, mainly at lower midslopes or flat plateaus. Thereby, extremely eroded soils occur at hilltops, ridges, and slope shoulders, while strongly eroded soils are found from slope shoulders to upper midslopes. Groundwater-influenced colluvial soils have developed at footslopes of closed depressions, which are often covering fossil peat (see more details in Öttl et al., 2021).

Within the Quillow study area, two agricultural fields (Fig. 1) were chosen to test the plausibility of the modelling results (i.e. current estimates of SOC stocks and patterns). They were selected because of existing SOC data from previous studies (Wehrhan and Sommer, 2021; Wilken et al., 2020). Test site A is located approximately in the centre of the study area, belonging to the village Christianenhof (53.3550° N, 13.6643° E), and has a size of ca. 4.4 ha and a mean slope of 8.7 % ± 3.9 %. Test site B is in the northeast of the study area, close to the village Holzendorf (53.3836° N, 13.7818° E), and has an area of ca. 20.5 ha and a mean slope of 5.5 % ± 2.9 %.

2.2 Modelling approach

The spatially explicit soil redistribution and C turnover model SPEROS-C (Dlugoß et al., 2012; Fiener et al., 2015; Van Oost et al., 2005) was applied for modelling tillage- (TIL) and water-induced (WAT) soil redistribution in the mesoscale study catchment over the past millennium. Thereby, lateral soil and SOC redistribution, SOC turnover, and vertical mixing within the profile (spatial and vertical resolutions of 5 m × 5 m and 10 increments of 0.1 m soil depth increments, respectively) were simulated. To isolate C fluxes that occur solely due to total soil redistribution (TOT is the sum of TIL and WAT), a reference run simulating C fluxes without lateral soil redistribution was calculated, i.e. vertical C fluxes solely due to C input and decomposition. Modelling soil redistribution and C dynamics required estimating and calibrating model input parameters and their uncertainty, as well as evaluating the model outputs. The single steps are described in detail in the following section.

2.2.1 Modelling soil redistribution and SOC dynamics

Tillage-induced (TIL) soil redistribution is calculated based on a diffusion-type equation developed by Govers et al. (1994) (Eq. 1). The net soil flux due to tillage Q_til (kg m⁻¹ yr⁻¹) can be written as

$$\begin{array}{}\text{(1)}& Q_{\mathrm{til}}=-k_{\mathrm{til}}\cdot s=-k_{\mathrm{til}}\cdot {\displaystyle \frac{\mathit{\delta }h}{\mathit{\delta }x}},\end{array}$$

whereby k_til is the tillage transport coefficient (kg m⁻¹ yr⁻¹), s is the local slope (%), h is the height at a given point of the hillslope (m), and x is the soil translocation distance in horizontal direction (m). The local tillage-induced soil redistribution rate E_til (kg m⁻¹ yr⁻¹) is calculated as

$$\begin{array}{}\text{(2)}& E_{\mathrm{til}}=-{\displaystyle \frac{\mathit{\delta }{Q}_{\mathrm{til}}}{\mathit{\delta }x}}=-k_{\mathrm{til}}\cdot {\displaystyle \frac{{\mathit{\delta }}^{\mathrm{2}}h}{\mathit{\delta }{x}^{\mathrm{2}}}}.\end{array}$$

Thereby, the intensity of the calculated erosion rates is determined by the k_til, and the change in slope gradient determines the spatial pattern of tillage-induced soil redistribution.

Water-induced (WAT) soil redistribution is calculated according to a slightly modified approach of the Revised Universal Soil Loss Equation (RUSLE; Renard et al., 1997), which is described in detail in Van Oost et al. (2000). A local transport capacity T_c (kg m⁻¹ yr⁻¹; Eq. 3) determines whether erosion, sediment transport, or deposition occurs. If the sediment inflow is higher than T_c, the excess is deposited, while the T_c is further routed downstream.

$$\begin{array}{}\text{(3)}& T_{\mathrm{c}}=k_{\mathrm{tc}}\cdot P\cdot C\cdot K\cdot R\cdot {\mathrm{LS}}_{\mathrm{2}\mathrm{D}}\end{array}$$

In the above, k_tc is the transport capacity coefficient (m); P, C, K, and R are the RUSLE factors; and LS_2D is a grid-cell-specific topographic factor calculated following Desmet and Govers (1996).

With regard to the SOC turnover model, SOC stocks are modelled for a soil profile with 10 soil layers of 0.1 m. The model equations describing the SOC depth profile and SOC decay are based on the Introductory Carbon Balance Model (ICBM; Andrén and Kätterer, 1997; Kätterer and Andrén, 1999). The ICBM considers a young (Y) and old (O) C pool with different turnover rates (k_Y=0.8 yr⁻¹, k_O=0.006 yr⁻¹). The fraction of the annual flux from Y to O is determined by the humification coefficient h. External environmental factors relating to climate and soils are combined in the factor r, and the mean annual C input to the soil is represented by the parameter i (Andrén and Kätterer, 1997). The dynamics of the two SOC pools are described by the following differential equations (Andrén and Kätterer, 1997):

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathit{\delta }\mathrm{Y}}{\mathit{\delta }t}}=i-k_{\mathrm{Y}}\cdot r\cdot \mathrm{Y},\text{(5)}& {\displaystyle \frac{\mathit{\delta }\mathrm{O}}{\mathit{\delta }t}}=h\cdot k_{\mathrm{Y}}\cdot r\cdot \mathrm{Y}-k_{\mathrm{O}}\cdot r\cdot \mathrm{O}.\end{array}$$

SOC turnover rates are assumed to decrease exponentially with depth due to a decreasing influence of environmental conditions (Eq. 6; Rosenbloom et al., 2001).

$$\begin{array}{}\text{(6)}& k_{\mathrm{Y}/\mathrm{O}z}=k_{\mathrm{Y}/\mathrm{O}z}\cdot e^{(-u\cdot z)}\end{array}$$

Annual C input i (g C m⁻² yr⁻¹) is derived from crops (i_c) and manure (i_m; Eq. 7). Thereby, i_c is made up of an aboveground and a belowground component. Crop residues are determined by the residue-to-aboveground-biomass (AGBM) ratio (Res). The fraction of C input from roots and rhizodeposition (p_z) at a given soil depth z (m) is defined by the root-to-AGBM ratio (RS). For i_c, a C content (C_cont) of 0.45 is used (Eq. 8; Tum and Günther, 2011).

$$\begin{array}{}\text{(7)}& {\displaystyle }i=i_{\mathrm{c}}+i_{\mathrm{m}}\text{(8)}& {\displaystyle }{i}_{\mathrm{c}}={\mathrm{C}}_{\mathrm{cont}}\cdot \left[\left(\mathrm{Res}\cdot \mathrm{AGBM}\right)+\left(p_z\cdot \mathrm{RS}\cdot \mathrm{AGBM}\right)\right]\end{array}$$

The C input into the soil is modelled by assuming an exponential root density profile (Gerwitz and Page, 1974; Van Oost et al., 2005), while manure input is only assigned to the plough layer (or layers). The allocation of total root dry matter to each soil layer z (m) was calculated according to a reference soil depth z_r=0.25 m (Van Oost et al., 2005) and a constant c that determines the proportion of the roots per soil layer (p_z; Eq. 9).

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \text{for}z\le z_{\mathrm{r}}:p_z={\displaystyle \frac{z}{z_{\mathrm{r}}+\frac{\mathrm{1}-e^{-c\left(\mathrm{1}-z_r\right)}}{c}}}\\ & \text{for}z>z_{\mathrm{r}}:p_z={\displaystyle \frac{z_{\mathrm{r}}+\left(\mathrm{1}-e^{-c\left(z-z_{\mathrm{r}}\right)}\right)/c}{z_{\mathrm{r}}+\left(\mathrm{1}-e^{-c\left(\mathrm{1}-z_{\mathrm{r}}\right)}\right)/c}}\end{array}\end{array}$$

The humification coefficient h is weighted according to the proportion of the source of i and depends on clay content cp (%) (Eq. 10; Kätterer and Andrén, 1999).

$$\begin{array}{}\text{(10)}& h={\displaystyle \frac{i_{\mathrm{c}}\cdot h_{\mathrm{c}}+i_{\mathrm{m}}\cdot h_{\mathrm{m}}}{i}}\cdot e^{\mathrm{0.0112}\cdot (\mathrm{cp}-\mathrm{36.5})}\end{array}$$

The temperature response factor r that accounts for the environmental influence on SOC decay is calculated with the following exponential Q₁₀ function (Kätterer et al., 1998; Van Oost et al., 2005):

$$\begin{array}{}\text{(11)}& r=Q_{\mathrm{10}}{\displaystyle \frac{T-\mathrm{5.4}}{\mathrm{10}}}.\end{array}$$

Thereby, r is estimated with a Q₁₀ value of 2.07 (Kätterer et al., 1998) and a temperature T (°C) calibrated for this study (as described below) and by correcting temperature by the annual mean temperature of central Sweden (+5.4 °C; Andrén and Kätterer, 1997).

With regard to the soil profile update, after every time step, the SOC profile is updated considering yearly soil loss and gain due to tillage and water erosion. At eroding sites, a fraction of SOC from the first subsoil layer equal to the thickness of the eroded layer is incorporated into the plough layer. Hence, erosion also leads to an uplift of soil into the deepest layer. At depositional sites, a fraction of the SOC from the plough layer is shifted downwards into a buried plough layer. The underlying subsoil layers are further buried according to the depth of the soil deposition in that time step (Dlugoß et al., 2012; Van Oost et al., 2005). Topographic change corresponding to soil redistribution was not taken into account to avoid blurring the mass balance of SOC. As this overstates the amount of buried C, we created two model runs of vertical C fluxes, either taking deeply buried C (> 1 m soil depth) into account or not to show the effect of deep C burial on the C balance of the whole study region.

To account for the development of tillage implements and practices (Fig. 2a, Table  in the Appendix), plough depth was updated with time but kept constant through periods without significant changes in historical plough development. Based on a literature review, we changed plough depth from 0.1 m for the first 800 years of the model simulations to 0.2 m for 1800–1900 CE and to 0.3 m for 1900–2000 CE (Fig. 2c). The yearly vertical C fluxes are then calculated following the profile update.

2.2.2 Model implementation

One of the major challenges in performing a model-based analysis of the impacts of 1000 years of soil erosion upon C fluxes in an area of ca. 200 km² is to estimate reasonable model inputs and to determine appropriate model parameters. Obviously, this is associated with large uncertainties and requires substantial simplifications. It is important to note that the model allows a reasonable analysis of the importance of soil redistribution for the C balance of the entire study area, but it is not expected to exactly mimic the current observational data.

In terms of model realisations, due to the uncertainties in the main model input parameters for the erosion modelling and to account for the varying importance of TIL and WAT, we created nine model realisations (R1–R9, Table 1). The realisations were simulated by a combination of low, medium, and high water erosion pathways, indicated by the minimum, the mean, and the maximum values of the C, K, and R factors, as shown in Fig. 2b, with a low, medium, and high tillage erosion pathway using the different k_til values from Fig. 2a. The theoretical background that led to the erosion pathways is explained in detail in the next paragraphs. It is important to note that the variation in TIL and WAT is set to the relative importance of tillage and water erosion in the region as determined in earlier studies (Öttl et al., 2021; Wilken et al., 2020). Due to the large computing requirements in simulating 1000 years for roughly 8 × 10⁶ raster cells of the entire study area, we only modelled the different realisations for the test fields (Fig. 1; together, roughly 10 × 10³ raster cells). The most plausible realisation (as defined below) was later used to model the entire study area.

Table 1 Combination of three scenarios of the tillage-induced (TIL) soil redistribution parameter k_til with three scenarios of the water-induced (WAT) soil redistribution parameters (C, K, and R factor) results in nine realisations.

Figure 2 Range of reasonable input parameters for modelling tillage- (a) and water-induced soil redistribution (b) and SOC dynamics (c) for the model period of 1000 years. The range of parameters in (a) and (b) (dashed-dotted for the lower range and dashed lines for the upper range; solid line represents the mean) is used in the different model realisations. Please notice the different scales of the y axes. Abbreviations: k_til is the tillage transport coefficient; AGBM refers to aboveground biomass; and C, R, and K factors are the factors of the RUSLE (Eq. 3). Information on data sources and an explanation of the parameters can be found in the text.

With regard to tillage-induced soil redistribution, a comprehensive literature review (comprising 47 original publications representing 137 k_til values; Table ) was performed to assess the tillage erosion intensity of different soil cultivation techniques. According to the land use history of the study region, the model period was subdivided into five periods representing different soil cultivation techniques.

For the first period (1000–1100 CE), the median k_til of 98 kg m⁻¹ (min. 9, max. 300 kg m⁻¹) was calculated from 23 k_til values for manual hoeing or the use of a simple ard (Table ). Although the medieval mouldboard plough was invented around 200–900 CE (Van der Meij et al., 2019), it was assumed that not every farmer practised mouldboard ploughing and that manual hoeing or simple ard ploughs were still widely used in the first period (Behre, 2008; Herrmann, 1985).

For the second period (1100–1800 CE), it was assumed that an increasing number of farmers used a rudimentary chisel or mouldboard plough drawn by an animal as the turning plough was introduced around 1000 CE (Behre, 2008; Herrmann, 1985). As not much further information is available until the end of the 18th century, we used a set of 30 k_til literature values representing an ard, chisel, or mouldboard plough drawn by a single animal. The median k_til of these studies is 88 kg m⁻¹ (min. 14, max. 300 kg m⁻¹; Table ).

The fourth period (1800–1900 CE) was characterised by the industrial revolution that tremendously changed the way land was managed. From 1800 onwards, the so-called “Ruchadlo”, a steep-turning tipping plough (Herrmann, 1985), and the “Mecklenburgischer Haken”, for seedbed preparation, were used (Behre, 2008). Both implements were pulled by animals (oxen or horses). A median k_til of 100 kg m⁻¹ (min. 14, max. 300 kg m⁻¹) was calculated from 15 k_til values for an ard, chisel, or mouldboard plough pulled by one or two animals (Table ).

The last period (1900–2000 CE) is characterised by the introduction of automotive tractors that were able to pull heavy implements; in consequence, the ploughing depths increased to 20–40 cm (Behre, 2008; Bork et al., 1998; Van der Meij et al., 2019). The median k_til of 234 kg m⁻¹ (min. 13, max. 900 kg m⁻¹) was calculated from 69 k_til values for tractor-pulled heavy machinery (early and recent chisel and mouldboard plough, harrow, cultivator, tandem disc, etc.; Table ).

With regard to water-induced soil redistribution, a range of C factor values were estimated to represent the changes in crop cover and/or management throughout the simulation period (Fig. 2b). As such, two different conditions were assumed: for the upper limit of the parameter space, it is assumed that the crop cover was low (i.e. high C factor) at the beginning of the simulation period due to relatively lower yields and high row spacing. For the lower limit, it is assumed that a much lower historic C factor might be reasonable due to a high vegetation cover related to a high proportion of weeds and grasses between the crops, which decreased over time due to improved weeding methods. To account for this uncertainty over time, we assumed that, at 1000 CE, the C factor might be either 50 % higher or lower than the current mean value. This range decreased according to a polynomial function (degree =3) until reaching ± 10 % of the current value in 2000 CE. The current mean C factor of 0.1 was calculated assuming a small-grain crop rotation (e.g. winter wheat – winter wheat – winter barley – winter rapeseed), typically applied under today’s conditions (Deumlich et al., 2002; Öttl et al., 2021; Schwertmann et al., 1987).

The soil erodibility factor K was assumed to remain constant throughout the simulation period and was calculated based on a soil group map (Bundesministerium der Finanzen, 2007; Rust, 2006), following the approach as described in DIN ISO (2017). The area-weighted mean K factor of 0.021 Mg ha h ha⁻¹ MJ⁻¹ mm⁻¹ was used as medium realisation (Fig. 2b). The lower and upper parameter values used for creating the model realisations are the area-weighted mean plus or minus the standard deviation of the K factor, respectively (0.021 ± 0.007 Mg ha h ha⁻¹ MJ⁻¹ mm⁻¹).

The rainfall erosivity factor R was calculated based on a long-term precipitation reconstruction for Europe (1500-2000 CE; Pauling et al., 2005) and the approach of Diodato et al. (2017) developed to estimate long-term erosion changes from historic precipitation data. As no precipitation data were available for the period 1000–1500 CE, the mean R factor of the available data (362 MJ mm ha⁻¹ h⁻¹ yr⁻¹ for 1500–2000 CE) was used as the mean for the whole modelling period (Fig. 2b). To address a potential range in the R factor, we used the mean ± 95 % confidence interval (362 ± 8.3 MJ mm ha⁻¹ h⁻¹ yr⁻¹; Fig. 2b).

For this study, a constant transport capacity coefficient k_tc of 150 m was used as this value was found to be suitable for cropland and a grid resolution of 5 m × 5 m (Dlugoß et al., 2012; Van Oost et al., 2003). The grid-cell-specific topographic factor LS_2D was calculated based on a digital elevation model (DEM; derived from airborne laser scanning; original spatial resolution of 1 m resampled to 5 m; Landesamt für Umwelt, Gesundheit und Verbraucherschutz Brandenburg und Landesvermessung und Geobasisinformation Brandenburg, 2012). The support practice factor P is 1.0 for all realisations for the whole modelling period as no erosion control practices are assumed.

Both water and tillage erosion are sensitive to field sizes and layouts, which, according to historic maps and later aerial photographs, substantially changed over time. As we could not reconstruct field layout over 1 millennium for the entire test area, it was decided to use recent field layouts. However, as the recent fields are very large, this leads to an underestimation of potential field border effects.

With regard to SOC turnover, to model SOC dynamics over 1000 years, SPEROS-C needs yearly estimates of C inputs based on AGBM and estimates of ploughing depths (Fig. 2c) as these variables change the C incorporation into the soil. To calculate the temporal evolution of AGBM (Fig. 2c), yield data for the federal state of Brandenburg from 1950 to 2018 CE (Federal Statistical Office, 1990–2018; Staatliche Zentralverwaltung für Statistik, 1956–1990) were combined with a long-term winter wheat yield dataset of the UK (1270–2014 CE; Ritchie and Roser, 2013). Yield was converted to AGBM by multiplying with the harvest index (HI; Donald and Hamblin, 1976), which determines the proportion of yield to total biomass for specific crop species. The HI was calculated with winter wheat grain and straw data from Brandenburg (mean HI = 0.449; Kuratorium für Technik und Bauwesen in der Landwirtschaft e. V. (KTBL), 1951). We assumed that, for every third year of the simulation, the AGBM would not be harvested in order to account for the so-called “three-field economy” (Rösener, 1985; Volkert, 1991), i.e. a crop rotation regime that has been commonly used in Germany since medieval times and in which a field is left fallow every third year.

While changes in AGBM and ploughing depths can be reasonably estimated based on existing data, it is hardly possible to estimate the temporal (or even spatial) variability of other model parameters used in ICBM, e.g. root:shoot ratio and manure application. Therefore, we used SOC depth distributions from four standard soil profiles representing undisturbed (i.e. non-eroded) arable soils in the study area (soil database of ZALF e.V. and Sommer et al. (2020); Fig. 1) and values from the literature as initial model parameters (Table 2). We assumed that, at the beginning of the modelling period, the soils had higher SOC stocks due to the conversion from forest to cropland at the onset of agricultural use. A mean SOC depth profile from three undisturbed soil profiles located in a forest in close proximity to the study area (Calitri et al., 2021) was used for the calibration of the first year of the model period (green line in Fig. 3).

Sommer et al. (2020)Van Oost et al. (2005)Van Oost et al. (2005)Verch (2020)Herbrich et al. (2018)Van Oost et al. (2005)Dlugoß et al. (2012)DWD Climate Data Center (CDC) (2018)Behre (2008)Herrmann (1985)

Table 2 Model input parameters for modelling SOC dynamics in agricultural soils as determined by Monte Carlo simulations (n=1000). The initial values are varied by ± 10 % for sampling, and the final value is the parameter set that yielded the highest Nash–Sutcliffe model efficiency. The references prove that the initial values and their ranges are valid assumptions.

The initial model parameters were later optimised to derive a model parameter set used for the entire modelling period (second year onwards). That is, first we varied the initial parameter values one at a time until they matched the observation data (i.e. combination of the four non-eroded SOC depth profiles; orange line in Fig. 3). Second, the obtained representative initial values for the observed SOC profile were used in a Monte Carlo simulation (n=1000). Each parameter was sampled from a uniform distribution in a range of ± 10 % around its initial value, which resulted in 1000 different modelled SOC–depth profiles. Hence, we considered not only the direct influence of each parameter on the model output but also the joint influence due to interactions between the parameters (Pianosi et al., 2016). The parameter set which yielded the highest Nash–Sutcliffe model efficiency (Nash and Sutcliffe, 1970) was selected for the final modelling (black stars in Fig. 3).

Figure 3 Depth profile of the mean observed SOC stocks for the forest (green) and agricultural soils (orange) with error bars of ± 1 standard deviation. The forest soils are used as the initial soil condition, while the calibration of the agricultural soils is used as the parameterisation for modelling the 1000 years (black stars).

In terms of model evaluation, a straightforward, traditional model-testing approach to the correspondence between observational data and model outputs after simulating 1000 years of soil redistribution and C turnover in a study area of about 200 km² is obviously hardly possible. There are no commensurate quantitative measurements of erosion available at this spatiotemporal scale, and it is not appropriate to directly compare soil truncation or SOC patterns of individual fields with a model output based on a parameterisation for the entire study area. As such, we focused on an investigative model evaluation approach (Baker, 2017), in which two independent datasets were used to evaluate the model’s capability to consistently represent the long-term erosion-induced C balance for the study area.

The first independent data used for model evaluation were derived from a remote sensing approach for identifying spatial patterns of severe soil erosion and soil truncation. Typical features in the study area are signs of soil truncation at hilltops, which most likely result from prolonged tillage erosion (Deumlich et al., 2010; Sommer et al., 2008). We qualitatively defined heavily eroded areas as the locations where bright subsoil material could be identified at the land surface by remote sensing images, which is indicative of the partial incorporation of glacial till into the plough layer due to extreme soil truncation (Fig. 4). The exposure of such subsoil material implies that ca. 1 m of soil was removed by erosion (Van der Meij et al., 2017).

Figure 4 Exemplary aerial photos of the study area showing eroded hilltops as indicated by the lighter soil colours. Notice that the aerial photo to the left was taken in 1953 (© ZALF e.V.), while the one to the right is from 6 September 2016 (© Google).

Heavily eroded areas can be straightforwardly detected by using remote sensing data for the entire catchment area. Therefore, 24 multispectral Sentinel-2 satellite images (ESA, 2015) acquired during bare-soil conditions were classified (support vector machine tool; ArcGIS version 10.7.1, ESRI, 2019). As the classification can only be performed for fields with bare-soil conditions at the time of satellite image acquisition, about 21 % of the study area (ca. 1.5 × 10⁶ raster cells) was classified, whereby 6 % of the study area was classified as heavily eroded (ca. 4.2 × 10⁵ raster cells). As hilltop erosion might also lead to a movement of the surface-exposed subsoil into the surrounding areas not affected by erosion, a buffer of −5 m was created on the inside of the area, resulting in 5.2 km² or 2.1 × 10⁵ raster cells being classified as heavily eroded. These raster cells are used to evaluate the consistency of the modelled erosion patterns, which have been shown to be dominated by tillage erosion in a previous study (Wilken et al., 2020).

The second source of independent model evaluation data was derived from measured SOC stocks for two different test sites in the study area. For test site A (Fig. 1b), plough layer SOC stocks are available from a nested sampling design (20 m × 20 m; see Wilken et al., 2020) carried out in 2018. The data were geostatistically interpolated using a kriging approach to a regular grid with 5 m × 5 m resolution. At test site B (Fig. 1c), the topsoil SOC stocks were derived from a regression analysis of ground truth SOC measurements (first 0.15 m) against multispectral images taken by a remotely piloted aircraft system (Wehrhan and Sommer, 2021). Both observed SOC patterns were compared to model outputs.

3 Results

The modelling approach resulted in a millennium of annual vertical C fluxes due to soil redistribution by tillage and water leading to a change in SOC stocks, as well as lateral C export by water erosion. According to the model evaluation (Sect. 3.2), realisation R4 (medium TIL, low WAT; Table 1) was used for the model analysis of the entire study area.

3.1 Results of modelling erosion-induced C flux dynamics for 1000 years

The modelled C fluxes without soil redistribution indicated a C loss to the atmosphere following conversion to arable land for about the first 800 years of the simulation (Fig. 5a, b; without (w/o) soil redistribution), with some interannual variability of vertical C fluxes due to the three-field economy (i.e. crops left on the field every third year). The resulting decrease in SOC stocks (Fig. 5e) was more pronounced for the first 500 years, nearly reaching a new equilibrium around 1700 CE. Soils turned into a slight C sink in the beginning of the 19th century after an abrupt change in modelled plough depth from 0.1 to 0.2 m. This changed again at the beginning of the 20th century, after the modelled plough depth was increased to 0.3 m, and especially after the end of the three-field economy, which substantially reduced the modelled soil C input (Fig. 5a, b; w/o soil redistribution). Finally, soils turned into a C sink again after 1950 due to the extremely increased yields (associated with a substantial increase in soil C input) following the end of the Second World War (Fig. 2c).

Figure 5 Temporal variation (1000 years) of annual vertical C fluxes, lateral C export, C balance, and SOC stocks modelled for the study region (R4) following conversion from forest to agricultural land (grey boxes). Plough depth was increased from 0.1 to 0.2 m and to 0.3 m in the years 1800 and 1900, respectively (vertical dotted lines). Until 1900, AGBM was left on the field every third year. Vertical C fluxes at erosional (a) and depositional sites (b), total lateral C export (c), soil-redistribution-induced C balance of all modelled fluxes (d), and mean soil SOC stocks of the entire study area (e; log-scaled y axis). Notice that negative vertical C fluxes indicate a loss of C to the atmosphere, while positive C fluxes indicate a gain in soil C.

Based on the model simulations with the representation of lateral soil redistribution processes, we found that, at erosional sites (Fig. 5a), the C loss to the atmosphere was less pronounced compared to sites without soil redistribution, and from about 1550 CE onwards, eroded soils became a C sink. From this time onwards, the C sink function steadily increased until 1900 CE, when it dropped due to changes in soil C input (i.e. end of the three-field economy). The C sink function at eroding positions increased more pronouncedly again in the 20th century compared to sites without soil redistribution until the end of the modelling period.

Moreover, the simulations considering lateral soil redistribution processes revealed a decrease in vertical C fluxes in depositional sites compared to the simulations without soil redistribution (Fig. 5b). Such differences became more pronounced over time as increasingly more C was stored in colluvial soils. It is important to note that deep C burial (> 1 m soil depth) did successively become more important for the vertical C fluxes over time (Fig. 5b), whereas the increase in C mineralisation in the upper 1 m of depositional sites (first-order kinetics) was reduced by deep burial. This is especially important in the case of prolonged severe deposition as more and more C-rich former topsoil is moving to depths below 1 m.

The simulated lateral C export due to water erosion (Fig. 5c; exp_tot) in the test region was very small, especially since a new SOC equilibrium in arable land is reached. For approximately the first 200 years of the simulation, the C export steadily decreased as the SOC content of the topsoil being eroded substantially declined following land use conversion. For about the next 600 years, the modelled export steadily increased due to rising water erosion rates as we assumed an increase in the C factor for the period in the low-water-erosion pathway (Fig. 2b; R4) used to analyse the entire study area. At the beginning of the 19th and 20th centuries, respectively, the increase in ploughing depth led to a reduction in the topsoil SOC content and hence less C export via water erosion.

Based on the conservative assumption that all lateral C fluxes leaving arable land were lost to the atmosphere (Fig. 5d; vf−exp _tot), soil redistribution resulted in a steadily increasing C sink of about 3 g C m⁻² yr⁻¹ at the end of the modelling period. This sink function is only slightly more pronounced if we assume that most C exported from arable land into neighbouring land uses or into kettle holes was stored and not mineralised, while only the proportion entering the stream network was lost to the atmosphere (Fig. 5d; vf−exp _rest). Especially in the last 2 centuries of the simulation, deep C burial became more important for the entire soil-redistribution-induced C balance of the study area (Fig. 5d; $vfw/o$ deep C burial).

The sum of all C fluxes with and without lateral soil redistribution is also mirrored in the changes of the mean SOC stocks of the study area (Fig. 5e). Here, it is interesting to note that, for about the first 300 years after conversion to arable land, soil redistribution led to a faster decline in SOC stocks compared to the system without soil redistribution. After about 500 years, the reverse was simulated, leading to a mean difference between mean SOC stocks with and without redistribution of about 0.42 kg C m⁻² (8.7 %).

3.2 Model evaluation

A comparison between modelled and remotely sensed soil redistribution patterns (Table 3) indicated that the most severely eroded sites were associated with tillage-induced (TIL) redistribution and total soil (TOT) redistribution. Overall, about 81 % of the areas classified as heavily eroded according to the remote sensing approach correspond to the modelled erosion class. On average, those areas show a modelled soil loss of −0.23 mm yr⁻¹ (R4), most of which was caused by tillage erosion (Table 3).

Table 3 Agreement between modelled WAT, TIL, and TOT erosion classes and remote-sensing-derived erosion classification. Note that a threshold of −0.05 m erosion per year was used to exclude areas with minimal erosion after modelling 1000 years of soil redistribution. The area classified from the remote sensing data represents about 21 % of the entire study area. Within the classified area, about 28 % is heavily eroded (about 4.2 × 10⁵ raster cells).

A comparison of modelled topsoil SOC stocks against observed topsoil SOC patterns based on the nine model realisations shows that the quality of the results is mostly determined by the differences in tillage erosion intensity (Fig. 6). The best results with respect to the used goodness-of-fit parameters can be reached for the medium- (R4–6) and high-tillage-erosion realisations (R7–9; Fig. 6). In contrast, WAT plays only a minor role in explaining the spatial distribution of SOC. It is important to note that, especially in the case of test site B, where topsoil SOC stocks are estimated with a remote sensing approach, the model substantially underestimates the SOC contents. Taking this into consideration while also trying to perform a somewhat conservative estimate of the extent of tillage erosion, we used realisation R4 (medium TIL, low WAT) for the model analysis of the entire study area. This realisation resulted in a Spearman's rank correlation coefficient of 0.83 and 0.82 for test sites A and B, respectively.

Figure 6 Modelled versus observed topsoil (first 0.1 m) SOC stocks for the two test sites, A (circles) and B (triangles), and the nine realisations (in panels; Table 1). Data are grouped into classes of total soil redistribution ranging from extreme erosion (≤ −1 m, red) to high deposition (≥ 1 m, blue). Error bars indicate the 95 % confidence interval of the mean per class. Black lines show the regression of the classified data (solid for A, dashed for B) with the respective adjusted coefficient of determination (R²; ns= p value ≥ 0.05, ^* p value < 0.05, and ≥ 0.01, ${}^{**}$ p value < 0.01).

4 Discussion

4.1 Challenge of long-term soil redistribution and C turnover modelling

Understanding current agricultural soil–landscape relations requires a consideration of the long-term soil change as today's soil and SOC patterns cannot be explained by the short-term soil redistribution history. Our results demonstrated that long-term soil redistribution processes in agricultural landscapes are particularly important in the Quillow catchment. Although soil redistribution in the study area increased with the intensive agricultural mechanisation that has been taking place since the 1960s (Frielinghaus and Vahrson, 1998), this does not explain the observed erosion rates and patterns in the area (Wilken et al., 2020), especially at slope shoulders, where signs of tillage erosion are clearly visible in aerial photographs from the 1950s (Fig. 4, left). A comparison between our results with typical soil truncation and accumulation rates for the study area (Van der Meij et al., 2017) shows that it is necessary to consider the past millennium (i.e. since the beginning of agricultural management) to understand the landscape C dynamics.

However, any long-term and, in particular, landscape-scale modelling approaches are subject to considerable uncertainties. Here, we did not intend to mimic detailed observational data of lateral soil fluxes (which are, in any case, not available at a commensurate temporal resolution in relation to our model outputs) from individual sites of the 200 km² study area with a high degree of accuracy and precision. On the contrary, our investigative model evaluation approach was focused on testing the model’s consistency in simulating long-term, landscape-scale spatial patterns of soil truncation and SOC stocks while partially representing the uncertainties associated with parameter estimation in such an ambitious modelling experiment. As such, a set of model realisations (Figs. 2, 6) that combined different soil redistribution assumptions was considered. The entire study area was ultimately analysed following the model realisation R4 (i.e. medium tillage and low water erosion), which could explain 69 % and 43 % (see R² in Fig. 6) of the current spatial pattern of SOC stocks in test sites A and B (Fig. 1). This leads to an underestimation of the mean SOC stocks by 40 % and 20 % in the topsoil of test sites A (50 cm soil depth) and B (plough layer), respectively. Importantly, the model outputs displayed a high agreement (81 %) with independent data used for estimating areas of severe soil truncation.

Overall, these results are encouraging, considering that (i) we only calibrated C turnover parameters, while the tillage and water erosion components of the model were applied “blindly” to derive a set of plausible realisations for the whole study area; (ii) the model was parameterised to represent the average conditions in the entire study area, not accounting for the anyway unknown specific land use and management history of the individual test sites; and (iii) in the relatively rare cases in which soil erosion models have been tested against independent spatial data, results have generally shown a poor agreement with observational data (Batista et al., 2019). As such, our modelling outputs are consistent with independent lines of evidence of related phenomena and with our current understanding of long-term soil and SOC redistribution processes at the landscape scale. This corroborates the usefulness of the employed modelling approach for elucidating soil redistribution and C dynamics in the study area over the last 1000 years.

4.2 Long-term soil redistribution and C dynamics

This model-based analysis of the long-term, landscape-scale effects of soil redistribution following land conversion from forest to arable land upon C dynamics extends previous studies that mostly combine soil redistribution with SOC turnover over shorter time periods and smaller areas and that were based on soils that are already in C equilibrium due to long-term arable use (e.g. Dlugoß et al., 2012; Nadeu et al., 2015; Wilken et al., 2017b). Taking the conversion from forest to arable land into account clearly indicates that time since conversion is essential for the understanding of soil-redistribution-induced C fluxes. This was not included in previous long-term modelling studies of larger areas (e.g. Bouchoms et al., 2017; Wang et al., 2017) but only applied to an artificial topographic setting (2.25 ha; van der Meij et al., 2020). Our results demonstrate that there is no dynamic replacement at erosional sites as long as topsoil soils still lose C following conversion from SOC-rich forest to SOC-depleted arable soils. This is particularly important as dynamic replacement is assumed to be one of the key processes for a potential C sink function of soil erosion (Doetterl et al., 2016; Harden et al., 1999). Within our simulation, it took about 500 years until eroded soils in the study region started to act as a C sink (Fig. 5a). This period would be substantially shorter in smaller, more erosion-prone catchments where SOC-rich topsoil from former forested areas is lost faster (Dlugoß et al., 2012; Juřicová et al., 2023; Wilken et al., 2017b). This result underlines that it is essential to model entire landscapes instead of upscaling conclusions from small-scale studies.

Erosion-induced SOC loss and its partial deposition are most pronounced shortly after land conversion as the topsoil is still rich in SOC. Therefore, results from studies in regions where arable land was established centuries ago (e.g. Dlugoß et al., 2012; Juřicová et al., 2023; Nadeu et al., 2015) might not allow us to draw general conclusions for regions where land conversion happened recently. This corroborates the argument from Van Oost and Six (2023) that our understanding of coupled erosion and C turnover processes is strongly biased towards humid and/or temperate settings, where land conversion mostly occurred centuries ago, while little is known for regions with ongoing land conversion, often located in tropical regions (Song et al., 2018).

4.3 Tillage-induced soil redistribution and C dynamics

4.3.1 Tillage as the main driver of the erosion-induced C pump

Within our study area, tillage erosion was demonstrated to be a critically important process dominating the catchment's C balance and the C sink function induced by soil redistribution. Water erosion cannot be neglected due to extreme events that are responsible for crop losses, high sedimentation rates, and off-site damage (Frielinghaus et al., 1992; Frielinghaus and Schmidt, 1993). However, as illustrated by the historical aerial photograph in Fig. 4, tillage-induced soil redistribution in this area is dominating and not only important since the introduction of heavy machinery 70 years ago (Van der Meij et al., 2017; Wilken et al., 2020; Winnige et al., 2003). In addition, tillage is known to further increase the susceptibility of arable soils to water erosion (Lobb et al., 1995) due to its effect on soil microporosity and changes in surface roughness (Poesen and Govers, 1985). The impact of tillage outcompetes soil redistribution by water due to specific conditions in the study area. The farming structures in the post-socialistic study area (large field sizes) perform tillage on a high optimisation level. Within a typical 5-year crop rotation (rapeseed – winter wheat – maize – winter wheat – winter barley), inversion tillage is only applied once between two small-grain cereals, while chisel ploughing is applied in all other years. The rationale is that chisel ploughing can be applied faster and requires less mechanical force (i.e. kinetic energy) compared to inversion ploughing, which reduces the time, effort, and fuel consumption (Dumanski et al., 2006; Helsel, 2007). Furthermore, the hummocky topography of the young morainic study area shows a short summit–footslope distance that benefits tillage erosion, which does not increase with slope length, such as in the case of water erosion. This characteristic topography also leads to large depositional areas (41 % of the study area in comparison to 25 ± 7 % in a global estimate of Van Oost et al., 2007) that favours C burial and sequestration. In addition, there is a low hydrological and sedimentological connectivity to the river system in the study area. Only 5 × 10⁻⁴ t ha⁻¹ of eroded soil is exported by water to the river system, while ca. 200 times more C (0.1 t ha⁻¹) is buried in kettle holes. Even if this limited hydrological connectivity benefits the water-redistribution-induced C sink function, tillage erosion is still dominant as erosive rainfall appears only at a relatively low frequency in this region (Deumlich, 1999; Wilken et al., 2018).

4.3.2 Recent developments in tillage-induced soil redistribution

Tillage erosivity partly decreased due to the introduction of pesticides for weed control, which reduced the relevance of inversion tillage over the past decades (Lobb et al., 2007). Also, in our study area, non-inversion conservation tillage receives more attention and is already applied to 47 % of the cropland area (44 % conventional tillage and 0.06 % no-till; Statistisches Bundesamt (Destatis), 2017). Nevertheless, it needs to be mentioned that conservation tillage focuses on water erosion mitigation, while tillage erosion is not accounted for. A recent study demonstrated that soil tillage by chisel ploughing leads to substantially more tillage erosion than inversion implements (Öttl et al., 2022). Hence, tillage practices show high spatiotemporal variation, which is fused to a large sensitivity of tillage erosion predictions. The sensitivity is reflected in the erosion pathways, whereby the difference between the C balance of the low- and high-tillage-erosion pathways (C balance R8–R2) is higher than the difference in the C balance without and with consideration of deep C burial, respectively. Hence, the C balance and corresponding sequestration potential of agricultural soil systems is mainly driven by individual farmers’ decisions.

4.4 The way ahead for long-term and large-scale soil redistribution and C dynamics modelling

It is evident that long-term and large-scale simulations are needed to gain an understanding of C dynamics, not only for scientific purposes but also to find adapted management strategies to increase soil C sequestration. From our perspective, the implementation of the following three processes would substantially increase the simulation quality of coupled soil redistribution and C turnover models.

4.4.1 Keeping track of topographic change by soil redistribution

The model does not account for topographic change related to soil redistribution (i.e. DEM update). For shorter temporal scales (ca. 50–100 years; e.g. Dlugoß et al., 2012; Nadeu et al., 2015; Wilken et al., 2017b), the topographic change has a limited impact, but for a modelling period of 1000 years, neglecting DEM updates affects lateral and vertical C dynamics. In a tillage-erosion-dominated study area like the Quillow catchment, both erosion and deposition processes will be substantially overestimated at individual raster cells (erosion: slope shoulders; deposition: footslopes and field borders). This is due to a constant erosion and deposition pattern, which becomes more relevant towards the end of the simulation period. This means that severe erosion is simulated for a smaller spatial area than that in which it would take place in reality. As a result, at erosional sites, substantial dynamic replacement is calculated for a limited number of raster cells, and SOC is more likely to be buried below 1 m at severe depositional sites. The latter is especially critical if the modelled deposition is large enough that deposited C-rich topsoil reaches soil layers below 1 m, where it is assumed that SOC is stable in time (Rumpel and Kögel-Knabner, 2011). Hence, taking the topographic change corresponding to soil redistribution into account would be an important step forward to improve the quality of soil patterns.

4.4.2 Plant feedback on soil degradation

Coupling the impact of soil redistribution against plant growth would be a great step towards a better representation of C dynamics in disturbed landscapes. A cornerstone for a landscape to function as a C sink is dynamic replacement of eroded C by fresh biomass C due to the uplift of unsaturated reactive minerals (Doetterl et al., 2016; Harden et al., 1999). However, this calls for constantly high yields and corresponding C input at eroding landscape positions (Doetterl et al., 2016; Van Oost and Six, 2023). As severe long-term soil erosion typically causes declining yields (e.g. Bakker et al., 2004; Den Biggelaar et al., 2001; Herbrich et al., 2018), which was also demonstrated in the study area (Öttl et al., 2021), C input is overestimated at erosional areas. On the other hand, C input is underestimated at depositional areas due to more favourable growing conditions (Öttl et al., 2021; Papiernik et al., 2005; Heckrath et al., 2005), which attenuates overstating the C sink term (Öttl et al., 2021; Quinton et al., 2022).

4.4.3 SOC burial in deeper soil layers (< 1 m)

Long-term soil redistribution following land conversion from natural forest to arable land leads to deep burial of SOC (< 1 m; Hoffmann et al., 2013). In our modelling approach, the assumptions regarding the stability of SOC buried below 1 m are of tremendous importance in the range of soil-redistribution-induced C fluxes (Fig. 5d). Assuming that all SOC allocated below 1 m is stabilised, the overall soil-redistribution-induced current-day C sequestration potential would lead to an increase in SOC stocks of 0.66 ‰ yr⁻¹. However, long-term modelling of SOC turnover in these deep layers is challenging due to the generally limited knowledge of SOC turnover in deep soils (Rumpel and Kögel-Knabner, 2011) and the fluctuating stagnic soil conditions partly associated with landscape positions where soil is deposited.

5 Conclusions

In this study, the long-term (1000 years) effect of soil redistribution upon C fluxes and SOC stocks was modelled in a study region of 200 km² in northeastern Germany. Differently to earlier studies focussing on erosion-induced C fluxes, we included the change in SOC stocks following conversion from a natural forest to arable land and accounted for changes in agricultural practices and production over time.

The modelling results from a study area representing ground moraine landscapes as typically found in northern Europe, Asia, and North America indicate that soil redistribution in such regions is resulting in a slight C sink, increasing the landscape-scale SOC stocks by 0.66 ‰ yr⁻¹ as compared to an area without erosion. This sink function mostly results from tillage-induced soil redistribution, while soil redistribution by water only plays a minor role, which is also quite typical for more continental climatic conditions. Modelling a representative segment of a larger landscape instead of focussing on a small (water-) erosion-prone area with steeper slopes indicates that the C sink function is less pronounced at this scale. The study also underlines the importance of addressing the soil-redistribution-induced C fluxes starting with forest-related SOC stocks before conversion to arable land because focussing only on the phase of arable soil use alone overestimates the erosion-induced sink function.

Appendix A

Dupin et al. (2009)Dupin et al. (2009)Kimaro et al. (2005)Kimaro et al. (2005)Su and Zhang (2010)Su and Zhang (2010)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Turkelboom et al. (1999)Van Muysen et al. (2000)Turkelboom et al. (1999)Turkelboom et al. (1999)Wang et al. (2021)Wang et al. (2021)Kimaro et al. (2005)Zhang et al. (2004a)Zhang et al. (2004a)Zhang et al. (2004b)Van Oost et al. (2006)Zhang et al. (2004b)Zhang et al. (2004b)Zhang et al. (2009)Zhang et al. (2009)Quine et al. (1999c)Van Oost et al. (2006)Quine et al. (1999c)Van Oost et al. (2006)Nyssen et al. (2000)Nyssen et al. (2000)Nyssen et al. (2000)Nyssen et al. (2000)Dercon et al. (2007)Dercon et al. (2007)Dercon et al. (2007)Dercon et al. (2007)Quine et al. (1999b)Quine et al. (1999b)Quine et al. (1999b)Quine et al. (1999b)Barneveld et al. (2009)Barneveld et al. (2009)Barneveld et al. (2009)Barneveld et al. (2009)Barneveld et al. (2009)Barneveld et al. (2009)Rymshaw et al. (1997)Van Oost et al. (2006)Thapa et al. (1999a)Van Oost et al. (2006)Thapa et al. (1999b)Van Oost et al. (2006)Thapa et al. (1999b)Van Oost et al. (2006)Quine et al. (1999b)Quine et al. (1999b)Govers et al. (1994)Govers et al. (1994)Lobb et al. (1999)Van Muysen et al. (2000)Marques da Silva and Alexandre (2004)Van Oost et al. (2006)Marques da Silva and Alexandre (2004)Van Oost et al. (2006)Poesen et al. (1997)Van Oost et al. (2006)Quine et al. (1999a)Quine et al. (1999a)Tiessen et al. (2007)Tiessen et al. (2007)Van Muysen and Govers (2002)Van Muysen and Govers (2002)Van Muysen et al. (2000)Van Muysen et al. (2000)Van Muysen et al. (2000)Van Muysen et al. (2000)Mech and Free (1942)Van Oost et al. (2006)Mech and Free (1942)Van Oost et al. (2006)Marques da Silva and Alexandre (2004)Marques da Silva and Alexandre (2004)Mech and Free (1942)Van Oost et al. (2006)Tiessen et al. (2010)Tiessen et al. (2010)Tiessen et al. (2010)Tiessen et al. (2010)Tiessen et al. (2010)Tiessen et al. (2010)De Alba (2001)Van Oost et al. (2006)De Alba (2001)Van Oost et al. (2006)Gerontidis et al. (2001)Van Oost et al. (2006)Gerontidis et al. (2001)Van Oost et al. (2006)Gerontidis et al. (2001)Van Oost et al. (2006)Gerontidis et al. (2001)Van Oost et al. (2006)Gerontidis et al. (2001)Van Oost et al. (2006)Govers et al. (1994)Govers et al. (1994)Heckrath et al. (2006)Van Oost et al. (2006)Heckrath et al. (2006)Van Oost et al. (2006)Heckrath et al. (2006)Van Oost et al. (2006)Heckrath et al. (2006)Van Oost et al. (2006)Heckrath et al. (2006)Van Oost et al. (2006)Heckrath et al. (2006)Van Oost et al. (2006)Heckrath et al. (2006)Van Oost et al. (2006)Kosmas et al. (2001)Van Oost et al. (2006)Kosmas et al. (2001)Van Oost et al. (2006)Lindstrom et al. (1992)Van Oost et al. (2006)Lindstrom et al. (1992)Van Oost et al. (2006)Lobb et al. (1995)Van Oost et al. (2006)Lobb et al. (1999)Van Muysen et al. (2000)Marques da Silva and Alexandre (2004)Marques da Silva and Alexandre (2004)Montgomery et al. (1999)Montgomery et al. (1999)Quine and Zhang (2004)Van Oost et al. (2006)Quine and Zhang (2004)Quine and Zhang (2004)Quine et al. (2003)Van Oost et al. (2006)Revel and Guiresse (1995)Van Oost et al. (2006)Tiessen et al. (2007)Tiessen et al. (2007)Tsara et al. (2001)Tsara et al. (2001)Van Muysen et al. (1999)Van Muysen et al. (1999)Van Muysen and Govers (2002)Van Muysen and Govers (2002)Van Muysen and Govers (2002)Van Muysen and Govers (2002)Mech and Free (1942)Van Oost et al. (2006)Van Oost et al. (2006)Govers et al. (1994)Van Muysen et al. (2006)Govers et al. (1994)Govers et al. (1994)Govers et al. (1994)Govers et al. (1994)Govers et al. (1996)Govers et al. (1996)Govers et al. (1996)Govers et al. (1996)Quine et al. (1994)Quine et al. (1994)Van Muysen et al. (2006)Van Muysen et al. (2006)Wilken et al. (2020)Wilken et al. (2020)Van Muysen et al. (2006)Van Muysen et al. (2006)Van Oost et al. (2000)Van Oost et al. (2000)Lobb et al. (1999)Van Muysen et al. (2000)

Table A1 Literature review as basis for the temporal variation of the tillage transport coefficient (k_til) values shown in Fig. 2a, sorted by period. “Original reference” is the publication where the experiment is described, which does not necessarily include the calculated k_til value. “Reference k_til value” is the publication that refers to the original reference but includes the k_til value. Abbreviations: CP refers to chisel plough, MP refers to mouldboard plough, TE refers to tillage implements, D refers to downslope, UD refers to alternating up- and downslope, C refers to contour-parallel, EXP refers to experiment, and SIM refers to simulation.