The NSB is a long and thin belt across northern China (Figure 1A), located in the transition zone between arid and semi-arid areas, divided into three sections from west to east, the Tarim Basin sand-stabilization belt (TB), the Hexi Corridor sand-stabilization belt (HC), and the Inner Mongolia sand-stabilization belt (IM). With a total area of 869,558.5 km². The slope of the study area ranges from 0 to 55 degrees, with the larger slopes concentrated in the northern part of the TB, the HC, and the partial eastern part of the IM. The mean annual precipitation of around 300 mm, with a significant precipitation gradient (Zhang, 2017). The mean annual temperature ranges from −1.9 to 13.5°C. The average wind speed of sand-raising is about 6.5–8.0 m/s. The average wind speed of sand-raising is larger in the middle and smaller in the west and east. FVC increases gradually from west to east. The vegetation types are desert grassland, steppe, forest grassland, forest meadow grassland, scrub meadow, alpine meadow, coniferous forest, deciduous broadleaved forest, hyper-arid shrubs, semi-shrubs, small shrubs, and semi-trees. The distribution of precipitation is generally uneven and highly variable.

We refer to some conventions and classification standards to divide FVC and precipitation gradient. For example, the classification of FVC refers to the Standard for “Soil Erosion Classification and Grading Standards” (SL190-2007) by the Ministry of Water Resources of the People’s Republic of China in 2008. And combined with the actual situation of the study area, the FVC was divided into three classes: FVC<30% (low vegetation cover), 30% ≤ FVC < 60% (medium vegetation cover), FVC ≥ 60% (high vegetation cover; Gao et al., 2015; Figure 1B). According to the convention of annual precipitation division, above 800 mm is the humid zone, 800–400 mm belongs to the semi-humid zone, 400–200 mm belongs to the semi-arid zone, and below 200 mm belongs to the arid zone (Zhang et al., 2011; Gao et al., 2015). The maximum annual precipitation was 685.5 mm in 2020, and most of the study area belongs to the arid zone. For convenience of analysis, this study divides the grades according to the gradient of 0–200, 200–400, and 400 mm<, which accounted for 46.31, 22.76, and 30.93%, respectively, (Figure 1C). The main land use and land cover (LULC) of NSB are grassland and desert (Figure 1D). Demonstrating the spatial distribution pattern of LULC in the belt-shaped study area helps to describe the distribution of ES related to LULC.

Loess plateau science data center, National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China,¹ including average annual precipitation and temperature (Peng et al., 2019), relative humidity, and hours of sunshine. DEM data with 90 m resolution came from the geospatial data cloud platform of the Computer Network Information Center of the Chinese Academy of Sciences.² Land cover and FVC data were obtained from MODIS-based land use product extraction.³ The soil dataset is derived from the Harmony World Soil Database_China subset of the National Qinghai-Tibet Plateau/Third Pole Environmental Data Center, which includes soil texture, sandy soils, chalky soils, clays, and organic carbon. The above data came from different platforms, so this study resampled all the raster data to the same image size with Arcgis, and then computed the raster data in Python.

In this study, sand-stabilization service is represented by the sand stabilization amount, which is equal to the difference between the potential wind erosion amount and the actual wind erosion amount. This project uses the revised wind erosion equation (RWEQ) model (Fryrear et al., 2000) to estimate the sand-stabilization service in NSB. The calculation formula is shown below:

where, G is the amount of sand fixation, kg/m²; S L 1 and S L 2 are the potential amount of wind erosion and the actual amount of wind erosion, kg/m²; Q max 1 , Q max 2 are the maximum transport capacity of potential soil erosion and actual soil erosion, respectively, kg/m²; S 1 , S 2 are the critical field lengths for potential and actual soil erosion, respectively, m; x is the distance of maximum wind erosion occurrence in the downwind direction, m; WF for meteorological factors; EF is soil erodibility factor; S C F as soil crust factor; K′ is the surface roughness factor; and COG is the vegetation cover factor.

Soil conservation service is represented by the soil retention amount, which is equal to the difference between the potential soil retention and the actual soil retention. The Universal Soil Loss Equation (USLE) is the most widely used remote sensing quantitative model of soil loss with good practicality and has been applied and practiced in large regional soil conservation studies in China. Therefore, USLE (Pandey et al., 2007) was selected to assess soil conservation services in NSB ecosystems for this project. The calculation formula is shown below:

