Understanding the role of water and tillage erosion from ^239+240Pu tracer measurements using inverse modelling

Soil redistribution on arable land is a major threat for a sustainable use of soil resources. The majority of soil redistribution studies focus on water erosion, while wind and tillage erosion also induce pronounced redistribution of soil materials. Tillage erosion especially is understudied, as it does not lead to visible off-site damages. The analysis of on-site/in-field soil redistribution is mostly based on tracer studies, where radionuclide tracers (e.g. ¹³⁷Cs, ^239+240Pu) from nuclear weapon tests are commonly used to derive the erosion history over the past 50–60 years. Tracer studies allow us to determine soil redistribution patterns but integrate all types of soil redistribution processes and hence do not allow us to unravel the contribution of individual erosion processes. The aim of this study is to understand the contribution of water and tillage erosion leading to soil patterns found in a small hummocky ground moraine kettle hole catchment under intensive agricultural use. Therefore, ^239+240Pu-derived soil redistribution patterns were analysed using an inverse modelling approach accounting for water and tillage erosion processes. The results of this analysis clearly point out that tillage erosion is the dominant process of soil redistribution in the study catchment, which also affects the hydrological and sedimentological connectivity between arable land and the kettle hole. A topographic change up to 17 cm (53 yr)⁻¹ in the eroded parts of the catchment is not able to explain the current soil profile truncation that exceeds the ^239+240Pu-derived topographic change substantially. Hence, tillage erosion already started before the onset of intense mechanisation since the 1960s. In general, the study stresses the urgent need to consider tillage erosion as a major soil degradation process that can be the dominant soil redistribution process in sloped arable landscapes.

1 Introduction

Soil erosion is a major threat to the supply of soil-related ecosystem services (Montanarella et al., 2016). Over the past decades, the off-site effects associated with water erosion, like nutrient inputs from arable lands into inland waters (Pimentel and Burgess, 2013) or siltation of reservoirs (Krasa et al., 2019), have been in scientific and political focus. Within the European Union, the focus on off-site erosion effects is partly caused by the definition of the goals of the EU Water Framework Directive (EU 2000/60/ES) that focuses mainly on water bodies and floodplains but not on a fully integrated catchment management that would call for complex shared responsibilities between different administrative units. Thereby, other soil erosion drivers like tillage and wind are somewhat out of the scope of most studies. Tillage erosion is a mostly ignored soil erosion process (Fiener et al., 2018) that, however, substantially contributes to on-site effects on soil properties and hence agricultural productivity (Winnige, 2004; Nie et al., 2019). Van Oost et al. (2006) pointed out that tillage erosion rates are globally at least of the same order of magnitude as water erosion rates. In the Parisian Basin of France, Chartin et al. (2013) demonstrated the dominating role of tillage erosion. Particularly areas of a hummocky topography with short summit–footslope distances, such as young morainic areas, can be subject to pronounced in-field soil degradation patterns caused by tillage erosion (Winnige, 2004; Deumlich et al., 2017). Young morainic areas that are under intense arable cultivation and associated tillage erosion are widespread in northern Europe, Canada, the northern USA, Russia and eastern Argentina.

Measuring or monitoring water and tillage erosion is challenging as both processes are interlinked and are strongly controlled by topography (Van Oost et al., 2005b, 2006). The quantification of water erosion requires a sufficiently long monitoring time (typically decades) to cover a statistically representative variation of rainfall events occurring in different land cover conditions (Fiener et al., 2019). Therefore, thousands of plot experiments, driven by either natural or artificial rainfall simulations, were carried out in different environments and different land cover conditions (Cerdan et al., 2010; Auerswald et al., 2014). Furthermore, a large number of small catchment studies were performed to quantify both erosion and depositional processes (for an overview, see Fiener et al., 2019). However, soil erosion monitoring is mostly based on sediment delivery monitoring, which cannot address catchment-internal redistribution. In contrast, tillage erosion can only be measured based on the movement of previously applied tracers or by morphological change monitoring (for an overview of tillage erosion measuring techniques, see Fiener et al., 2018). However, these previously applied tracers and change monitoring methods cannot provide a reconstruction of soil redistribution of the past. Natural or anthropogenic tracers in soils can be used to understand soil redistribution (Fiener et al., 2018). Especially anthropogenic radionuclides (e.g. ¹³⁷Cs, ^239+240Pu, ²¹⁰Pb, ⁷Be) have demonstrated their ability to determine changes in topography (Mabit et al., 2014; Alewell et al., 2017; Deumlich et al., 2017). The force of atmospheric nuclear weapon tests transported radioisotopes outside the troposphere, where circulation led to a (regionally) homogeneous spatial distribution and subsequent fallout on soils via precipitation (Meusburger et al., 2016). The main period of atmospheric nuclear weapon tests was from 1953 to 1964 (Schimmack et al., 2001), while the Test Ban Treaty caused a rapid decrease in atmospheric bomb tests in 1963–1964 (Wallbrink and Murray, 1993; Evrard et al., 2020). This rapid decrease led to a distinct peak in the activity of radioisotopes in soils, which enables the use of radioisotopes as redistribution tracers in soils (Alewell et al., 2017). The radioisotope ¹³⁷Cs has been used as a soil redistribution tracer in a large number of studies (e.g. Porto and Walling, 2012; Chartin et al., 2013; Evrard et al., 2020) and has become a widely used method in soil erosion science. However, the Chernobyl disaster in 1986 supplied additional radioactive fallout to soils across large areas of Europe (Evangeliou et al., 2016). For some years after the Chernobyl disaster, an unmixing of the Chernobyl fallout from the original 1960s bomb peak signal was possible by the use of the ¹³⁴Cs∕¹³⁷Cs ratio (Lust and Realo, 2012). However, due to the short ¹³⁴Cs half-life of 2 years (Schimmack et al., 2001), this method cannot be applied anymore. Hence, the use of ¹³⁷Cs as a soil redistribution tracer in Europe is associated with uncertainties and requires special attention concerning a potential Chernobyl contamination (Evangeliou et al., 2016). Furthermore, due to the ¹³⁷Cs half-life of about 30 years, decay has already led to a pronounced reduction (73 % in 2020) of activity until today (Alewell et al., 2017). Over the past decade, ^239+240Pu has been discussed and tested as an alternative radioisotopic tracer for soil erosion studies. Decay is not an issue as the half-life of ²³⁹Pu and ²⁴⁰Pu is long (²³⁹Pu = 24 000 years; ²⁴⁰Pu = 6563 years) and the ^239+240Pu contamination by the Chernobyl accident was spatially very limited (<100 km; Kashparov et al., 2004; Matsunaga and Nagao, 2009) and can be determined by the ²³⁹Pu∕²⁴⁰Pu ratio (Alewell et al., 2014, 2017).

