The Cross-Verification of Different Methods for Soil Erosion Assessment of Natural and Agricultural Low Slopes in the Southern Cis-Ural Region of Russia

Abstract

The conventional measuring methods (runoff plots and soil morphological comparison) and models (WaTEM/SEDEM and regional model of Russian State Hydrological Institute (SHI)) were tested with regard to the Southern Cis-Ural region of Russia, along with data from rainfall simulation for assessing soil erosion. Compared with conventional methods, which require long-running field observations, using erosion models and rainfall simulation is less time-consuming and is found to be fairly accurate for assessing long-term average rates of soil erosion and deposition. In this context, ¹³⁷Cs can also be used as a marker of soil redistribution on the slope. The data of soil loss and sedimentation rates obtained by using conventional measuring methods were in agreement with the data based on the used contemporary modeling approaches. According to the erosion model calculations and data on the fallout of radionuclides in the Southern Cis-Ural (54°50–25′ N and 55°44–50′ E), the average long-term annual soil losses were ~1.3 t·ha⁻¹ yr⁻¹ in moderate (5°) arable slopes and ~0.2 t·ha⁻¹ yr⁻¹ in meadows. In forests, surface erosion is negligible, or its rates are similar to the rate of soil formation of clay–illuvial chernozems. The rates of soil erosion and sediment deposition on the arable land obtained using different methods were found to be very close. All the methods, including the WaTEM/SEDEM, allowed us to measure both soil erosion and intra-slope sedimentation. The regional SHI model fairly accurately assesses soil erosion in the years when erosion events occurred; however, soil erosion as a result of snowmelt did not occur every year, which should be taken into account when modeling. The concentrations of ¹³⁷Cs in the topsoil layer (0–20 cm) varied from 0.9 to 9.8 Bq·kg⁻¹, and the ¹³⁷Cs inventories were 1.6–5.1 kBq·m⁻², with the highest values found under the forest. The air dose rate in the forest was higher than in open areas and above the average of 0.12 μSv·h⁻¹ on the slope (0.1 μSv·h⁻¹ in the meadow and 0.08 μSv·h⁻¹ on the arable land), with the value increasing from the watershed to the lower part of the slope in all the areas. The γ-background level in the studied ecosystems did not exceed the maximum permissible levels.

1. Introduction

Soil erosion is a natural process of the destruction of rocks and soils by surface water and wind, including the separation and removal of fragments of material, as well as their further transport and redeposition. About 11% of the total global land area is subject to erosion, and the annual loss of fertile soil is about 24 billion tons [1,2]. In the European Union alone, soil loss due to water erosion is 970 million tons per year [3,4], and in Russia, more than 500 million tons of the upper, most fertile soil horizon is washed away from arable slopes every year (even for entirely plowed areas including flat territories, it amounts to ~1 mm) [5]. Erosion results not only in direct soil loss and a decrease in agricultural land productivity, including crop yields, but it also leads to siltation in aquatic ecosystems and deterioration of water quality in aquabiocenoses [6].

Soil losses from agricultural fields are observed in many regions of Russia [7,8,9,10], including the study area, namely, the Republic of Bashkortostan (RB). The part of eroded agricultural land in RB comprises about 64% [11], and the main causes of soil erosion in the RB are (i) the rugged relief with long and steep slopes [11]; (ii) the continental climate with harsh and cold winters, leading to freezing soils and erosion processes during snowmelt [12]; and (iii) the high frequency of erosive rainfall events, amounting to 37.5% of the annual precipitation [13]. It is worth noting that land resources are very important for the RB because it is a major agricultural region, producing ~3.1% of all crop and livestock products in the country and ranking 6th among 89 regions of Russia, and it is therefore crucial to monitor/assess erosion losses for the subsequent adoption of soil-conservation measures.

Previous large-scale studies of erosion processes in the RB were carried out in the 1950–1980s [14]. The obtained erosion rates were mainly based on traditional field methods of accounting for erosion losses: soil morphological and pin methods or those based on the turbidity of water flowing from slope catchments or runoff plots [15]. The need for conducting additional studies on the erosion processes in the RB, in particular in the Southern Cis-Ural as a region of active manifestation of soil erosion, is also driven by the contemporary climatic changes that have occurred on planet Earth over the past several decades. Numerous studies have shown that global warming has led to a reduction in erosion losses in Russia [16,17], Europe [4], and China [18] but creates a small increase in such losses in the USA [19]. According to the data provided by the Intergovernmental Panel on Climate Change (IPCC) and other researchers, the average temperature in the world over the past 50 years has risen by 0.6 ± 0.2 °C [20,21]. In the Cis-Ural region, over the past decade, climate changes also occurred: the average monthly and average annual air temperatures increased, and the amount of precipitation and the average annual wind speed decreased [22]. This means that the likelihood of developing erosion due to rainfall and wind has decreased in the region. In winter, the depth of soil freezing has significantly decreased, which, according to Barabanov et al. [16], leads to an increase in the infiltration of water during snow/ice melt and, accordingly, to a decrease in surface runoff and soil erosion. All of the above changes are supposed to lead to a general decrease in the intensity of erosion processes in the region.

Erosion and accumulation processes are very variable across time. The total long-term soil erosion consists of individual erosion events, the intensity of which varies greatly and depends on both relatively stable and dynamic natural factors. The assessment of changes in erosion and accumulation processes over time is a difficult task and poorly studied in some regions. Moreover, simultaneous assessments are necessary to understand the accuracy of measurements of erosion-accumulative processes and their prediction adequacy. Currently, there are many different methods for quantifying soil erosion. The most accurate of these are field measurement methods. Such measurements, however, cover very short periods and can be used with some assumptions to calculate the average long-term rates of soil erosion in the past or forecast the future. On the other hand, methods that evaluate soil erosion over long periods of time (for example, soil morphological) are not very accurate. Actually, each of the field methods has its drawbacks. Zhidkin et al. [23] found that soil erosion estimates carried out using different methods for representative slopes in the southern part of the European territory of Russia can differ by more than 10 times. Recently, methods of modeling soil erosion have become increasingly popular, their advantage being the rapidity in obtaining results, low cost, and reduced labor intensity. Nevertheless, in order to avoid gross errors, it may be worthwhile to use several methods at the same time and identify the most suitable one for specific natural and climatic conditions.