where, S C for soil retention, [t / ( hm 2 × a)] ; S E p and S E a are the potential soil erosion and actual soil erosion, respectively, [t / ( hm 2 × a)] ; R is the rainfall erosion force factor, M J × mm / ( hm 2 × h × a ) ; K is the soil erodibility factor, t × hm 2 × h / ( hm 2 × MJ × mm ) ; LS and COG are topography factor and vegetation cover factor, respectively, and are dimensionless. The calculation process for the different factors is shown below:

The water conservation service was calculated using the balance equation of water quantity (Zhang et al., 2001), which treats the water service as precipitation minus evapotranspiration and storm runoff. The calculation index includes annual precipitation, annual evapotranspiration, and annual stormwater production volume. The water balance equation is calculated as follows:

where, WR is the water-bearing capacity, mm; PET is the annual precipitation, mm; QF is the storm water runoff, mm. ET is the actual evapotranspiration, mm.

The calculation of the main indicator factors required for the calculation of the water conservation services is shown below:

Carbon sequestration service is calculated using aboveground biomass multiplied by the biomass-carbon conversion coefficient (Fang et al., 2001; Piao et al., 2009). The main calculation formula is as follows:

where 𝐶O𝑆 is the above-ground carbon stock of terrestrial ecosystems. i is the i-th type of ecosystem; j is the total number of ecosystem types; 𝐴𝐺𝐵𝑖 is the aboveground biomass of the ith ecosystem type; and 𝐶_𝑖 is the biomass-carbon conversion coefficient for ecosystem type i.

Reflected by the biological habitat quality index, the regional biodiversity was evaluated in this project mainly in terms of regional habitat quality and habitat scarcity, which were obtained using the INVEST model calculation (Aneseyee et al., 2020; Wu et al., 2021). The calculation formula is as follows:

Synergies between two ES occur when an enhancement of one ES leads to an increase in another ES (Bennett et al., 2009). While, Trade-offs between two ES occur when one ES is promoted at the expense of the other (Rodriguez et al., 2006). In addition to trade-offs and synergies, there is also a neutral category, which means that an increase in one ES does not result in an increase or decrease in the other. A negative correlation coefficient for two ES that passes the significance test is considered a trade-off, while a positive correlation coefficient that passes the significance test is considered a synergy relationship (Jopke et al., 2015). NSB is composed of three sub-barriers. A certain ES has spatial heterogeneity in each sub-barrier, and this spatial distribution pattern may produce trade-offs and synergies in the whole area. Therefore, this paper used pixel-by-pixel (resampling to align the raster) spatial correlation analysis to calculate the correlation coefficients between ES pairs for five periods from 2000 to 2020 in the study area and conducts significance tests (Figure 2).

Based on the pixel-by-pixel spatial correlation analysis is method-Pearson product–moment correlation coefficient method (Mukaka, 2012). The correlation coefficients between the two groups of ES were calculated separately, and the trade-offs and synergies between the ES were measured based on the positive and negative correlation coefficients and the absolute magnitude of the relationship. The formula is as follows:

where, R is the correlation coefficient, if R is positive, the relationship between two services is synergistic, and vice versa is a trade-off. If R is zero, there is no relationship, and a larger absolute value indicates a stronger correlation, i.e., a greater degree of synergy or trade-off. x and y are the two ES variables; i is the i-th year.

Significance of trade-offs and synergies between ES determined by t-test (Niu et al., 2022). The formula is as follows:

The significance of the interrelationships between ES was judged based on the null hypothesis test t-test of the correlation coefficient. When ∣T∣ < T_{0.05, 3}, i.e., p > 0.05, the original hypothesis is valid and the correlation result is not significant. When T_{0.05, 3} ≤ ∣T∣<T_{0.01, 3}, i.e., 0.01<p ≤ 0.05, rejection of the original hypothesis and more significant correlation results. When ∣T∣ ≥ T_{0.01, 3}, i.e., p ≤ 0.01, the original hypothesis was rejected and the correlation result was highly significant.

Figure 3 shows the spatial distribution of ES. High values of soil conservation service are clustered in the southern part of IM, while low values of soil conservation service are mainly distributed in the west and the plains between HC and IM. Not only the forests and grasslands in eastern IM and eastern HC, but also parts of the cropland in TB show a high carbon sequestration supply. The lower carbon sequestration service supply is mainly found in the TB and desert areas at both ends of the HC. The high sand-stabilization service is distributed in the northern part of IM and the grassland between IM and HC, and the rest of the area had low sand-stabilization service values. The high habitat provision values are mainly distributed in the eastern part of HC where the elevation changes are dramatic, and low habitat provision is generally distributed in grassland and desert areas. The water conservation service in the study area is generally small, but the values of the water conservation service are high in lakes and high-altitude mountain areas.