Radionuclide tracers integrate soil erosion processes over time (e.g. since the bomb peak of the 1960s in the case of ¹³⁷Cs and ^239+240Pu), which somewhat averages out the large temporal variability of water (episodic nature) and tillage (mechanisation) erosion. However, the use of radioisotope tracers integrates all types of soil redistribution processes and does not provide information on the relative contribution of the driving processes at play (e.g. water, tillage, wind). Unravelling the respective contributions of these different processes requires the use of an inverse modelling approach carrying out model runs with different parameterisations to alter the contribution and mechanisms of different soil redistribution drivers. There are only a few models that take both water and tillage erosion processes into account. Physically oriented models like MCST-C (Wilken et al., 2017b) and LandSoil (Ciampalini et al., 2012) simulate individual erosion events and are developed to enhance process understanding, while conceptual USLE-based models (WaTEM/SEDEM: Van Oost et al., 2000; Van Rompaey et al., 2001) aim at a robust prediction of long-term soil erosion rates (Alewell et al., 2019). For an inverse modelling approach to unravel tillage and water erosion based on a radionuclide tracer, it is necessary to use a parsimonious approach with a limited parameter space covered by available input data over 5 to 6 decades, which suggests the use of conceptual models.

In this study, we will determine the soil redistribution patterns in a small 4.2 ha catchment based on high-resolution ^239+240Pu measurements and analyse the contribution of water and tillage erosion processes based on an inverse modelling approach using a combined water and tillage erosion model. The general objective is to unravel the importance of water and tillage erosion driving the current variability of soil properties in an intensively used arable landscape of north-eastern Germany.

Figure 1 (a) Topography and nested soil sampling scheme in the young morainic study area in north-eastern Germany. P1 to P10 indicate the locations for high-resolution depth profile sampling (see Fig. 3). (b) Idealised catena for soil landscapes in hummocky ground moraines; 1: Calcaric Regosol (extremely eroded), 2: Nudiargic Luvisol (strongly eroded), 3: Calcic Luvisol (non-eroded), 4: Gleyic-Colluvic Regosol (colluvial).

2 Methods

2.1 Study area

The study area (53^∘21^′2^′′ N, 13^∘39^′5^′′ E) is situated in the hummocky ground morainic landscape of the Weichselian glacial belt (“young morainic area”) of north-eastern Germany (Fig. 1). Characteristics of these landscapes include widespread closed depressions, so-called kettle holes, which result from a delayed melting of dead ice blocks. They are nowadays filled with mineral soil, (degraded) peat or water. The study area is part of a kettle hole catchment (4.2 ha; Fig. 1) showing a high morphological variability covering convex hilltops, steep slopes and flat areas. The recent crop rotation is rape (Brassica napus L.)–winter wheat (Triticum aestivum L.)–winter barley (Hordeum vulgare L.)–winter barley, cultivated without cover crops, which is a typical conventional crop rotation that is adapted for the highly fertile soils of the Uckermark region. The mean arable land of a farm in the region is 352 ha, which is much larger compared to the mean of the State of Brandenburg (250 ha) and Germany (60 ha; Troegel and Schulz, 2018). The associated large field sizes are explained by land consolidation programmes implemented in the 1960s during the socialist period. However, also before, agriculture in the region was already characterised by large-scale farming structures and corresponding high agricultural mechanisation. In 1939, large farms that manage more than 100 ha of arable land cultivated 7 % of the total arable area of Germany. In contrast, within the present-day federal states of Mecklenburg–Pomerania, Brandenburg (study area location) and Saxony-Anhalt, large farms cultivated 30 % of corresponding arable land (Wolz, 2013). The catchment is part of a single large field (54 ha), which is a size that can be frequently found in the region. The soils are developed from glacial till and vary with respect to their location in the landscape. Convex hilltops and steep slopes are dominated by extremely eroded A–C profiles (Calcaric Regosols, soil classification according to IUSS, 2015), while Luvisols showing different degrees of erosion that are typically situated at the up and mid slopes, the footslopes and depressions are dominated by Gleyic-Colluvic Regosols (Fig. 1b; Sommer et al., 2008; Gerke et al., 2010). These soils regularly reveal fossil surface horizons below 1 m depth (fAh, fH). The closed kettle hole depression itself is built up by degraded Histosols and covered by a thin colluvial layer of mineral soil (40±8 cm mean, 25 % and 75 % quantile, n=20; based on soil auger prospection). The study area has a continental climate (Köppen: Dfb) with low annual precipitation (500 mm; Fiener et al., 2018) and high temperature amplitude (July = 18 ^∘C, January = 0 ^∘C; mean 8.9 ^∘C WMO-CLINO 1981–2010 for the meteorological stations Gruenow and Angermuende). Between 7 and 11 erosive rainfall events take place each year (Deumlich, 1999). In the region, maximum intensities up to 162 mm h⁻¹ (10 min interval) were recorded during an extreme event in 2016 (Wilken et al., 2018).

2.2 Soil sampling design and preparation

The soil sampling design was organised according to a regular 20 m × 20 m grid with at least one sampling point of the transect line exceeding the spatial extent of the catchment under study (Fig. 1) to avoid boundary effects. To assess small-scale spatial variability (for distances of sampling points between 5 and 20 m), a nested sampling approach was applied (Fig. 1; Hengl and MacMillan, 2019). Therefore, five densified sub-grids, located at different topographical locations (hilltop, ridge shoulder, moderately steep mid slope, steep mid slope, footslope/valley; Fig. 1), were selected. In total, 209 locations were included in the sampling design.

