Improved Understanding of Trade-Offs and Synergies in Ecosystem Services via Fine Land-Use Classification and Multi-Scale Analysis in the Arid Region of Northwest China

Abstract

Ecosystem services (ESs) serve as a fundamental cornerstone for upholding global biodiversity and promoting human well-being. ESs trade-off and synergy are supposed to be significantly affected by climate change (CC) and land use/cover change (LULC). However, the limited availability of finely classified future land-use data and integrated landscape change models incorporating climate change scenarios has hindered our understanding of the trade-off and synergistic patterns and controls of ESs at multiple scales, particularly in arid areas. Here, a future multi-scenario ESs trade-off/collaborative assessment framework (SD-PLUS-InVEST model) for multi-scale conversion and refined land-use classification was developed by coupling the patch-generated land-use simulation (PLUS) model, system dynamics (SD) model, InVEST model, geographically weighted regression (GWR) model, optimal parameter geographical detector (OPGD) model, and structural equation model (SEM). The four ESs, namely carbon storage (CS), habitat quality (HQ), water conservation (WC), and soil conservation (SC), were assessed. Further, multi-scale ESs were evaluated under different climate change and development scenarios (i.e., the SSP1-2.6 and ecological protection scenario, SSP1-2.6-EP; SSP2-4.5 and natural development scenario, SSP2-4.5-ND; SSP5-8.5 and economic growth scenario, SSP5-8.5-EG). The results demonstrated that the arid region of northwest China (ANWC) was experiencing a significant and continuous warming trend accompanied by increased humidity. There will be a significant decrease in the areas occupied by paddy fields, natural forests, and permanent glaciers among the 24 LULC types. Conversely, there will be a substantial increase in dry land, high-coverage grassland, and urban construction land areas. According to the SSP1-2.6-EP, SSP2-4.5-ND, and SSP5-8.5-EG scenarios, the comprehensive land-use dynamic degrees were estimated to reach 2.58%, 4.08%, and 4.74%, respectively. The LULC resulting from CC exacerbates the differences in the four ESs of ANWC. In particular, CS and HQ experience significant reductions in 2100. Conversely, WC and SC show notable increases during the same period. The changes in CS, HQ, WC, and SC reach 11.36 × 10⁸ m³, 1735.25 × 10⁸ t, −1.29 × 10⁸ t, and −0.009, respectively. The four ESs of CS, HQ, WC, and SC in ANWC display a synergistic relationship. This synergy is influenced by the heterogeneous spatial distribution of CS, HQ, WC, and SC, with the strongest synergy observed between CS and HQ and the weakest between CS and WC. Interestingly, the distribution differences in ESs synergy were amplified at watershed, county, and grid scales in mountainous areas, with the most significant detection differentiation occurring at the grid scale. Furthermore, the detection of spatial heterogeneity in the four ESs can be attributed to various factors. These factors include the drought index (q = 0.378), annual average precipitation (q = 0.375), economic density (q = 0.095), vegetation coverage (q = 0.262), and soil bulk density (q = 0.077). Our results highlight the importance of CC in influencing ESs. The spatial variations in ESs trade-offs and coordination at different scales, particularly the pronounced differences observed in mountainous areas, underscore the need to prioritize the conservation of arid mountainous regions in terms of future policy making.

1. Introduction

The ecosystem, as a vital pillar of the global social ecosystem, constitutes an important foundation for human existence and progress and provides a variety of products and services for mankind [1,2]. In recent decades, the global landscape has witnessed substantial transformations due to climate change (CC) and land use/cover change (LULC). These alterations have had a far-reaching consequences on the delivery of ecosystem services (ESs) worldwide, as they have altered both the structure and functioning of Earth’s ecosystems [3]. In the 21st century, the acceleration of global change has been primarily attributed to rapid urbanization and industrialization. These factors have resulted in the widespread depletion of both natural resources and non-renewable resources, while simultaneously exerting mounting pressures on ecosystems. Consequently, this poses a substantial risk to the quality of human living environments and the sustainable development of socio-economic systems [4]. The UN Millennium Ecosystem Assessment report indicated that over 60% of global ecosystem services (ESs) have experienced or are currently undergoing degradation [5,6]. It was urgent to deeply understand the process and mechanism of ecosystem change under the dual disturbance of global CC and human activities and to evaluate and predict the change trend of ESs [7,8,9].

Over the past few decades, the research community has made numerous endeavors in the categorization of ESs, exploration of underlying drivers, spatial quantitative assessment, identification of trade-offs, and understanding of synergistic interactions [10,11,12,13]. However, there are still some problems worth exploring in the existing research. Firstly, there was a lack of finely classified future land-use data [14]. Secondly, the scale effect of ESs was not fully understood [15]. Thirdly, less attention was paid to arid areas with minimal human disturbance, particularly natural lands such as grasslands, forests, and sandy lands [16]. Moreover, the combined effects of CC and human activities had introduced substantial uncertainties in comprehending the underlying mechanisms of evolution and spatial dynamics, which were crucial for accurate prediction of future ESs [17,18,19,20].

The IPCC (Intergovernmental Panel on Climate Change) proposed the ScenarioMIP (Scenario Model Intercomparison Project) climate prediction scenario in the latest sixth International Coupled Model Intercomparison Project (CMIP6) [21,22]. This provides strong support for multi-scenario simulation of ESs that take into account CC and socio-economic development [17]. Furthermore, the advancement of land-use simulation technology has also played a vital role in enabling multi-scenario simulations of ESs. Future land-use simulation endeavors primarily focus on anticipating the spatial distribution patterns and quantitative makeup of land. These efforts aim to forecast LULC across different classes, taking into account both the quantitative and spatial aspects of land allocation [23,24,25,26,27]. Some recent studies have used a combination of spatial simulation of land utilization (CA model, FLUS model) and land quantity structure models (Markov model, SD model) to simulate future LULC and corresponding ecological environment changes [17,28]. Nonetheless, current models exhibit a deficiency in their ability to flexibly handle changes in multi-class land-use patches, particularly for fine-scale LULC simulations. These models lack a flexible mechanism to effectively capture the intricate dynamics of land-use transitions across various classes.