The distribution of five ES in the study area is inhomogeneous, showing obvious differences in the precipitation gradient and the FVC gradient. We made statistics on the distribution of precipitation gradient, as shown in Figure 4. 33.80–78.10% of the total amount of the five ES are mainly distributed in areas with precipitation greater than 400 mm, followed by areas with precipitation of 200–400 mm. The areas with precipitation in the range of 0–200 mm are the least, accounting for 4.43–33.80% of the total amount of ES. It shows that ES increases with increased precipitation. Figure 5 shows the distribution of five ES across the FVC gradient, which is not obvious compared to the distribution characteristics of the precipitation gradient. The total distribution of habitat provision and carbon sequestration service in the three gradients is basically the same. The soil conservation service is mainly distributed in areas of FVC greater than 0.6. On the contrary, sand-stabilization service and water conservation service are mainly distributed in areas with FVC less than 0.3. From the perspective of ES per unit area, except soil conservation service and water conservation service, other ES increases with the increase of FVC. When the FVC range is 0.3–0.6, the soil conservation service and water conservation service per unit area reach the maximum.

Figure 6 shows the areas where synergies and trade-offs occur between the five ecosystem services. As shown in Table 3. The synergy between sand-stabilization service and habitat provision accounts for 86.67% of the study area, and the significant synergy accounts for 49.08% of the study area. Synergy, especially significant synergy, is widely distributed in the study area, but the trade-off area only accounts for 2.73%. The significant synergy/trade-off between sand-stabilization service and soil conservation service is staggered in the east of IM, and the significant synergy area (13.54%) is larger than the significant trade-off area (3.91%). The area of trade-off between soil conservation service and habitat provision is the largest among the 10 relationships (34.52%), mainly distributed in HC and TB areas with weak trade-off relationship. The synergy between sand-stabilization service and carbon sequestration service accounts for 43.52% of the study area, and the remarkable synergy is clustered and distributed in the east of the study area. The distribution of habitat provision and carbon sequestration service is similar to that of sand-stabilization service and carbon sequestration service. The unconnected relationship with carbon sequestration service and soil conservation service (43.97%) is distributed in the southern part of TB, the western part of HC, the plain between HC and IM, and the land use type is desert. The significant synergy/trade-off relationship are staggered in the eastern part of the study area. The four relationships between sand-stabilization service and water conservation service, carbon sequestration service and water conservation service, soil conservation service and water conservation service, habitat provision, and water conservation service are the relationships between water conservation service and other four ES. Obviously, the relationships including water conservation service services, unrelated relationships all account for more than 75% of the study area. Apart from the forests and grasslands in the middle of HC, the forests and grasslands in the southeast of IM, and the forests and grasslands in the north of TB, unrelated relationships are widely distributed in grasslands and deserts in the study area.

Since the unrelated area accounts for most of the total study area, the relationship between the two ES is better captured in terms of the area ratio of synergies to trade-offs when exploring the trade-offs between ES. The ratio of synergy to trade-off between sand-stabilization service and habitat provision was 31.74: 1, and the area of significant synergy represented 49.08% of the total area of the study area. The ratio of synergy to trade-off between carbon sequestration service and water conservation service was 21.53: 1. The area of significant synergy accounted for 9.45% of the total area of the study area, while the area of significant trade-off accounted for only 0.01% of the total area of the study area. In terms of the area ratio of synergy to trade-off, the synergy area was greater than the trade-off area in all 10 ES pairs, with an average ratio of 8.85:1.