The sampling points were located using a differential GPS (AgGPS™ 132; Trimble Navigation Ltd, Sunnyvale CA, USA) applying the SAPOS (LGB) correction signal. With respect to drill penetrability, soil sampling was carried out under moderately wet soil moisture conditions (ranging from 18 % to 32 % and 16 % to 33 % for topsoil and subsoil, respectively) in December 2015. Closed soil cores, using a steel drill containing a plastic liner (4.6 cm inner diameter), were driven by a percussion corer (Cobra TTe; Atlas Copco Power Techniques GmbH, Stockholm, Sweden) into the ground down to a depth of 50 cm at 209 sampling points.

To minimise physical and biogeochemical disturbance, the soil cores were stored in a freezer until sample preparation. After complete thawing, the soil cores were separated into topsoil (Ap horizon) and subsoil. The separation was done by visual interpretation of soil horizon characteristics (colour, structural and density differences), showing a variation in topsoil thickness between 16 and 30 cm (mean of 23.5 cm) depending on the topographic position. Aliquot (n=3) samples for each topsoil and subsoil location were taken for gravimetric water content (weighted before and after drying at 105 ^∘C) and dry bulk density (known sample volume) measurements. The soil samples were air dried and subsequently sieved with a 2 mm mesh to separate stones from the fine soil.

2.3 ^239+240Pu measurements

All topsoil samples (n=209) were analysed, while 145 (∼ 70 %) of the subsoil samples of highest topsoil activities were measured on their ^239+240Pu activity. This was done to reduce the number of samples with Pu activities below the detection limit. Furthermore, at 10 locations higher-resolution (5 cm) depth increments were measured to assess the depth distribution of ^239+240Pu activity below the mixed plough layer.

The fallout radionuclides ^239+240Pu were used to estimate effective soil redistribution since the 1960s. Plutonium isotope measurements were conducted following Calitri et al. (2019) based on the procedure of Ketterer et al. (2004). Before the mass spectrometry analysis, 10 g of milled fine earth was dry-ashed for at least 8 h at 600 ^∘C to remove organic matter. Subsequently, the samples were spiked using 30 pg (ca. 0.0044 Bq) of a ²⁴²Pu tracer solution (NIST 4334). The sample leaching was applied using 16 M nitric acid (HNO₃) overnight at 80 ^∘C and subsequently filtered and adjusted to a concentration of 8 M HNO₃. Plutonium species were adjusted to the Pu (IV) oxidation state using first an acidified FeSO₄•7H₂O solution (2 mg mL⁻¹ of leached solution) and subsequently a sodium nitrite (NaNO₂) solution (20 mg mL⁻¹ of leached solution). The samples were heated at 75 ^∘C for 2 h. Tetravalent Pu was separated from the leached solution using a Pu-selective TEVA resin (2 mg of TEVA per millilitre of leached solution). Following occasional agitation for 2 h, the resin was collected in a pipette tip equipped with a glass wool plug. This disposable column was rinsed with 2 M aqueous HNO₃ to remove unretained matrix elements (e.g. uranium (U)), then rinsed with 8 M HCl to elute thorium (Th) and finally rinsed again with 2 M aqueous HNO₃ (rinse volume = 1 mL per 30 mg of TEVA). Plutonium was eluted using 0.05 M aqueous ammonium oxalate. Data quality was evaluated through the analysis of blanks (soils or rocks devoid of Pu), duplicates and control samples of known ^239+240Pu activities (Standard Reference material 4350b – River sediment for radioactivity measurements from NIST). Activities of ^239+240Pu were measured using a Thermo X Series II quadrupole ICPMS, located at Northern Arizona University. The ICPMS instrument is equipped with an APEX HF high-efficiency sample introduction system. The masses of ²³⁹Pu and ²⁴⁰Pu in the samples were converted into the summed activity ^239+240Pu.

2.4 ^239+240Pu-based soil erosion assessment

We applied the proportional conversion approach of Walling et al. (2011). Erosion is calculated following Eqs. (1) and (2):

$$\begin{array}{}\text{(1)}& {\mathrm{SL}}_i=\mathrm{10}\cdot {\mathrm{BD}}_i\cdot {\mathrm{TD}}_i\cdot {\mathrm{RR}}_i\cdot T^{-\mathrm{1}},\end{array}$$

where SL_i is the mean annual soil loss in Mg ha⁻¹ yr⁻¹, BD_i is the soil bulk density in kg m⁻³, TD_i is the vertical depth of the Ap horizon in metres (tillage depth), RR_i is the relative reduction of the reference inventory of the ^239+240Pu inventory and T is the years that have elapsed since the end of atmospheric nuclear weapon tests (mainly since 1964).

$$\begin{array}{}\text{(2)}& {\mathrm{RR}}_i={\displaystyle \frac{({\mathrm{Pu}}_{\mathrm{ref}}-{\mathrm{Pu}}_i)}{{\mathrm{Pu}}_{\mathrm{ref}}}},\end{array}$$

where Pu_i is the inventory at sampling point i and Pu_ref is the reference inventory of undisturbed sites in Bq m⁻² (see implementation Sect. 2.5).

For a three-dimensional representation of the ^239+240Pu redistribution by water and tillage erosion, SPEROS-Pu was developed that is based on a modified version of the SPEROS-C model (Van Oost et al., 2005a; Fiener et al., 2015; Nadeu et al., 2015). SPEROS-Pu is a spatially explicit water and tillage erosion model. Water erosion is simulated based on a gridded application of the Revised Universal Soil Loss Equation (RUSLE: Renard et al., 1996) coupled with a sediment transport and deposition approach. Erosion is calculated according to a slightly modified RUSLE approach and transport and deposition are based on the grid-cell-specific local transport capacity TC (kg m⁻¹ yr⁻¹), which multiplies RUSLE factors by a transport capacity coefficient (k_tc; in metres):

