SPYR (N32°02′–36°13′, E95°42′–102°17′) comprises an area of about 130,798 km², and the area records for 35% of the total runoff of the Yellow River (Zheng et al., 2018). The elevation of the SPYR ranges from 2,052 to 6,227 m and descends from the Southwest to the Northeast (Figure 1). Topographic features and geographical situations pose an impact on the natural environment, with characteristics such as frost temperatures, intense diurnal temperatures, frequent seasonal rainfall, recurrent strong winds, and powerful solar radiation (Ge et al., 2017). Most of the study area has a typical continental plateau climate, with an average annual temperature ranging −4 and 2°C rising from Western to Eastern. Average annual rainfall ranges from 350 to 750 mm and is concentrated from June to September. The soil types are mainly alpine meadow soil and steppe soil. Alpine grassland accounts for 95.12% of the total area of the SPYR, and it includes alpine meadow at 83.12% and alpine steppe at 12% (Liang et al., 2016).

Among the 33 environmental variables, the number of variables selected by the GA, RFE, BOR, UF, and SA algorithms was 29, 26, 22, 22, and 13 variables, respectively (Appendix Table 2). The five algorithms selected precipitation, wind velocity, rainfall erosivity, frozen earth period, slope, aspect, water flow direction, the NDVI, land cover, and vegetation type as variables in common that may play significant roles in explaining the soil erosion mechanism.

Figure 2 compared five variable search methods across four modeling techniques (GLM, cubist, RFM, and BST) in terms of both prediction accuracy and model complexity. The cross-validation results of GA BST model that selected 29 variables showed slightly inferior performance to the UF cubist model. This indicated that the 22 variables selected by the UF cubist model contained better predictive power to that of 29 variable selected by the GA BST model. Comparing with the UF cubist model, the SA cubist model had fairly comparable performance in similar model accuracy, but dramatically decreased model complexity. The SA cubist model struck a balance between model complexity and performance. It greatly reduced the number of variables to only 13 while preserving most of the predicting power to infer soil erosion. The cross-validation identified and optimized model parameters. The 13 variables selected by the SA cubist model captured the major factors – soil organic carbon, soil property and soil erodibility besides nine common variables mentioned above involving climate (e.g., precipitation, wind velocity, and rainfall erosivity), soil (e.g., frozen earth period, soil organic carbon, soil property, and soil erodibility), topography (e.g., slope, aspect, and water flow direction), vegetation (e.g., vegetation type and NDVI) and management (e.g., land cover) properties (Figure 2 and Appendix Table 2).

The goodness of fit for the RUSLE models and machine learning algorithms combined with variable selection models showed that the GLM (R = 0.63) and cubist model (R = 0.67) performed the best while RUSLE (R = 0.49) performed the worst (Figure 3). The Taylor diagram show that (1) the fitting result was poorer for the RUSLE model than the 24 machine learning models (the location labeled with ‘RUSLE’ were the outermost from the observation point). (2) Compared with the other five variable selection models, the SA variable set combined with the GLM and cubist model was more consistent with the actual measurements, and the SA variable set combined with BST was the least consistent with the actual measurements. (3) Compared with the BST and RF models, the GLM and cubist models were slightly more consistent with the field measurements (the points labeled with ‘GLM’ and ‘cubist’ are nearer to the observation point than the points labeled with ‘BST’ and ‘RFM’). The cross combination of the BST model and all variable sets had the lowest agreement (the points labeled with ‘BST’ is the farthest to the observation point than the other points). (4) Contrast to the AR, BOR, GA, UF, and RFE variables, the GLM did not present a significant difference in fitting results (the gaps of the points labeled with ‘RFM’ from the observation point were nearly the same). The results show that in terms of classification accuracy, the RFE performs best among the five feature-selection algorithms, and the Cubist model performs best among the four machine learning algorithms. In these 24 combinations, the optimal model is the SA cubist model (R = 0.67, RMSD = 1368).

The SA cubist model estimated the total erosion loss was about 8.45 10⁸ t⋅a^–1 in 2014, which is 12.64 times the allowable soil erosion loss. The soil erosion modulus showed substantial spatial heterogeneity in the SPYR (Figure 4). The soil erosion modulus ranged from 4841.68 to 8054.09 t⋅km^–2⋅a^–1, and the average soil erosion modulus was 6460.95 t⋅km^–2⋅a^–1. In general, 99% of the area in the SPYR belongs to the severe soil erosion modulus (5,000–8,000 t⋅km^–2⋅a^–1).