1.1. The Review and Description of Methods and Models for Soil Loss Estimation Used in Study

1.1.1. Short-Term Sediment Yield Measurements

Usually, these methods include field observations using runoff plots and sediment tanks and are aimed at the calculation of soil loss after erosional events. At the same time, obtaining reliable long-term average data on soil losses or accumulation caused by erosion during snowmelt, erosive rainfall events, and winds of different intensities requires conducting observations for 20 to 25 years [24]. Since the values of soil erosion vary from year to year, and field monitoring studies are required to estimate long-term average values, an express method for assessing soil erosion or deposition can possibly give better results. We are referring, in particular, to rainfall simulation (see Section 2.3.1) and methods/models described below.

1.1.2. Medium-Term Radiocesium (¹³⁷Cs) Technique

The advantages of the radioactive tracer method include the ability to reliably estimate the long-term average annual rate of erosion-accumulation processes over a known time interval. This is especially important when studying the dynamics of the rates of loss and accumulation of soil sediments due to climate fluctuations and changing land use conditions. New research approaches will help to obtain a better understanding of the mechanisms of erosion processes due to changes of climatic, landscape, and soil conditions. Furthermore, updated data on the dynamics of soil erosion at different time periods will allow us to know the current erosion level and adequately predict the loss of soil mass, humus, and nutrients depending on fluctuations in climatic characteristics, the change of crop rotation, and soil cultivation methods and, accordingly, to develop more effective systems of soil-protection measures [7].

Radiocesium (¹³⁷Cs) for investigating soil erosion and the accumulation level was first used in the 1970s in the USA, Australia, and Canada [25,26]. It was noticed that when cesium falls on the Earth’s surface, it is quickly and firmly adsorbed and fixed by soil particles, mainly by the clay fraction [27,28]. During erosion-accumulation processes, it is fine particles of soil that are primarily transported along the slopes [29,30], with ¹³⁷Cs embedded in suspended sediments [31]. Thus, ¹³⁷Cs can be used as a tracer of sediment redistribution. At the same time, the vertical migration of ¹³⁷Cs in the soil is very slow [32,33]; therefore, based on the data on its concentration in different layers of the soil profile, one can also judge the dating of the accumulation of erosion layers—that is, use ¹³⁷Cs as a chronomarker [34]. At the same time, it is possible to assess the rates of denudation in different elementary landscape units of the relief based on data about the inventory of radioactive isotopes of ¹³⁷Cs in the soil. Sedimentation of eroded soil particles from the upper slopes is evidenced by higher ¹³⁷Cs inventories than on the reference plot. Conversely, lower ¹³⁷Cs inventories suggest erosion. Studies conducted in different countries and physico-geographical conditions have shown the high accuracy and effectiveness of the radiocesium method in determining erosion losses [35,36,37,38,39]. In Russia, ¹³⁷Cs is widely used as a chronomarker for the evaluation of sedimentation rates in different sediment sinks for the two time windows (1963–1986 and 1986–the time of sampling) ([17,40], etc.). However, it is more difficult to use the ¹³⁷Cs technique for the evaluation of erosion rates in most parts of European Russia because of the comparable fallout of weapon testing (before 1986) and Chernobyl-derived ¹³⁷Cs [41].

Radiocesium has entered the soil, and biogeocenoses have entered the soil mainly as a result of tests of nuclear weapons, peaceful nuclear explosions, and accidents at nuclear power plants. The main sources of soil contamination with artificial radionuclides in the Southern Cis-Ural region and RB is associated with the following events:

(i) The formation of global radioactive contamination of the Earth’s biosphere by artificial radionuclides began in the middle of the last century as a result of massive tests of nuclear weapons (501 tests from 1954 to 1963). The significant deposition of bomb-derived ¹³⁷Cs can be attributed to 1954, where two peaks are observed: a smaller one in 1957 to 1959 and the maximum one in 1963 [42,43]. After 1964, the intensity of the fallout quickly decreased due to the signing in Moscow of the International Treaty banning the testing of nuclear weapons in the atmosphere, outer space, and underwater on 5 August 1963—Limited Test Ban Treaty [44];

(ii) An additional contribution to the formation of global contamination of biogeocenoses by radioisotopes was made by local technological disasters and accidents at nuclear facilities. A catastrophic accident at the Chernobyl nuclear power plant (NPP) (Figure 1a) occurred in the USSR in 1986. During the explosion, 1850 Pbq of radionuclides were released from the destroyed reactor, including 270 Pbq of radioactive cesium. The spread of radionuclides occurred on a global scale. The major part of radionuclides deposited in the USSR fell out in the Ukrainian SSR, Byelorussian SSR, and the Central Economic Region of the RSFSR [45];