$$\begin{array}{}\text{(3)}& \mathrm{TC}=k_{\mathrm{tc}}\cdot R\cdot C\cdot K\cdot L\cdot S\cdot P,\end{array}$$

where R, C, K, L, S and P are the RUSLE factors (see Renard et al., 1996).

The tillage erosion module of SPEROS-Pu follows a diffusion-type equation adopted from Govers et al. (1994) that derives tillage erosion based on change in topography and management-specific coefficients:

$$\begin{array}{}\text{(4)}& Q_{\mathrm{til}}=-k_{\mathrm{til}}{\displaystyle \frac{\mathrm{\Delta }h}{\mathrm{\Delta }x}},\end{array}$$

where Q_til is the soil flux in kg m⁻² yr⁻¹, Δh is the elevation difference in metres, Δx is the horizontal distance in metres, and k_til is the tillage transport coefficient in kg m⁻¹ yr⁻¹:

$$\begin{array}{}\text{(5)}& k_{\mathrm{til}}={\mathrm{BD}}_i\cdot {\mathrm{TD}}_i\cdot x_{\mathrm{til}},\end{array}$$

where x_til is the tillage translocation distance in metres (for BD_i and TD_i, see the description of Eq. 1).

The representation of the ^239+240Pu redistribution in the SPEROS-Pu model is three dimensional and accounts for ^239+240Pu source area depletion and corresponding redistribution of depleted sediments. The horizontal distribution of ^239+240Pu is grid based, while the vertical distribution is represented by ten 10 cm depth layers. The two uppermost layers are assumed to be homogeneously mixed by tillage operations and have the average ^239+240Pu activity of the upper two soil layers. At the beginning of the simulation, the ^239+240Pu reference activity is homogeneously distributed within the mixed plough layer and over the entire catchment. Subsequently, the local ^239+240Pu inventory is altered by soil redistribution processes. Soil erosion processes lead to a reduction of the ^239+240Pu inventory per m² due to soil and associated ^239+240Pu loss, which causes mixing of non-contaminated subsoil. Deposition adds contaminated material on top and increases the ^239+240Pu inventory.

2.5 Implementation and inverse modelling

2.5.1 Soil redistribution based on ^239+240Pu measurements

Spatially distributed topsoil and subsoil bulk density and tillage depth information for each individual sampling location was applied to the proportional conversion approach. The reference inventory of undisturbed sites follows the value determined by Calitri et al. (2019), who found a ^239+240Pu inventory of 43±3 Bq m² based on four sites that did not show a profile morphological or geochemical indication of soil redistribution at a location 8.5 km from the study area. To address the uncertainty inherent to the reference measurements, a reference range from 40 to 46 Bq m² was accounted for in the simulations.

The point data of the ^239+240Pu inventory for the depth 0–50 cm were geostatistically (block kriging) interpolated for a gridded spatial representation that matches the spatial resolution of the soil redistribution model (5 m × 5 m). Block kriging was used to reduce small-scale scattering that is naturally inherent to soil cores of 4.6 cm diameter that are supposed to represent a decametre scale. Different block sizes were tested for the kriging approach. A block size of 20 m was selected that matches the sampling resolution and did not to cause over-smoothening of the interpolation result. The theoretical semivariogram model (Fig. 2) was fitted using all 209 sampling points, including the nested samples to account for variations over short distances. However, the input data of the interpolation itself solely use the regular 20 m × 20 m grid points. The interpolation and geostatistical analysis were carried out using the statistical software GNU R (R Core Team, 2019), version 3.5.3, and the add-on package gstat (Pebesma, 2004).

Figure 2 Semivariogram and semivariogram model fit of the ^239+240Pu block kriging interpolation.

2.5.2 Inverse modelling of water and tillage erosion

An inverse modelling approach was used to understand the proportion of water and tillage erosion that is inherent to the ^239+240Pu-based soil erosion map. The inverse modelling iterates three parameter sets from low to extreme settings over the 53-year modelling period (1964–2016): (i) the k_til tillage translocation coefficient (as given for illustration in Eq. 5), (ii) the product of all RUSLE factors (as given in Eq. 3, in the following referred to as water erosion strength) and (iii) the water erosion transport capacity coefficient k_tc (see Table 1) that controls the transport distance and is the standard calibration parameter of the model. While changes in the tillage translocation coefficient and water erosion strength only alter the quantity of soil redistribution, the k_tc has a pronounced impact on spatial patterns of modelled soil redistribution. The parameter range covers very low to extreme soil redistribution rates ($k_{\mathrm{til}}max.=\mathrm{1000}$ kg m⁻¹; $k_{\mathrm{tc}}max.=\mathrm{500}$ m; RUSLE factors product deviation = 100 %; see Table 1). The water erosion reference parameterisation for k_tc is 150 m (Van Oost et al., 2003). The interplay of parameter combinations was assessed in 35 722 different model runs (Table 1).

Table 1 Parametrisation of the inverse modelling approach.

To determine the model match, different goodness-of-fit parameters were calculated that compare the interpolated ^239+240Pu raster map against the results calculated by the inverse modelling approach. To address the high spatial variability of the 5 m × 5 m raster-by-raster comparison, a classification of results was carried out. Therefore, mean values were calculated based on the cells that fall into a specific class of the interpolated map of ^239+240Pu inventories. The classification covers 20 classes with 2.5 Bq m⁻² steps from 17.5 to 65 (and a class >65) Bq m⁻² of ^239+240Pu. First, the spatial correlation was calculated for both the raster-by-raster comparison and the classified results. Second, the classified results of the inverse modelling were tested using goodness-of-fit parameters that take absolute differences between observed and predicted values into account (RMSE, model efficiency coefficient: MEF according to Nash and Sutcliffe, 1970). As a last step, the results of the measured and modelled ^239+240Pu inventories were transferred into tillage and water erosion maps (given in Mg ha⁻¹ and topographical change in centimetres) applying the proportional conversion approach (see Sect. 2.4).

3 Results

3.1 ^239+240Pu activities and inventories