According to the Standards for Classification and Gradation of Soil Erosion, the region with the most severe soil erosion modulus (7,000–8,000 t⋅km^–2⋅a^–1) occupied 12.38% of the total soil erosion regions. It mainly appeared in the Northern part of the SPYR. The moderately severe region of the soil erosion modulus (6,000–7,000 t⋅km^–2⋅a^–1) occupied 71.19% of the total soil erosion regions and was distributed in the middle of the study area. The slightly severe region of the soil erosion modulus (5,000–6,000 t⋅km^–2⋅a^–1) occupied 16.43% of the total soil erosion regions and was distributed South of the SPYR (Figure 4).

The simulated soil erosion modulus values were more significant than that in two previous report (276 and 3,208 t⋅km^–2⋅a^–1), which used the RUSLE model (Teng et al., 2018) and RWEQ model (Teng et al., 2021) to calculate the soil water and wind erosion modulus on the Qinghai-Tibetan Plateau, respectively. At the same time, the RUSLE model emphasizes soil erosion loss by rainfall, slope, aspect, vegetation, and soil erodibility but ignores soil erosion loss from winds, freeze-thaw cycles, and human activity. Although few research studies have provided a comprehensive prediction of soil erosion (Sadeghi et al., 2018; Shen et al., 2018), high levels of precipitation have been reported to lead to severe erosion, while low temperatures combined with adequate precipitation may cause even greater soil losses in the study area (Yang et al., 2018). The same as the RWEQ model.

Under freeze-thaw conditions, soil moisture is transported, and the soil structure is damaged; moreover, the soil porosity, bulk density, shear strength, aggregate stability, and organic matter are all changed, which leads to high soil erodibility and erosion intensity increases (Wang L. et al., 2014; Sadeghi et al., 2018). In addition, soil erosion modulus caused by wind erosion is also an important aspect of soil erosion because the wind speed is a vital factor that affects soil erosion (Teng et al., 2021). Soil particles erode as the sheer pressure exerted by wind-force leads to soil particles being hard to grip on the soil surface (Jiang et al., 2018; Zhang et al., 2018). Furthermore, grazing livestock and human activity are both key factors in the intensification of soil erosion. Intensive grazing will alter and destroy soil surfaces, causing whole sod to erode and expose mineral soil, especially in areas with low vegetation cover (Lin et al., 2008; Evans et al., 2017; Zhang et al., 2018).

Besides, frequent human activities will cause soil erosion, such as collecting Chinese Caterpillar Fungus (Wang C. G. et al., 2018) and overexploitation of tourist attractions, and so on (Yang Y. et al., 2016). The comprehensive set of variables based on the inclusion of freeze-thaw erosion, wind erosion, water erosion, and human-caused erosion was more objective than that used in the RUSLE or RWEQ model. The inclusion of these variables may explain why the soil erosion results derived from this set of variables were more excellent than the results from the RUSLE or RWEQ model because this comprehensive set of variables involved more soil erosion factors, including meteorological, soil, topographic, vegetation, and management factors, which enhance our understanding of soil erosion processes using a data-based model.

Evaluation results of model accuracy comparisons among five variable selection methods and four prediction models are illustrated in Figure 3. The model validation included machine learning models and the RUSLE model. The results showed that variable screening technology could further reduce model complexity. The SA method selected 13 out of the 33 relevant variables, i.e., climatic factors, including precipitation, wind velocity, and rainfall; soil factors, including SOC, soil properties, soil erodibility, and frozen earth period; topography factors, including slope, aspect, and water flow direction; and vegetation factors, including vegetation type and the NDVI. The UF method largely reduces the number of variables that were most able to explicate soil erosion. The Boruta algorithm screening out the majority of climatic and topographic variables, thus abated model collinearity. The feature-selection method can simplify model parameters and reduce operational costs; in addition, this method can also improve the estimation accuracy of a model by eliminating redundant interference parameters.

Our results showed that in terms of classification accuracy, the RFE performs the best among the five feature-selection algorithms and the Cubist model performs the best among the four machine learning algorithms. The projected results of the RUSLE model were limited by fixed parameters, which included only rainfall, slope, slope direction, vegetation, and soil texture. The accuracy of the RUSLE model is the lowest among the models because some factors that may affect soil erosion were ignored, such as diurnal temperature, strong winds, freeze-thaw, and so on.

In order to gain much more insight through statistical analysis and modeling compared to biased linear sampling along roads, the conditional Latin Hypercube design for soil sampling to even its spatial distribution would be adopted in our near future studies (Ma et al., 2020; Tang et al., 2021).