(iii) Another contamination by radionuclides in the South Urals is linked with the Production Association (PA) Mayak, which was involved in the storage, processing, and disposal of radioactive waste. PA Mayak dumped low- and medium-level radioactive waste into the Techa River for several years in the middle of the 20th century. A substantial discharge of radioactive waste was carried out in the 1950s, amounting to 102 Pbq, including 12.4 Pbq of ¹³⁷Cs. Later, in 1957, at the same place, a thermal explosion occurred in the radioactive waste storage, resulting in the release of radionuclides into the environment outside the industrial site with a total activity of 74 Pbq, including 0.2 Pbq of ¹³⁷Cs [46]. In the years 1948 to 1998, as a result of the facility operation and emergency situations, PA Mayak released more than 1.8 × 10¹⁷ Bq of radionuclides into the air and water bodies, and contamination occurred over 25,000 km² [46,47]. Since the RB is located in the immediate vicinity of the East-Ural Radioactive Trace (about 100 km west of PA Mayak, Figure 1a), there is a possibility of contamination of the territory with man-made radionuclides via wind transport;

(iv) Another additional factor of ecosystem radioactive pollution in the region may be associated with the fact that in the RB in the years 1960 to 1970, a series of peaceful underground nuclear explosions were carried out. The first series of explosions (years 1965, 1970) was intended to study the possibility of intensifying oil production (the “Butan” object—40 km east of the city of Meleuz (Figure 1b), wells No. 617, 618, 622, 401, and 403). The second series (the years 1973 and 1974) were pilot works on the disposal of biohazard industrial wastewater from petrochemical plants into deep horizons—2000–2100 m (Kama-1 facilities—28 km west of Sterlitamak; Kama-2—22 km west of Salavat) [48]. In recent years, due to the tightening of environmental regulations, the number of enterprises using atomic energy has significantly decreased in the RB (in 1999, there were 149; in 2002, there were 69) [49].

1.1.3. Soil Profile Truncation Method

The depth of erosion/deposition is estimated based on a comparison between the composition of a typical undisturbed soil profile and soil profiles on cultivated land with the same soil type and position on the slope [50,51]. The mean annual soil loss/gain is determined using the following equation [52]:

$$R={\displaystyle \frac{Mr-Ma}{T}}=k{\displaystyle \frac{Mr}{T}}$$

(1)

where R = the mean annual erosion/deposition rates, mm·yr⁻¹; T = the duration of land cultivation, years; Mr = the thickness of the A + AB + illuvial B1 soil horizons at an undisturbed location, mm; Ma = the thickness of soil horizons A + AB + B1 on cultivated fields, mm; and k = the coefficient representing the eroded/depositional part of Mr.

Difficulties in the practical application of the method are associated with the natural variability of the upper horizons of soils due to heterogeneity of the composition of surface deposits, moisture, and other factors of soil formation [53]. It is necessary to note that the soil truncation method provides a basis for identifying the total reduction/increase in soil horizon depths ∆h = M_r − M_a due to the combined influence of different erosive processes (water, wind, and tillage erosion and soil losses with root crop harvesting) for the entire period of agriculture [52].

1.1.4. The WaTEM/SEDEM Model

In recent years, a large number of soil erosion models were developed [54,55,56,57]. In this study, we applied WaTEM/SEDEM model. WaTEM/SEDEM is a spatially distributed erosion and sediment transport model based on the RUSLE equation [58]. The model was developed by the Physical and Regional Geography Research Group in K.U. Leuven, Belgium [59,60,61]. WaTEM/SEDEM has been applied to catchments representing a wide range of environmental conditions. In Europe, these include: Belgium [60,61], Czech Republic [62,63], Slovenia [64,65], Slovakia [66], Italy [67], Spain [68,69], Turkey [70], Hungary [71], etc. Model applications have also been performed for areas in other parts of the world, e.g., in China [72,73,74], Ethiopia [75], Brazil [76], Australia [77], etc. WaTEM/SEDEM has also showed reliable results in soil erosion modeling for different soil types in the Central Russian Upland [78,79]. The model structure, wide geography of its application, and the similarity of natural conditions in some of the studied areas and Cis-Ural provide a basis for the application of the model in this region. A moderate amount of input data required to run the WaTEM/SEDEM is an important aspect because of limited field data often available in practice.

1.1.5. The Snowmelt Erosion Models and the Regional Russian Model of State Hydrological Institute

At present, the development of rain erosion models is far ahead of that for soil erosion during snow melting. One of the main reasons for this is the complexity of processes during soil thawing—in particular, the detachment of soil aggregates and formation of the rill network. There are models that describe the processes of snow melting [80] or calculate runoff volume [81,82], but not soil loss. Despite this, there were attempts from some researchers to devise snowmelt erosion models. For example, Ollesch et al. [83] used the integrated winter erosion and nutrient load model (IWAN) for a small low-mountain catchment in Germany, and IWAN showed good performance for modeling snowmelt erosion and calculating sediment yield. Weigert et al. [84] simulated snowmelt erosion by using the EROSION 3D model. It is worth noting that EROSION 3D model initially was intended for determination of soil loss mainly during rainfall events [85]. The snowmelt-induced sediment yield is possible to estimate by using the EUROSEM-GRIDSEM, because this model could differentiate soil losses caused by snowmelt and rainstorms [86].

In Russia, there have also been attempts to develop snowmelt erosion models, since this type of soil degradation is very common for the climatic conditions of the country. The first mathematical model of snowmelt erosion was proposed by the State Hydrological Institute [87]. In addition, a physically based equation for soil particle detachment from thawing soil has been suggested by Kuznetsov et al. [88,89], and promising results have been obtained from its application to different soil types [90]. However, the proposed equation includes parameters that require special investigations, which has so far been limiting its practical application. Subsequently, a detailed snowmelt erosion model for arable slopes (<6°) was proposed by Sukhanovskii [91]. In our study, the snowmelt erosion was estimated using the regional Russian model of the State Hydrological Institute (SHI) or the SHI model modified by Larionov [92] and Litvin et al. [8].