The patch-generating land-use simulation (PLUS) model, built upon enhancements made to the FLUS model, tackles the challenge of effectively mimicking transformations in types of land utilization at the micro-scale for future circumstances [18,19]. Still, the model lacks the ability to predict the quantity demand for land use under prospective CC situations [19]. This is where the SD model comes in, as it can address the aforementioned issues and complement the technical path offered by the PLUS model. Based on this premise, the present study puts forward a comprehensive framework for trade-off and co-prediction of ESs, taking into account the impacts of CC and socio-economic development. The primary objective was to unveil the intricate fine-scale variations of ESs in future multi-scenarios within representative arid regions [18,29].

The CC and human disturbance exert a substantial influence on ESs within the arid region of northwest China (ANWC). Whether in the past or in the future, the ANWC is experiencing a process of warming and wetting, which brings more challenges to the changes of ESs in the region. Thus, this study places emphasis on the management and restoration of diverse mountain–oasis–desert systems in the ANWC under the influence of CC and human disturbances.

The specific objectives of the present study were (1) to demonstrate the spatiotemporal pattern traits of LULC in the past–present–future scenario, (2) to discloses the evolution characteristics of ESs profit and loss and the trade-off and synergy pattern of ESs under different scales and scenarios, and (3) to identify the main driving forces of ESs changes in different scenarios. Capitalizing on these research objectives, the framework proposed in this study offers a comprehensive reference for regional ecological management decisions. It provides valuable insights for complementing status quo assessment results, proactively identifying ecologically fragile areas, and ensuring regional ecological security.

2. Materials and Methods

2.1. Study Area

The ANWC is located in the Eurasian hinterland, specifically in the northwest region of China, spanning geographical coordinates ranging from 73°~107°E and 35°~50°N. With an approximate extent of 2.35 × 10⁶ km², it encompasses the entire Xinjiang Uygur Autonomous Region, Gansu Province, Inner Mongolia Autonomous Region, and a fraction of Ningxia Hui Autonomous Region (Figure 1a,b). The ANWC is characterized by a diverse topography, predominantly consisting of mountains, plateaus, and basins, with elevations ranging from −160 to 7311 m (Figure 1c). Notably, the region encompasses vast desert and Gobi landscapes, including China’s major deserts such as the Taklimakan, Badain Jaran, Gurbantunggut, and Tengger. Major mountain ranges within the region include the Altai, Tianshan, Kunlun, Qilian, and Helan. The primary basins comprise the Junggar and the Tarim. In terms of climate, the ANWC experiences a typical temperate continental climate, with most areas receiving an average annual precipitation of less than 200 mm. The region is primarily drained by inland rivers, such as the Tarim, Shiyang, and Heihe, with the Irtysh River being an outflow river (Figure 1d). The ANWC stands out as a highly impacted area in the face of global CC [30]. In recent years, the severity of ESs issues, such as oasis shrinkage, glacier retreat, and water scarcity, has been intensified by the progression of CC and socio-economic development.

2.2. Data Sources

2.2.1. LULC Datasets

The LULC dataset spanning from 1980 to 2020 was acquired from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn, accessed on 14 October 2022). Specifically, the remote-sensing image data for the years 1980, 1990, 1995, 2000, 2005, and 2010 were sourced from NASA Landsat 5 remote-sensing data. Additionally, the remote-sensing image data for 2015 and 2020 were obtained from Landsat 8 remote-sensing data, providing a spatial resolution of 30 m. To generate these LULC datasets, a series of processing steps were implemented, including band synthesis, geometric correction, and image enhancement. These steps were carried out through a combination of human–computer interactions. To ensure accuracy, the overall classification accuracy was verified to be 90.1% by comparing random sampling with field survey data.

Types of land utilization in the research area were categorized into 24 classes, namely paddy field (PF), dry land (DL), forest land (FL), shrub forest (SB), sparse land forest (SF), garden land (GL), high-cover grassland (HG), medium-cover grassland (MG), low-cover grassland (LG), river (RV), lake (LK), artificial reservoir (AR), permanent glacier snow (GS), wetland (WL), urban construction land (UL), rural residential land (RL), industrial and mining construction land (ML), sandy land (SL), Gobi (GB), saline–alkali land (SA), marsh land (MS), bare land (BL), bare rock land (RL), and alpine desert (AD). It was important to note that the classification of high, medium, and low coverage grasslands was based on vegetation coverage. High-coverage grasslands corresponded to a coverage greater than 50%, medium-coverage grasslands ranged from 20% to 50% coverage, while low-coverage grasslands featured a vegetation coverage between 5% and 20%.

2.2.2. Meteorological and Remote-Sensing Datasets