See Figures 7, 8. As precipitation increased, trade-offs and synergies between 10 paired ES increased, while their trade-offs decreased. Specifically, the synergies between soil conservation service and habitat provision, sand-stabilization service and habitat provision, and sand-stabilization service and soil conservation service all exceeded 80% of the study area. And their distribution pattern was consistent with the precipitation gradient. Trade-offs and synergies between pairs of ES, including water conservation service (habitat provision and water conservation service, soil conservation service and water conservation service, carbon sequestration service and water conservation service, and sand-stabilization service and water conservation service) were mainly distributed in the northern end of TB, central HC, and central and eastern IM, accounting for only 16.21–24.10% of the study area. However, note that the uncorrelated relationships between these paired ES were widely distributed throughout the study area. Nonetheless, the synergies still account for the vast majority in different precipitation gradients, and the synergies between paired ES also increase with the increase of precipitation. For sand-stabilization service and carbon sequestration service, carbon sequestration service and soil conservation service, and carbon sequestration service and habitat provision, the trade-offs and synergies between these ES pairs accounted for about 50–60% of the study area, and again, synergy effects dominate. Similarly, synergies accounted for the majority, which were mainly distributed in IM and scattered in TB and HC. Trade-offs and synergies did not exist in about 40–50% of the areas, which were mainly desert areas. Among the 10 paired ES, sand-stabilization service and water conservation service changed most obviously with the gradient of precipitation. In the annual precipitation range of 0–200 mm, the area ratio of trade-offs and synergies was 4.23, while in the annual precipitation range of 400 mm <, it reached 24.59. Compared with other paired ES, the synergies of sand-stabilization service and habitat provision accounted for the largest proportion of each gradient area, indicating that in arid and semi-arid regions, increasing biodiversity can promote more sand-stabilization service in ecosystems. In addition, the ES pairs containing soil conservation service had the smallest synergies than other paired ES, indicating that it was more difficult to play the synergy between ES pairs containing soil conservation service than other paired ES, but with the increase of precipitation, this gap will be reduced. Although less than 20% of the regions have trade-offs and synergies between carbon sequestration service and water conservation service, the ratio of synergy and trade-off is close to sand-stabilization service and habitat provision in each precipitation gradient. The study showed that the synergy between ES was greater in areas with more precipitation, and this ratio increased with increasing precipitation.

We compared the trade-offs and synergies of ES along the FVC, See Figures 8, 9. Similarly to the change in trade-offs and synergies between paired ES with precipitation gradients, as the FVC increased, the synergies between other paired ES increased except sand-stabilization service and habitat provision. Furthermore, the distribution law of soil conservation service and habitat provision, sand-stabilization service and habitat provision, and sand-stabilization service and soil conservation service, as well as the various characteristics of the trade-offs and synergies with FVC of paired ES containing water conservation service and carbon sequestration service, respectively, were also similar to those between pairs of ES and precipitation gradient. Nonetheless, along the FVC, there were some different regularities in the trade-offs and synergies relationship between paired ES. Compared with other paired ES, the synergy of sand-stabilization service and habitat provision accounted for the largest proportion of each gradient area, indicating that in arid and semi-arid regions, increasing biodiversity can promote more sand-stabilization service in ecosystems. However, FVC had little effect on this synergy, and the synergy of sand-stabilization service and habitat provision did not increase with the increase of FVC. Although less than 20% of the regions had trade-off and synergy between carbon sequestration service and water conservation service, the ratio of synergy and trade-off is close to sand-stabilization service and habitat provision in each FVC gradient. It showed that in areas with low FVC, the synergy between ES was less, with the increase of FVC, its proportion was increasing. However, with the change of FVC, paired ES has a different trend. As FVC increases, the synergy of carbon sequestration service and water conservation service increases. This shows that increasing FVC in arid and semi-arid regions can first increase carbon sequestration, and in some regions can increase water conservation service through synergy effects.

Previous studies suggest that land use/land cover affect trade-offs and synergies between paired ecosystems (Hasan et al., 2020; Liu et al., 2022). We also investigated the characteristics of the trade-offs and synergies of paired ES across different land uses and land covers (Figure 10). Among the eight land use/land cover, except that the trade-off values of sand-stabilization service and water conservation service and soil conservation service and water conservation service in the impervious is larger, paired ES has been synergy dominant in other land use types. However, the same paired ES differed among different land use types. Contrary to Shao et al. (2020), research that different spaces in the same LULC will have different effects on ES there are also differences in the performance of different ES in the same land use/land cover. The research of Vigl et al. (2017) also shows that land management types and biophysical conditions make the provision of ES unstable. In forests, grasslands, and barrens, sand-stabilization service and habitat provision have greater synergies than other paired ES. This finding has practical implications, suggesting that through anthropogenically increased vegetation cover (i.e., afforestation), similar benefits may be expected for other paired ES, while the synergy of sand-stabilization service and habitat provision may increase more in multiple paired ES. And this result also shows that with the emphasis on biodiversity, sand-stabilization service gains the most synergy from it. This is of great significance to protect wildlife and reduce soil wind erosion. In shrub areas, water conservation service and carbon sequestration service have greater synergy than other paired ES. In semi-arid and arid regions, these types of land use have better hydrothermal conditions, and increasing carbon sequestration service can promote the synergy growth of other ES. The finding is unexpected, but has practical implications, as low shrubs produce more ecological benefits than tall trees in semi-arid and arid regions. Adding low shrubs can first increase carbon sequestration service, and it is beneficial to reduce surface wind speed, reduce wind erosion, and can also increase habitat provision growth through synergy. However, it is more significant that in precipitation-limited arid and semi-arid regions, the transpiration and actual evapotranspiration of low shrubs are lower than those of tall trees (Zhang and Huang, 2013), and low shrubs can produce more water conservation service. Therefore, in arid and semi-arid regions, in an attempt to improve the ecological environment through artificially increased vegetation coverage, appropriate vegetation types should be selected according to the target (the type of ES expected to be improved). For example, in desert-dominated NSB, by improving the sand-stabilization service of the ecosystem to control desert expansion, the choice of planting trees and grasses may reap more benefits. In smaller areas, such as cities, adding shrubs may be more appropriate. Therefore, in order to make better use of the synergistic effect of paired ES and improve the ecological environment, the trade-off between trees/grass and shrubs is very important.