Given the limitations of the existing soil loss data, advancements in erosion modeling techniques, climate change concerns, and potential radionuclide contamination in the RB, our study aimed to (i) assess soil erosion rates in arable, forest, and meadow watersheds using conventional (runoff plots, soil truncation, and rainfall simulation) and modern (modeling and radionuclide marker) methods and compare the obtained results; (ii) evaluate the WaTEM/SEDEM and SHI models, and for radionuclides, employ the ¹³⁷Cs technique and assess radiological conditions in the Cis-Ural region; and (iii) investigate changes in soil properties based on slope position and land use.

2. Materials and Methods

2.1. Study Site and Main Geomorphological and Climatic Characteristics

The studies were carried out on the territory of a Water-Balance Station located 20 km northwest of Ufa (54°50′25″ N; 55°44′50″ E, 165–170 m above sea level; Figure 1). The study areas were natural catchments of arable land, meadows, and forests formed on a south-facing moderate (5°) slope. The main soil type is clay–illuvial chernozem (Luvic Chernozems (Aric, Pachic)), which is the most common soil in the Southern Cis-Ural, covering more than 70% of moderate slopes. The south-facing slopes are more susceptible to erosion than others [93], while during snowmelt period, they are subject to a high level of insolation and have less-developed vegetation and dry soil in summer.

The study area is characterized by a temperate continental, rather humid climate with an average annual air temperature of +3.8 °C. The average annual precipitation is 589 mm, with about a third falling in solid form (from November to March) [11]. The hydrothermal coefficient of Selyaninov [94], representing the level of natural moisturization during the growing season, is 1.0–1.2, which means that the studied territory is slightly arid. The study area, and the South Cis-Urals in general, is considered an erosion-prone region: the potential and maximum possible rates of erosion vary from 10 t·ha⁻¹ yr⁻¹ [93] and 20 t·ha⁻¹ yr⁻¹ [95] to 60 t·ha⁻¹ yr⁻¹ [11], or between about 0.6 and 5 mm per year, respectively.

Fodder crops prevail in the crop rotation system of the studied arable catchment. Land management is carried out in a conventional way with plowing of arable land to a depth of 20–25 cm; the plowing has been conducted annually in late April or early May since 1960. Weed vegetation growing on the meadow catchment consists of field chamomile (Anthemis arvensis), field sow thistle (Sonchus arvensis), absinth sage (Artemisia absinthium L.), and corn lily (Convolvulus arvensis L.). The total vegetation projective cover of the slope is 85–95%. Sometimes, cattle from nearby farms graze on the meadows (approximately 10–30 heads per ha, 5–20 days during the vegetation period), which can lead to overgrazing and erosion. The vegetation of the forest catchment consists of broad-leaved trees: small-leaved limetree (Tilia cordata), common elm (Ulmus laevis), and white poplar (Populus alba). The average height of woody vegetation is 10 to 20 m, the thickness of trunks is 0.2 to 0.4 m, and the density is 15 to 20 trunks per 10 m².

2.2. Field Work and Sampling and Laboratory Analysis

On each catchment, soil profile pits (2 m long, 0.5 m wide, and 0.5 m deep) were dug in three relief landscape elements (upper, middle, and lower part of the slope (accumulation zone)) (Figure 1c). A total of 9 profiles were excavated to determine sampling and morphological soil property description (and determination of soil type); the thicknesses of the humus-accumulative (A + AB) horizons were also measured. Additionally, to ensure statistically reliable measurements, near each soil pit, the thicknesses of A + AB horizons were measured (3 times; 27 in total) using a core soil sampler (probe). Samples from the soil profile were taken from layers at 0–2.5, 2.5–5, 5–10, and 10–20 cm and further with an interval of 10 cm to a depth of 50 cm. The soil samples, each about 1 kg in weight, were taken from each layer and from 3 different parts of the profile using a small shovel and stored in a plastic package, with total of 3 replications from one soil pit. Specifically for ¹³⁷Cs technique of soil erosion estimation (see Section 1.1.2), additional 3 soil reference pits (Figure 1c) were excavated on flat parts in unforested sites, where no erosive processes occur; the soil sampling procedure was similar to that for catchment pits.

Soil samples were dried in an oven (at 90 °C) to constant weight, ground in a mortar, and sieved through a 2 mm sieve. Then, the activity concentration of ¹³⁷Cs was determined via gamma (γ) spectrometry using an MKS-01A “MULTIRAD” (Scientific production company “Doza”, Moscow, Russia). Measuring times (usually 30 min per sample) were chosen to be such that a statistical uncertainty of less than 5% could be assured. All measured activity concentrations were corrected for radioactive decay to each sampling date. In the field, air dose rate was measured (1 m above the ground) in different parts of the slope with a search dosimeter–radiometer MKS/SRP-08A (Scientific production company “Amplitude”, Zelenograd, Russia). This γ survey meter continuously displays the air dose rate. The meter was turned on in the field, and once the display stabilized (usually after about 30 s), three measurements were recorded within 1 min. The average of the three measurements was calculated.

Since it is the topsoil layer that is susceptible to water-erosion processes, particle size distribution (using the standard pipette/sedimentation method of Kachinsky–Robinson–Kohl [96]) and main agrochemical parameters (obtained using the methods described by Sokolov [97]) were determined. These chemical soil properties included organic carbon (C_org) content, pH (H₂O), exchange cations (Ca²⁺ and Mg²⁺), alkaline hydrolyzable nitrogen (N_alc), plant-available (P_av), and total (P_t) phosphorus. In particular, C_org was determined using Tyurin titrimetric (wet combustion) method with Nikitin’s modification with spectrophotometric termination according to Orlov and Grindel; soil pH using glass electrode and 1:1 (weight:volume) solution of soil to water; and Ca²⁺ and Mg²⁺ using the trilonometric method. N_alc was extracted at 1 mol/L KCl (1:2.5 soil to solution) according to the Cornfield method; P_av and P_t were extracted at 0.5 mol/L CH₃COOH at a 1:2.5 soil-to-solution ratio using the Chirikov method.

