Evaluation of Qinghai-Tibet Plateau Wind Erosion Prevention Service Based on RWEQ Model

Abstract

Ecosystem service research is essential to identify the contribution of the ecosystem to human welfare. As an important ecological barrier zone, the Qinghai-Tibet Plateau (QTP) supports the use of a crucial wind erosion prevention service (WEPS) to improve the ecological environment quality. This study simulated the spatiotemporal patterns of the WEPS based on the Revised Wind Erosion Equation (RWEQ) and its driving factors. From 2000 to 2015, the total WEPS provided in the QTP ranged from 1.75 × 10⁹ kg to 2.52 × 10⁹ kg, showing an increasing and then decreasing trend. The average WEPS service per unit area was between 0.72 kg m⁻² and 1.06 kg m⁻². The high-value areas were concentrated in the northwest and north of the QTP, and the total WEPS in different areas varied significantly from year to year. The average retention rate of the WEPS in the QTP was estimated to be 57.24–62.10%, and high-value areas were mainly located in the southeast of the QTP. The total monetary value of the WEPS in the QTP was calculated to be between 223.56 × 10⁹ CNY and 321.73 × 10⁹ CNY, and the average WEPS per unit area was between 0.08 CNY m⁻² and 0.13 CNY m⁻², showing a declining–rising–declining trend. The high-value areas gradually expanded to the west and east of the QTP. The slope was the most important factor controlling the spatial differentiation of the WEPS, followed by the landform type, average annual precipitation, and average annual wind speed, and human activities such as land-use change could improve the WEPS by returning farmland to grassland and desertification control in the QTP.

1. Introduction

The Qinghai-Tibet Plateau (QTP) is the most unique geographical unit in the world because of its high altitude and mountainous landforms. The QTP is considered an important ecological barrier in south-western China, as it supports various ecosystem services in the region and in neighboring countries [1,2]. Among these ecosystem services, wind erosion prevention service (WEPS) is very important because of the extreme climate and low vegetation coverage in most parts of the plateau [3,4,5,6,7]. The average wind erosion rate for the plateau as a whole reached the medium erosion standard in 2001 [8], and the area of desertification in the QTP accounted for 15.1% of the total region in 2015 [9]. Therefore, it is critical to assess the regional WEPS and analyze the influencing factors in order to forecast wind erosion and provide a scientific basis for desertification control [10].

The WEPS can be estimated by the difference between potential and actual sand erosion by ecosystems. Several models have been adopted to simulate sand erosion by ecosystems, including the sediment transport equation [11], the wind erosion equation (WEQ) [12], the Texas erosion analysis model (TEAM) [13], the description model [14], the revised wind erosion equation (RWEQ) [15], the wind erosion prediction system [16], and the wind erosion stochastic simulator (WESS) [17]. Among these models, the RWEQ model is most commonly used for sand erosion simulations in the QTP. For example, Jiang [18] used the RWEQ model to evaluated wind erosion modules in Qinghai and showed that the total area of wind erosion in Qinghai Province accounts for more than half of the regional area, and the level of wind erosion damage in the Qaidam Basin is relatively serious. Huang [19] analyzed the WEPS of ecosystems in the Tibet Plateau from 1990 to 2010 by RWEQ, showing that the annual average soil wind erosion modulus, the average sand fixation service quantity, and the retention rate of the average annual WEPS in the Tibet Plateau were 1.58 kg m², 18.99 × 10⁸ t, and 66.5%, respectively. Teng [10] adopted the RWEQ model to simulate the WEPS of the QTP from 2000 to 2015. However, studies on WEPS in this area have been limited, with most studies only involving local areas [20,21]. Accordingly, few studies have comprehensively described the WEPS and its dynamic characteristics over the whole QTP. As the terrain of the QTP is complex, the spatial pattern of sand erosion and WEPS may differ significantly between the east and west parts [9]. Few studies have examined wind erosion throughout the QTP or investigated the changes among different areas and ecosystems.

In recent years, many scholars have used multi-source data and methods such as the linear correlation model, multiple regression model, structural equation model, and constraint effect to investigate the factors (vegetation coverage, wind speed, precipitation, animal husbandry development, land use, etc.) driving the WEPS [10,22,23,24,25,26,27]. Latocha et al. showed that changes in land use have affected soil erosion to a larger extent than climate change in the Sudetes Mts. within the last 150 years [28]. Garbrecht et al. showed that warmer temperatures and reduced rainfall could exacerbate soil wind erosion in farmland areas in the southern Great Plains [29]. Li argued that the increase in the afforestation area is the main factor driving the improvement of WEPS in Inner Mongolia [30]. Sharratt showed that increase in the crop yield in the Colombian Plateau could reduce damage from soil wind erosion [31]. Few studies have conducted a driving factor analysis of the WEPS in the QTP. Teng [10] simulated the influences of climate changes and human activities on sand erosion and the WEPS in the QTP from 2000 to 2015 through comparative analysis of spatiotemporal dynamics. The above studies have helped us to understand the driving mechanism responsible for the spatial distribution pattern of WEPS, but there are still gaps in our knowledge on the interaction of different driving factors.