Through the above analysis of LULC, increasing the greenness of land is the most effective way to increase the synergy relationship among ecosystem services. For managers, it is the most effective and scientific management measure to continuously promote desertification control and increase vegetation cover in the Northern Sand-stabilization Belt area.

From the precipitation distribution, 46.31% of the study area in 2020 has less than 200 mm of annual precipitation, which is an arid region. The precipitation distribution is similar to FVC, which gradually decreases from east to west. However, there are higher precipitation and FVC values in the middle of the HC. In TB region, precipitation increases with increasing altitude. Paired ES have different sensitivities to vegetation cover and precipitation due to differences in distribution patterns. Compared with precipitation, there are more pairwise synergistic changes in the effect of FVC. For example, among the 10 ES pairs, the synergy between sand-stabilization service and water conservation service, sand-stabilization service and carbon sequestration service, carbon sequestration service and water conservation service, carbon sequestration service and habitat provision increased significantly with increasing FVC. Because in arid and semi-arid regions, high evaporation and low precipitation are the norms (Yin et al., 2005). Vegetation can increase surface roughness and shade and reduce surface evaporation, wind erosion, and water erosion. Synergy of ES pairs that increase with precipitation is only sand-stabilization service and soil conservation service. The study by Zhang et al. (2022) shows the high contribution of climatic factors to vegetation recovery in northern China. Wang X. et al. (2022) showed that a 30% increase in precipitation in desert areas stimulated the growth of sophora alopecuroides, stipa breviflora, pennisetum centrasiaticum, and other species. Thus, increased precipitation indirectly contributes to the synergy of sand-stabilization service and soil conservation service through increase vegetation. Due to the high overlap between the distribution of precipitation and FVC, the other four ES pairs were not differently affected by precipitation and FVC.

There are some limitations to this study. Uncorrelated relationships between pairs of ES account for a large proportion, especially for ES pairs containing water conservation service. We analyze the reasons, and the study area is located in northwest China, dominated by deserts and grasslands. While most of the water conservation service values in these areas are 0. When exploring the trade-offs and synergies between water conservation service and other ES, the correlation coefficient will be 0, resulting in a large area of irrelevant relationships. When the value of a certain ES of a grid is zero in five periods, then the relationship between this ES and other ES is irrelevant. There is also a situation that when the value of a certain ES does not change for five periods, it will also get an irrelevant result with other ES calculations. In addition, the interpolation method and resolution of the underlying data are also important influencing factors. For example, the uncorrelated relationships of sand-stabilization service and soil conservation service, carbon sequestration service and soil conservation service, soil conservation service and water conservation service, and soil conservation service and habitat provision are distributed in the eastern part of the study area in a grid-like form. The reason for this phenomenon may be that the grid meteorological data obtained by spatial interpolation is used in estimating soil conservation service. As far as the research content is concerned, this study only evaluated 5 ES, and did not evaluate other important ES (such as crop production and climate control). The estimation model involves a limited number of factors, and the model itself has some inaccuracies. In future studies, in addition to further investigating and addressing these limitations, we will also explore the driving mechanism of ES trade-off and synergy and the model of the relationship between its sustainability, so as to better provide information for ecosystem management measures, so as to achieve the goal of harmonious development between man and nature.