2.3. Soil Loss Estimation

2.3.1. Field Soil Erosion Measurements/Observations and Simulation Experiments

The measurements of snowmelt erosion were carried out in the investigated area from 2009 to 2019. In 2008, long runoff plots (200 × 20 m) were constructed on each catchment (forest, meadow, and arable). The runoff volume and soil loss were monitored using Thomson sediment weirs, which were installed in the lower parts of slope on each plot. Soil loss was assessed by measuring the content of suspended sediments using turbidity calculations [15] and a sediment tank installed at the outlet of the sediment weir. Studies have also been conducted to assess erosion due to rainfall. For this purpose, small runoff (10 × 5 m) plots were used in the field and in the same investigating area. The plots have a similar inclination (~5°) and include a meadow and arable land. The plots were subjected to rainfall simulation using different types of sprinklers. In studies with rainfall simulators, the duration (time) of irrigation varied from 16 to 1245 min, and rainfall intensity was from 0.11 to 2.1 mm·min⁻¹. Also, experiments were conducted in laboratory while using a rainfall simulator. Different slope gradients (3 and 7°) and intensities of artificial rainstorms (2, 4, and 6 mm·min⁻¹) were used for plowed and meadow sites. For both experiments (field and laboratory), 40–42 mm was the daily absolute maximum precipitation in the region. A more detailed description of the methodology and measuring procedures of snowmelt erosion and rainfall simulations can be found in the previous works [12,13,30].

2.3.2. Modeling

• WaTEM/SEDEM. Calculation of soil erosion rates was carried out using mathematical models separately for rainfall and snowmelt erosion only for cultivated slopes. Rainfall erosion was estimated using the WaTEM/SEDEM v. 2004 model [59,60]. We used a computer program developed at the Laboratory for Experimental Geomorphology (LEG), K.U. Leuven, Belgium [61]. The calculation algorithm is based on the RUSLE [58]. The following input parameters were used:

R = the rain erosivity factor was taken from the Global Rainfall Erosivity database [98]. On the studied site, R = 260 MJ∙mm∙ha⁻¹∙h⁻¹∙yr⁻¹.

K = the soil erodibility factor was calculated on the basis of the analytical properties of the soil in accordance with the equation [58]. The clay fraction content (<0.002 mm) of the studied soil varied from 24.8 to 29.9%, the silt fraction content (0.002–0.05 mm) accounted for 62.7–68.4%, the very fine sand fraction content (0.05–0.1 mm) was 3.6–4.8%, and the C_org comprised 4.3–5.8%. Thus, the K-factor varied from 37 to 45 t∙ha∙h∙ha⁻¹∙MJ⁻¹∙mm⁻¹ and averaged = 41 t∙ha∙h∙ha⁻¹∙MJ⁻¹∙mm⁻¹.

A digital elevation model was compiled based on a topographic survey with a 1:500 scale with the recommended (for WaTEM/SEDEM) cell sizes of 20 × 20 m.

C = the crop and management factor for the arable plot was 1 due to the fact that this plot was kept without vegetation during the experiment.

• SHI model. As was noted in the section above, the snowmelt-induced sediment yield in this study was estimated using the SHI model. Water reserves in snow were obtained by field measurements for the period 2009–2019. The snow–water equivalent was measured by calculation of snow height and density. It varied from 68 to 211 mm in different years. Soil erodibility factor was calculated based on C_org content and particle size distribution of the soil. LS factor was calculated on the basis of digital elevation model constructed from the topographic map. C-factor = 1 due to the fact that the arable plot was kept without vegetation.

2.4. Statistical Analysis

Statistical data processing was performed in Microsoft Excel 2010 (Microsoft Corporation, Redmond, WA, USA), with the error for the confidence probability p = 0.95, not exceeding 10%.

3. Results and Discussion

Overall, no visible areas with noticeable water erosion were detected in different years of the study, either in rainfall or snowmelt periods, during the reconnaissance survey of the slope of the forest and meadow catchments. Also, at the forest site, no surface runoff was observed from 2008 to 2019; runoff only occurred at the meadow and arable site, with a low runoff intensity and concentration of suspended sediments in the meadow. In contrast to the meadow site, small rills were identified (1 to 4 cm wide and 3 cm deep) after snowmelt and in the fall on the arable land in the middle part of the slope, and brown particles with a transitional (AB) and illuvial (B) horizon were sometimes found on the soil surface. Occasionally, visible sediment accumulations were found in the lower part of the arable land site, with the soil having a darker color and higher moisture there.

3.1. Results of Field Soil Erosion Measurements/Observations and Simulation Experiments

Detailed results about soil loss caused by snowmelt erosion and rainfall simulations can be found in previous works [12,13,30]. The results of soil erosion rates in those studies are summarized, averaged, and given in (

Table 1. The average values of soil erosion and deposition in the studied areas, t·ha⁻¹·yr⁻¹.

Location	Assessment Method				Mean Value of Soil Loss for Slope (All Methods)
	WaTEM/SEDEM (Rainfall Erosion), SE ± 0.01	SHI Model (Snowmelt Erosion), SE ± 0.01	Soil Profile Truncation, SE ± 0.1	Field Monitoring/Simulation, SE ± 0.01
Forest			U (±0.2), M (±0.2), L (±0.2)	U (±0.01), M (±0.2), L (±0.01)	~0
Arable land	U (0); M (−2.7); L (+1.9)	M (−0.14)	U (±0.2); M (−1.4); L (+1.5)	U (±0.01); M (−0.6); L (+0.7)	~1.3
Meadow			U (±0.2); M (−0.3); L (+0.3)	U (±0.01); M (−0.2); L (+0.1)	~0.2

Notes: U—the upper, M—the middle, L—the lower part of the slope; (+)—sediment accumulation; (−)—soil loss; (±)—soil loss/deposition/pedogenesis; SE—standard error of measurement.

).