The geographic detector is a spatial statistical method based on spatial differentiation that is used to reveal the driving mechanism of a variable [32]. The geographic detector can analyze the interaction between two independent variables and determine the specific effect type while detecting the influences of variables with attributes of different types and values on the spatial distribution of the dependent variable. At present, the geographic detector is widely used in the analysis of the change process of spatial patterns in many fields such as meteorology, environmental pollution, ecology, human health, and regional planning [33,34,35,36,37,38,39]. Few studies have adopted the geographic detector to analyze the driving factors related to the WEPS. Zhang [40] conducted an exploratory study on the factors driving changes in the WEPS in Xilingol League using the geographic detector, and proved that the application of the geo-detector analysis on the WEPS is rational and scientific. The application of geographic detectors to analyze the factors driving changes in the spatial pattern of WEPS can account for the influences of type variables and numerical variables on service changes, analyze the interactions between various factors, and make up for the lack of analysis of factors driving the spatial distribution of the WEPS.

Therefore, this study adopted the RWEQ model to simulate the WEPS in the QTP from 2000 to 2015, estimated its monetary value using the restoration cost method, evaluated the differences in the material and monetary value of the WEPS among different cities and ecosystems, and analyzed the contributions and interactions of natural and human factors in the formation and evolution of the spatial pattern of the WEPS. This study can be used to improve the understanding of the spatial patterns of the WEPS in the QTP and provides crucial knowledge that can be used to make decisions about ecological restoration and protection to improve the WEPS in different ecosystems.

2. Methodology and Data

2.1. Study Area

The QTP is located in southwestern China, and involves five provinces (autonomous regions): Qinghai, Tibet, Sichuan, Xinjiang, Gansu, and Yunnan (Figure 1). In the western part of the QTP is the Pamir Plateau and the Karakoram Mountains, bordering Tajikistan and Pakistan, and in the southern part there are the Himalayas, bordering Nepal, Bhutan, and other countries. The QTP has an average elevation of more than 4000 m. The highest altitude is found on Mount Everest, 8844.43 m, whereas the Jinsha River is only 1503 m above sea level. The special geographical location of the QTP means that it has a unique plateau climate with a cold, hypoxic, and arid climate with strong radiation. The temperature is 15–20 °C lower than that of other areas at the same latitude. The complex terrain in this area also causes the climate to vary. The temperature decreases from southeast to northwest, and the average annual temperature can reach below −6 °C. Due to the blocking effect of the multiple high mountains, the annual precipitation decreases from south to north, and the annual precipitation is less than 50 mm. The total annual solar radiation on the QTP can be as high as 5400–7900 MJ·m⁻²·year⁻¹, and the total annual sunshine hours are 2500–3200 h.

2.2. Method and Data

2.2.1. Calculation of WEPS

The WEPS represents the reduction of wind erosion caused by vegetation and is calculated as the difference between the potential wind erosion under bare soil conditions and the actual wind erosion under vegetation conditions. The potential and actual wind erosion are calculated with the RWEQ model [15,41]. The actual wind erosion is the level of soil erosion under the condition of vegetation coverage in a real scenario, whereas the potential wind erosion is the level of soil erosion under bare land conditions without vegetation. The WEPS was calculated as follows:

$$SL=\frac{2z}{s^2}\times Q_{max}\times e^{-{\left(\frac{z}{s}\right)}^2}$$

(1)

$$Q_{max}=109.8\times \left(WF\times EF\times SCF\times K^{\prime }\times C\right)$$

(2)

$$s=109.8\times {\left(WF\times EF\times SCF\times K^{\prime }\times C\right)}^{-0.3711}$$

(3)

$$SLR=\frac{2z}{s_r{}^2}\times Q_{rmax}\times {\mathrm{e}}^{-{\left(\frac{z}{s_r}\right)}^2}$$

(4)

$$Q_{rmax}=109.8\times \left(WF\times EF\times SCF\times K^{\prime }\right)$$

(5)

$$s_r=109.8\times {\left(WF\times EF\times SCF\times K^{\prime }\right)}^{-0.3711}$$

(6)

$$G=SLR-SL$$

(7)

where SLR represents the potential wind erosion (kg m⁻²); Q_rmax represents the maximum potential sediment transport capacity (kg m⁻¹); S_r represents the potential critical plot length (m); SL represents the actual wind erosion (kg m⁻²); Q_max represents the maximum sediment transport capacity (kg m⁻¹); S represents the critical plot length (m); G represents the WEPS (kg m⁻²); z represents the calculated downwind distance (m); and WF, EF, SCF, K′, and C correspond to the factors of weather (kg m⁻¹), the soil erodibility fraction (%), soil crust factor (dimensionless), the soil roughness factor (dimensionless), and vegetation (dimensionless), respectively. The WF was calculated as follows:

$$WF=Wf\times \frac{\rho }{g}\times SW\times SD$$

