Modeling Dynamics of Soil Erosion by Water Due to Soil Organic Matter Change (1980–2020) in the Steppe Zone of Russia

: This research aims to evaluate the dynamics of soil loss through soil erosion by water in agricultural lands in steppe areas using a modiﬁcation of the RUSLE2 model from the 1980s to the 2010s. The calculation was performed using a raster model of data that included a model of the slope angle, slope length, soil erodibility, rainfall and snowmelt erosivity factors, types of land use, and cover management factor. All data were taken from open sources. The average soil erosion in the territory studied amounted to 1.48 t ha − 1 year − 1 in the 1980s and 1.72 t ha − 1 year − 1 in the 2010s. The discrepancy with other studies was 12% for the level of the 1980s and 2–7% for the level of the 2010s. The main factor leading to an increase in soil loss was soil erodibility due to the loss of soil organic matter, which affected about 52% of the studied lands. The increase in the amount of soil loss occurred against a background of compensating processes: reduction in precipitation and climate change (getting drier), as well as the overgrowth of agricultural lands with natural steppe vegetation. Modiﬁed model RUSLE2 has shown good results correlated with other studies for the research area.


Introduction
Water is an essential natural resource for sustaining life on Earth, but its unchecked movement and excessive flow can also have devastating consequences.One of the most significant impacts of water on the environment is soil erosion by water, which is the process of soil and rock being worn away by the force of water [1].Soil erosion by water has been a natural process throughout the Earth's history, but human activities such as deforestation, intensive agriculture, and urbanization have accelerated this process, causing severe environmental problems [2].The adverse effects of soil erosion by water on the environment include a loss of soil fertility, reduced crop yields, sedimentation in water bodies, and the destruction of habitats [3][4][5][6][7].Furthermore, soil erosion by water constitutes one of the most potent processes of soil and land degradation, resulting in the degradation of the quality of ecosystem services provided.These services hold significance not only for human populations but also for the entire biocenosis [2,8].Therefore, understanding the mechanisms and impacts of soil erosion by water is critical for developing effective strategies to mitigate its effects and protect our environment [9].
With the increasing effects of climate change, soil erosion by water has become an even more pressing issue.Rising temperatures, altered precipitation patterns, and more frequent extreme weather events can exacerbate erosion rates and intensify the damage caused by erosion [7].Climate change can affect erosion processes in various ways, such as by altering Agronomy 2023, 13 the amount, intensity, and frequency of rainfall and modifying the vegetation cover and the topography of landscapes.These changes can increase erosion rates, cause soil loss, and reduce agricultural land productivity [10].Furthermore, the impacts of climate changeinduced erosion are not limited to specific regions or countries but have global implications, as the sediment and pollutants carried by eroding water can travel long distances and affect ecosystems far from the source [10].To mitigate the impacts of soil erosion by water in the context of climate change, it is crucial to understand the underlying processes, assess the vulnerability of different regions and landscapes, and develop appropriate management strategies [2,11].Arid lands are vulnerable zones in climate change [12,13].The combination of low vegetation cover, high soil erodibility and erratic rainfall makes arid lands particularly prone to erosion.Even small amounts of rainfall can cause significant soil loss and degradation [14,15].It is fascinating to understand how climate changes affect soil erosion.On the one hand, the precipitation decreases, leading to decreased soil losses due to rain.On the other hand, the number of extreme weather events increases, which can lead to an increase in soil losses due to short rainfall but high kinetic energy phenomena.Under the latest Intergovernmental Panel on Climate Change (IPCC) Report [16], there has been a discernible upward trend in the frequency of intense precipitation events observed since the 1950s across numerous geographical regions, including the specific research site under consideration.Furthermore, the report posits a likelihood of further exacerbation in the incidence of heavy precipitation events in the foreseeable future [16][17][18].
This research represents a pioneering effort in the high-resolution spatial modeling of soil erosion by water in steppe areas considering national dynamics (the 1980s to 2010s) in Russia.To our knowledge, this study is the first to comprehensively investigate the dynamics of soil erosion by water processes, utilizing advanced modeling techniques and incorporating detailed spatial data.Other papers focused on a large-scale current level of soil erosion by water [3,[19][20][21][22] or provided data for the last century without spatial distribution [19,23,24].By addressing this critical knowledge gap, our work aims to significantly contribute to understanding erosion patterns and distribution and provide valuable insights for sustainable land management and conservation efforts in the steppe regions of Russia.
The primary objective of this study was to comprehensively assess the dynamics of soil erosion by water over the past four decades (from 1980 to 2020) within the steppe zone of the Volgograd region, Russia.Specifically, our research sought to investigate the influence of climate variability and change on the regional-scale water erosion patterns.This endeavor was motivated by the imperative to gain a deeper understanding of the intricate interplay between climate fluctuations and their impact on the levels of waterinduced erosion, which is vital for informed land management and conservation strategies in the region.

Research Area
The Volgograd region, Russia is located in the southeast of the East European Plain and is the most significant region in the south of Russia.The region's relief is diverse-from a lowland plain in the Volga region in the east to an elevated dismembered territory in the north and west.The DEM and Slope maps are shown in Figures 1 and 2.  The climate of the region is arid, with pronounced continentality.The north-western part is located in the forest-steppe zone, the central and south are in the steppe zone, and the eastern part is in the semi-desert zone, approaching real deserts [25].The map of land use according to Chen and co-authors [26] is shown in Figure 3.The climate of the region is arid, with pronounced continentality.The north-wes part is located in the forest-steppe zone, the central and south are in the steppe zone, the eastern part is in the semi-desert zone, approaching real deserts [25].The map of l use according to Chen and co-authors [26] is shown in Figure 3.The climate of the region is arid, with pronounced continentality.The north-western part is located in the forest-steppe zone, the central and south are in the steppe zone, and the eastern part is in the semi-desert zone, approaching real deserts [25].The map of land use according to Chen and co-authors [26] is shown in Figure 3.The study area exhibits distinct climatic gradients, with average annual precipitation ranging from approximately 500 mm in the northwest to less than 300 mm in the southeast.Moreover, extreme climatic conditions are observed, including absolute maximum temperatures of +42 °C to +44 °C typically occurring during July and August, while absolute minimum air temperatures plummet to −36 °C to −42 °C in the winter months of January and February [25].
The region is located within two soil zones-chernozem and kastanozems.The soil and texture maps are shown in Figures 4 and 5.The study area exhibits distinct climatic gradients, with average annual precipitation ranging from approximately 500 mm in the northwest to less than 300 mm in the southeast.Moreover, extreme climatic conditions are observed, including absolute maximum temperatures of +42 • C to +44 • C typically occurring during July and August, while absolute minimum air temperatures plummet to −36 • C to −42 • C in the winter months of January and February [25].
The region is located within two soil zones-chernozem and kastanozems.The soil and texture maps are shown in Figures 4 and 5.The study area exhibits distinct climatic gradients, with average annual precipitation ranging from approximately 500 mm in the northwest to less than 300 mm in the southeast.Moreover, extreme climatic conditions are observed, including absolute maximum temperatures of +42 °C to +44 °C typically occurring during July and August, while absolute minimum air temperatures plummet to −36 °C to −42 °C in the winter months of January and February [25].
The region is located within two soil zones-chernozem and kastanozems.The soil and texture maps are shown in Figures 4 and 5.The total square area of agricultural land is 76,100 km 2 .The total population of people-2,492,808.The share of urban residents is 77%.The people density is 21.87 people/km 2 [25].

The Technique for Calculating Potential Soil Erosion
We used the RUSLE2 [27] model as a basis to carry out the simulation.However, we approached the choice of the formula for calculating its indicators by considering the available data and the characteristics of the object of study.For example, the RUSLE-2 model needs to work better on soils with a high organic matter content (>4%), which are widespread in the Volgograd region [25,27,28].Because of this, we applied a different algorithm for calculating the K-factor, taking into account the features of the object.In the case of calculating the R-factor, the problem was the absence of pluviographs at most meteorological stations in the region.Furthermore, standard precipitation data are published as the sum of precipitation for 3 h (180 min), which makes it impossible to apply the formula for calculating the R-factor for precipitation that fell in 60 min or less according to the formula laid down in the model.
Two periods were chosen for modeling to assess the dynamics of soil erosion by water-the level of the 1980s and the level of the 2010s.All model calculations for the research area were performed in raster layers with 30 m resolution (Figure 6).
(1) Model calculations were performed as follows: (2) Collecting data for calculating factors of soil erosion.
(3) Choosing optimal model parameters for performing calculations.The total square area of agricultural land is 76,100 km 2 .The total population of people-2,492,808.The share of urban residents is 77%.The people density is 21.87 people/km 2 [25].

The Technique for Calculating Potential Soil Erosion
We used the RUSLE2 [27] model as a basis to carry out the simulation.However, we approached the choice of the formula for calculating its indicators by considering the available data and the characteristics of the object of study.For example, the RUSLE-2 model needs to work better on soils with a high organic matter content (>4%), which are widespread in the Volgograd region [25,27,28].Because of this, we applied a different algorithm for calculating the K-factor, taking into account the features of the object.In the case of calculating the R-factor, the problem was the absence of pluviographs at most meteorological stations in the region.Furthermore, standard precipitation data are published as the sum of precipitation for 3 h (180 min), which makes it impossible to apply the formula for calculating the R-factor for precipitation that fell in 60 min or less according to the formula laid down in the model.
Two periods were chosen for modeling to assess the dynamics of soil erosion by water-the level of the 1980s and the level of the 2010s.All model calculations for the research area were performed in raster layers with 30 m resolution (Figure 6).

Rainfall and Snowmelt Erosivity Factors
Data from the meteorological station of the Roshydromet network located in the Volgograd region were used to calculate the R-factor.The data were provided by RIHMI-WDC [29].Data were separated into two ranges: 1970-1998 (Rainfall erosivity in the 1980s) and 1999-2019 (Rainfall erosivity in the 2010s).
The R-factor was calculated according to the approach proposed by Naipal [30].The formula for calculating the R-factor varies depending on the type of climate (according to Köppen).Data on the distribution of climate types by Köppen in the Volgograd region were taken from the article by Peel et al. [31].
For area with «Hot-summer humid continental climate (Dfa)», we used Formula (1) ln R = −1.99+ 0.737 × ln P + 2.033 × ln SDII, For area with «Warm-summer humid continental climate (Dfb)», we used Formula (2): Values of variables: ln R-the natural logarithm of the R-factor; ln P-the natural logarithm of the average annual precipitation; and ln SDII-the natural logarithm of the precipitation intensity index.SDII is calculated as the ratio of precipitation per year to the number of days per year with precipitation >1 mm; ln Z is the natural logarithm of the height of the location of the meteorological station.
The amount of snowmelt runoff was calculated by multiplying the precipitation that fell from December to March by a coefficient of 0.0591 [32].After the calculation, this value was added to the amount of stormwater runoff.

Soil Erodibility Factor
As is well known, the RUSLE model works poorly with soils with content of SOM > 4% [27,28,33].In our case, 65% of the research area had chernozem and kastanozem soils

Rainfall and Snowmelt Erosivity Factors
Data from the meteorological station of the Roshydromet network located in the Volgograd region were used to calculate the R-factor.The data were provided by RIHMI-WDC [29].Data were separated into two ranges: 1970-1998 (Rainfall erosivity in the 1980s) and 1999-2019 (Rainfall erosivity in the 2010s).
The R-factor was calculated according to the approach proposed by Naipal [30].The formula for calculating the R-factor varies depending on the type of climate (according to Köppen).Data on the distribution of climate types by Köppen in the Volgograd region were taken from the article by Peel et al. [31].
For area with «Hot-summer humid continental climate (Dfa)», we used Formula (1) ln R = −1.99+ 0.737 For area with «Warm-summer humid continental climate (Dfb)», we used Formula (2): Values of variables: ln R-the natural logarithm of the R-factor; ln P-the natural logarithm of the average annual precipitation; and ln SDII-the natural logarithm of the precipitation intensity index.SDII is calculated as the ratio of precipitation per year to the number of days per year with precipitation >1 mm; ln Z is the natural logarithm of the height of the location of the meteorological station.
The amount of snowmelt runoff was calculated by multiplying the precipitation that fell from December to March by a coefficient of 0.0591 [32].After the calculation, this value was added to the amount of stormwater runoff.

Soil Erodibility Factor
As is well known, the RUSLE model works poorly with soils with content of SOM > 4% [27,28,33].In our case, 65% of the research area had chernozem and kastanozem soils with a high SOM content.Because of this, we used a modern step-by-step approach for soil erodibility calculation presented in an article by Auerswald and co-authors [28].
Soil structure and soil permeability were defined according to Panagos and coauthors [34].The algorithm is presented in Tables 1 and 2. For the calculation of soil erodibility in the 2010s, we used open-source data.The list of sources is below: The Red Book of soils of the Volgograd region [25].
To calculate soil erodibility in the 1980s, we used the 'Soil map of Volgograd region' in scale 1:400,000 [37].We calculated the average K-factor based on the data used for preparing this map for each kind of soil type on this map [38].

Length Slope Factor
The LS-factor was estimated based on the Digital Elevation Model (DEM) with a 30 m resolution [39].This DEM is a data product of the Shuttle Radar Topography Mission (SRTM) [40].The DEM was processed according to the method proposed by L. Wang and H. Liu [41].This method efficiently identifies the topographic lows, overland flow paths, and watershed subdivision levels, which is convenient for the LS calculation.The calculation was performed using SAGA [42].The LS-factor of the area under study was mapped using ArcGIS 10.4.1 software.
After the processing, the LS-factor was calculated using the formula proposed by Desmet and Govers [43], Equation ( 11): where A ij-in is the contributing area at the inlet of the grid cell (i,j) measured in m 2 , D is the grid cell size (meters), X i,j = sin aaii,jj + cos aaii,jj and the aaii,jj is the aspect direction of the grid cell (i,j).m is related to the ratio β of the rill to interill erosion (Equation ( 12)): where β (Equation ( 13)): β= sinθ 0.0896 0.56 + 3 * (sin θ) 0.8 (13) The S-factor is calculated using the Formulas ( 14) and (15) proposed by Renard [27]:

Land Use and Crop Factors
The GLOBAL LAND 30 model with a spatial resolution of 30 m was used to determine agricultural land [26].Due to the lack of data on the contour of agricultural land during the USSR, we used the same model for data from the 1980s.
To calculate the C-factor, we took Landsat data.Landsat scenes were taken to the territory of the Volgograd region during the main vegetation period (from May to September).In order to avoid the impact of atmospheric phenomena on the quality of images, only images subjected to atmospheric correction (L2A level product) were used.The images were taken annually from 1984 to 1988 (5 years) and from 2017 to 2021 (5 years).For each year, the scenes were stitched into a single raster, and the value of the C-factor was calculated.Next, the average value of the C-factor was calculated for each of the 5-year periods.
The C-factor was calculated according to the Equation ( 16) suggested by Almagro et al. [44]: The values of the P-factor were assumed to be equal to 1.

Data Interpolation
In this paper, Empirical Bayesian Regression Kriging (EBRK) was used to set up local spatial models of the R-factor (1980s and 2010s) and K-factor (only 2010s).Two proxy variables, DEM and LS-factor, were used for the R-factor, and for the K-factor four proxy variables were used: LS-factor, DEM, soil types and a texture map of the Volgograd region.
EBRK is an interpolation technique that combines Ordinary Least Square regression and simple kriging to accurately predict moderately non-stationary data at a local scale [45].To improve the accuracy of ordinary kriging and account for uncertainty introduced in the variogram, EBRK estimates the semivariogram through subsetting and repeated simulations.This approach is particularly beneficial for small datasets and provides better accuracy than other kriging methods [46].EBRK also offers the advantage of reducing the smoothing effect of spatial interpolation, allowing for more precise predictions of moderately non-stationary data.Basically, the EBKR is decomposed as follows Equation ( 17): where µ (s) is the mean structure to model a large-scale spatial trend in the random field (RF), w(s) is a Gaussian process GP (0, C (s 1 , s 2 ; σ 2 , ϕ)) and ε(s) is Gaussian white noise N (0, τ 2 ).The mean structure µ (x) is explained by a set of covariates; µ (s) = x t (s) β, where x is a vector of covariates values at location s and β is the coefficients for the covariates.Denote the model parameters as θ = {β, σ 2 , ϕ, τ 2 }, with σ 2 representing spatial variance and ϕ controlling spatial decay effects.C (s 1 , s 2 ) is correlation function, which provides a way to characterize spatial dependence between pairs of locations (s 1 , s 2 ) [47].The interpolation's quality control was performed using cross-validation, which involves removing one data point to evaluate the interpolation model's fit.This method uses the remaining data points to interpolate at the removed point's location and then compares the interpolated value with the measured value to determine the interpolation's accuracy.This approach is advantageous compared to other cross-validation methods because it evaluates the data against the full distribution rather than only predicting a specific point.Statistics like the root mean standard error (RMSE) and average CRPS (continuously ranked probability score) were used to assess accuracy.The average CRPS measures the deviation from all data points and predicts cumulative distribution function, and a smaller value indicates better accuracy [48].The CRPS is decomposed as follows Equation (18): X is a variable.F is the cumulative distribution function (CDF) of X, such as F(y) = P[X ≤ y].x is the observation, and F is the CDF associated with the probabilistic forecast [49].
All calculations were performed with the ArcGIS PRO 2.9.2 program.The ArcMap 10.8 program was used to prepare the cartographic material.

R-Factor
The Figures show the distribution maps of the R-factor on the territory of the Volgograd region in the 1980s (Figure 7) and 2010s (Figure 8).The results of the R-factor for two range computations are presented in Table 3.The percentage distribution by classes of R-factor values is presented in Table 4.The results of the R-factor for two range computations are presented in Table 3.The percentage distribution by classes of R-factor values is presented in Table 4.The results of the R-factor for two range computations are presented in Table 3.The percentage distribution by classes of R-factor values is presented in Table 4.The average value of the R-factor in the territory decreased by 9% compared to the level of the 1980s (Table 5).The share of precipitation with a low erosion potential increased, and with a high one decreased (Table 6).This fact is evident in the southern part of the region: in the territories where R-factor values were previously recorded in the range of 60-100 Mj•mm/ha•h•g, today values with R-factor <60 Mj•mm/ha•h•g are recorded (Figures 2 and 3).In the north of the region, there was also a "displacement" of the erosion zone with an R-factor >190 Mj•mm/ha•h•g to a less dangerous zone with a value of 140-190 Mj•mm/ha•h•g (Figures 2 and 3, Table 6).However, the standard deviation of the R-factor increased from 64.70 in the 1980s to 74.71 in the 2010s, indicating more significant variability in the R-factor during the more recent period.
The RMSE and CRPS are measures of model accuracy, with lower values indicating a better fit to the data.The lower values in the 1980s suggest that the model performed better in that period than in the 2010s.

LS-Factor
The distribution map of the LS-factor in the Volgograd Region is shown in Figure 9. Statistical characteristics of LS-factor values on the territory of the Volgograd region, as well as the percentage distribution by class, are given in Tables 5 and 6.
The minimum and maximum values of the LS-factor were significantly lower or higher than average, but their total percentage value was the same (2.68% and 1.93%).
The average LS value of 0.52 is a measure of the central tendency of the dataset (66%), indicating that most points in the dataset have an LS-factor value around this value.However, as the standard deviation of 2.20 indicates, the dataset also has considerable variation.5 and 6.The minimum and maximum values of the LS-factor were significantly lower or higher than average, but their total percentage value was the same (2.68% and 1.93%).
The average LS value of 0.52 is a measure of the central tendency of the dataset (66%), indicating that most points in the dataset have an LS-factor value around this value.However, as the standard deviation of 2.20 indicates, the dataset also has considerable variation.7 and 8.   7 and 8.

K-Factor
The figures show the distribution maps of the K-factor on the territory of the Volgograd region in the 1980s (Figure 12) and 2010s (Figure 13).The statistical characteristics of the K-factor values and the percentage distribution by class are given in Tables 9 and 10

K-Factor
The figures show the distribution maps of the K-factor on the territory of the Volgograd region in the 1980s (Figure 12) and 2010s (Figure 13).The statistical characteristics of the K-factor values and the percentage distribution by class are given in Tables 9 and 10.            11 and 12.   11 and 12.The obtained data on soil loss by water erosion represent typical values for this region, averaging approximately 2 tons per hectare per year.This can be attributed to the continental climate of the area, characterized by limited precipitation, as well as the flat topography of the terrain [3,20].

Accuracy of the Results of the Calculation of Soil Erosion
To verify the results obtained, they were compared with the results of another independent research work.
We compared our results with the results from applying the soil erosion model developed under the guidance of G. A. Larionov [50].This model was used for different regions of the European part of the Russian Federation in the work of L.F. Litvin and co-authors [23] and for the assessment of soil loss by soil erosion by water in small river basins in the work of Maltsev and Yermolaev [20].In addition, the results were compared with those presented by Golosov with co-authors (Golosov et al., 2018 [19]).In this work, modified versions of the USLE and SHI erosion models were used for the evaluation mean annual erosion rate for 1986-2012.Only in the work of Litvin and co-authors were results for the regions presented which we could thoroughly compare with our data.The paper by Yermolaev and Maltsev provides data for the current level but by biomes.In the work of Golosov and co-authors, data are given for both the level of the 80s and 2012, but by natural zones.Therefore, we compared the results presented in the work of Maltsev and co-authors with the current value of the soil loss and the results of Golosov and co-authors to the value of the washout of the 1980s to the level of the 2010s.Results are presented in Table 13.According to Litvin's data [23] soil loss increased by 23.4% compared to the 1980s.We calculated the annual soil loss from the Litvin data for our area.The difference between the total area was because no USSR land use spatial data exist.Therefore, we used actual land use data for the USSR time (Methodology).The discrepancy with Litvin's results for 1980s level data was 12%, which is extremely small and indicates a good reproducibility and reliability of the results.
The discrepancy with Litvin's data was only 7.2%, demonstrating the high reproducibility and reliability of the data obtained.The discrepancy for the current level with Golosov and co-author's data was less than Litvin's and only 2.5%.But these data were not only produced for the steppe area in the Volgograd region, and represented the entire steppe area in the European plane.
Additionally, we compared our results with the results of local field erosion stations [51,52].The results are presented in Table 14.
Table 14.Comparison of the results with the results of field studies.

Source
Soil Loss, t ha −1 year −1 Soil Loss, t ha −1 year −1 (Kriuchkov, Makarov 2023) The Difference in the Results Obtained, % Babayan, 2014 [52] 12.3 10.5 15% Pynzar, 1980 [51] 0.51 0.48 6.3% As can be seen from the data obtained, the difference between the field measurements and our results was not critical and ranges from 6 to 15%.Also, a comparison was made between the Volgograd region's erosive soil loss map and European Russia's soil erosion map, created using the same calculation methods [20].According to the map, erosive soil losses for the Volgograd region were 0-5 t ha −1 year −1 .These values are not very different from ours, obtained here (2 t ha −1 year −1 ).Unfortunately, the authors did not published results by regions, only data according to natural zones.The Steppe zone is much broader than the Volgograd region and includes some parts of different areas.However, the results match if compared in the range between several zones presented in the Volgograd region.Semi-desert and steppe areas had mean 2 t ha −1 year −1 soil losses, and in our research we had 2 t ha −1 year −1 .
In addition, we compared our results with the soil erosion map of the European Union [21].A close country in landscape and climate conditions is Poland.Poland's mean soil loss in arable land is 1.61 t ha −1 year −1 .The difference between our data was 7%.

Discussion and Perspectives
The average R-factor in the Volgograd region has decreased in recent decades, and there has been more significant variability in the R-factor during the more recent period.
The R-factor depends directly on the amount of precipitation.Therefore, its decrease can be explained by a reduction in precipitation and the development of the aridification process.The results of other independent studies confirm the development of this process; for example, the work of Zolotokrylin and co-authors [53].The authors of this work analyzed the dynamics of the Standardized Precipitation Index (SPI) for several southern regions of Russia, including the Volgograd region.It was shown that from 1970 to 2000, a wetter humidification period was observed in the region's territory.Since 2001, it has been replaced by more arid conditions.
Another paper that has confirmed that fact is the article by Eltsov and co-authors [54].In this paper, based on the analysis of climatic and meteorological data for the last 60 years, the authors demonstrated an increase in the aridity of the climate of the Lower Volga region and a shift in the isolines of the aridity index values (according to De Morton-IDM) of the region by about 100 km to the northwest, which also correlates with our data.
The high values of the LS-factor in the Volgograd region were associated with elevated forms of relief and, accordingly, with a higher steepness of the slopes.The main erosionhazardous areas are the slopes of the right bank of the Volga and Don rivers.A similar pattern was observed with low values of the LS-factor, which were recorded in the region's territories with low relief forms (the territory of the left banks of the above rivers).The highest values of the LS-factor were observed in the territory of districts, a significant part of the area of which are hills: the southeastern part of the Central Russian Upland (Alekseevsky (No. 34) and Kumylzhensky (No. 18) municipal districts) and the southern part of the Volga Upland (Kamyshinsky (No. 23) and Danilovsky (No. 29) municipal districts) (Figure 4).The lowest value of the LS-factor was recorded in territories with predominantly low-lying terrain-Sredneakhtubinsky (No. 11), Pallasovsky (No. 7), and Nikolaevsky (No. 15) municipal districts.Approximately 8% of the region's total agricultural land is dangerous for agriculture (LS-factor > 1).
The values of the LS-factor calculated for the Volgograd Region were in good agreement with the results of similar studies performed for EU countries with predominantly flat territory-Poland, Hungary, Estonia, Netherlands, Latvia, Lithuania, and Denmark [55].The average values of the LS-factor for the territories of these countries were 0.52, 0.59, 0.32, 0.19, 0.39, 0.35, and 0.32, respectively.
The area of contours with a high protective potential has increased, while that with a low one has fallen (Table 10).The indicator of the C-factor in the region's territory decreased by approximately 6% compared to the level of the 1980s (Table 9).That means that the protective capacity of vegetation in the region increased by 6%.Most likely, this is due to the overgrowth of agricultural land with natural vegetation, which increased the protective degree of vegetation cover.Natural steppe vegetation has a higher erosion resistance than field crops [50].It is worth noting that all this is happening against the background of a decrease in the amount of incoming moisture, both from natural sources (Section 3.1) and artificial ones.According to the data from Gorokhova and co-authors and Pankova and Novikova [56,57] the amount of irrigated land in the Volgograd region has decreased relative to the level of the 1980s.The value of the C-factor in the region increased from north to south and from west to east, which is associated with a decrease in the total biomass of vegetation and an increase in the sparsity of vegetation cover [50].Low values of the C-factor corresponded to areas with a high value of the R-factor, and vice versa.This is because territories with a high R-factor value correspond to territories with a large amount of precipitation.Furthermore, moisture is one of the limiting factors in steppe biomes [57].
Looking at the data, we can see that the range of K-factor values increased from 0.156 to 0.973 in the 1980s to 0.182 to 0.994 in the 2010s, indicating a wider variety of soil erodibility in the latter period.The minimum K-factor value also increased from 0.156 to 0.182, suggesting that even the least erodible soils have become slightly more susceptible to erosion.The maximum K-factor value increased significantly from 0.973 to 0.994, indicating that the most erodible soils have become even more susceptible to erosion.This trend is also reflected in the average K-factor value; the average value of the K-factor increased by 19%, and the erosion resistance of the soil decreased.
The K-factor's standard deviation (SD) decreased slightly from 0.129 to 0.088, suggesting that the variability of soil erodibility was reduced in the 2010s compared to the 1980s.The RMSE and CRPS, which are measures of the accuracy of a model's predictions, were only available for the 2010s data, with values of 0.098 and 0.050, respectively.These values indicate that the model used to estimate the K-factor in the 2010s is reasonably accurate in predicting soil erodibility.RMSE and average CRPS for the data of the 1980s were not determined because they were calculated using a direct method based on soil map data.
According to the government data, 52% of agricultural land's territory is subject to soil organic matter losses [58].This is also due to the deterioration of the soil's basic physical and chemical properties, which affect high erosion resistance, for example a decrease in organic matter content and a change in the soil texture to a lighter one.It is logical to assume that the main factor in reducing erosion resistance is the loss of soil organic matter.
The reasons for soil organic matter losses can be varied.Climate change can change the ratio of humification and mineralization processes in soil erosion.Climate change leads to changes in the ratio erosion level.On the one hand, rainfall erosivity is decreasing.On the other hand, the level of wind erosion should have been increasing.As is known, chernozems soils lost 11-25 t/ha over the last decades in the Volgograd region, or 0.3-0.8%Agronomy 2023, 13, 2527 20 of 23 organic matter loss through soil erosion [59].Currently, there are no studies to assess the magnitude of wind erosion in the research area, so it is difficult to determine the wind and soil erosion by the water ratio in this area.
The number of soil losses in the territory of the Volgograd region had increased by 16.2% compared to 1980.That demonstrates the trend of an increase in the erosion load on the region's territory.It is especially noticeable through the change in the structure of the distribution of the amount of soil loss: the area of contours with minimal values of soil erosion by water level (up to 1 t/ha/year) decreased, and the areas of all other shapes increased.The fact that the increase in the amount of soil erosion by water occurred against the background of two compensating phenomena-a decrease in the erosion potential of snowmelt and rainfall runoff and an increase in the protective potential of vegetation, because of changes in land use-is exciting.These processes should lead to a decrease in the territory's soil erosion by water level.The main factor that led to this is the deterioration of the anti-erosion ability of the soil as a result of losses of SOM.About 53% of the Volgograd region's agricultural lands are subject to SOM losses [58].Such a wide spread of this process has led to an increase in the value of the K-factor by 19% and, as a result, an increase in the amount of soil loss.
Soil erosion by water in the region's territory decreased from west to east and north to south.That is due to a decrease in the values of the leading actors of erosion processesthe R-factor and LS-factor.In the very northwest of the region, soil erosion by water level is also decreasing.This is due to the widespread chernozem soils, which have high erosion resistance.

Conclusions
This study provides quantitative and spatial data on soil erosion intensity across time in agricultural lands on a regional scale for the first time.Previous studies in this field have been restricted to specific regions in Russia or conducted on a small scale, resulting in high spatial generalization.Our research is the first attempt at the high-resolution evaluation of soil erosion by water across time.On the one hand, we suggest a new algorithm for assessing soil erosion by water levels, with good reliable results.On the other hand, certain limitations related to the quality and level of detail of the original data and unresolved methodological objectives are apparent: (1) The problem with soil data.The problem with open-source soil data is familiar in Russia and has been discussed many times in scientific papers [6,20].For K-factor calculations in 2010s, we used open-source data already published in books and our own data.However, more is needed to cover the entire region.For the K-factor in the 1980s, we used a soil types map of the Volgograd region, which had a scale of 1:400 000 [37].The simplification of the soil property's spatial distribution, and consequently soil erosion, was significantly facilitated by this model.(2) The problem with the USSR's actual areas of agricultural land.We used a present land use model for both time ranges.However, the extent of farmlands in the USSR was different from the extent nowadays.Furthermore, we need to find out how current areas match with the 1980s extent.
In order to advance research on soil erosion, obtaining more accurate assessments is crucial.To clarify the direction of future investigations, we propose the following areas of study.
(1) Updating the rainfall erosivity (R-factor) equation for research area.In our study, we used a regression equation of the magnitude of the R-factor from the type of climate, obtained mainly for the territory of the United States, and tested in other territories.
It is logical to assume that despite similar natural areas, it can work poorly in other territories.We need to collect data from all pluviographs on the territory of Russia and calculate similar equations for that area.(2) Calculation K-factor for 1980s based on USSR data.In our study, we assessed the K-factor based on average properties for each soil type based on the soil map of the Volgograd region.This method is ineffective because a spatial variation in the soil parameters needs to be considered.It can improve our results significantly.To do this, we must collect data from all agricultural land surveys conducted from 1980 to 1991.These data exist, but it will take time to collect, digitalize, and prepare the dataset because it is on paper sources.In our assessment, the dataset's total points were approximately 10,000.(3) Define actual areas agricultural lands in the USSR time period.For defining actual USSR agricultural lands, we can use remote sensing data.The Landsat mission started in the 1970s, which can give us enough high-resolution data to define areas we need.Modern Russian land-use models have data from 2000 to the present day.Nevertheless, the areas of agricultural lands and fields were changed after the USSR collapsed.Some farms still exist, but many were closed.(4) Estimating wind erosion.The Volgograd region is a region where water and wind erosion exist together.The level of wind erosion increases from northwest to southeast, and soil erosion by water vice versa.Nowadays, we do not have actual data about the level of wind erosion in that region.Researchers highlight the problem of modeling soil erosion by water when wind erosion also exists in modeling.We need to set the reason for soil losses in our area precisely.For a more precise evaluation, we should also know the value of wind erosion.We also have plans for modeling wind erosion levels based on international or national models.
Within the climatic changes observed in the last four decades, including a decrease in rainfall and snowmelt runoff, and farmlands being overgrown by natural steppe vegetation, there has been a notable increase in soil losses within the Volgograd region.On average, the intensity of soil losses from the 1980s to 2010s has increased by approximately 16%, increasing from 1.48 to 1.72 t ha −1 year −1 .This has happened because the soil organic matter's massive losses are likely caused by wind erosion and soil erosion by water.Our modeling has produced similar results to other studies.
It will be possible to refine the cartographic model for predicting soil erosion in the future by incorporating more precise spatial data on soil erodibility (K-factor) and land use factor (C-factor).This will help in better understanding the reason for the decreasing soil organic matter in soils.Furthermore, wind erosion can be evaluated.

( 4 )
Modelling of soil erosion.(5) Checking the quality of results with other research.

( 1 )
Model calculations were performed as follows: (2) Collecting data for calculating factors of soil erosion.(3) Choosing optimal model parameters for performing calculations.(4) Modelling of soil erosion.(5) Checking the quality of results with other research.

Figure 6 .
Figure 6.The design of the experiment.

Figure 6 .
Figure 6.The design of the experiment.

Figure 7 .
Figure 7. R-factor spatial distribution in the 1980s.

Figure 8 .
Figure 8. R-factor spatial distribution in the 2010s.

Figure 8 .
Figure 8. R-factor spatial distribution in the 2010s.

Figure 8 .
Figure 8. R-factor spatial distribution in the 2010s.

Figure 9 .
Figure 9. LS-factor spatial distribution.Statistical characteristics of LS-factor values on the territory of the Volgograd region, as well as the percentage distribution by class, are given in Tables5 and 6.

Figures 10
Figures 10 and 11 show maps of the distribution of the C-factor in the Volgograd region.

Figures 10 and 11 13 Figure 10 .
Figures 10 and 11 show maps of the distribution of the C-factor in the Volgograd region.Agronomy 2023, 13, x FOR PEER REVIEW 13

Figure 10 .
Figure 10.C-factor spatial distribution in the 1980s.Figure 10.C-factor spatial distribution in the 1980s.

Figure 10 .
Figure 10.C-factor spatial distribution in the 1980s.

Figure 11 .
Figure 11.C-factor spatial distribution in the 2010s.Statistical characteristics of the values of the C-factor in the Volgograd region, as well as the percentage distribution by class, are given in Tables7 and 8.

Figure 11 .
Figure 11.C-factor spatial distribution in the 2010s.Statistical characteristics of the values of the C-factor in the Volgograd region, as well as the percentage distribution by class, are given in Tables7 and 8.

Figures 14 and 15
Figures 14 and 15 show maps of the distribution of soil erosion levels in the territory of the Volgograd region.Figure 16 shows a map of the difference in soil loss between the 1980s and 2010s.Agronomy 2023, 13, x FOR PEER REVIEW 16 of 24

Figure 14 .
Figure 14.Soil loss by soil erosion by water: spatial distribution in the 1980s.Figure 14.Soil loss by soil erosion by water: spatial distribution in the 1980s.

Figure 14 .
Figure 14.Soil loss by soil erosion by water: spatial distribution in the 1980s.Figure 14.Soil loss by soil erosion by water: spatial distribution in the 1980s.

Figure 14 .
Figure 14.Soil loss by soil erosion by water: spatial distribution in the 1980s.

Figure 15 .
Figure 15.Soil loss by soil erosion by water: spatial distribution in the 2010s.Figure 15.Soil loss by soil erosion by water: spatial distribution in the 2010s.

Figure 15 . 24 Figure 16 .
Figure 15.Soil loss by soil erosion by water: spatial distribution in the 2010s.Figure 15.Soil loss by soil erosion by water: spatial distribution in the 2010s.Agronomy 2023, 13, x FOR PEER REVIEW 17 of 24

Figure 16 .
Figure 16.The difference in soil loss between the 1980s and 2010s.

Table 1 .
Definition of soil structure index.

Table 2 .
Definition of soil permeability index.

Table 3 .
Statistical characteristics of the R-factor in (Mj mm)/(ha h g).

Table 4 .
Percentage distribution by classes of R-factor values in %.

Table 5 .
Statistical characteristics of the LS-factor.

Table 6 .
Percentage distribution by classes of LS-factor values.

Table 5 .
Statistical characteristics of the LS-factor.

Table 6 .
Percentage distribution by classes of LS-factor values.

Table 7 .
Statistical characteristics of the C-factor.

Table 7 .
Statistical characteristics of the C-factor.

Table 8 .
Percentage distribution by classes of C-factor values. .

Table 9 .
Statistical characteristics of the K-factor.

Table 9 .
Statistical characteristics of the K-factor.

Table 10 .
Percentage distribution by classes of K-factor values.

Table 11 .
Statistical characteristics of soil losses by soil erosion by water.

Table 11 .
Statistical characteristics of soil losses by soil erosion by water.

Table 12 .
Percentage distribution by classes of soil losses by soil erosion by water.

Table 13 .
Comparison of the results obtained with other studies.
* From the 1980s to 2010s.