3.2. Results of Application of Soil Truncation Method

The calculation of soil erosion rates using this method shows no soil loss at the forest site, which is in agreement with the granulometric (Figure 2) and agrochemical (

Table 2. Morphological and agrochemical properties of the slope soils (0–20 cm).

Location	Slope Part	Humus-Accumulative Horizon Thickness (A + AB), cm	Corg, %	pH H2O	Exchange Cations, cmol(+) kg−1		Nalc	Pav	Pt
					Ca2+	Mg2+	mg·kg−1
Forest	U	64 ± 4	6.0	6.0	45	9	306	3.4	271.1
	M	62 ± 4	5.9	5.9	43	10	295	3.2	255.4
	L	66 ± 5	6.0	6.0	44	10	301	3.4	258.7
Arable land	U	42 ± 3	5.0	5.9	35	12	175	2.5	188.1
	M	35 ± 3	4.3	6.4	29	13	168	2.0	135.6
	L	51 ± 4	5.8	6.2	38	17	205	2.9	230.3
Meadow	U	45 ± 3	5.3	5.9	33	11	231	2.7	219.3
	M	43 ± 2	5.1	6.1	29	10	219	2.5	204.4
	L	47 ± 4	5.4	6.0	31	12	240	2.8	228.8

Note: U—upper, M—middle, L—lower parts of the slope.

) soil composition there. The content of clay and sand in the forest soil was comparable; regardless of the position on the slope, the content of C_org and nutrients was approximately the same. This is consistent with the results of other authors suggesting that in the presence of erosion on the slope, along with the water–erosion flow and suspended sediments from the slope soil, predominantly silt and clay fractions with a high content of nutrients are washed out [99,100,101]. Even at the meadow site, some asynchrony was seen in the particle size distribution, specifically, a slight accumulation of the silt and clay fractions in the lower part of the slope and their removal from the middle part. A more pronounced trend was observed for the arable land, with the proportion of sand increasing in the middle part of the slope and that of silt and clay increasing in the lower part. This is indicative of the ongoing erosion processes on the arable land and transport of fine soil particles enriched with C_org and nutrients along the slopes and deposition at the slope toes.

In the arable plot, the average rate of soil erosion was estimated to be ~1.3 t·ha⁻¹·yr⁻¹, and in the meadow, it was estimated to be ~0.2 (Table 1). It is important to note that according to the Russian scale of assessment of the soil degradation degree [87] with respect to the morphological properties of the soil, the slopes with forests and meadow soil are characterized as non-eroded. Soils on the arable catchment are not subject to erosion either, being slightly eroded in the middle part of the slope and poorly inwashed in the lower part. The thickness of the humus-accumulative horizons in the upper and middle part of the arable land slope is less than on the meadow site (Table 2), the average being 7 ± 3 cm along the slope. If we assume that the meadow site is not prone to erosion and that plowing of the same meadow began in 1960 (the Virgin Lands campaign), the soil erosion from the arable land is estimated to be ~1.2 mm or 1.4 t·ha⁻¹·yr⁻¹, which correspondents to the soil morphological method results.

3.3. Results of WaTEM/SEDEM and SHI Modeling

According to the mathematical model (WaTEM/SEDEM) calculations, the average rainfall soil erosion rate on the arable plot is 2.7 t·ha⁻¹·yr⁻¹, and average sediment accumulation is 1.9 t·ha⁻¹·yr⁻¹ (Table 1). The accumulation of sediments is predominantly noted in the lower part of the plot (Figure 3). The soil erosion is manifested locally in small areas, especially with greater steepness. In addition to soil washout in the middle parts of the slope, a significant accumulation of sediments is noted in the lower (stock zone) part.

Snowmelt soil erosion was calculated over 10 years based on the data of standard input SHI model parameters (K-, LS-, and C-factors) and existing field measurements (the snow–water equivalent data). According to the SHI model, the snowmelt soil erosion rates varied up to five times in different years—from 0.05 to 0.28 t·ha⁻¹·yr⁻¹ (Table 1, Figure 4). On average, the snowmelt soil erosion according to WaTEM/SEDEM is estimated to be 0.14 kg·ha⁻¹·yr⁻¹. Thus, the total erosion losses (due to rainfall and snowmelt events) calculated on the basis of the mathematical models (WaTEM/SEDEM and SHI) for the arable land amount to 1.0 t·ha⁻¹·yr⁻¹.

According to the 10 years of field monitoring, the snowmelt erosion was registered only four times and varied on the arable plot from 0.1 to 0.3 t·ha⁻¹·yr⁻¹ (Figure 4) and 0.03 to 0.2 t·ha⁻¹·yr⁻¹ on meadows [12]. Moreover, the runoff was only formed in the years when the depth of soil freezing was more than 40 cm. In addition, no snowmelt-induced soil erosion was observed under the forest site.

The variation of snowmelt erosion is caused not only by varying water reserves in the snow but also by the depth of soil freezing. The absence of this parameter (freezing depth) during the mathematical modeling of snowmelt-induced soil loss causes significant errors in calculating the rates and volumes of runoff. As a result, the total (sum) snowmelt induced soil loss according to the data of field runoff measurements amounted to 0.7 t·ha⁻¹·yr⁻¹ for 10 years of observations, and mathematical estimates according to the SHI model yielded an almost 2.5-fold excess or 1.7 t·ha⁻¹·yr⁻¹. Under the conditions of global climate change, snowmelt erosion is less and less manifested in the forest–steppe zone of Russia [16]; thus, it is necessary to take into account this fact and soil depth freezing during snowmelt erosion modeling.