The topsoil and subsoil ^239+240Pu activities differ substantially from each other and show distinct spatial patterns according to their topographic position. All topsoil samples (n=209) showed a ^239+240Pu activity above the detection limit, while 7 of the 145 subsoil samples fall below the detection limit (<0.002 Bq kg⁻¹). Those seven samples are all located at positions with ^239+240Pu inventories below the lower reference boundary (40 Bq m⁻²). The average ^239+240Pu activity is 0.078±0.016 Bq kg⁻¹ and 0.035±0.038 Bq kg⁻¹ for topsoil and subsoil, respectively. All high-resolution depth profiles (5 cm increments) at erosional sites show a sharp reduction of the ^239+240Pu activity below the plough layer (Fig. 3), while depositional sites show more complex depth distributions. Location P4 (Fig. 3) does not show a drop in activity until a depth of 0.5 m, while P1 (and partly P3) shows an increase in the ^239+240Pu activity with depth that is potentially caused by Pu enrichment processes during lateral transport or the deposition of already Pu-depleted source material into the topsoil. The topsoil ^239+240Pu/reference (43 Bq m⁻²) ratio indicates soil erosion-related ^239+240Pu depletion or enrichment according to a ratio lower and higher than one, respectively. The highest depletion (min. = 0.28; 5 % quantile = 0.37) can be found at the hilltops that are most affected by tillage, while enrichment (max. = 1.18; 95 % quantile = 0.92) is spatially limited to the flat surroundings of the kettle hole (Fig. 1), where topsoil material was potentially deposited by both water and tillage soil redistribution processes. At 14 (of 209) sampling locations, a higher subsoil than topsoil ^239+240Pu activity was found, which points to deposition of ^239+240Pu-depleted sediments. The majority of these locations show enriched ^239+240Pu activities in the subsoil (11 of 14). These locations are all, except for one, located at the kettle hole surrounding plateau where both water and tillage erosion cause deposition. The ^239+240Pu depletion of sampling points (including locations outside the study catchment) goes down to 12 Bq m⁻², while the highest five locations exceed 86 Bq m⁻², which means that the ^239+240Pu inventory has more than doubled compared to the reference inventory of 43 Bq m⁻², which has to be attributed to enrichment processes. As the enrichment processes inherent to these five locations cannot be corrected, the locations were excluded from the analysis.

Figure 3 Depth distribution of ^239+240Pu at different geomorphological positions. Locations P1 to P10 are given in Fig. 1.

The distribution of the interpolated maps covering the study catchment shows that substantially more locations fall below the reference inventory than exceed it (Fig. 4). This indicates that a larger area is subject to erosion processes compared to depositional processes.

Figure 4 Distribution histogram of 5 m × 5 m interpolated ^239+240Pu measurements in 20 classes with descriptive statistics.

3.2 ^239+240Pu measurements vs. inverse water and tillage erosion modelling

To understand the drivers of current ^239+240Pu and associated soil degradation patterns, an inverse modelling was carried out that was quantitatively analysed by goodness-of-fit parameters. The spatial correlation between the ^239+240Pu-derived patterns and the modelled best knowledge soil redistribution, including both water and tillage erosion, is only moderate (R²=0.45, ρ = 0.73) on a raster-by-raster comparison (n=1699, 5 m × 5 m grid points; Fig. 5a). To reduce small-scale variability and understand the average goodness-of-fit, the inverse modelling results were classified according to the measured ^239+240Pu activity. The classified results average out the spatio-temporal dynamics and show a very high correlation (R²=0.95, Rho = 0.99; Fig. 5b), which illustrates the great agreement of the spatial soil redistribution patterns between the ^239+240Pu measurements and the model results.

Figure 5 Linear correlation between measured and modelled ^239+240Pu inventories redistributed by water (k_tc: 150, P factor: 1) and tillage erosion (k_til: 350 kg m⁻¹; ^* = p value <0.001). (a) Point-by-point correlation at 5 m × 5 m resolution (n: 1699); (b) class aggregation according to ^239+240Pu-derived soil redistribution. Minimum and maximum class n is 27 and 184, respectively. While the points and classes are calculated for a reference of 43 Bq m⁻², the trend lines display the offset sensitivity of different reference ^239+240Pu activities.

While the analysis of the spatial correlation is a relative comparison, the absolute deviation is considered according to the MEF (Nash and Sutcliffe, 1970; 1 = perfect prediction, 0 = as good as mean of all measurements, <0 = worse than mean). The quality of model predictions shows hardly any sensitivity to water erosion-related parameterisations (k_tc and erosion strength; Fig. 6c). In contrast, the tillage erosion strength, represented by k_til parameter iterations, showed a substantial impact on the MEF (Fig. 6a and b). A MEF better than 0.8 and RMSE below 6.5 Bq m⁻² were found for a k_til range from 225 to 475 kg m⁻¹, while the best model fit was found for a k_til of 350 kg m⁻¹, achieving a MEF of 0.87 and a corresponding RMSE of 5.2 Bq m⁻². The best model fit was found without the contribution of water erosion. The highest impact on the best-fit model run was a 0.31 MEF reduction by an extreme water erosion parameterisation (k_tc=500, water erosion strength = 200 %).

Figure 6 Inverse modelling of tillage and water erosion compared to ^239+240Pu-derived soil redistribution. Three parameter combinations (tillage transport coefficient, k_til; water transport capacity coefficient, k_tc; deviation in water erosion strength compared to reference parameterisation) are tested on their effect on the goodness-of-fit, represented by the MEF (model efficiency coefficient: perfect model fit = 1; model prediction as good as the mean = 0; model prediction worse than mean = <0).