The climatic information utilized was obtained from two sources. The first set of data covers the period ranging from 1980 to 2020 and was obtained from the National Meteorological Station Scientific Data Center (http://data.cma.cn, accessed on 28 May 2022). This collection of data was acquired based on monitoring data from meteorological stations within the studied region using spatial interpolation, coordinate transformation, resampling, and accuracy verification, among other techniques. The second set of data pertains to the year 2100 and comprises monthly average temperature and precipitation datasets from multiple scenarios and models, with a resolution of 1 km, spanning from 2021 to 2100, within China [31]. These data were obtained through diverse climate models, and the dataset is publicly accessible at http://data.tpdc.ac.cn, accessed on 14 October 2022. Additional information about the data sources is presented in Figures S3–S6.

The evapotranspiration data and drought index data used were obtained from two separate datasets. The evapotranspiration data were sourced from the Global Aridity Index (Global-AI) dataset, while the drought index data were acquired from the Global Reference Evapo-Transpiration (Global ET0) dataset. Both datasets have a spatial resolution of 30 arc-seconds. These datasets are available online through the CGIAR-CSI GeoPortal at the following address: https://cgiarcsi.community/2019/01/24/globalaridity-index-and-potential-evapotranspiration-climate-database-v3/, accessed on 24 January 2019. By utilizing these datasets, we were able to access comprehensive and reliable information on evapotranspiration and drought conditions for our analysis.

The NDVI remote-sensing image data utilized was obtained from the MOD13Q1 dataset available on the Geospatial Data Cloud platform (www.gscloud.cn, accessed on 22 November 2022). The spatial resolution of this dataset was 250 m. For each research period, the imaging time spanned from July to August, and the cloud cover was maintained below 5%. The data underwent several preprocessing steps, including mosaic, projection conversion, cutting, and unit conversion.

Furthermore, the initial version of the VIIRS global nighttime light annual data was obtained from the Earth Observation Group (EOG) and can be accessed at https://eogdata.mines.edu/, accessed on 14 October 2022. This dataset provides valuable information regarding global nighttime light emissions for our analysis.

2.2.3. Soil and DEM Datasets

The soil data (i.e., soil organic carbon density, soil depth, soil bulk density, soil type, soil available water content, soil pH) utilized was obtained from the basic attribute dataset of China’s high-resolution national soil information grid for the period 2010–2018. This dataset, with a resolution of 90 m, was sourced from the National Earth System Science Data Center, which is part of the National Science & Technology Infrastructure of China (https://www.geodata.cn, accessed on 23 September 2022). The soil erosion data used had a resolution of 1 km and was acquired from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn, accessed on 23 September 2022). The soil erosion encompasses primarily three main types: hydraulic erosion, wind erosion, and freeze–thaw erosion.

The digital elevation model (DEM) employed had a resolution of 30 m and was obtained from the Geospatial Data Cloud platform (http://www.gscloud.cn, accessed on 12 October 2022). It was utilized to calculate various topographic parameters, such as aspects, topographic relief, and slopes. The inclusion of topographic relief was crucial for macroscopic landform characterization and analysis [32].

2.2.4. Socioeconomic and Other Environmental Datasets

The traffic network and river network data used in this study were obtained from the 1:1 million National Geodatabase, which can be accessed at https://mulu.tianditu.gov.cn/, accessed on 12 August 2020.

The gross domestic product (GDP) and population (POP) grid dataset with a resolution of 1 km originated from the Data Registration and Release System of the Resource and Environmental Science Data Center of the Chinese Academy of Sciences. The dataset is available at http://www.resdc.cn/, accessed on 14 October 2022.

Additionally, it utilized the socioeconomic database titled “Shared grid dataset of population and economy along socioeconomic paths” published by Chen et al. (2020) [33]. The dataset can be found at https://cstr.cn/31253.11.sciencedb.01683, accessed on 3 November 2022.

Furthermore, the oasis dataset was obtained from the spatial and temporal distribution dataset of oases in arid areas of China (2000–2019). This dataset, with a resolution of 30 m, was downloaded from the National Glacier Permafrost Desert Science Data Center. It can be accessed at http://www.ncdc.ac.cn, accessed on 9 November 2022.

2.2.5. Environment Variables

To investigate the influencing factors of LULC and ESs, a total of 24 environmental variables were selected. These variables were carefully chosen by referring to relevant studies that focused on the environmental characteristics of the study area [28,34,35].

These factors (i.e., X1: drought index, X2: soil available water content, X3: slope aspect, X4: distance from oasis, X5: distance from railway, X6: distance from residential area, X7: distance from river, X8: digital elevation model (DEM), X9: soil erodibility, X10: potential evapotranspiration, X11: vegetation coverage, X12: economic density, X13: nighttime light index, X14: soil pH, X15: annual average precipitation, X16: relief intensity, X17: terrain slope, X18: soil depth, X19: soil organic carbon (0–5 mm), X20: soil type, X21: perennial mean temperature, X22: volume weight of soil, X23: population density, and X24: distance from road). For a visual representation of the spatial distribution characteristics of these variables, Supplementary Materials Figure S1 is provided.

Based on the principle that the lowest resolution determines the highest resolution in grid spatial analysis and is limited by the spatial operation simulation of multi-layer superposition, a spatial resampling process of 1 km × 1 km was employed for all grid layers to obtain the final results of ESs.

2.3. Methods

The research framework, as depicted in Figure 2, includes several key steps. Firstly, three future development scenarios were constructed by integrating the CMIP6 climate model scenario with potential future development models. These scenarios were named as follows: SSP1-2.6 and ecological protection scenario (SSP1-2.6-EP), SSP2-4.5 and natural development scenario (SSP2-4.5-ND), and SSP5-8.5 and economic growth scenario (SSP5-8.5-EG).

Next, the future land-use quantity demand under each scenario was simulated using a system dynamics (SD) model. Subsequently, the allocation of land utilization requirements to various dispersion patterns was accomplished using the PLUS model. This allows for the examination and evaluation of the spatial distribution attributes and temporal dynamics of land use across diverse scenarios.

Ultimately, the InVEST model was utilized to assess the trade-offs and synergies among four major ESs: carbon storage (CS), water conservation (WC), soil conservation (SC), and habitat quality (HQ). The evaluation was conducted across historical, present, and future time periods of the ANWC, thereby providing an in-depth analysis of these ESs over time. This analysis allows for the identification of break-even evolution characteristics and ESs trade-offs/synergies under varying CC and socio-economic development scenarios. Furthermore, it helps in identifying the primary driving forces responsible for changes in ESs within different scenarios.

2.3.1. Measurement of LULC

The land-use dynamic degree consists of two components: the single land-use dynamic degree (K) and the extensive land-use dynamic degree (LC). These components quantify the degree of change in a specific or comprehensive land-use type over a given period. A higher value indicates a stronger level of LULC. The formula was as follows [34,35,36].

$$K=\frac{U_b-U_a}{U_a}\times \frac{1}{T}\times 100\%$$

(1)

$$LC=\left(\frac{{\displaystyle \sum _{i=1}^n\Delta L{U}_{i-j}}}{2{\displaystyle \sum _{i=1}^{n}L{U}_i}}\right)\times \frac{1}{T}\times 100\%$$

(2)

where U_a and U_b denote the spatial extent of specific land-use types from the commencement (time point a) to the conclusion (time point b) of the designated study period. T corresponds to the duration of the study period, measured in years. LU_i represents the initial area occupied by the land-use type i under investigation, while ΔLU_i−j denotes the absolute magnitude of land-use conversion between type i and another land-use type j during the study period. Here, n symbolizes the total number of distinct land types considered in the analysis (n = 1, 2, 3, …, 24).

The land-use transition matrix was utilized to depict the dynamic process of mutual transformation among different land-use types within the ANWC from the beginning to the end of the study period. The formula for calculating the transition matrix was as follows [37,38].

$$S_{ij}=\left[\begin{array}{cccc}{S}_{11}& S_{12}& \cdots & S_{1n}\\ S_{21}& S_{22}& \cdots & S_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ S_{n1}& S_{n2}& \cdots & S_{nn}\end{array}\right]$$

(3)

where S was the area (km²), n was land-use type in ANWC number, and i and j represent the land-use types at the starting and ending points of the study period, respectively.

2.3.2. Simulation of Future Land-Use Quantity Demand

When constructing a system dynamics (SD) model, the simulation area was considered as a relatively autonomous system. To develop the model effectively, it was divided into six distinct subsystems, namely population, economy, society, resources, environment, and policy. Each of these subsystems represents a specific aspect or dimension that contributes to the overall dynamics of the system. By structuring the model in this manner, it allows for a comprehensive analysis and understanding of the interdependencies and feedback loops within the system [37]. The construction of the system dynamics (SD) model involves considering various factors within each subsystem, including population growth, social and economic development, energy consumption, environmental protection, policy adjustment, and other relevant influences. These factors collectively drive the changes in LULC patterns, ensuring a balance between the supply and demand dynamics of the overall system [32].

Within the ANWC, the land-use types were categorized into 24 distinct classifications. The interactions between each subsystem and variable were analyzed based on the respective driving factors. By conducting multiple simulation tests, the relationships and changes between variables were determined. The SD model was subsequently established using Vensim software (Figure S2). For more detailed information and guidelines regarding the Vensim software, it is possible to download them from https://vensim.com/, accessed on 12 June 2022.

To validate the accuracy of the simulation results, this study employed a simulation period spanning from 1980 to 2100, with an annual time step. Specifically, the simulation was divided into two stages: the simulation stage (1980–2020) and the prediction stage (2020–2100). During the simulation stage, model parameters were established, evaluated, and tested using historical data. The historical LULCs were matched with the period from 1980 to 2020, and the accuracy of the model was verified against the historical stage data. The absolute error for all land-use types was found to be less than 5%, indicating the effectiveness of the model. Considering the base year of scenario simulation as 2020, the future changes in land-use demand were simulated by inputting parameters corresponding to different scenarios.

2.3.3. Simulation of Future Land-Use Spatial Distribution

The PLUS model was an LULC simulation model that combines the land expansion analysis strategy (LEAS) rule mining framework with the multi-class random patch (CARS) cellular automaton (CA) model. This integration allows for a comprehensive analysis of land expansion patterns and their impacts on LULC dynamics. The LEAS framework provides a systematic approach for identifying rules governing land expansion based on historical data, while the CARS CA model simulates the spatial patterns of land conversion and allocation across multiple land-use classes. By leveraging these two components, the PLUS model offers an enhanced understanding of LULC changes and supports informed decision making regarding land management and resource allocation [23,38,39]. To analyze land-use expansion between two periods, the spatial patches of land-use expansion were extracted using the LULC data. The random forest (RF) algorithm was employed to identify the types of land expansion and their associated drivers. Additionally, during the base period, the development probability of each land-use type and the contribution of each driver to the growth of specific land-use categories were assessed. By integrating multi-objective planning algorithms, the structure of land use under different scenarios was determined, considering factors such as land-use types, development probabilities, and driver contributions [38]. This approach allows for a comprehensive understanding of land-use dynamics and aids in formulating effective land-use strategies [40,41]. The specifics of the PLUS model can be found in Supplementary Material S3.1, while the simulation verification is presented in Supplementary Materials S3.2, providing a comprehensive overview of its functionality and performance.

2.3.4. Future Land-Use Scenario Setting and Goal

To simulate and predict LULC under future CC conditions, key parameters such as temperature, precipitation, population, and GDP need to be incorporated. The Coupled Model Intercomparison Project Phase6 (CMIP6) model offers researchers a range of future development scenarios for global CC. These scenarios were created by integrating shared socioeconomic pathways (SSPs) with representative greenhouse gas concentration emission pathways (RCPs). By combining these datasets, researchers can generate comprehensive and realistic projections of LULC changes under future CC scenarios [8,42]. Among the available scenarios, researchers can consider utilizing the SSP1-2.6 (representing a low-coercive socio-economic development scenario), SSP2-4.5 (representing a medium-coercive socio-economic development scenario), and SSP5-8.5 (representing a high-coercive socio-economic development scenario). For further detailed information on these scenarios, visit the website http://esgf-node.llnl.gov/search/CMIP6/, accessed on 12 June 2022. This resource provides comprehensive information and materials related to the CMIP6 project.

The LULC within the ANWC were driven by both CC and various socio-economic factors. In addition, the local land management policies in China exerted a profound influence on LULC patterns. It is important to note that the driving forces behind these land policies also play a significant role in shaping the dynamics of LULC [43]. The observed changes can be attributed to the cumulative impacts of regional development policies, urbanization policies, and ecological protection policies. Notably, China has implemented several national ecological restoration projects in recent years, including the Three-North Shelter Forest Project, the Natural Forest Protection Project, the Grain for Green Project, and the Grazing and Grassland Project. These initiatives have significantly contributed to the enhancement of the ecological environment and the restoration of degraded ecosystems [44]. In light of various CC scenarios, researchers have employed simulation and prediction techniques to generate multiple potential LULC scenarios, considering diverse land management policies. This approach allows for the comparison of potential outcomes under different future environmental conditions. By providing multiple options, it enables decision makers in current land management to assess the potential consequences and make informed choices.

We specifically selected three representative shared socioeconomic pathway (SSP) scenarios from the CMIP6 dataset. These SSP scenarios were combined with potential future development models to formulate the following scenarios: SSP1-2.6 in conjunction with the ecological protection scenario (SSP1-2.6-EP), SSP2-4.5 paired with the natural development scenario (SSP2-4.5-ND), and SSP5-8.5 coupled with the economic growth scenario (SSP5-8.5-EG). The design of these scenarios was guided by specific principles and objectives, which are outlined as follows:

Scenario 1: SSP1-2.6-EP signifies a focus on green sustainable development, with a specific emphasis on ecological protection by the local government. Under this scenario, efforts were made to enhance the preservation of forests, shrubs, medium/high-coverage grasslands, and waters, aiming to attain an interdependent relationship and balanced growth between human engagements and the environment. In line with this objective, strict regulations were implemented to limit the expansion of construction land [45,46]. Based on the aforementioned protection policy, the transfer matrix for the scenario is presented in Table S3. Notably, 72 designated conservation areas and significant open water bodies, including Qilian Mountain National Park, Ningxia Helan Mountain National Nature Reserve, and Xinjiang Altun Mountain National Nature Reserve, were designated as spatial constraints within this framework.

Scenario 2: SSP2-4.5-ND symbolizes a future development scenario that closely resembles the prevailing circumstances. It indicates the continuation of the historical trend in social economy and LULC. To assess this scenario, an initial LULC transfer matrix was computed based on historical data from 1980 to 2020. The resulting transition matrix for this scenario is displayed in Table S4.

Scenario 3: SSP5-8.5-EG embodies a high-speed development paradigm driven by the extensive utilization of fossil fuels. This model was characterized by rapid growth in population and economy, accompanied by a substantial increase in the rate of expansion for construction land [42,47]. Table S5 presents the transition matrix corresponding to the SSP5-8.5-EG scenario.

2.3.5. Future ESs Simulation

Given the crucial role of the ANWC in water and biodiversity conservation, it serves as a significant ecological security barrier in China. However, the region has been faced with severe soil erosion due to both CC and anthropogenic activities. To evaluate the ESs provided by four representative types (CS, WC, SC, and HQ) within these areas, the InVEST software (version 3.9.2) was utilized.

To conduct the assessment, the InVEST model, which can be obtained from the website www.naturalcapitalproject.org, was employed in conjunction with LULC data. This allowed for the evaluation of the ecological functions and values associated with the identified ESs in the ANWC [48]. To gain a comprehensive understanding of the model settings, visit the website https://naturalcapitalproject.stanford.edu/software/invest, accessed on 2 May 2021. This offers detailed information regarding the configuration of the model and the parameters employed in the assessment process.

2.3.6. Trade-Off and Synergy Relationship of ESs

Pearson correlation coefficients (PCCs) were utilized to assess the linear relationship between two variables, namely X and Y, with values ranging from −1 to 1. A value closer to 1 indicates a strong connection, while a value closer to −1 indicates an indirect association. Pearson correlation analysis was conducted to scrutinize the trade-offs and synergies among the four ES functions (CS, WC, SC, and HQ) over time. A higher correlation coefficient value close to 1 indicates a stronger degree of synergy between the variables, while a value closer to −1 suggests a greater degree of trade-off. To perform this analysis, the “corrplot” package in R (install.packages (“corrplot”)) was utilized. This package facilitated the visualization and interpretation of the correlation matrix for further analysis and discussion [49]. To access detailed instructions on visualizing correlation matrices using correlograms, it is recommended to visit the website titled “Visualize correlation matrix using correlogram—Easy Guides—Wiki—STHDA”.

Geographically weighted regression (GWR) is a modeling approach that relaxes the assumptions inherent in conventional regression models. It enables localized regression modeling of parameters at individual locations, accommodating spatial variations and yielding regression outcomes that incorporate spatial dynamics. This technique accounts for spatial heterogeneity by estimating regression results that capture spatially varying relationships. By considering the spatial changes within a dataset, GWR provides a more nuanced understanding of the relationships between variables [50]. The current study utilized the GWR model to scrutinize the territorial correlation among ESs within the designated study area. The GWR model is expressed as follows [51,52]:

$$y_i={\beta }_0\left(u_i,v_i\right)+{\displaystyle \sum _{j=1}^k{\beta }_k\left(u_i,v_i\right)}{x}_{ij}+{\epsilon }_i$$

(4)

$$\widehat{\beta }\left(u_i,v_i\right)={\left[X^{T}W\left(u_i,v_i\right)X\right]}^{-1}{X}^{T}W\left(u_i,v_i\right)Y$$

(5)

$$\left\{\begin{array}{cc}{w}_{ij}={\left[1-{\left(d_{ij}/h\right)}^2\right]}^2,& d_{ij}<h\\ w_{ij}=0,& d_{ij}\ge h\end{array}\right.$$

(6)

$$AICc=-n\ln \left({\sigma }^{\prime }\right)+n\ln \left(2\pi \right)+n\left\{\left[n+tr\left(S\right)\right]/\left[n-2-tr\left(S\right)\right]\right\}$$

(7)

where the dependent variable is represented by y, while the independent variable is denoted by x. The coordinates of the i sample point are given by (u_i, v_i). The intercept term is denoted as β₀ (u_i, v_i), and the k regression coefficient on the i sample point, which varies with geographical location, is represented as β_k (u_i, v_i). The regression residual was denoted as ε_i. The regression coefficient in the GWR model was estimated using the weighted least squares method. Here, $\widehat{\beta }\left(u_i,v_i\right)$ refers to the spatial weight diagonal matrix, X represents the matrix of independent variables, and Y represents the vector of dependent variables. The spatial weight was determined using the Bi-square function, where W_ij denotes the weight assigned to known point j for estimating the unknown point i. The Euclidean distance between points i and j is represented by d_ij, and h signifies the bandwidth. The bandwidth was selected based on the minimum AICc (corrected Akaike Information Criterion). Additionally, n denotes the number of data points, σ′represents the error term estimator, and tr (S) indicates the trace of hat matrix S.

2.3.7. ESs Scale Setting, Spatial Heterogeneity Detection and Attribution

To analyze the trade-off and synergy effects of ESs across various scales, watershed, county, and grid scale units were employed. Moreover, the optimal parameter geodetector (OPGD) model was utilized to uncover the spatial heterogeneity of these effects. Additionally, the evolution of ESs was analyzed using a structural equation model (SEM). For detailed information on the methodology, refer to Supplementary Materials Sections S3.3–S3.5, where specific methodological details are provided.

3. Results

3.1. LULC from 1980 to 2100

The ANWC exhibits significant spatial heterogeneity in the distribution of LULC (Figure 3, Figures S8 and S9). Within the study area, the Tarim Basin, the core region of the Junggar Basin, and the northeastern part of the Hexi Corridor are characterized by the presence of unused land, such as sandy land and Gobi. This particular region also represents the primary distribution area of the four major deserts within the ANWC. Woodland areas were primarily found along the northern and southern slopes of the Tianshan Mountains, Altai Mountains, and Kunlun Mountains at elevations ranging from 1800 to 2900 m. Additionally, woodland was observed on the northern slope of the Qilian Mountains at elevations between 1800 and 4200 m. Grassland regions can be found in the Tianshan Mountains, Altai Mountains, and Qilian Mountains at elevations ranging from 1050 to 1800 m. In certain areas of the Qinghai-Tibet Plateau, grasslands were also present on the southern slope of the Kunlun Mountains at elevations between 2900 and 4200 m. Waters were predominantly concentrated in the oasis areas of piedmont plains. Notably, permanent glaciers can be found above 4200 m in mountainous regions, contributing to the overall water resources. Construction land was mostly situated in transitional oasis plain areas of the Tarim Basin, peripheries of the Junggar Basin, and mountainous regions. Among these areas, the largest extent of construction land was observed in the Tianshan Mountains, Kunlun Mountains, and the northern foot of the Qilian Mountains, accounting for approximately 97% of the total construction land area within the ANWC.

During the period from 1980 to 2020, the primary land-use categories in the ANWC included sandy land, Gobi, bare rock texture, and low-coverage grassland. These land types collectively accounted for approximately 73.71% to 74.49% of the total area, which amounts to approximately 15.73 × 10⁵ km². The simulated land-use pattern of different scenarios in the ANWC in 2100 still continues this structure. In the past 40 years, paddy fields, natural forest land, and permanent glacier decreased, while dry land, high-coverage grassland, and urban construction land increased (Figures S10 and S11). The land transfer area reached 57.1 × 10⁵ km², and the comprehensive land-use dynamic degree was 19.47% (Figure S11). The land-use types that exhibited the most rapid growth in terms of dynamic degree were plantation land and urban construction land, with growth rates of 2.66% and 2.35%, respectively. Furthermore, under the SSP1-2.6-EP, SSP2-4.5-ND, and SSP5-8.5-EG scenarios, the comprehensive land-use dynamic degrees were projected to reach 2.58%, 4.08%, and 4.74%, respectively. Moreover, the extent of land transfer was anticipated to reach 14.28 × 10⁵ km², 17.36 × 10⁵ km², and 33.67 × 10⁵ km², respectively, under these scenarios (Figures S10 and S11).

3.2. Temporal and Spatial Variations of ESs

Between 1980 and 2100, the spatial configuration of four land-use types, namely CS, HQ, WC, and SC, exhibited similarities within the ANWC. The functions of CS, HQ, WC, and SC in mountainous regions surpassed those in non-mountainous areas. Specifically, the areas with higher functional values were concentrated in the Tianshan Mountains, Altai Mountains, Kunlun Mountains, and Qilian Mountains (Figures S12–S15). Compared with 1980, the four ESs in the ANWC showed two types of evolution trends. Among them, CS and HQ decreased significantly in 2020 and 2100, and SC and WC increased significantly in 2020 and 2100 (Figure 4). During the 40-year period from the past, notable changes have occurred in the quantities of WC, SC, CS, and HQ within the ANWC. Specifically, WC has experienced an increase of 11.36 × 10⁸ m³, SC has witnessed a substantial change of 1735.25 × 10⁸ t, CS has decreased by 1.29 × 10⁸ t, and HQ has shown a decrease of 0.009. Additionally, it was projected that the overall status of the four ESs within the ANWC will generally exhibit improvement by the year 2100. Under the SSP1-2.6-EP scenario, these improvements were expected to be slightly higher compared with the other two scenarios.

3.3. Trade-Off and Synergy of ESs

Throughout the time span from 1980 to 2100, the ESs of CS, HQ, WC, and SC within the ANWC exhibited a synergistic relationship. These ESs functioned synergistically, enhancing their combined effects on the overall ecosystem (Figure 5). The strongest synergistic relationships were CS and HQ, and the weakest were CS and WC. The synergistic relationship of the four ESs decreased with the interannual variation, and the synergistic relationship under the SSP1-2.6-EP scenario was higher than that under other scenarios.

Over the course of the previous four decades, the spatial arrangement of the interrelationship among the ESs of CS, HQ, WC, and SC within the ANWC displayed a synergistic pattern (Figures S16–S18). The distinctions in the spatial distribution patterns concerning the correlations among the four ESs at various scales, such as watershed, county, and grid scales, primarily manifest in mountainous regions. At the watershed scale, the correlations between the four ESs tend to be relatively low. The areas exhibiting a synergistic relationship were predominantly concentrated in the Tianshan, Altai, Kunlun, and Qilian Mountains. On the other hand, regions around the Taklimakan, Badain Jaran, Gurbantunggut, and Tengger deserts were characterized by a trade-off relationship. It is worth emphasizing that the Tianshan and Qilian Mountains were particularly significant locations with a high concentration of strong trade-offs among the four ESs. Under the SSP1-2.6-EP, SSP2-4.5-ND, and SSP5-8.5-EG scenarios, the spatial distribution patterns of the four ESs (CS, HQ, WC, and SC) within the ANWC demonstrate a continuation of historical trends (Figures S19–S27). Still, the distribution area of the synergistic relationship was shrinking; in particular, the strong trade-offs of the four ESs of the Tianshan Mountains and Qilian Mountains were weakening.

3.4. Spatial Heterogeneity Detection of ESs

The top five main factors detected were X1 drought index (q = 0.378), X15 average annual precipitation (q = 0.375), X12 economic density (q = 0.095), X11 vegetation coverage (q = 0.262), and X22 soil bulk density (q = 0.077) (Figure 6a). Furthermore, the interactions between the aforementioned driving factors were further explored (Figure 6a). The findings additionally indicate that the combined influence of the aforementioned driving factors had a more significant impact compared to individual factors alone. This signifies that the spatial variations in ESs within the ANWC were not solely attributed to a sole determinant but rather result from the collective interplay of diverse essential factors (Figure 6b).

Furthermore, a structural equation model (SEM) was employed to provide further evidence regarding the contributions of CC, socio-economic population (SEP), eco-environmental background (EB), and LULC in influencing alterations within ESs processes (Figure 7). Specifically, CC in the ANWC has an important impact on ESs. In particular, the average annual precipitation has a positive driving effect on ESs (0.83), and the average annual temperature (−0.61), drought index (−0.90), and potential evapotranspiration (−0.81) have a negative driving effect. In addition, EB was also an important factor driving ESs in the ANWC (0.89). The driving effect of SEP and LULC was smaller than that of the previous two factors, indicating that CC and EB play the greatest role in ESs in the ANWC.

4. Discussion

4.1. Causes of LULC

This study made improvements in the classification of LULC within the ANWC compared to previous research. We evaluated the LULC in the ANWC over the past 40 years. We found that the area of paddy fields, natural forest land, and permanent glaciers in 24 types of LULC decreased significantly, while the area of dry land, high-coverage grassland, and urban construction land increased significantly (Figures S10 and S11). In the past 40 years, the ANWC has experienced a warm and humid process. In particular, there had been a notable increase in the annual average temperature and precipitation within certain mountainous regions, including the Tianshan, Altai, and Qilian Mountains. The maximum observed increase in temperature was 0.044 °C·yr⁻¹, while the maximum increase in precipitation reached 2.70 mm·yr⁻¹ (Figure S3). Recent studies have shown that glaciers in China have generally experienced rapid melting and shrinkage under global warming [19,51,52]. Our study provides evidence of the reduction in permanent glacier coverage within the ANWC. Our simulations for the year 2100 revealed a persistent trend of warm and humid conditions, with the most significant glacier ablation observed under the SSP5-8.5-EG scenario. The changes in glaciers has lead to an augmented discharge in mountainous watersheds, which subsequently facilitates the process of ‘greening’ the Gobi and desert regions. This phenomenon had the potential to enhance water availability for agricultural cultivation and ecological rehabilitation, ultimately supporting the expansion of cultivated areas and mitigating the extent of sandy land. From the land-use transfer, we found that, in the past 40 years, more Gobi has turned to low-coverage grassland, bare rock texture has turned to medium-coverage grassland, and desert has turned to medium-coverage grassland (Figures S10 and S11). Part of the increase in these grasslands is the conversion of mountain bare rock texture to alpine grassland under glacier ablation and part of the expansion of more runoff oasis areas downstream. Furthermore, it is important to highlight that the expansion of cultivated and construction land predominantly takes place within oasis regions characterized by densely populated areas, well-developed transportation infrastructures, abundant water resources, and flat topography. This expansion has been facilitated by the government’s implementation of a series of ecological initiatives (i.e., the protection of grasslands, preservation of natural forests, and safeguarding of cultivated areas) [53,54]. Moreover, the adoption of water-saving agriculture practices has led to the elimination of more water-intensive paddy fields. Simultaneously, the migration and concentration of population has resulted in an increase in urban land area and a decrease in rural residential land.

4.2. Driving Factors of ESs Change

CC, SEP, EB, and LULC are important factors in the change of ESs in the study area (Figure S28). Our research findings indicate that CC and EB play a significant role in driving the alterations observed in ESs within the ANWC. In essence, it can be concluded that natural environmental factors primarily drive the changes in ESs within the ANWC. Specifically, we found that the average annual precipitation in CC plays an essential positive driving role in ESs in the ANWC. Some studies have also shown that the warm and humid climate in the ANWC, especially the increase in precipitation in some areas, increases vegetation coverage and is conducive to the improvement of ESs (HQ, CS, SC, and WC) [55]. In the meantime, it was observed that the drought index, potential evaporation, and annual average temperature negatively influence CC, consequently leading to indirect negative effects on ESs. The aforementioned findings are in accordance with the research outcomes presented by Dou et al. (2020) [56] in select regions of Inner Mongolia, which are also encompassed within the geographical scope of the ANWC. Furthermore, EB is also an important factor that dominates the change of ESs in the study area. Vegetation coverage, topography, soil, and other factors explain more than 50% of ESs. In particular, vegetation coverage has had the greatest positive direct impact on ESs [57,58]. On the contrary, altitude has had an indirect negative impact on ESs. SEP and LULC are also the other two main factors driving the ANWC. We find that both of them had a negative effect on ESs. This concurs with the prevailing body of research indicating that LULC exerts an adverse influence on ESs as a consequence of persistent demographic expansion, economic progress, and the proliferation of both agricultural land and urban development [3,4]. In addition, it is worth noting that these factors have served as the primary drivers of spatial variations in ESs within the ANWC over the past four decades (Figure 6).

4.3. Exploring the Significance of ESs in Ecological Protection and Management

This study emphasizes the management and restoration of diverse mountain–oasis–desert systems in the ANWC under the influence of CC and human disturbances. In the face of mountain ecosystems, we consider the role of mountains in regulating WC and HQ in the ANWC. It was found that the area of permanent glaciers decreased significantly, and this may supply more water downstream in a short time, which also means that WC in some areas, especially in mountainous valleys, will be improved [32]. However, the loss of glaciers is irreversible and difficult to estimate [59]. It is necessary to achieve the goal of “carbon neutrality and carbon peak” in the world to alleviate climate warming, especially to strengthen greenhouse gas emission reduction. At present, studies have also shown that engineering measures, such as laying sunshade blankets on permanent glaciers, can delay the ablation of glaciers [59]. Furthermore, the mountain ecosystem maintains a higher level of HQ, providing wildlife habitats for nationally protected species such as Panthera uncia, Lynx, Snow Lotus Herb (Saussurea involucrata (Kar. & Kir.) Sch. Bip.), and Tangut Rhodiola (Rhodiola tangutica (Maxim.) S. H. Fu) [60]. Nonetheless, the activities of mining, overgrazing, and tourism development in mountainous regions has caused severe degradation to the habitats of various fauna and flora. Consequently, there is an urgent need to enhance environmental management strategies in order to ameliorate the vulnerable ecosystems found in these mountains [47].

The local government in the oasis region should enhance the management of water resources, take proactive measures to combat desertification, accurately assess agricultural water consumption, and refrain from encroaching upon ecologically vital water sources [32]. In addition, controlling the size of oasis cities is also important to avoid unnecessary waste of water resources. Urban management measures should improve and restore the quality of green spaces.

In the context of desert ecosystems, we acknowledge that ecological restoration in the ANWC entails higher temporal and economic costs. Consequently, we recommend prioritizing measures such as enclosure protection, minimizing human disturbance, and allocating adequate resources for the ecological restoration of nature reserves, key ecological function areas, and national parks.

4.4. Advantages and Limits

Based on the PLUS model, SD model, InVEST model, GWR model, Pearson correlation coefficient (PCCs), OPGD model, and structural equation model, we evaluated and attributed the ESs in different climate and socio-economic scenarios in the ANWC from 1980 to 2100. In contrast to prior research, this study has enhanced the precision of land-use classification, enabling a more comprehensive monitoring of the nuanced alterations in LULC types within the examined region. This refined classification methodology offers enhanced insights into the intricate dynamics of ESs evolution, facilitating a more detailed analysis of their changes. In addition, we provide more choices for the watershed, administrative region, and grid scale establishing of ESs trade-offs and synergies, which provides more specific solutions for decision makers. However, recent studies have focused on the scale effect of ESs [52,53]. How to choose the appropriate scale and what the relationship is between scales is still something we need to study in depth. Moreover, we were limited by more high-resolution and field-measured data and need to consider the changes in surface runoff and groundwater level in the selection of simulation covariates [32]. At the same time, we only selected four ESs: CS, HQ, WC, and SC. In further research, we need to detect more changes in ESs. We also note that recent studies have begun focusing on the supply and demand relationship, flow direction, and flow of ESs. Therefore, we need to further study such issues in ANWC.

5. Conclusions

During the period from 1980 to 2100, the ANWC has been undergoing a notable shift towards warmer and more humid conditions. There has been a substantial reduction in the permanent snow glacier area, particularly under the SSP5-8.5-EG scenario, which had the most pronounced impact. Consequently, these glacier transformations have led to amplified runoff in mountainous basins, subsequently fostering an accelerated process of vegetation expansion in downstream Gobi and desert areas, commonly referred to as “greening”. The degrees of change in land-use dynamics were considerable, reaching 2.58%, 4.08%, and 4.74% under scenarios SSP1-2.6-EP, SSP2-4.5-ND, and SSP5-8.5-EG, respectively. Additionally, there was a transfer area of 14.28 × 10⁵ km², 17.36 × 10⁵ km², and 33.67 × 10⁵ km² under the three scenarios, respectively, indicating substantial shifts in land uses across the ANWC. At the same time, the alterations in LULC resulting from CC intensified the disparities in the four ESs of the ANWC. Specifically, the quantities of CS and HQ experienced significant declines by the year 2100, while SC and WC witnessed substantial increases during the same period. The changes in WC, SC, CS, and HQ were measured at 11.36 × 10⁸ m³, 1735.25 × 10⁸ t, −1.29 × 10⁸ t, and −0.009, respectively. The four ESs of CS, HQ, WC, SC in the ANWC exhibit a synergistic relationship. The spatial distribution of these ESs displayed heterogeneity, resulting in the strongest synergy between CS and HQ, while the weakest synergy was observed between CS and WC. The distribution variations of ES synergies were amplified at watershed, county, and grid scales in mountainous areas, with more significant detection differentiation at the grid scale. Furthermore, the detection of spatial heterogeneity in the four Ess was primarily attributed to various factors: drought index (q = 0.378), annual average precipitation (q = 0.375), economic density (q = 0.095), vegetation coverage (q = 0.262), and soil bulk density (q = 0.077).

The findings of our study underscored the critical role of CC and EB in driving the fluctuations witnessed in the ESs within the ANWC. In future studies, there is a need to further enhance the precision of land-use classification, ensuring a more nuanced identification of land-use patterns and facilitating the detection of additional changes in ESs. Moreover, it was crucial to investigate the intricate interplay between the supply and demand of ESs, as well as conduct an in-depth analysis of the flow and distribution of these vital services. Such endeavors will contribute to a deeper understanding of the complexities involved in ESs dynamics and provide valuable insights for sustainable land management and conservation efforts.