(8)

$$Wf=u_2\times {\left(u_2-u_1\right)}^2\times N_d$$

(9)

where Wf is the wind factor (m³ s⁻¹); g represents the gravitational acceleration (m² s⁻¹); ρ represents the air density (kg m⁻³); SW is the soil moisture factor (dimensionless); and SD represents the snow cover factor (dimensionless), which is the ratio of days with no snow cover to the total number of days studied. A snow cover depth less than 25.4 mm indicates no snow cover. The expression u₁ represents the threshold wind velocity at a height of 2 m (m s⁻¹). The threshold wind velocity for sand lands and sandy grasslands is 5 m s⁻¹; u₂ represents the monthly average wind velocity at a height of 2 m (m s⁻¹); and N_d represents the number of days with a monthly average wind velocity that exceeds the threshold wind velocity.

In this study, EF and SCF were assumed to be unchanged over time and were calculated using the soil raster data that were converted from a 1:1,000,000 soil data shapefile provided by the Harmonized World Soil Database (HWSD) developed by the International Institute for Applied System Analysis (IIASA) of the FAO. For the calculation, the classification criteria for soil texture taken from international data in the database were converted to American criteria to meet the factor calculation requirements by using the logistic growth model proposed by Skaggs et al. [42]. The values of EF and SCF were calculated as follows:

$$EF=29.09+0.31\times Sa+0.17\times Si+0.33\times \frac{Sa}{Cl}-2.59\times OM-0.95\times CaC{O}_3$$

(10)

$$SCF=\frac{1}{1+0.0066\times c{l}^2+0.021\times O{M}^2}$$

(11)

where Sa, Si, Cl, OM, and CaCO₃ represent the content of sand, silt, clay, organic matter, and calcium carbonate (%), respectively.

The value of K′ was calculated with the following equation:

$$K^{\prime }=\cos \alpha$$

(12)

where α represents the slope gradient and was calculated by a digital elevation model (DEM) in ArcGIS. The value of C was computed as follows:

$$C=e^{-0.0483\times SC}$$

(13)

$$SC=\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}$$

(14)

where SC represents the vegetation coverage (%); NDVI_max represents the normalized difference vegetation index (NDVI) value of a highly vegetated grid; and NDVI_min represents the NDVI value of a bare land grid. The values of NDVI_max and NDVI_min correspond to the NDVI values at cumulative frequencies of 95% and 5%, respectively.

2.2.2. Calculation of Retention Rate of WEPS

The retention rate of WEPS more accurately quantified the contribution of ecosystems on WEPS and avoided the impact of climatic factors on WEPS. The retention rate of WEPS was computed as follows:

$$F=\frac{G}{SLR}$$

(15)

where F is the WEPS retention rate.

2.2.3. Calculation of the WEPS Monetary Value

The monetary value of the WEPS can be measured as the cost required to restore sandy land. It was calculated with the following formula:

$$V_{sf}=\frac{SFQ}{\rho ·h}\times c$$

(16)

where V_sf is the WEPS value in CNY (CNY is the Chinese Currency); ρ is the soil bulk density (g cm⁻³), which was 1.25 g cm⁻³ [43]; h is the thickness of the soil, which was calculated as 0.71 m [43]; and c is the average cost of the sand control project (CNY m⁻²), based on the costs of afforestation and grass grid desertification control in The QTP (1.135 CNY m⁻²).

2.2.4. Geographical Detector

The geographical detector (geo-detector) was mainly used to analyze correlations of the WEPS with the five selected impact factors and multiple impact factor interactions. This is because the geo-detector q values have clear physical meaning with no assumption of linearity, allowing us to objectively show that the dependent variable explains 100 × q% of the difference [44]. This study applied factor detectors and interaction detectors in geo-detectors to analyze the correlations between the WEPS and selected impact factors more comprehensively. The formula used to calculate for the value of q was

$$q=1-\frac{SSW}{SST}={\displaystyle \sum }_{h=1}^L{N}_h{\delta }_h{}^2$$

(17)

where h = 1,…, L is the stratification of variable Y or factor X, that is, the classification or partition; N_h and N are the number of units in the layer and the whole area, respectively; δ_h² and δ² are the variance of Y with in the layer and in the whole area, respectively; and SSW and SST are sum of the variances within the layer and the total variance of the whole area, respectively.

Risk area detection was used to judge whether there was a significant difference in the average WEPS between the sub-areas of two influencing factors. This was tested using the statistic

$$t=\frac{{\overline{Y}}_{h=1}-{\overline{Y}}_{h=2}}{{\left[\frac{Var\left({\overline{Y}}_{h=1}\right)}{n_{h=1}}+\frac{Var\left({\overline{Y}}_{h=2}\right)}{n_{h=2}}\right]}^{1/2}}$$

(18)

where ${\overline{Y}}_h$ is the average WEPS of the sub-areas h, n_h is the number of samples in the sub-areas (h), and Var represents the variance.