Soil redistribution determined by the proportional conversion approach using ^239+240Pu measurements indicates substantial geomorphological dynamics in the study catchment over the past decades. Soil erosion at hilltop locations is shown to reach up to 14.9 cm (43 Bq m⁻² reference; 40 Bq m⁻² reference = 14.1 cm; 46 Bq m⁻² reference = 15.6 cm), while deposition can build a colluvium layer with a maximum thickness of 21.5 cm (43 Bq m⁻² reference; 40 Bq m⁻² reference = 24.9 cm; 46 Bq m⁻² reference = 18.6 cm) over the past 53 years. The inverse modelling stresses that substantial soil erosion, which takes place over large areas, is almost exclusively attributed to tillage translocation (modelled max. water erosion = 3.8 cm (53 yr)⁻¹ vs. max. tillage erosion = 13.5 cm (53 yr)⁻¹; Fig. 7c, d). In turn, both processes contribute to deposition in the kettle hole surrounding flats (max. water deposition = 27.1 cm (53 yr)⁻¹ vs. max. tillage deposition = 22.4 cm (53 yr)⁻¹; Fig. 7c, d).

Figure 7 The figure consists of four parts. (a) Soil redistribution derived from ^239+240Pu top and subsoil measurements using 43 Bq m⁻² as the reference inventory; (b) geostatistically interpolated soil redistribution based on ^239+240Pu point measurements; (c) modelled tillage erosion with a tillage transport coefficient (k_til) of 350 kg m⁻¹; (d) modelled water erosion according to reference parameterisation (k_tc: 150; also see Table 1).

4 Discussions

4.1 ^239+240Pu methodological benefits and limitations

The use of fallout radionuclides to determine soil redistribution patterns and rates over the past decades has been used in many studies in various study areas around the world (see reviews: Mabit et al., 2014; Alewell et al., 2017; Evrard et al., 2020) and contributed substantially to understanding soil degradation processes. However, the most frequently used fallout radionuclide, ¹³⁷Cs, faces upcoming limitations (Chernobyl fallout that adds to the global fallout over large parts of Europe and ongoing decay below the detection limit of standard measuring devices; also see Sect. 1) in the use as a soil redistribution tracer (Evrard et al., 2020). The fallout radionuclide ^239+240Pu has demonstrated its suitability for determining the recent soil redistribution history (since the 1960s; see review in Alewell et al., 2017) and is a potential alternative for ¹³⁷Cs as a soil redistribution tracer (Mabit et al., 2013; Alewell et al., 2017). In Europe, where large parts were re-contaminated by ¹³⁷Cs fallout of the Chernobyl accident (Evangeliou et al., 2016), additional information on the spatial change in the inventory is needed to derive accurate soil redistribution rates. Particularly in the area of the former GDR, almost no information that can be used for a correction on the ¹³⁷Cs Chernobyl re-contamination is available (Evangeliou et al., 2016). The ^239+240Pu fallout caused by the Chernobyl disaster was very local (radius <100 km) and has a distinct fingerprint based on the ²³⁹Pu∕²⁴⁰Pu ratio. While the ²⁴⁰Pu∕²³⁹Pu ratio of global fallout in the Northern Hemisphere is 0.180±0.014 (Kelley et al., 1999), the ²⁴⁰Pu∕²³⁹Pu ratio in soils that received high Chernobyl fallout is about twice as high (0.408±0.003, determined for soils within the 30 km exclusion zone of the Chernobyl reactor; Muramatsu et al., 2000; Boulyga and Becker, 2002). The 95 % interval of confidence and average of the ²⁴⁰Pu∕²³⁹Pu ratio found in the soil samples of this study were 0.281 and 0.199, respectively. Hence, a relevant ^239+240Pu re-contamination by Chernobyl fallout can be ruled out for the study area. Another limitation for the use of ¹³⁷Cs as a soil redistribution tracer is the ongoing decay due to short half-life times that has already caused a substantial reduction of the inventory. Due to lower activities, measuring devices of much higher complexity are needed in the future (Evrard et al., 2020). Decay is not an issue for ²³⁹Pu and ²⁴⁰Pu as both nuclides have long half-life times that allow for a quasi-unlimited use; however, it needs to be mentioned that sample preparation for ^239+240Pu ICP-MS measurements is much more laborious compared to the standard procedure of physical ¹³⁷Cs measurements.

Enrichment processes, due to selective transport of soil constituents that fallout radionuclides are preferentially associated with, are a critical issue for the use of most (e.g. ^239+240Pu, ¹³⁷Cs, ²¹⁰Pb) radionuclide tracers (Parsons and Foster, 2011; Mabit et al., 2014; Alewell et al., 2017). While ¹³⁷Cs is mainly associated with clay particles, ^239+240Pu binds to soil organic matter and oxides (Alewell et al., 2017) that are less affected by selective water transport and corresponding ^239+240Pu enrichment (Meusburger et al., 2016; Xu et al., 2017). However, it needs to be mentioned that radionuclide-associated particles are typically not transported as primary particles, but in soil aggregate complexes (Hu and Kuhn, 2014; Hu et al., 2016), which has a pronounced effect on enrichment processes (Wilken et al., 2017b). Nevertheless, the ^239+240Pu activity at depositional sites, which are redistributed by water (transport by tillage is typically assumed to be non-grain size selective; Fiener et al., 2018), can be higher in relation to the activity of the source material. A soil profile that shows a distinct indicator of enrichment processes in this study is sampling profile P1 (25–45 cm; Fig. 3) that is situated in the kettle hole surrounding flat. Hence, enrichment in fine particles of relatively high ^239+240Pu activity is to some extent also an issue within this study that causes an overestimation of deposition. A particle size correction factor was not applied as topsoil enrichment (topsoil Bq m⁻² > ref. 43 Bq m⁻²) was exclusively found at very few sampling locations (<6 %) in the kettle hole surrounding flats. Furthermore, the mean topsoil ratio of enriched sediments is moderate (1.2) and supports the general assumption that ^239+240Pu is less affected by selective transport compared to ¹³⁷Cs (Alewell et al., 2017) and that transport by tillage is non-grain size specific. The counteracting process of enrichment is the deposition of ^239+240Pu-depleted sediments that are transported from highly eroded locations. Such highly depleted locations can be found at the hilltops of the study area (Fig. 7b). Hence, the hilltops are the main source of highly depleted sediments that are deposited in kettle hole surrounding flats. However, the minimum horizontal distance from the hilltops to the kettle hole surrounding flat is roughly about 70 m and the approximate tillage translocation distance 0.5 to 1 m per pass (Fiener et al., 2018). Hence, deposition of depleted ^239+240Pu material has to be mainly attributed to surface runoff that can flow across longer transport distances. SPEROS-Pu takes depletion of deposited sediments into account but does not address enrichment processes. Furthermore, a maximum soil sampling depth down to 50 cm was carried out within this study that technically allows us to derive a maximum depositional depth of 25 cm using the proportional conversion approach of Walling et al. (2011), which was exceeded at four sampling locations. Nevertheless, also with a deeper soil sampling, it would be arguable whether these potentially enriched or depleted sampling locations should be excluded from the statistical analysis, as was done for extreme depositional locations (four sampling locations) within this study.