According to the Standards for Classification and Gradation of Soil Erosion, 99% of the SPYR area experiences severe soil erosion. (Ministry of Water Resources of the People’s Republic of China, 2007), soil erosion modulus was further reclassified into three grades: level I (severe, soil erosion modulus in 5,000–6,000 t⋅km^–2⋅a^–1), level II (moderately severe, the soil erosion modulus in 6,000–7,000 t⋅km^–2⋅a^–1), and level III (most severe, the soil erosion modulus 7,000–8,000 t⋅km^–2⋅a^–1). The soil erosion modulus of each county was reclassified into three grades (Figure 5). Level III erosion occurred mainly in the triangle zone of Qumalai County and Maduo County. Level II erosion occurred in Xinghai County, Dari County, Chengdu County, Western Zeku County, Western Maqin County, Tongde County, and parts of Jiuzhi County, Henan County, Gande County, and Yushu County. Level I erosion primarily occurred in Eastern Maqin County, Eastern Zeku County, and most areas of Jiuzhi County, Henan County, Gande County, and Yushu County (Figure 5).

Severe soil erosion causes significant ecological degradation as well as substantial ecological, economic losses in the SPYR. For example, the economic monetization of SOC losses caused by soil erosion total $47.90 billion in the SPYR. The economic losses were approximately fivefold of the output value of animal husbandry ($10.10 billion) in 2014 (Qinghai Bureau Statistics, 2014). We can infer that the economic monetization of soil erosion would be far huger than the SOC loss. Actually natural resources are priceless. Monetization of soil erosion loss is only a caution light for taking a partial solution in eco-compensation. Therefore, a new ecological policy recommendation based on soil erosion reclassification should be proposed to balance grassland conservation and animal husbandry production for the prevention and control of soil erosion loss (Figure 6).

A conservative management strategy was suggested for some fragile areas (soil erosion at level III) since grassland’s economic products have led to serious soil erosion (Figure 6). This grassland sector should not be used for grazing (Yin et al., 2019). At this level, the government should provide ecological subsidies to compensate for the economic losses caused by not allowing grazing. Additionally, the government should relocate local herder households to new countryside neighbors to cities. At the same time management department should also supply training projects to sustain transitions to alternative livelihoods for these herdsmen, such as making Tibetan typical handicrafts and building Tibetan distinctive hotels for developing Tourism. Notably, over the long term, simply banning grazing may result in a snowballing population increase of deleterious wild animal species that are short of plenty natural predators, such as rabbits and pika, thus probably resulting in new ecological problems (Wang Y. et al., 2018). Therefore, after 5 years of banning grazing, grasslands should be managed under a grass-animal balance policy, or grazing should be continued to be banned. The implemented policy will depend on the condition of the grassland.

The moderately severe sector (level II) should be managed intensively to improve the basic conditions (Figure 6). Even now the balance of livestock and pasture production is an important factor at level II soil erosion regions (Yin et al., 2019). The study conducted by Fan et al. (2010) propose that the determination of appropriate carrying capacity should follow the principle of the “law of the minimum limiting factors” in order to minimize degradation of the grasslands, i.e., stocking rates should be primarily based on herbage yields of low-productivity years, there is also the opportunity to conserve forage in high-yield years and carry this over to alleviate grazing pressure in low-yield years. Further, there is the possibility of reducing grazing pressure on more degraded parts of the grassland in high-yield years to aid their recovery. For sustainable development, it is insufficient to manage grazing systems solely within the framework of pasture productivity (Lin et al., 2011) especially for distinctive erosion areas (level II). The ‘optimum’ stocking rate should be obtained from grazing experiments not only considering the forage-livestock balance, but also maintaining the soil erosion at a tolerable level. Sun et al. (2015) conducted a 3-year study on Tibetan sheep on the Eastern Qinghai-Tibetan Plateau. They found that the ‘optimum’ stocking rate that ensures economic sustainability is 24 sheep months ha^–1 (Sun et al., 2015). Considering moderately severe soil erosion conditions in this area, the stocking rate should not exceed one-third of the ‘optimum’ stocking rate in order to prevent and control the soil erosion loss at level II soil erosion regions. If herder households comply with this policy, then they should be able to receive a subsidy.

The severe sector (level I) should be managed intensively to improve grassland productivity via the application of certain measures (fertilization, tillage, and resowing) for the sustainable utilization of grassland (Figure 6). When soil erosion is at the level I, the stocking amount should not exceed two-thirds of the ‘optimum’ stocking rate in order to obtain a certain production income, meanwhile prevent and control the soil erosion. In addition, the stocking rate in the cold-period pastures should be less than the stocking rate in the warm-season pastures (Sun et al., 2015).