3.4. Results of Application of ¹³⁷Cs for Evaluation of Relative Soil Losses Along the Slope and Current Radiological Conditions in the Southern Cis-Ural

According to the latest data, the average density of ¹³⁷Cs contamination in soils due to global atmospheric fallout in the zone between 50° and 60° N is 2.2 kBq·m⁻² and 1.3 kBq·m⁻² for ⁹⁰Sr [102,103]. At the same time, the soil and vegetation cover in the Ural region is characterized by an increased content of ¹³⁷Cs compared to its average level at the same latitudes of the Northern Hemisphere [104]. According to Mikhailovskaya et al. [105], the background level of soil contamination of ¹³⁷Cs in the Urals between 55° and 57° N is 1.0–7.6 kBq·m⁻², which means that the upper boundary is 3.5 times higher than the average value for these latitudes. The density of ¹³⁷Cs contamination within the studied slope varied in the range of 1.6 to 5.1 kBq·m⁻² (Table S1). These results are consistent with the data on the radiation situation in the territory of the RB available from the national reports on the state of the environment and other literature sources. For example, Gorelov [49] found that the density of radioactive contamination of the soil in the RB due to the isotopes of ¹³⁷Cs and ⁹⁰Sr in the early 2000s varied from 1.1 to 3.9 kBq·m⁻². Based on the latest data from 2016 [106], the average levels of ¹³⁷Cs deposition for 2013–2016 in the RB were within 0.9–1.2 kBq·m⁻². According to regulatory documents [107,108], the degree of contamination of the territory, both within the studied catchments and the RB as a whole, is assessed as “permissible” for the “radioactive substances” (contamination level is less than 1 Ci·km⁻² = 37 kBq·m⁻²), and for the level of ¹³⁷Cs concentration in the soil—it is believed that “the soil is clean”. Thus, the cultivation of crops and conduct of any other economic activity in these territories can be permitted without restrictions.

The air dose rate on the studied slopes did not exceed the permissible level (0.20 μSv·h⁻¹) and varied from 0.07 to 0.12 μSv·h⁻¹ (Figure 5) and was comparable with regional values. According to [49], the radiation background in the RB ranges from 6 to 17 µR·h⁻¹ (0.06–0.17 μSv·h⁻¹). It is worth noting that the γ-background was different depending on the position on the slope and the land use type. In all catchments, the maximum was determined in the lower part of the slope and the minimum in the middle of the slope, which is associated with washout from steeper slopes and the accumulation of sediments enriched with radiocesium in the lowlands. The dose rate in the forest was higher than on meadow soil by ~9 and on arable land by ~23%. Similar results were obtained for territories exposed to radionuclide contamination as a result of accidents at the Chernobyl and Fukushima nuclear power plants. For example, Komissarov and Ogura [109] found a significant (R² = 0.81, p = 0.95) air dose rate dependence on the concentration of radiocesium in the surface soil layers: the higher the activity concentration, the higher the level of the γ-radiation exposure dose rate. Yasumiishi et al. [110] indicated that the air dose rate in elevated areas was lower than at the foot of the slopes. Ramzaev et al. [111] showed that the dose rate in the Bryansk region (Russia) decreased in the ecosystems as follows: forest–meadow–arable land. In Miyagi Prefecture (Japan), the γ-background in the forest was 10% higher than in the adjacent meadow [31], and remediation of contaminated meadows (plowing) reduced the air dose rate by 30–50% [109]. The increased gamma-background in the forest and lower on the arable land is due to the fact that forests are known for their filtering and adsorbing properties. They primarily accumulate radioactive aerosols from the atmosphere [112,113]. In turn, during plowing, the “burial” of a radiation-contaminated layer occurs in the lower horizons, which results in the “dilution” of the radiocesium concentration due to its mixing throughout the entire arable layer.

A visual analysis of the diagrams of radiocesium distribution in soil profiles (Figure 6) shows that in forest ecosystems, the concentrations of ¹³⁷Cs at the corresponding depths are similar for all parts of the slope, with a distinctive common trend: the maximum values occur in the upper layers, and they gradually decrease with depth. Many authors have shown that in undisturbed ecosystems, radiocesium is predominantly located on the upper horizon [32,33,109,114]. The radiocesium inventories in the studied forest soils turned out to be almost the same and varied in the range of 4.9 to 5.1 kBq·m⁻² (Table S1) but were larger than in the open areas. Similar results were obtained in other regions, such as Scotland [115] and Spain [116]. It was also found that the difference in inventories appeared with an increase in slope and distance from the open area. Thus, we can conclude that forests under the moderate slopes of the Cis-Urals are geomorphologically stable units, which is confirmed by calculations and observations showing that there is no erosion under the forest (Table 1). The insignificant soil losses on the meadow site (0.1–0.2 t·ha⁻¹ per year) are indirectly confirmed by the presence of a ¹³⁷Cs peak at a depth of 2.5–5 cm in the lower part of the slope. About 60 years have passed since the time of the maximum ¹³⁷Cs deposition (1954–1963) to the moment of sampling (2019), and during this time, there was a gradual “burial” of the upper radiation-contaminated layer with an average rate of 0.4 mm per year. The deepening of ¹³⁷Cs could also occur due to its migration with a filtration flow or gravitational migration, as well as bioturbation.

3.5. Limitations of Applied Methods and Techniques

It should be noted that each of the methods and approaches used in this study has its own shortcomings and limitations. The accuracy of calculating the washout rates using erosion models is determined by the ability to obtain precise values of the parameters. Thus, the parameters of the topography and soil properties were obtained by us through a detailed survey of the relief and analyses of samples taken from the studied area. However, the accuracy of determining the erosion index of precipitation depends on the location of the weather station relative to the site where the research is carried out. This is due to the high spatial variability of rainfall that forms surface runoff and washout. In our case, we were forced to use data from the global estimates of the erosion index of precipitation. In addition, any calculations using erosion models must be verified using independent data. In particular, the WaTEM/SEDEM model was verified based on field data in the western part of the European territory of the Russia [117,118]. The results obtained showed satisfactory convergence of the calculated data with the estimates obtained using field methods.

