1. Introduction
The Xilingol League of Inner Mongolia, located in the middle of the farming–pastoral ecotone in northern China, is an important part of the “Two Barriers and Three Belts”, and has been given the main function of “a windbreak and sand-fixing belt of northern China”. This functional positioning is not only related to the sustainable development of local grassland animal husbandry and farmland planting but is also of great significance for ensuring ecological safety and air quality of northern China and even the entire East Asia region. Due to global climate change, intensifying human land reclamation and grazing activities, the region has experienced increasingly serious ecological and environmental problems, such as grassland degradation, soil desertification, and frequent sandstorms, since the 1980s. Soil erosion with wind erosion are important causes of land degradation, desertification, and dusty weather in this area [
1,
2]. Since 2000, the Chinese government has carried out ecological management and restoration projects, namely, the Returning Grazing Land to Grassland project, the North Shelter Forest, and the Natural Forest Protection project, achieving major results. The quality of the local and downwind ecological environment has been significantly improved [
3,
4]. Therefore, accurately assessing the soil erosion situation in the farming–pastoral ecotone in northern China, and quantitatively measuring the dynamic changes in soil erosion moduli, are of great significance to the restoration of degraded land and the construction of ecological civilization.
At present, the global research on soil erosion has changed from qualitative to quantitative. Relevant research methods are becoming more and more abundant; for example, remote sensing combined with USLE and RUSLE model methods, to achieve soil erosion moduli measurements using high-tech, although these are still subject to the limitations of remote sensing image resolution and scale, and thus field sampling data are required to verify the model results [
5,
6,
7]. The rare earth element tracing method (REE) is used to measure the concentration of REE and the amount of erosion or deposition of the slope section after a certain erosion time by releasing rare earth elements, which is more used for indoor simulation and small-scale field research [
7]. Although the runoff plot method can calculate the law of soil and water loss on a slope, the cost of time, funds, and energy required are relatively high [
8]. The geomorphological research method mainly judges the degree of soil erosion qualitatively by observing the geomorphological phenomena related to soil erosion in the field [
7]. Although the artificial rainfall simulation method can save manpower, material resources, and time costs, it is difficult to apply the artificial rainfall simulation device in large-scale field experiments [
5,
9]. The isotope-tracing method is mainly used to collect soil samples at the site to obtain information on the erosion and deposition status of soil particles in the soil, and then estimate the soil erosion moduli. The choice of research methods varies according to the research scale (over a large spatial scale or small-scale research), the erosion properties (wind erosion, water erosion or farming erosion, etc.), and applicability.
The application of the isotope-tracing method to soil erosion research is relatively mature in theory and technology. Menzel [
10] pointed out that some radioisotopes (such as
137Cs,
210Pb
ex,
7Be, etc.) are closely combined with surface fine particles after settling from the atmosphere, and their content in the soil would change accordingly with the transport process of soil materials; therefore, these radioisotopes can be used as ideal tracers for soil erosion or accumulation processes. For
137Cs isotope-tracing techniques, scientists have conducted more in-depth research and established many quantitative soil erosion models (the empirical model, the proportion model, mass balance model, profile distribution model, diffusion and migration model, etc.), which have been applied in studies related to hydraulic erosion, tillage erosion, and wind erosion worldwide [
11,
12,
13,
14,
15,
16]. In China, researchers have conducted many soil erosion studies based on
137Cs tracer techniques in the Loess Plateau [
17], the Yangtze River Basin [
18,
19], the Northeast Plain [
20,
21], the Tibet Plateau [
22], and the southern hilly areas [
23]. Most of these studies have concentrated on soil erosion in agricultural cultivation regions caused by hydraulic erosion, whereas in the farming–pastoral ecotone of northern China, which is dominated by wind erosion, researchers have conducted relatively few studies on soil erosion changes caused by natural processes as well as ecological restoration processes [
24,
25,
26].
On the other hand, since
137Cs is derived from atmospheric nuclear blasts, its characteristics, including a single source, a lack of supplies, and a short half-life (30.17 a), make
137Cs increasingly difficult to detect in soil samples. Furthermore, various spatial distribution of tracer isotopes in sandy soils lead to greater uncertainty in the measurement results [
25]. Therefore, researchers have been trying to develop alternative isotope technology to
137Cs.
210Pb
ex soil tracing technology or combined
137Cs–
210Pb
ex tracing technology is an important technology that reflects this development trend. There have been some successful cases of soil tracing studies using
210Pb
ex worldwide [
27,
28,
29,
30,
31,
32,
33]. For example, Walling and He [
33] applied
210Pb
ex to measure the erosion moduli of cultivated land in the UK, confirming the potential of using
210Pb
ex to estimate hydraulic erosion on cultivated land. Benmansour [
27] applied
137Cs and
210Pb
ex in a Mediterranean agricultural area in Morocco, demonstrating that soil erosion processes have not changed significantly over the last 100 years. In general, the number of
210Pb
ex soil erosion study cases are much less than that based on
137Cs, and its application is not as extensive as the latter, indicating that there are still many technical details worth exploring, such as the
210Pb
ex reference inventory calculation and the soil erosion measurement models improvement.
Researchers have applied a single-isotope (
137Cs) tracing technique combined with a multi-isotope (
137Cs–
210Pb
ex) tracing technique to quantify soil erosion in the farming–pastoral ecotone of northern China, particularly, in a typical temperate grassland with a wind-sand area, Inner Mongolia, China. Hu et al. [
34] were the first to establish high-precision
137Cs vertical distribution profiles in grassland and cultivated land in the farming–pastoral ecotone of northern China (Taipusi Banner, Inner Mongolia). Furthermore, Hu measured the soil erosion moduli and pointed out that the application of the
137Cs soil tracing technique for quantifying soil erosion has important advantages due to the vast land area and invading wind erosion on the Mongolian plateau. Based on the
137Cs reference inventory estimation model developed by Walling and He [
35], Liu et al. [
26] established the
137Cs reference inventory for different areas in the southern and northern parts of the Tariat–Xilingol transect on the Mongolian plateau, and further noted that the wind erosion moduli at two typical grassland sample sites in the southern part of the transect (Xilinhot and Zhengxiangbai, Inner Mongolia) were nearly three times that of the typical grassland sample sites in the northern part of the transect (Bayannur, Mongolia), which may be caused by the differences in population density and livestock carrying capacity. However, directly using the theoretical model to estimate
137Cs reference inventory without any modification may lead to uncertainty in the results. Qi et al. [
36] compared the theoretical model-estimated reference inventory with the measured reference inventory at 66 sample sites in China and found that the theoretical model underestimated the actual reference inventory. Hu et al. [
24] further showed that the actual reference inventory was 120–150% of the theoretical model-estimated values in the central-eastern region of Inner Mongolia. This result provides a basis for identifying regional
137Cs reference points and discerning a reasonable
137Cs reference inventory interval.
The deposition process of
137Cs particles mainly occurred in the 1950s and 1970s when the world’s major nuclear powers conducted atmospheric nuclear tests, while the deposition of
210Pb
ex series particles have been occurring for a long time in nature, potentially being continuously replenished; thus, the erosion rates measured by
137Cs and
210Pb
ex as soil tracers reflect the soil erosion process on a short time scale (since the 1970s; that is, a 50-year scale) and long-term scale (a scale of more than 100 years). Therefore, the composite tracer of the above two isotopes can clearly analyze the dynamic characteristics of soil erosion. Hu and Zhang [
25] also developed a combined
137Cs–
210Pb
ex tracing method and estimated soil erosion moduli in the southern margin of Hunsandake, Inner Mongolia. The results showed that wind and sand activity in the southern margin of the Hunsandake sandy land has significantly weakened over the past 50 years. In addition, based on the information of the
137Cs and
210Pb
ex reference inventory provided by the established literature, Wang et al. [
5] carried out an analysis of the effectiveness of soil conservation by regional ecological management and construction. The results concluded that soil erosion moduli in cultivated and fallow lands were significantly reduced after the implementation of ecological projects, even converting some lands from being in a soil erosion state to a soil-accumulation state.
From a comprehensive analysis of the global and northern China farming–pastoral ecotone, it can be concluded that the application of the 137Cs–210Pbex combined tracing method for quantitative soil erosion assessment has become the frontier of soil isotope-tracing research. Using the characteristics of the different decay cycles of the two tracer isotopes, to measure and compare changes in soil erosion moduli on different periods, is crucial to analyze the evolution of the characteristics of soil erosion processes, especially to evaluate the effectiveness of ecological restoration and governance projects. In this paper, two soil tracer isotopes, 137Cs and 210Pbex, were applied to quantify soil erosion moduli and conduct comparative analysis in the farming–pastoral ecotone of northern China. Our specific objectives were to answer the following questions: (1) How do one select background plots and determine a reference inventory based on multiple characteristics? (2) How does soil erosion or accumulation vary on different time scales? (3) What are the main factors contributing to the uncertainty of the study results?
2. Study Area and Method
2.1. Study Area
Xilingol League is located at the southeastern margin of the Mongolian Plateau and the east-central part of the Inner Mongolia Autonomous Region, China. The sampling plots were selected from Taipusi Banner and Zhengxiangbai Banner (
Figure 1), which are adjacent to Zhenglan Banner to the east, bordering Xianghuang Banner and Huade County to the west, Guyuan County to the south, and Sunitezuo Banner and Abaga Banner to the north. The study region lies between 114.2° E–115.9° E and 41.5° N–43.1° N.
The study area has a temperate continental monsoon climate. Springs are windy and dry, summers are warm and short, autumns are cool, and winters are cold and long. The annual average precipitation is about 343 mm. From south to north, the dryness of the climate gradually increases, and the soil type gradually changes from grassland chestnut-calcium soil to grassland wind-sand soil. The southern part is mainly covered by high plain grassland and dry cultivated land, while the northern part mainly presents fixed, semi-fixed, semi-shifting sandy and dune landform types. The study area also has a small amount of forest, water bodies, and other surface cover types scattered throughout.
2.2. Sampling and Testing
We collected 14 soil samples from 21 to 23 October 2018, with 7 sample plots each in Taipusi Banner and Zhengxiangbai Banner (
Figure 1). There were 5 sample plots and 9 sample plots in cultivated land and grassland, respectively. On the other hand, 2 sample plots were in semi-shifting sandy land and 3 sample plots were in fixed sandy land.
The sample plots are generally set in flat and open areas, and the specific locations are determined according to different topographical conditions (semi-shifting sandy land, shifting sandy land, semi-fixed sandy land, and fixed sandy land, etc.) and surface cover (high coverage grassland, medium coverage grassland, low coverage grassland, cultivated land, and fallow land, etc.).
After clearing the ground vegetation and litter, 4 groups of soil samples were collected at each plot according to the triangular distribution method: 3 full-profile samples were excavated at vertexes of the triangle by a soil driller with an internal diameter of 35 mm and a maximum sampling depth of 30 cm; a suite of layer samples were collected at the center via the stratified stripping method. For the layer samples, a layer thickness of 3 cm was used for the top 0–15 cm column, a layer thickness of 5 cm was used for the 15–20 cm column, and a thickness of 10 cm was used for the lowermost 20–30 cm column. The depth was set at 30 cm for the whole sample collection, and the mean values of the 3 sets of full-profile sample data at each plot were calculated statistically in the subsequent data processing.
A low-background photon (γ-ray) multi-channel (214) energy spectrometer from EURISYS, France, was used for soil isotope activity detection, equipped with high purity germanium (HPGe) well probe, with a germanium crystal activity volume of 114 cm3, a diameter of 5.5 cm, and height of 6.0 cm. The effective detection depth was 4.1 cm, the inner diameter was 1.5 cm, and the detection γ-ray energy covered 15–2084 keV. A total of 8–10 g of the dried sample was packed into a plastic test tube with an inner diameter of 1.5 cm before testing. After sealing for 20 days, we measured it using the HPGe probe for more than 24 h. The minimum detection limit activity was 0.5 dpm, at the 99% confidence level. According to the national standard of the γ-ray spectrum analysis for soil radionuclides (GB 08304), the activities of 210Pbex and 137Cs were obtained by counting the γ-rays with energies of 46.5 keV and 662 keV (137Ba–m) in the energy spectrum, respectively.
The activity information obtained from the instrument measurements is the activity value per unit mass (Bq∙kg
−1). Combined with soil bulk density and sampling depth, this can be converted to an area-based isotopic activity value (Bq∙m
−2). For the layer samples, the area-based cumulative isotopic activity of
137Cs and
210Pb
ex can be calculated from the following formula:
where
(
137Cs point inventory) or
(
210Pb
ex point inventory) is the total activity of
137Cs and
210Pb
ex at the sample point (Bq∙m
−2);
is the layer number of the sampling layer;
is the number of sampling layers;
is the isotope activity of
137Cs and
210Pb
ex at layer
(Bq∙kg
−1);
is the bulk density of layer
(t∙m
−3); and
is the sampling depth of layer
(m).
For the full-profile samples, the total activity of
137Cs and
210Pb
ex per unit area can be calculated by the following equation:
where
is the isotope activity of
137Cs and
210Pb
ex in the whole sample (Bq∙kg
−1);
is the total weight of each soil sample after sieving (kg); and
is the cross-sectional area of the soil driller (m
2).
2.3. Reference Inventory Determining
The reference inventory (CRI: 137Cs Reference Inventory; PRI: 210Pbex Reference Inventory), also known as background values, refers to the activities of 137Cs and 210Pbex in the sample without any soil erosion, accumulation/deposition, or any kind of manual disturbances. It is crucial to measure reliable regional CRI or PRI to estimate the soil erosion modulus, to analyze the soil erosion dynamics.
In an earlier-proposed method, a small number of plots are selected through field investigations and on-site observations, and then the mean CPI (or PPI) of them is regarded as the regional CRI (or PRI). Considering the relationship model between
90Sr atmospheric deposition and precipitation, Walling analyzed the relationship between CRI and precipitation in different spots around the world and proposed a model to estimate CRI globally. However, Qi et al. [
36] systematically compared the estimated values of this model with the measured values in China and concluded that the model normally underestimates the actual CRI. Hu et al. [
24] further pointed out that the measured CRI in Inner Mongolia was about 120–150% of the model estimates. The above studies provide an important basis for discerning the accuracy and reasonableness of the measured CRI; they also provide support for modifying the model estimates and determining the regional CRI in the absence of accurate and reliable background points.
In this study, we consciously selected high coverage grasslands as alternative background points based on actual land-cover/land-use information. Then we compared the difference between the CPI values from the alternative background points and the model-estimated CRI value to determine whether the CPI values fall within the potential CRI interval (120–150% of the model values). Finally, we also investigated the 137Cs distribution characteristics in the soil profile. Combining the information from the above three aspects, the background sample plots and regional 137Cs reference inventory can be determined. Furthermore, we also took the 137Cs background sample plots as the 210Pbex background sample plots and determined the corresponding 210Pbex reference inventory.
2.4. Erosion/Accumulation Moduli Measurement
Soil erosion or accumulation moduli, also known as the soil erosion or accumulation rate, is mainly quantified by soil bulk density and annual average soil erosion thickness. Due to the difference in the self-decay rate between 137Cs and 210Pbex radioisotopes, based on previous research results, this paper uses a profile distribution model (Formulas 3 and 4) and mass balance model (Formula 8) to calculate the soil erosion or accumulation thickness at different time scales at each sampling point, and further calculate the soil erosion or accumulation moduli (soil erosion or accumulation rate).
The calculation of the
137Cs soil erosion modulus in grassland was performed using a profile distribution model [
16,
34,
37]. The model is as follows:
where
is the modified
137Cs reference inventory (Bq∙m
−2), which is the result of performing the wind and snow distribution influence (
k = 0.95) on the
CRI;
is the
137Cs activity at the current sampling plot (Bq∙m
−2);
is the sampling year, which is 2018 in this study;
is the morphological parameter of the
137Cs distribution curve, with depth in the background plot profile (Qi et al., 2008), and which can be obtained from the layer-by-layer
137Cs activity values and fitted by the least-square method; and
is the annual average soil erosion thickness (cm∙a
−1) at the sampling plots since 1963 (the peak year of
137Cs deposition).
After obtaining the average annual soil erosion thickness, the soil erosion modulus can be calculated by the following equation:
where
is the soil erosion modulus (t∙km
−2∙a
−1); and
is the soil bulk density (t∙m
−3).
Currently, several models apply
210Pb
ex to measure soil erosion moduli [
38,
39,
40]. Considering the continuous sedimentation process and the decay and loss processes of
210Pb
ex, Zhang et al. [
41] derived a mass balance model by measuring the steady distribution state of
210Pb
ex in uncultivated soils. Sun et al. [
42] simplified the model and established a new formula to depict the relationship between the soil erosion moduli and
210Pb
ex content in uncultivated soil. The model form is as follows:
where
is the
210Pb
ex reference inventory (Bq∙m
−2);
is the total area activity of
210Pb
ex in the current sampling plot (Bq∙m
−2);
is the
210Pb
ex sedimentation flux (Bq∙m
−2∙a
−1);
is the
210Pb
ex decay coefficient (0.031∙a
−1);
is the relaxation depth (kg∙m
−2); and
is the multi-year average erosion mass thickness (kg∙m
−2∙a
−1).
From the above two equations, the multi-year average erosion mass thickness (kg∙m
−2∙a
−1) of
210Pb
ex is derived as follows:
Therefore, applying a simple inversion formula, the soil erosion modulus can be roughly calculated for two different periods—1920s–1970s and 1970s–present—by following Equation (9):
where
is the soil erosion modulus at the 100-year scale, which is measured by the
210Pb
ex tracer method;
is the soil erosion modulus at the 50-year scale (1970s–present), which is measured by the
137Cs tracer method; and
is the soil erosion modulus for the 1920s–1970s period.
3. Results
3.1. 137Cs Distribution Characteristics and Reference Inventory
Based on the CRI estimation software, the theoretical value of CRI in the study area is 1404 Bq∙m
−2. Further considering the amplification coefficient of the measured CRI in the Inner Mongolia Plateau relative to the model simulation value, we determined that the reasonable range of CRI values in the study area is 1684~2167 Bq∙m
−2 (
Table 1).
Among all 14 sample plots, only the CPI value of sample plot ZXB-3 was within the above value range. Further investigation of the
137Cs distribution characteristics at this plot showed that the
137Cs activity of sample plot ZXB-3 was the highest in the surface soil (3–6 cm) and decreased regularly with increasing depth (
Figure 2). These values exhibited the typical distribution type, characterized by “a single peak + a negative exponential curve”, indicating ZXB-3 is a typical
137Cs background sample plot. Therefore, by combining four types of information, including soil type (chestnut calcium soil), land cover (grassland),
137Cs activity value, and
137Cs distribution characteristics, sample plot ZXB-3 can be identified as a background sample plot, and the CPI (1927.9 Bq∙m
−2) of sample plot ZXB-3 can be seen as the CRI.
In addition to the ZXB-3 sample plot, we also investigated the CPI and 137Cs distribution characteristics of the remaining sample plots.
The analysis showed that, for sample plots ZXB-1, ZXB-2, ZXB-4, ZXB-5, ZXB-6, ZXB-7, TPS-1, and TPS-3, the measured CPI values were lower than the inferred CRI interval from the theoretical model. Although the distribution of the 137Cs profiles showed a “negative exponential” or “single peak + a negative exponential curve” pattern, the maximum occurrence depth of 137Cs was very shallow (less than 20 cm), which was much shallower than the normal occurrence depth of 25–30 cm. Therefore, by combining land cover characteristics, total CPI, and 137Cs distribution characteristics, the above eight sample plots can be regarded as eroded sandy land/grassland plots.
For sample plots TPS-2 and TPS-4, the measured CPI values (991.9 Bq∙m−2 and 304.5 Bq∙m−2) were significantly smaller than the inferred CRI value interval. In deeper soils, below 21 cm or 15 cm, the 137Cs activity was extremely small, while the 137Cs activity of each layer showed a uniform distribution in the upper soils. Therefore, according to the land cover/use characteristics, total CPI, and 137Cs profile distribution characteristics, it can be determined that sample plots TPS-2 and TPS-4 were regarded as a plot that had inherited the characteristics of a former cropland.
The measured CPI values in the TPS-5, TPS-6, and TPS-7 sample plots (735.6 Bq∙m−2, 417.1 Bq∙m−2, and 629.7 Bq∙m−2, respectively) were significantly smaller than the inferred CRI value interval, and the profile distribution pattern of the 137Cs specific activity was complex and differed from the common distribution pattern of a “negative index”, or “single peak + negative exponential curve”, or “uniform distribution of each layer”. Therefore, by combining total the CPI and 137Cs profile distribution, it is assumed that the above three sample plots may had experienced more complex soil erosion and soil accumulation processes due to anthropogenic activities, but the overall process is dominated by soil erosion.
Overall, the distribution patterns of the 137Cs profiles at the sampling plots were closely related to their land-use characteristics. The 137Cs profile distribution at the sample plots located in grassland and sandy grassland (ZXB-1, ZXB-2, ZXB-4, ZXB-5, ZXB-6, ZXB-7, TPS-1, and TPS-3) usually showed a “single peak + negative exponential curve” profile distribution pattern. The peak 137Cs activity was generally found in the surface or subsurface layer (0–6 cm) of the soil, and the 137Cs specific activity decreased gradually with increasing soil depth. In the sample plots located on cultivated land (TPS-2 and TPS-4), 137Cs showed a uniform distribution pattern in the plow layer (0–20 cm), and the 137Cs specific activity did not change significantly with soil depth. At the sample plots located in fallow land (TPS-5, TPS-6 and TPS-7), the 137Cs specific activity was often influenced by a combination of anthropogenic and natural factors, and its profile distribution pattern was more complex, without any uniform and obvious pattern.
3.2. 210Pbex Distribution Characteristics and Reference Inventory
We investigated the profile distribution patterns of
210Pb
ex in the sample plots (
Figure 3). In general, the maximum depth of
210Pb
ex in soil was generally higher than that of
137Cs in soil.
The results showed that 210Pbex was still present in the soil layers below 30 cm in grassland (ZXB-3), semi-shifting sandy land (ZXB-4), as well as in cultivated land (TPS-2, TPS-4) and fallow land (TPS-5, TPS-7).
In the grassland (ZXB-1, ZXB-3, TPS-1, TPS-3), fixed, or semi-fixed sandy land (ZXB-2, ZXB-4, ZXB-5, ZXB-6, ZXB-7), 210Pbex basically showed a complete “negative exponential curve” distribution. The peak activity of 210Pbex was usually found at the 0–3 cm depth of the profile, and the specific activity of 210Pbex decreased gradually with increasing depth of the profile. Unlike the “single peak + negative exponential curve” pattern of 137Cs in grassland, the profile distribution of 210Pbex only had a “negative exponential curve” distribution pattern.
In the two sample plots (TPS-2 and TPS-4) located on cultivated land, 210Pbex showed approximately equal activity from the surface layer to 20 cm depth, while, with increasing soil depth, 210Pbex showed less activity, even close to 0 in soils deeper than 20 cm. This reflects the “mixing” process caused by plowing.
In the three sample plots (TPS-5, TPS-6, and TPS-7) located in fallow land, the distribution of 210Pbex was still homogenized by plowing in the deeper soil layer (9–30 cm), but in the upper soil layer (0–9 cm), the spatial distribution of 210Pbex clearly showed a “negative exponential” distribution pattern.
In the previous section, sample plot ZXB-3 was identified as a background plot based on 137Cs activity, 137Cs profile distribution, and land-cover information. Therefore, the PPI of sample plot ZXB-3 (10,041.2 Bq∙m−2) is naturally regarded as the PRI.
3.3. Soil Erosion Modulus
Based on the qualitative analysis in
Table 2 and
Figure 2, the CPI of all sample points (except reference point ZXB-3) was less than the CRI value, indicating different intensities of soil erosion in the above sample plots. According to the
137Cs soil erosion modulus measurement model (Equations (3)–(5)), the average annual erosion thickness at each plot ranged from 0.1 to 0.7 mm∙a
−1 on a 50-year scale, and the erosion modulus ranged from 139.9 to 1029.5 t∙km
−2∙a
−1 (
Table 2).
Specifically, sample plot ZXB-6, located in fixed sandy land with high vegetation cover, had the highest 137Cs area activity (1353.4 Bq∙m−2) and the lowest annual average erosion modulus (139.9 t∙km−2∙a−1). In contrast, sample plot ZXB-2, located in a semi-shifting sandy land with a wind-sand soil type, had the lowest 137Cs area activity (238.3 Bq∙m−2) and the highest annual average erosion modulus (1029.5 t∙km−2∙a−1).∙The soil erosion modulus in this area is closely related to land cover, land-use history, and vegetation coverage.
It can also be seen from
Table 2 that the sample plots of the semi-shifting sandy land and fixed sandy land (ZXB-2, ZXB-4, ZXB-5, ZXB-6, ZXB-7) usually had the highest soil erosion modulus, with an average soil erosion modulus of 605.5 t∙km
−2∙a
−1. The second was the samples (TPS-2, TPS-4, TPS-5, TPS-6, TPS-7) from the cultivated and fallow land, which are more influenced by anthropogenic disturbance, showing an average soil erosion modulus of 576.1 t∙km
−2∙a
−1. Finally, the lowest soil erosion modulus occurred in the medium- and high-cover grassland areas (ZXB-1, TPS-1, TPS-3), where the surface vegetation can resist wind and sand erosion to a certain extent, with an average value of 379.5 t∙km
−2∙a
−1.
Based on the same principle, the PPI of all sample points (except the background point ZXB-3) was smaller than the PRI (
Table 3), indicating different intensities of soil erosion in the above sample points at the century scale. According to the
210Pb
ex soil erosion modulus measurement model (Equations (6)–(8)), the average annual erosion thickness at each plot ranged from 0.02 to 1.9 mm∙a
−1 and the erosion modulus ranged from 34.9 to 2637.9 t∙km
−2∙a
−1 at the century scale (
Table 3).
Among them, the sample plot (TPS-4) on the cultivated land with chestnut soil and lower vegetation coverage had the lowest 210Pbex area activity (3396.3 Bq∙m−2) and the highest average annual erosion modulus (2637.9 t∙km−2∙a−1). In contrast, the sample plot (TPS-1) with medium vegetation coverage had the highest 210Pbex area activity (9776.6 Bq∙m−2) and the lowest average annual erosion modulus (34.9 t∙km−2∙a−1).
In contrast to the soil erosion modulus regulations at the 50-year scale obtained from the 137Cs measurements, the sample plots located on the cultivated and fallow lands (TPS-2, TPS-4, TPS-5, TPS-6, TPS-7) had the largest average soil erosion modulus (1182.5 Bq∙m−2) at the past 100-year scale, indicating that complex human land reclamation activities lead to more intense soil erosion processes. Furthermore, the sample plots located on semi-shifting sandy land and fixed sandy land (ZXB-2, ZXB-4, ZXB-5, ZXB-6, ZXB-7) had an average soil erosion modulus of 523.3 Bq∙m−2. Similar to the above findings, the sample plots located in medium and high grassland coverage areas (ZXB-1, TPS-1, TPS-3) had the lowest average soil erosion modulus (454 Bq∙m−2).
3.4. Changes in Soil Erosion and Accumulation Modulus
According to Equation (9), after obtaining the average soil erosion modulus at the 100-year scale (i.e., since the 1920s) based on the
210Pb
ex soil tracer, and the soil erosion modulus at the 50-year scale (i.e., since the 1970s) based on the
137Cs soil tracer, the soil erosion modulus from the 1920s–1970s in the 20th century can be further inferred (
Figure 4). According to the changes in the soil erosion modulus, all 13 erosion sample plots were divided into three categories.
There are six plots (ZXB-1, ZXB-4, ZXB-6, TPS-2, TPS-4, and TPS-5) in the first category. In the past 100 years, there has been a continuous erosion trend, but the erosion intensity has decreased in the past 50 years. The soil erosion modulus since the 1970s has been 40–92% that of the previous period. In terms of spatial distribution, the above sample plots are mainly distributed in the central grassland area and southern cultivated area of the study area.
There are four plots (ZXB-5, ZXB-7, TPS-6 and TPS-7) in the second category. They showed continuous erosion over the past 100 years, but the intensity of erosion has intensified over the past 50 years, mainly in the cultivated land of the southern area and the semi-shifting sandy land of the northern area. In particular, sample plots ZXB-5 and ZXB-7 changed from a state of basically no erosion or no accumulation to a mildly eroded state. However, in general, the erosion amount of these two plots was not high.
There are three plots (ZXB-2, TPS-1, and TPS-3) in the third category, where the soil process is converted from soil accumulation to erosion, and the erosion modulus was higher than that of the accumulation process. In terms of spatial distribution, the above sample plots are mainly distributed in the semi-shifting sandy areas in the northern region, cultivated land, and grassland of the southern area. This indicates that these sample plots have undergone a more complex soil erosion or accumulation process.
In summary, the southern part of the study area (cultivated land and grassland) was generally dominated by erosion processes. With the 1970s as the boundary, the soil erosion modulus in the former period is smaller than that in the latter period, which indicates that land reclamation and grassland grazing practices in the southern area may have intensified the soil erosion process. In the northern part of the study area (sandy grassland dominated by fixed and semi-fixed sandy land), both soil accumulation and soil erosion were present. The soil erosion or accumulation modulus after 1970 was generally lower than the previous period, indicating that wind and sand activities in the northern area had weakened.
4. Discussion
Accurate determination of regional
137Cs and
210Pb
ex background values and a reference inventory are the scientific basis for the quantitative assessment of soil erosion by applying
137Cs and
210Pb
ex composite tracer techniques. However, there are uncertainties in this process, which involves the following aspects. The study area, which is in the farming–pastoral ecotone in northern China, experienced complicated land-cover changes and human activities, such as cultivation–abandonment and grassland grazing–abandonment alternate states, leading to uncertainties in the backgrounds of the sample plots. Frequent dune consolidation activities lead to the formation of either blown-out or wind-sand soil accumulation. In the wind-sand area, the distribution of tracer isotopes vary greatly in plane and profile. In some sites, the isotope cannot be detected at all, while in other sites the depth of the isotope distribution exceeds 30 cm. The slope may also affect the isotopic content of the soil surface; for example, the
137Cs content in the upper or steep slope will be lower than that in the middle and lower slopes [
43]. The geomorphology of the study area in this paper mainly includes a high plain grassland or sand dune geomorphology, and the gradient fluctuation is relatively small. In addition,
137Cs in soils become increasingly difficult to detect over time due to the single source and a short half-life, and
210Pb
ex comes through multiple decay processes of
226Ra,
222Rn, and
238U, resulting in the complex measurement process of
210Pb
ex [
44]. It should be noted that it is difficult to separate the two isotope particles after they are combined with fine particles. The activity of
137Cs is mainly related to the local rainfall, and the activity of
210Pb
ex is mainly affected by the local geological structure and rock stratum type. Both particles are mainly distributed in the soil layer above 30 cm in the soil profile, and their content in the soil is very low. Therefore, due to the above characteristics, isotope particles such as
137Cs and
210Pb
ex will not affect the ecological environment of the soil, nor will they threaten the safety of humans and animals.
In this study, various characteristics, such as the land cover, theoretical simulation values, and profile distribution pattern of
137Cs, were considered to determine ZXB-3 as the background plot and the corresponding CRI value (1927.9 Bq∙m
−2). Compared with the previous results in adjacent areas (CRI values for different years were uniformly modified to 2018: 1740.1 Bq∙m
−2 in Miyun of Beijing, 1690.3 Bq∙m
−2 in Dehui of Jilin, 1712.8 Bq∙m
−2 in Jiutai of Jilin, 1947 Bq∙m
−2 in Keshan of Heilongjiang, 1938.7 Bq∙m
−2 in Hunsandake of Inner Mongolia, 1994.8 Bq∙m
−2 in Xingan of Inner Mongolia, 1930.4 Bq∙m
−2 in Xilingol of Inner Mongolia, and 1980.9 Bq∙m
−2 in Tongliao and Chifeng of Inner Mongolia) [
24,
25,
43,
45,
46], the CRI value obtained in this study is within the ranges of geographically adjacent CRI values and has considerable credibility. In this paper, the
137Cs background plot (ZXB-3) was also used as the
210Pb
ex background plot to obtain the PRI value (10,041.2 Bq∙m
−2), which is higher compared to the results from adjacent areas (6600 Bq∙m
−2 in Keshan of Heilongjiang, 5730 Bq∙m
−2 in the Loess Plateau, and 8112 Bq∙m
−2 in Hunsandake, Inner Mongolia) [
25,
41]. In fact, considering the fact that traces were detected in the soil layer below 30 cm in the ZXB-3 site (
Figure 3), the PRI in this study may be higher. The reason may be that the continuous sedimentation of
210Pb
ex resulted in more accumulation of
210Pb
ex.
The soil erosion modulus ranged from 139.9 t∙km
−2∙a
−1 to 1029.5 t∙km
−2∙a
−1 for grasslands and sandy grasslands in the study area based on
137Cs isotope tracing. It is higher than the previous studies conducted in Tariat–Bayannur–Xilingol of the Mongolian Plateau (64.58 t∙km
−2∙a
−1 to 419.6 t∙km
−2∙a
−1) [
26,
47]. This is related to the fact that the study area is located at the southern margin of the Hunsandake sandy area, where there are many fixed and semi-fixed sandy lands. Compared with the results of the abandoned land (Harahalin: 6723.06 t∙km
−2∙a
−1) in the northern part of the Tariat–Bayannur–Xilingol sample zone [
26,
47], the soil erosion modulus on the cultivated and fallow land in this study was lower (348.2–890.2 t∙km
−2∙a
−1) due to weaker winds and the presence of conservation tillage and shelterbelts in the study area.
In this study, changes in the soil erosion modulus were estimated for the past 100- and 50-year scales. In the southern part of the study area (cultivated land and grassland), the land reclamation and grassland grazing activities promoted the rate of soil erosion, while in the northern part of the study area (marginal area of Hunsandake sandy land), there is a slowing down trend in the soil erosion modulus due to climate change (weaker wind and more precipitation) [
48]. This is consistent with previous findings, showing there is a slowdown in the soil erosion modulus and even a shift from erosion to accumulation processes at most sample plots in the Hunsandake sandy land [
25]. These new quantitative findings provide new evidence for understanding the impacts of climate change, physical-geographical changes, and human activities in the farming–pastoral ecotone in northern China.
Considering that the Chinese government has carried out large-scale ecological restoration and treatment projects in the region over the past 20 years (since 2000) (e.g., the Beijing–Tianjin Sand Source Control project, Three-North Shelter Forest program, Returning Grazing Land, and the Natural Forest Protection project), the measurement of the soil erosion characteristics’ changes in the two periods before and after 2000 is an important basis for evaluating the effectiveness of ecosystem construction projects. However, current soil erosion measurement techniques based on
137Cs and
210Pb
ex cannot support this goal. In fact, after determining the background plots, we can qualitatively and roughly analyze soil erosion and accumulation processes in recent years based on the isotope distribution characteristics in soil profiles. A more reliable method is to develop soil tracing and erosion rate measurement techniques based on
7Be because of the characteristics of
7Be as a natural radionuclide with a shorter half-life (53.3d) [
49].