Ecological detection was used to evaluate whether the influences of two influencing factors, X1 and X2, on the spatial differentiation of the windbreak and sand fixation were significant. This was done using the F statistic:

$$\mathrm{F}=\frac{N_{X1}\left(N_{X2}-1\right)SS{W}_{X1}}{N_{X2}\left(N_{X1}-1\right)SS{W}_{X2}}=\frac{N_{X1}\left(N_{X2}-1\right){{\displaystyle \sum }}_{h=1}^{L1}{N}_h{\sigma }_h^2}{N_{X2}\left(N_{X1}-1\right){{\displaystyle \sum }}_{h=1}^{L2}{N}_h{\sigma }_h^2}$$

(19)

where N_X1 and N_X2 are the sample numbers of impact factors X1 and X2, respectively; SSW_X1 and SSW_X2 represent the sum of the variance within the layers corresponding to the impact factor stratification; and L₁ and L₂ represent the number of layers in impact factors X1 and X2.

Interactions occur between different factors. To evaluate whether the explanatory power of the dependent variable increases when the influencing factors work together, it was necessary to select the dominant factors. Interaction detection can identify interactions between different influencing factors, as showed by the following table (

Table 1. Detecting factor indicators and class break values.

Factor	Class Break Value
	1	2	3	4	5	6	7	8	9	10	11	12	13	14
Y	0.67	2.24	4.26	7.95	28.56
X1	PL	PLA	H	SRM	MRM	BRM	EBRM	-	-	-	-	-	-	-
X2	Ls	SLs	Cs	Drs	Des	Ps	Sas	Wfs	Ses	Ms	As	Fa	Ur	Os
X3	2.06	4.5	7.13	9.76	12.58	15.58	18.76	22.34	27.03	47.68	-	-	-	-
X4	0.1	0.17	0.24	0.33	0.43	0.53	0.62	0.7	0.77	0.91	-	-	-	-
X5	2126.31	3506.49	4523.47	5685.72	6847.98	8010.24	9390.42	11,279.09	13,894.16	18,979.04	-	-	-	-
X6	−0.95	0.68	2.01	3.34	4.67	6.07	7.77	9.54	11.24	14.57	-	-	-	-
X7	0.64	1.07	1.41	1.70	2.00	2.33	2.64	2.94	3.31	5.04	-	-	-	-
X8	2050	2624	3110	3597	4038	4426	4761	5074	5458	8313	-	-	-	-
X9	1.00	2.00	3.00	4.00	>5.00	-	-	-	-	-	-	-	-	-
X10	1.00	2.00	3.00	4.00	>5.00	-	-	-	-	-	-	-	-	-

Y: grade of average WEPS (g km⁻²), X1: landform types (PL: plains, PLF: platforms, H: hills, SRM: small rolling mountains, MRM: medium rolling mountains, BRM: big rolling mountains, EBRM: extremely big rolling mountains); X2: soil type (Ls: leaching soil, SLs: semi leaching soil, Cs: calcareous soil, Drs: dry soil, Des: desert soil, Ps: primary soil, Sas: semi-aqueous soil, Wfs: water-forming soil, Ses: saline earth, Ms: man-made soil, As: alpine soil, Fa: ferro-alumina, Ur: urban area, Os: other soil); X3: slope (°); X4: average annual NDVI, X5: average annual precipitation (mm), X6: average annual temperature (°C); X7: average annual wind speed (m·s⁻¹); X8: DEM (m); X9: average annual GDP (10⁴ CNY km⁻²); X10: population density (human km⁻²).

, Figure 2). Achievement of the geographic detector was operated with geo-detector software (http://geo-detector.cn/, accessed on 15 January 2022). The data needed to be separated and decentralized for the calculation due to their different spatial dimensions. ArcGIS software was used to extract the values of different data.

Changes in the WEPS are simultaneously affected by natural factors such as climate, soil, and vegetation types, as well as by human activities such as grazing withdrawal and grazing prohibition (see previous research [39,40]). In this study, eight natural factors and two socioeconomic factors were selected as detection factors and used to analyze the factors driving change in the average annual WEPS in the QTP from 2000 to 2015 (Table 1). Based on data types and operability, the landform types and soil types of the QTP were divided into 7 and 14 categories, respectively. The numerical variables were classified by the natural discontinuity method in ArcGIS. For example, the population density and GDP per unit area were divided into 5 categories, and the remaining variables were divided into 10 categories (Table 1).

2.2.5. Data Source

Meteorological datasets, including data on the daily temperature, precipitation, wind speed, and solar radiation data from 2000 to 2015, were acquired from meteorological stations. Snow cover data with a 25-km resolution were obtained from a long-term snow depth dataset [45], which was accessed from the Cold and Arid Region Science Data Center (http://westdc.westgis.ac.cn, accessed on 15 January 2022). NDVI values (spatial resolution 500 m) were derived from the international scientific data mirror website of the computer network information center of the Chinese Academy of Sciences (http://www.gscloud.cn, accessed on 15 January 2022). The DEM data (spatial resolution 30 m), GDP data (spatial resolution 1 km), 1:1,000,000 landform data, and land cover data (spatial resolution 1 km) were derived from the Resource and Environment Data Cloud Platform of the Chinese Academy of Sciences (CAS) (http://www.resdc.cn, accessed on 15 January 2022). The 1:1,000,000 soil data shapefile was provided by the Harmonized World Soil Database (HWSD) developed by the International Institute for Applied System Analysis (IIASA) of FAO. During the calculation, the international soil texture classification criteria in the database were converted to American criteria to meet the factor calculation requirements by using the logistic growth model proposed by Skaggs [42]. The 1:4 million calcium carbonate content data were accessed from the Data Sharing Infrastructure of Earth System Science (http://www.geodata.cn, accessed on 15 January 2022). All of the above data were resampled to a 1 km spatial resolution.

3. Results

3.1. WEPS and Its Monetary Value in the QTP

From 2000 to 2015, the total SLR (potential wind erosion) in the QTP was computed to range from 56.36 × 10⁸ t to 68.13 × 10⁸ t. There was an increase of 5.67 × 10⁸ t, and a trend of “rising–declining–declining” was shown. The SLR per unit area was calculated to be from 0 to 60.86 kg m⁻², and the average SLR per unit area was estimated to be from 2.40 kg m⁻² to 2.95 kg m⁻² (Figure 2). The SLR per unit area was higher in the western and northern parts, whereas it was lower in the southeastern parts. This shows that the western and northern areas of the QTP, which contain desert areas, are more threatened by wind erosion.

The total SL (actual wind erosion) in the QTP was calculated to be from 38.55 × 10⁸ t to 43.39 × 10⁸ t during the year (Figure 3). Compared to 2000, the total SL was 0.65 × 10⁸ t lower in 2015. The SL per unit area was the highest in 2005, reaching 60.83 kg m⁻², and a continuous decline occurred from 2005 to 2015. The average SL per unit area was estimated to be 1.67 kg m⁻², 1.90 kg m⁻², 1.75 kg m⁻², and 1.65 kg m⁻² in 2000, 2005, 2010, and 2015, decreasing by 0.02 kg m⁻². The SL per unit area was higher in the northern and western parts and lower in the eastern and southern parts. This shows that strong wind erosion occurred in the northern and western areas of the QTP. In addition, when comparing the SL and SLR in the western area, it can be seen that the vegetation cover significantly reduced the occurrence of wind erosion (Figure 1 and Figure 3).

The total WEPS in the QTP was computed as 17.48 × 10⁸ t in 2000 and 25.15 × 10⁸ t in 2015, and the total monetary value was estimated to be 22.36 × 10⁸ CNY in 2000 and 30.52 × 10⁸ CNY in 2015. The WEPS decreased by 6.39 × 10⁸ t in 2015 compared with the value in 2000, and its monetary value was reduced by 8.17 × 10⁸ CNY. The WEPS per unit area was calculated to range from 0 to 38.40 kg m⁻² in 2000 to 2015 (Figure 2), and the average WEPS per unit area was estimated to be from 0.09 kg m⁻² to 0.13 kg m⁻². The WEPS per unit area was higher in the central and southwestern parts of the QTP, where grassland is the main land type, whereas it was found to be lower in the northern parts (Figure 1 and Figure 3). It can be speculated that the grassland distribution area played a major role in preventing wind erosion in the QTP.

The average F value (retention rate of WEPS) in the QTP was between 57.24% and 62.10% from 2000 to 2015 (similar to that found in previous research [10]), showing an increasing trend. The F value was calculated to range from 0 to 99.94%. Different from the inter-annual variation in the WEPS, potential wind erosion, and actual wind erosion values, the retention rate of the WEPS in the QTP was highest in 2015 and lowest in 2000. Areas with higher WEPS retention rate were mainly located in the southeast where there is better vegetation coverage (Figure 1). A relatively obvious geographical distribution from northwest to southeast was shown.

3.2. WEPS and Its Monetary Value of Different Province and Cities in the QTP

Among the provinces in the QTP, Qinghai and the Tibet Autonomous Region were found to account for more than 90% of the WEPS. The total WEPS in Tibet showed the highest value, reaching 16.15 × 10⁸ t in 2005, and its monetary value was calculated as 20.66 × 10⁸ CNY. This was followed by Qinghai Province and Yunnan Province.

From 2000 to 2015, there were significant differences in the WEPS. The total WEPS in Xinjiang showed a gradual increasing trend over time, and the other four provinces showed a trend of “rising–declining”. The Tibet Autonomous Region and Yunnan Province reached the highest total WEPS values in 2005, whereas in 2010, the maximum values occurred in Qinghai and Sichuan (Figure 4). The total WEPS in all provinces increased from 2000 to 2015, with Qinghai Province showing the largest increase of 3.55 × 10⁸ t, followed by the Tibet Autonomous Region. It can be concluded that Tibet and Qinghai have greater WEPS concentrations.

In 2000, 2010, and 2015, the average WEPS per unit area was the highest for Qinghai province with values of 0.89 kg m⁻², 1.52 kg m⁻², and 1.42 kg m⁻², whereas in 2005, the Tibet Autonomous Region had the highest average WEPS per unit area, and the Yunnan province had the lowest values. From 2000 to 2015, the average WEPS per unit area in Sichuan province showed a gradual increasing trend over time, and that in Yunnan showed a trend of “rising-declining-rising”. The other three provinces showed a trend of “rising-declining”. However, the average WEPS per unit area increased for all provinces, with that in Qinghai increasing the most at 0.52 kg m⁻², followed by Xinjiang with 0.30 kg m⁻².

Among the 37 cities in the QTP, Nagqu was shown to have the highest total WEPS, accounting for 25–33% of the total. In 2015, its total WEPS amounted to 6.17 × 10⁸ t. This was followed by Haixi, accounting for 16–22%. From 2010 to 2015, all cities showed a trend of “rising–declining”. Compared with 2000, the total WEPS decreased in Shannan, Shigatse, Haibei, and Zhangye in 2015. Shigatse decreased by 0.54 × 10⁸ t. Zhangye decreased the least, by only 0.06 × 10⁸ t. The largest increase in the WEPS occurred in Haixi, where there was an increase of 1.91 × 10⁸ t, followed by Naqu and Yushu, which increased by 0.19 × 10⁸ t and 0.17 × 10⁸ t, respectively.

In 2000, the average WEPS per unit area was the highest for Hainan Region at 1.33 kg m⁻², whereas during 2005 and 2015, Naqu reached the highest average WEPS per unit area with a range from 1.82 to 2.43 Guangyuan city had the lowest WEPS per unit area. From 2000 to 2015, the average WEPS per unit area in cities such as Mianyang, Ganzi, Ya’an, Yushu, Chengdu, Deyang, and Guangyuan showed a gradual increasing trend over time. In Shannan, Shigatse, Haibei, and Zhangye there was a volatile decreasing trend, whereas in other cities, a volatile growing trend was shown. Compared with 2000, the average WEPS per unit area in Yushu increased the most (0.71 kg m⁻²), followed by Haixi with 0.64 kg m⁻². The average WEPS per unit area decreased the most in Shigatse (by 0.32 kg m⁻²). Overall, the WEPS in most of the cities in the QTP improved.

3.3. WEPS and Its Monetary Value of Different Ecosystems in the QTP

Among the different types of ecosystems in the QTP, the grassland ecosystem has the highest WEPS, accounting for more than 68% of the total. The WEPS of the grassland ecosystem was estimated to range from 12.7 × 10⁸ t to 18.29 × 10⁸ t from 2010 to 2015. The monetary value was calculated to be from 16.24 × 10⁸ CNY to 23.39 × 10⁸ CNY. The settlement ecosystem had the lowest WEPS. In 2015, the average WEPS per unit area was the highest for the grassland ecosystem, 8.06 kg m⁻², followed by the settlement ecosystem, and the value for the forest ecosystem was the lowest. From 2000 to 2015, the total WPES in all ecosystems showed an increasing trend, among which the forest ecosystem increased the most, reaching 3.61 × 10⁸ t, followed by the desert ecosystem (Figure 5).

3.4. Change of Spatial Pattern of WEPS in the QTP

From 2000 to 2015, the change in SLR per unit area in the QTP ranged from −27.74 kg m⁻² to 15.29 kg m⁻², with an average change of 0.24 kg m⁻². An increase in SLR per unit area mainly occurred in the western areas of northeastern QTP, on the eastern edge, and in the central and northern parts of the area. Both the central and western parts showed a downward trend, especially the western part. The change in SL per unit area in the QTP ranged from −27.04 kg m⁻² to 13.03 kg m⁻², with an average change of 0.02 kg m⁻². An increase in SL per unit area mainly occurred in northeastern QTP, on the eastern edge, and in the central and northern parts. Both the central and western parts showed a downward trend.

The change in WEPS per unit area in the QTP ranged from −14.85 kg m⁻² to 20.40 kg m⁻², with an average change of 0.26 kg m⁻². Increases in the WEPS per unit area mainly occurred in the western part of northeastern QTP, on the eastern edge, and in the central and northern parts. Both the central and western parts showed a downward trend. The WEPS per unit area increased most significantly in Huzuoqi, Xinbalhuyouqi, and Chenbalhu in the west of Hulunbuir. In Zhenglan Banner in the southeast of Xilin Gol League, West Uzhumuqin in the northwest, the central and northern parts of Chifeng City, and the southeast of Ulanchabu City, the decline in the WEPS per unit area was more significant. From 2000 to 2015, the change in F per unit area in the QTP ranged from −89.17% to 97.89%, with an average change of 5.00%. A decrease in F per unit area mainly occurred in the western parts of northeastern QTP, on the eastern edge, and in the central and northern parts. Both the central and western parts showed an upward trend (Figure 6).

3.5. WEPS in the QTP and the Detection of Its Driving Forces

This study set the grade of the average annual WEPS per unit area from 2000 to 2015 as the Y factor, and 23,853 samples were generated using geographic detectors to match and extract attribute values with other raster data layers (Figure 7). The landform type, soil type, slope, average annual NDVI, average annual precipitation, average annual temperature, average annual wind speed, DEM, average annual gross domestic product (GDP), and average annual population density were set as factors X1–X10 (Figure 8). The results show that the detected q values were as follows: slope (0.440) > landform type (0.275) > average annual precipitation (0.223) > average annual wind speed (0.186) > average annual NDVI (0.142) > average annual population density (0.124) > average annual GDP (0.076) > DEM (0.071) > average annual temperature (0.057) > soil type (0.053) (

Table 2. Result of factor detection.

Detected Factor	Landform Type	Soil Type	Slope	Average Annual NDVI	Average Annual Precipitation	Average Annual Temperature	Average Annual Wind Speed	DEM	Average Annual GDP	Average Annual Population Density
q statistic	0.276	0.053	0.440	0.142	0.223	0.057	0.186	0.071	0.076	0.124

). Slope was found to be the most important factor controlling the spatial differentiation of the WEPS, followed by landform type, average annual precipitation, and average annual wind speed. The WEPS grade was not significantly affected by the average annual temperature or soil type. The influences of human factors, such as the average annual population density and average annual GDP, on the grade of the average annual WEPS tended to be less than the influences of most natural factors, which indicates that natural factors have dominated the formation of wind erosion in the area, but human activities have still had an impact on the WEPS in the QTP.

The risk detector was used to detect the areas with the highest WEPS grade for each factor partition to reveal the internal mechanism responsible for the spatial pattern of the WEPS. The risk detection results show that, as the average annual wind increased, the average WEPS grade per unit area gradually increased. As the slope, average annual precipitation, average annual precipitation temperature, average annual GDP, and average annual population density increased, the average WEPS grade per unit area showed a trend of increasing and then decreasing. The highest WEPS grade occurred when the DEM ranged from 5074 m to 5458 m (Table 1, Figure 9). The average annual GDP and average annual population density trends show that appropriate levels of human activity can help to improve the WEPS grade, but excessive human disturbance in areas with moderate economic development and population density can lead to lower WEPS levels. However, the WEPS grade in the most economically developed areas was higher than that in the moderately developed areas, which could be due to the greater investment in, and influence of, environmental protection measures.

The results for the detected factors (

Table 3. Statistical significance between detecting factors.

Detecting Factor	Landform Type	Soil Type	Slope	Average Annual NDVI	Average Annual Precipitation	Average Annual Temperature	Average Annual Wind Speed	DEM	Average Annual GDP	Average Annual Population Density
Landform type
Soil type	N
Slope	N	N
Average annual NDVI	Y	N	N
Average annual precipitation	Y	Y	Y	Y
Average annual temperature	Y	Y	Y	Y	Y
Average annual wind speed	N	N	Y	Y	Y	N
DEM	N	N	N	N	Y	Y	N
Average annual GDP	Y	Y	Y	Y	Y	Y	Y	Y
Average annual population density	Y	Y	Y	Y	Y	Y	Y	Y	Y	Y

Y indicates a statistically significant difference (95% confidence level) between the two factors; N indicates no significant difference.

) show that the average annual precipitation, average annual GDP, and average annual population density were significantly different from all the other influencing factors. Landform type, slope, average annual NDVI, average annual wind speed, and average annual temperature showed significant differences with most of the other influencing factors, whereas the DEM showed no significant differences with most of the influencing factors.

The multi-year average interaction detector results (

Table 4. Interaction type and q value between different detecting factors.

Detecting Factor	Landform Type	Soil Type	Slope	Average Annual NDVI	Average Annual Precipitation	Average Annual Temperature	Average Annual Wind speed	DEM	Average Annual GDP	Average Annual Population Density
Landform type	0.28
Soil type	0.34	0.05
Slope	0.46	0.50	0.44
Average annual NDVI	0.39	0.19	0.57	0.14
Average annual precipitation	0.38	0.25	0.55	0.26	0.22
Average annual temperature	0.32	0.12	0.48	0.22	0.29	0.06
Average annual wind speed	0.37	0.22	0.53	0.22	0.27	0.23	0.19
DEM	0.32	0.13	0.46	0.22	0.28	0.14	0.26	0.07
Average annual GDP	0.31	0.14	0.46	0.17	0.24	0.14	0.20	0.14	0.08
Average annual population density	0.33	0.18	0.47	0.18	0.24	0.19	0.21	0.19	0.14	0.12

The shaded green indicates nonlinear enhancement interaction type between the two factors; the other indicates two-factor enhancement interaction type between the two factors.

) showed that each factor had an enhanced interaction effect on the WEPS of the QTP. Landform type, soil type, average annual NDVI, and average annual precipitation showed nonlinear enhancement relationships with other factors. Slope, average annual temperature, average annual wind speed, DEM, average annual GDP, and average annual population density showed nonlinear enhancement relationships with other factors. The two-factor enhancement between the slope and other factors was found to have a greater impact on the spatial variation of the WEPS grade, and the degree of influence (q value) on the WEPS grade exceeded 44%, which showed that the slope also strengthened the influences of other factors on the spatial distribution of the WEPS.

However, the nonlinear enhancement relationship between the annual average temperature and other factors was found to have a high degree of influence on the spatial variation of the WEPS with a nonlinear enhancement for all factors. It can be concluded that although the single-factor effect of the annual mean temperature and the soil type on the spatial variation of the WEPS were weak, the interaction between the two factors had a higher level of influence on the WEPS than the two factors alone by changing the regional natural environment or human activities. This suggests that the sum of the individual effects of two factors indirectly affects the spatial variation of the WEPS.

4. Discussion

The change pattern of the WEPS was similar to that of the SLR (Figure 6), mainly because the WEPS is more dependent on natural factors, such as the slope, precipitation, and wind. As shown in previous studies, soil type played an important role in the spatial distribution pattern of the sand-stabilization service function [40], soil wind erosion rate was intensified significantly with increasing wind velocities by a power function [46], and the increase in rainfall reduced wind erosion [10]. It seems to be proven that climate changes are the leading factors that affect the soil wind erosion in a region overall, which is consistent with the finding of Li et al. [47].

Increasing fractional vegetation cover can effectively reduce the amount of soil erosion. People can change factors such as the NDVI by improving vegetation coverage and reducing the SL. From 2000 to 2015, a series of ecological conservation and restoration programmes and policies such as Grazing Forbidden Policy and the Grassland Ecological Protection Subsidies and Rewards Program have been implemented in the QTP, which have proven to be effective in reducing soil wind erosion [10]. In terms of changes in land use from 2000 to 2015, the amount of area translated from grassland ecosystem to desert ecosystem was the highest, followed by that from desert to grassland. The increase in the forest ecosystem area was mainly due to conversion from grassland and settlement areas (Figure 10). The change area mainly occurred in eastern part and southern edge of the QTP, where human activities are concentrated [48,49].

Among the different types of land use change, the increase in the forest ecosystem area caused an increase in the WEPS, whereas the increase in the desert ecosystem area caused a decrease in the average WEPS in the QTP. The average WEPS per unit area only decreased in the area that was transformed from a water ecosystem to a settlement ecosystem (Figure 11). As same as the result shown in Chi’s study [4], the land use/cover change (such as returning farmland to forest, shrub, and grassland) played a key role in reducing the soil wind erosion modulus and enhancing regional WEPS. It was found that the increase in grassland area played a certain role in improving the WEPS due to the control of desertification and the implementation of the returning farmland to grassland project. However, due to unsuitable land use, the resulting significant increase in the desert area is the main reason for the decrease in the total WEPS.

Because land use change played an important role in WEPS change in QTP, reasonable land use management will help to promote regional ecological environment protection and sustainable development. In this study, grassland and desert changed significantly, over-grazing, logging, and crop planting have resulted in significant vegetation degradation and desertification, but the implementation of ecological protection projects has promoted the area of grassland and forest [50,51]. For further management, the grassland and desert in the eastern and south-eastern parts of the QTP could be restored to forests engineering practices, plantations, enclosures to limit human disturbances, etc., and the most widely adopted practice in the north and western part of the QTP should be natural restoration. In addition, water shortages in arid and semiarid regions were neglected and the conservation of water bodies should receive more attention.

5. Conclusions

This current study assessed spatial–temporal changes of the WEPS in the QTP from 2000 to 2015. The spatial dynamics of the WEPS were analyzed from the perspectives of various leagues and cities and ecosystems, and its driving forces such as contributions of influential factors to the WEPS and their interactions using geo-detector. The results indicate the following conclusions.

Areas with moderate and high soil erosion intensities were mainly distributed in the western and northern QTP, which contain distributed massive desert and sparse grassland. The annual mean WEPS showed an increasing trend from 2000 to 2015. Slope and weather factors, such as precipitation and wind speed, played leading roles in determining the spatial pattern and grade of the WEPS, and human activities showed an important factor in the variation of WEPS.

In addition, the implementation of projects, such as returning farmland to grassland and desertification control, increased the WEPS by increasing the grassland area. The increase in the desert ecosystem area due to unsuitable land use resulted in a significant reduction in the WEPS in the QTP.

However, there were still some limitations in our study. First, due to the lack of field observation data, this study did not consider the seasonal variation in the wind velocity threshold, thus the seasonal difference in the WEPS could not be analyzed. Second, further parameter localization in the determination of the erodibility boundary of the RWEQ model should be executed in the future to improve the accuracy of the simulation results. Finally, the study of the factors driving the WEPS using the geographic detector model showed the relationships among different factors, whereas the coupling effects of natural factors and human factors need to be studied further.