Observations of water and sediment runoff at runoff plots can provide objective data on average long-term washout rates if they are conducted for at least 10–12 years [119]. Then, the obtained washout values characterize the average long-term variation of erosion processes associated with the variability of meteorological conditions in different years. In our case, the data on observations of water and sediment runoff during the snowmelt period cover the required time window. However, observations at runoff plots characterize washout from slopes of limited length and one configuration. A much greater variety of slopes and slope catchments are used as arable land in the studied region.

The soil profile truncation method allows us to estimate the total soil losses from water and mechanical erosion, as well as soil losses during the harvesting of root crops [52]. The main disadvantages of this method are the difficulty in determining the total duration of plowing for the studied area and the fairly high natural variability of soil horizons’ A and A + B thickness. The area we studied began to be plowed relatively recently, and therefore, the plowing period is known. However, the accuracy of this method is not high enough due to the significant variability of the soil horizon thickness.

The use of the radiocesium technique also allows us to determine the total soil losses over a known time window [116]. However, in our study, we use the ¹³⁷Cs content only for the identification of erosion and sedimentation zones. This is due to the fact that both bomb-derived and Chernobyl-derived ¹³⁷Cs fallout occurred in the study area. In this case, it is almost impossible to correctly recalculate the ¹³⁷Cs content for the erosion/deposition rates based on the application of conversion models due to the superposition of Chernobyl-derived fallout on the initial field of bomb-derived ¹³⁷Cs fallout origin transformed by erosion–sedimentation processes for the time window from 1963 to 1986 [41].

The use of a set of different methods and techniques in our study allows us partly to overcome the limitations and disadvantages of each of them and to avoid gross errors in assessing the rates of erosion/sedimentation processes. The resulting values of soil loss from erosion on arable land, forests, and meadows, obtained using different methods, differ slightly from each other. However, they fairly objectively characterize the differences between the intensity of erosion processes on the most common types of land use in Southern Cis-Ural. Further research is needed to obtain more detailed quantitative estimates of the rates of washout on arable slopes and slope catchments of various configurations.

4. Conclusions

The degree of development of erosion processes in soils of the same type is known to depend on agroecological and geomorphological conditions, such as the land use type and position in the relief. Observations and calculations have shown that in the Southern Cis-Ural, on south-facing moderate slopes of 5° with forests, soil erosion is absent regardless of the position on the slope and the methods of soil erosion assessments. Erosion is also less pronounced on meadow soil, with soil losses averaging 0.2 t·ha⁻¹·yr⁻¹, while on arable land, the erosion rates were found to be five times higher (~1.3 t·ha⁻¹·yr⁻¹).

The rates of soil erosion and sediment deposition obtained using different methods were very close. In particular, on arable land in the middle part of the slope, the soil erosion rates obtained using different methods ranged from 0.8 to 2.7 t·ha⁻¹·yr⁻¹. In the Southern Cis-Ural region of Russia, intra-slope sediment accumulation processes have a significant impact on the sediment delivery ratio and sediment budget. On the studied arable slope, all the methods diagnosed sediment deposition in the lower part of the slope at close rates from 0.9 to 1.9 t·ha⁻¹·yr⁻¹. Importantly, the WaTEM/SEDEM calculations made it possible to quantify not only soil erosion but also intra-slope sediment deposition.

In recent decades, the depth of soil freezing has changed due to global climate change, which has led to a reduction in soil losses during snowmelt. Over the 10-year period, erosion events during snowmelt were only observed in half the cases. The regional SHI model is shown to fairly accurately assess snowmelt soil erosion in the years when surface runoff events occurred. This model, however, needs to be improved since it does not take into account the depth of soil freezing. Compared with field-based, long-term observation methods, the radiocesium technique, modeling (WaTEM/SEDEM and SHI), and rainfall simulation are all quite fast and convenient. It would be reasonable to test the models and radiocesium technique on other soils in the South Cis-Ural and evaluate the rate of sediment redistribution on steeper slopes.

In contrast to the forest, in the open areas, the mass transfer of fine fractions was noted in the surface soil layer along the slope, while on the arable land, this process was seen to be even more pronounced, which means that as a result of erosion, organic matter and nutrients are lost from the middle part of the slope and accumulate in the lower part.

The analysis showed that despite the relative proximity to the East Ural Radioactive Trace and underground nuclear explosions conducted in the Republic of Bashkortostan, the soils are suitable for economic activity. The activity concentration of ¹³⁷Cs in the surface soil layers ranges from 0.9 to 9.8 Bq·kg⁻¹, and inventories range from 1.8 to 5.1 kBq·m⁻², these quantities being higher within the forest. The air dose rate was also higher in the forest and equal to 0.12 μSv·h⁻¹, while on the meadow site, it was lower (about 0.10 μSv·h⁻¹), yet higher by 20% than on the arable land. For all the sites studied, the γ-background level showed a common pattern: a decrease in the middle and an increase in the bottom of the slope, but it was below the permissible level (0.2 μSv·h⁻¹) at all sites. To sum up, from the radioecological point of view, the studied ecosystems can be considered “clean”.

Thus, any or a combination of the applied methods can be used for areas with similar landscapes and geo-climatic conditions since the obtained results were comparable. Moreover, the SHI model was suggested to be used for predicting snowmelt erosion in the USA, Canada, China, etc. Also, radiocesium on contaminated areas can be used to calculate soil erosion rates in any place in the world.