The application of soil tracers like fallout radionuclides is based on the assumption that spatial and vertical patterns are exclusively caused by soil redistribution processes due to water and tillage erosion. Bioturbation or colloid migration can be drivers of vertical displacement, which can cause an incomplete representation of the inventory if the downward migration exceeds the soil sampling depth. Within this study, a rapid reduction of the ^239+240Pu activity was found below the plough and colluvial layer of the high-resolution depth profiles (see Fig. 3). Hence, there is no indication of a pronounced downward migration of ^239+240Pu below the sampling depth.

4.2 Using ^239+240Pu and inverse modelling to understand the recent soil erosion history

Within this study, an inverse modelling approach was carried out to understand the contribution of soil redistribution by water and tillage erosion. The model is subject to specific uncertainties that need to be mentioned. The diffusion-type equation for the tillage erosion simulation (Govers et al., 1994) follows the assumption of a spatially homogenous tillage translocation coefficient for the entire study catchment. However, as the tillage translocation coefficient is a function of tillage speed and depth, tillage translocation distances are likely to show spatial differences that may follow the topography of the study area. Furthermore, the model uses a static digital elevation model that does not account for topographical change during the simulation period, which causes static soil redistribution patterns that ignore feedback processes. This is an issue for the used digital elevation model that represents short-term topographical structures like agricultural tramlines at the acquisition time. To reduce this effect, a low-pass filter was applied on the digital elevation model to somewhat even out short-term structures. Numerical dispersion can affect the model results, which is related to the vertical and spatial representation of the model. The vertical soil profile representation in 10 cm increments is rather low for a fallout radionuclide tracer with a rapid decrease in ^239+240Pu between the plough layer/colluvium and the subsoil layer (see Fig. 3). However, the calculation of soil redistribution is based on relative changes in the ^239+240Pu inventory in relation to reference inventories. Therefore, the sensitivity of the soil redistribution calculations to the ^239+240Pu soil profile development is minor as the ^239+240Pu inventory is independent of the depth distribution. The relatively high spatial resolution of 5 m × 5 m was selected to adequately represent the high landscape variability in the study area and has reliably simulated spatial patterns in previous studies (e.g. Van Oost et al., 2003; Dlugoß et al., 2012). Overall, the model performs well and explains 95 % of the variance and achieves a MEF better than 0.8 in predicting ^239+240Pu patterns by simulating the combined effect of three model parameters (tillage translocation coefficient, water transport capacity coefficient, water erosion strength).

Within the intensively managed study catchment, substantial ^239+240Pu-derived soil redistribution was found with soil loss up to 45 Mg ha⁻¹ yr⁻¹ (ref. 43 Bq m⁻²; ref. 40 and 46 Bq m⁻² = 43 and 47 Mg ha⁻¹ yr⁻¹) and sediment deposition up to 65 Mg ha⁻¹ yr⁻¹ (ref. 43 Bq m⁻²; ref. 40 and 46 Bq m⁻² = 75 and 56 Mg ha⁻¹ yr⁻¹). Very high deposition can only be found in the spatially narrow area of the kettle hole surrounding flat where both water and tillage erosion processes lead to deposition (Fig. 1). The kettle hole surrounding flat is a spatially narrow area but the only zone where water erosion substantially contributes to pronounced geomorphological dynamics (Fig. 7d). As a result of the small spatial extent where this process takes place, the inverse modelling shows hardly any sensitivity to goodness-of-fit changes in reaction to the variation in model parameterisations (Fig. 6c). Nevertheless, sediment deposition and delivery by surface runoff is an important process in the study area. Evidence for runoff-based sediment delivery is a colluvial layer covering the peat in the kettle hole with an average depth of 40 cm. This sediment delivery into the kettle hole cannot be explained by the inverse modelling of water erosion applying a reasonable parameter range. Therefore, we assume the reference parameterisation for the region given by the State of Brandenburg to be the most appropriate (k_tc=150 m and RUSLE parameters according to Table 1). According to the model run using the reference parameterisation for water erosion, a colluvial layer of 1.7 cm (53 yr)⁻¹ would have been developed on top of the peat that has been exported from the arable part of the catchment (see Fig. 1) due to water transport over the past decades. This indicates a long water erosion history before the 1960s. This is not surprising as bare soil conditions and erosive rainfall events have taken place since the onset of arable use approximately 1 kyr before present (Van der Meij et al., 2017; Kappler et al., 2018). In contrast, tillage erosion is typically assumed to be a process that is linked to recent developments of increasing mechanical forces that have been applied to soils over the past century (Sommer et al., 2008; Calitri et al., 2019). Within this study, a maximum topographical change by hilltop erosion up to 17 cm (53 yr)⁻¹ was determined. In a review of tillage erosion by Van Oost and Govers (2006), tillage translocation coefficients of 44 experiments were reported for different tillage practices. This resulted in a mean k_til of 234 kg m⁻¹ (5 % percentile = 30 kg m⁻¹; 95 % percentile = 640 kg m⁻¹) for mouldboard and chisel plough. Within this study, a tillage translocation coefficient of 350 kg m⁻¹ per year was determined. The tillage translocation coefficients, determined by the inverse modelling approach, are rather high compared to other studies considering that fallout radioisotopic tracer approaches cover a phase of high mechanical development from low- to high-power farming machines (Sommer et al., 2008; Keller et al., 2019). Although recent tillage translocation rates are rather high, they cannot explain the soil depth patterns that are visible by augerings. In the study region, it can be observed that tillage erosion mainly affected hilltops. Calcaric glacial till is approaching the surface by soil profile truncation and is partially mixed into the plough layer. Within the study catchment, this is the case for 20 sampling locations (CaCO₃>0.5 %) that are also indicated as the most eroded sites by the ^239+240Pu measurements and the inverse modelling. Non-eroded reference profiles (n=210) in the region show the parent material (calcaric glacial till) is found at 102 cm depth on average (van der Meij et al., 2019). Hence, less than 17 % of soil depth reduction can be attributed to most recent process dynamics. This suggests that traditional hand- or cattle-based tillage systems, which have been used since the beginning of arable agriculture in the region (1 kyr BP; Kappler et al., 2018), must have caused extensive soil redistribution over long periods. This suggests that tillage erosion might be the dominant process even without mechanised soil tillage, which is the common practice in most developing countries that also partly cultivate very steep slopes. Therefore, the general assumption of tillage erosion being only an issue for highly mechanised agricultural systems (Van Oost et al., 2006) might need to be reconsidered across a range of contrasted agricultural environments.

4.3 Interplay of sediment redistribution by water and tillage

The inverse modelling has shown that soil redistribution by water has only a minor impact on erosion processes in the study area. However, sediment deposition by water has a complex interplay with tillage translocation (kettle hole surrounding flat; Fig. 7). Very high deposition by tillage translocation towards the field-kettle hole edge (typically >1 m known from soil augering; Kappler et al., 2018) builds up local hydrological depressions (Figs. 1 and 7). Only infrequent extreme events exceed the critical runoff quantity to connect the arable hillslopes with the inner peat area of the kettle hole (Fig. 1), while the majority of events lead to deposition of sediments in the kettle hole surrounding flats (see Figs. 7d; 3, P1 and P4). This statement is supported by surface runoff and sediment delivery monitoring in the study catchment (2015–2019) that has demonstrated that only very few rainfall events caused runoff and associated sediment delivery into the kettle hole (data not shown). Therefore, the study catchment shows a very limited hydrological and sedimentological connectivity between the cultivated area and kettle hole. Hence, tillage translocation in hummocky young morainic regions does also have a pronounced impact on hydrology and biogeochemistry.

4.4 Relevance of tillage erosion and scientific attention

Our results clearly indicate that soil erosion in the study area exceeds the tolerable soil loss rates (according to Schwertmann et al., 1990, 6 Mg ha⁻¹ yr⁻¹ in the study region) and is mainly attributed to tillage erosion (Figs. 6 and 7). During the socialist era (1949–1990), productivist agricultural management strategies were implemented that included land consolidation to merge large fields and the use of heavy farming machines (Forstner and Isermeyer, 2000; Wolz, 2013). For instance, annual ploughing was combined with a recommended practice of episodically using a paraplough (tillage depth ∼0.6 m; Fachbereichsstandard-DDR, 1985) to break the plough pan. The average field size in the region (Quillow catchment = 22 ha) is rather large; this has favoured big farming structures that utilise powerful machinery. However, tillage erosion does not receive reasonable scientific attention (Fiener et al., 2018), even if its effects on yields (Quine and Zhang, 2002; Papiernik et al., 2005), nutrient and carbon cycling (Wilken et al., 2017a; Zhao et al., 2018; Nie et al., 2019) and soil hydrology (Herbrich et al., 2017) are widely known. Globally, tillage erosion has been recognised as an environmental threat in the hummocky young morainic regions that have shallow soils that are subject to dropping yields at hilltop locations (Canada: Pennock, 2003; Tiessen et al., 2007a, b, North America: Li et al., 2007, 2008, Russia: Olson et al., 2002; Belyaev et al., 2005, and northern Europe: Quine et al., 1994; Heckrath et al., 2005; Wysocka-Czubaszek and Czubaszek, 2014). Most arable regions are subject to pronounced tillage erosion (e.g. illustrated in the landscape by tillage banks along downslope field borders; Chartin et al., 2013) but may not show a pronounced impact on yields (Lal et al., 2000). Loess-derived soils with a homogeneous grain size distribution for several metres of depth do not show major differences in soil structure (Blume et al., 2016), while nutrient losses are compensated by fertiliser applications. Another reason for not being a prominent soil degradation mechanism might be that the impacts of tillage erosion are not as visible as those caused by water erosion, which leads to rapid topographical dynamics (rills and gullies) and off-site damages (muddy floods, siltation). However, tillage erosion is a highly important soil redistribution process, taking place on the majority of sloped arable fields, that urgently needs scientific consideration and implementation in soil conservation management by policy makers.

5 Conclusions

In this study, ^239+240Pu was used as a tracer to reconstruct soil redistribution in a hummocky young morainic study catchment under intense arable use. To understand the role of water and tillage erosion in soil degradation patterns, an inverse modelling approach was carried out in the study catchment. The results clearly show that recent soil degradation in the study area is dominated by tillage translocation. Furthermore, tillage erosion has a substantial impact on surface runoff. Tillage forms hydrological depressions at the downslope border between the cultivated field and the kettle hole that limits the hydrological and sedimentological connectivity into the kettle hole and causes deposition of sediments that are transported by water. Soil redistribution by water has no major contribution to soil loss on the catchment hillslopes but causes pronounced deposition in the spatially narrow area of the kettle hole surrounding flat. Within this study, soil erosion up to 17 cm (53 yr)⁻¹ and deposition exceeding 25 cm (53 yr)⁻¹ of recent geomorphological dynamics (since the 1960s) were found. However, even these relatively high erosion rates cannot explain the current soil degradation patterns determined from soil prospection and chemical analysis that show both profile soil truncation and colluviation larger than 1 m. This indicates that tillage erosion might not be a process that exclusively takes place in highly mechanised agro-ecosystems but is potentially causing pronounced soil degradation in smallholder farming structures. Our results clearly underline that tillage erosion is a critically underrepresented soil degradation process that can be the main soil redistribution driver on catchment scale.