Trade-O ﬀ and Synergy Relationships and Spatial Bundle Analysis of Ecosystem Services in the Qilian Mountains

.


Introduction
Ecosystem services (ES) refer to the benefits and welfare that ecosystems provide to humans [1]. In 2005, the Millennium Ecosystem Assessment report, released by the United Nations, highlighted that 60% of the world's ecosystems are in a degraded state [2]. Consequently, ES have become a significant focus of research in ecology and related disciplines [3,4]. Given the complexity and diversity of ES, the evident spatial heterogeneity in their distribution, and the fact that humans utilize ES to improve their own well-being, there exit complex linear and nonlinear relationships among ES [5]. Previous studies have categorized the relationship between ecosystems as trade-offs, synergies, and non-response [6]. In recent years, exploring the trade-off and synergy between ES has become a core issue in ES research [7][8][9][10]. If the evaluation of services provided by land is insufficient, then the one-sided pursuit of service while ignoring another type of service will lead to trade-offs, eventually resulting in the degradation of the ecological environment on a large scale [11]. Differences in climate, upper soil, and vegetation cover can make landuse types spatially heterogeneous, thus affecting ES [12]. Water resources play a crucial role in connecting and controlling different ecosystems while ensuring a stable water supply is essential for maintaining the stability of ES in a region [13]. Climate change significantly impacts the global water cycle, leading to changes in precipitation, evapotranspiration, and runoff patterns, altering the spatial and temporal distribution of water resources, and increasing the frequency of extreme disasters such as droughts and floods [14], and thus exacerbating the conflict between supply and demand of ES. Scholars from many countries, at regional [15][16][17][18], watershed [19][20][21], and other scales, have studied the tradeoffs and synergies between ES and spatial-temporal differentiation characteristics [12,[16][17][18][21][22][23]. However, most previous research focused only on several important services [15,22], with limited exploration of the relationships between provisioning, regulating, habitat, and cultural services. Additionally, there is a lack of research on the spatial composition of ES.
In recent years, the concept of ES bundles (ESB) has emerged to enhance the management of multifunctional landscapes [24]. ESB involves clustering ES to identify dominant services in a specific region, thereby improving ecosystem management by promoting synergies and reducing trade-offs. ESB represents combinations of services within a particular time and space, offering insights into the interdependent relationships among ES. Based on the similarity principle of clustering analysis, it can spatially quantify the composition and structure of ES and is a valuable tool for the identification of several important services in a study area [25,26]. Common methods used to identify service bundles include principal component analysis [26], self-organizing maps [25,27,28], spectral clustering, and k-means clustering analysis [29]. Current research on service bundles is mostly focused on different scales, such as countries [30], regions [31], and basins [32]. While these scales facilitate zoning management, they may not be optimal for studying services based on ecological processes. For fine-scale land management and research on services, it is beneficial to partition the study area into appropriate grids. However, there is currently a lack of research on the dynamic partitioning of ESB at the grid scale [33]. Therefore, using service bundles to study the spatial dynamic evolution of ES on a grid scale is of great significance for maintaining the sustainable development of ecosystems and providing targeted scientific decision-making, considering local economic development.
In the global terrestrial ecosystem, mountainous areas form approximately onefourth of the total surface area; 12% of the global population lives in mountainous areas, and these areas are hot spots for more than 50% of the world's plant and animal species [34]. The Qilian Mountains are an important water-producing area and a priority biodiversity conservation area in China. In addition, the ecological environment of these mountains is fragile and has drawn the attention of the societal community. The ecological protection of the Qilian Mountains is related to the ecological security and sustainable social development of Western China.
In view of this, we selected the Qilian Mountains as a representative case to investigate spatiotemporal changes and hotspots of ES, trade-offs, and synergies among ES and identify four distinct types of ecosystem service clusters. By integrating these findings with a decision support system, our research provides a valuable reference for future decision-making optimization and institutional reform of Qilian Mountain.

Study Area
The Qilian Mountain range is located in the arid and semi-arid inland region of northwest China at the intersection of the Qinghai-Tibet Plateau, Inner Mongolia Plateau, and Loess Plateau. It spans the Gansu and Qinghai provinces and is the largest marginal mountain system in the northeastern Qinghai-Tibet Plateau. Tuanjie Peak, the highest peak of Mountain Shulenan, is 5808 m above sea level [35,36]. Starting from Wushaoling in the east and ending at the Dangjinshan Estuary in the west, the Qilian Mountains are an important water source conservation area for three inland rivers, namely the Shiyang River, Heihe River, and Shule River. This mountain range has a variety of land use and vegetation cover types, including forests, shrubs, grasslands, wetlands, deserts, and glaciers. It is a priority area for biodiversity conservation and an important eco-security barrier in northwest China. It has the important ecological function of preventing the Badain Jaran, Tengger, and Kumtag Deserts from encroaching to the south, maintaining the stability of the oasis in the Hexi Corridor, and safeguarding the supply of runoff from the Yellow River and inland rivers [37]. Figure 1 provides an overview of the study area.

Selection of ES
ES encompasses a broad range of topics. Although a great deal of work has been conducted by domestic and international scholars in the study of ES, there are still many controversies regarding the definition and classification of ES [1,2,[38][39][40][41]. In this study, our goal is to provide recommendations for ecological management and promote the sustainable development of ES for the benefit of human well-being. Therefore, we finally chose the classification system of The Economics of Ecosystems and Biodiversity. This system emphasizes the integration of ES into policy-making processes, provides a clearer distinction between ES and benefits, and has gained widespread international recognition [41][42][43][44][45]. To support our research, we conducted field studies on the ecological characteristics of the Qilian Mountains and conducted a comprehensive review of the relevant literature on ecosystem assessment [46][47][48][49][50][51][52]. We carefully considered data accessibility and the feasibility of calculation methods during the selection process, resulting in the identification of nine important ES for our study (Table 1). The Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model is a freely available, open-source software that simulates changes in ES under different land cover scenarios. It serves as a valuable tool for policymakers, providing a scientific foundation to assess the benefits and impacts of human activities. The InVEST model generates evaluation results that can be visualized cartographically, allowing for the easy depiction of spatial heterogeneity and dynamic changes in ES. In our study, we utilized the InVEST model to evaluate water resource supply, water conservation, soil conservation, water purification, and habitat quality in the Qilian Mountains over 19 years (2000-2018). To assess water resource supply and water conservation, we employed the water yield module of the InVEST model. For soil conservation, water purification, and habitat quality, we utilized the Sediment delivery ratio module, Nutrient delivery ratio module, and Habitat quality module, respectively. The specific modules and parameters of the InVEST model involved in this study can be found in Appendix A.
The water yield module in the InVEST model is based on the Budyko water-heat coupling equilibrium hypothesis [53]. It utilizes various parameters, such as precipitation, surface evapotranspiration, maximum root depth of plants in the soil, available water content of plants, and land use type, to calculate water yield. The water balance equation is employed to estimate water conservation [54][55][56], taking into account precipitation, net surface flow, and evapotranspiration as input data. To calculate soil retention, we used the Sediment delivery ratio module of the InVEST model coupled with the drainage basin wind erosion model developed by Dong et al. (1998) [57]. We used nitrogen export to reflect water purification. The principle of the nutrient delivery ratio module of the InVEST model is to calculate the long-term steady-state flow of nitrogen in the entire watershed based on the ability of different ecosystem types to transfer nitrogen into runoff and to reflect the water quality in the watershed according to nitrogen export in the water body. Higher nitrogen content indicates more severe water pollution or eutrophication in the watershed and indicates weaker water purification capacity of the ecosystem. For evaluating habitat quality, the habitat quality module of the InVEST model was applied. Research has demonstrated that the habitat quality simulated by the habitat quality module of the InVEST model can be a good reflection of the current state of habitat services [58].

Calculation Methods and Parameters of Other ES
The food supply service in the study area was assessed based on the output values of agriculture, animal husbandry, and fisheries. The normalized difference vegetation index (NDVI) has a significant linear relationship with crop production [59], enabling the spatialization of agriculture and animal husbandry output values using NDVI data [60]. The fishery output value is distributed equally across the area of the aquatic ecosystem. Raw material supply includes wood and non-wood products provided by forests, fresh grass, and refined feed provided by grasslands. The output value of forests and grass yield were used as indicators for quantifying the raw material supply provided by these ecosystems [61,62]. The remote sensing model, particularly the CASA model, was used to calculate the net primary productivity (NPP) of vegetation [62], and grassland yield was calculated according to the functional relationship between the dry weight of organic matter produced by the grassland ecosystem and the grass yield per unit area [63]. Climate regulation was assessed using the net primary productivity of vegetation [64,65] using the CASA [66,67]. For entertainment tourism services, four evaluation criteria were established [46] based on the landscape characteristics, natural attributes, and potential human activities of different ecosystem types in the Qilian Mountain ecosystem, including accessibility, tourism potential, unique natural landscape or cultural sites, and tourism comfort level. Finally, the results of the four evaluation criteria were weighted to obtain the evaluation results of entertainment tourism services. Detailed calculation methods and parameter settings are provided in Appendix A.

Data Requirement and Preparation
The required data for this study are presented in Table 2. The land use data for the years 2000, 2010, and 2018 used in this research were obtained from the Data Center for Resources and Environmental Sciences and the Chinese Academy of Sciences. These land use types were verified through field investigation, resulting in an overall accuracy of 85%. Digital elevation data were sourced from the China Geospatial Data Cloud Platform. Meteorological data, including temperature, precipitation, sunshine duration, and average wind speed, were obtained from the National Meteorological Information Center of China. By using elevation as a covariate, ANUSPLIN software was employed for spatial continuous interpolation of the meteorological data [68]. Soil data were acquired from the Harmonized World Soil Database version 1.2, developed by the Food and Agriculture Organization of the United Nations and the International Institute for Applied Systems of Vienna. The soil data for China originated from the 1:100,000 soil data provided by the Nanjing Soil Institute of the Second National Land Survey, and the data were interpolated in ArcGIS to a resolution of 30 m. Socio-economic data were sourced from the Statistical Yearbooks and Water Resources Bulletin of the research published by the National Bureau of Statistics of China, the Bureau of Statistics of Gansu Province, and the Bureau of Statistics of Qinghai Province. Additionally, a questionnaire survey was conducted to gather data on forage prices from herders in the study area. A total of 146 questionnaires were distributed, and after sorting, screening, and elimination, 133 valid questionnaires were obtained, resulting in an effective rate of 91.10%. These questionnaires provided data on feedstock prices and other relevant economic information for the assessment of raw material supply services.

Exploratory Spatial Data Analysis Methods
Exploratory spatial data analysis is a method used to thoroughly analyze and visualize data with a spatial distribution [69]. Its purpose is to identify geographic characteristics within the dataset and reveal spatial distribution patterns and evolutionary trends [70]. In this study, we extensively reviewed the scholarly literature on the construction of grids to study the spatial differentiation of ES [12,[69][70][71][72]. We considered various grid sizes, including 30 m × 30 m, 500 m × 500 m, 1 km × 1 km, and 5 km × 5 km as evaluation cells. Considering the topographic features and size of the study area, and the practicality of implementation of territorial spatial planning and future environmental protection policies, we finally chose a 1 km × 1 km grid as the evaluation cell for this study.
In this study, we used global and local spatial autocorrelation analyses to analyze the spatial association patterns of ES and measured and assessed the convergence or heterogeneity of different ES. ArcGIS software tools were used to create a 1 km × 1 km grid in the study area, extract different ES, and realize the different ES of microgrid scale refactoring.
Global spatial autocorrelation analysis is used to study the similarity and degree of correlation of regional grid attributes in the global space, which includes global Moran's I and Geary's C methods [73]. In this study, we used global Moran's I to analyze the global spatial autocorrelation of the study area with the following calculation formula: where n is the total number of grid cells in the study area; , are the measurement values of an attribute feature in regions and , respectively; ̅ is the average value of all observations of a certain attribute feature in the study area; is the standardized space-weight matrix; and is the variance. Local spatial autocorrelation analysis can measure the spatial difference or degree of correlation between each grid and the surrounding grids [73]. It provides insights into the spatial patterns and relationships within a study area. The main research methods include the Moran's scatter plot, local Moran's I, and statistics [69,73]. This study selected the local Moran's I to explore the spatial connections of ES and the similarities and differences of local spaces in the study area.

Analysis of Trade-Off and Synergy Relationships
In ecosystems, optimizing multiple ecosystem services simultaneously at the same spatial or temporal scale can be challenging and often results in tradeoff relationships [74]. In order to study the trade-off and synergistic relations among various ES in the Qilian Mountain area, this study used the Create Fishnet tool in ArcGIS to randomly create 1000 points in the years 2000, 2010, and 2018 and extracted different ES values corresponding to each point. Spearman correlation analysis was performed on the extraction values, and the results were validated using the t-test. If the correlation coefficient between ES was positive and the p-value was less than 0.05, the two services were considered synergistic; otherwise, they were considered trade-offs [75].

K-Means Clustering Analysis
The K-means clustering analysis adopted in this study is an unsupervised clustering method [76]. By providing the dataset and specifying the number of clusters (K), this algorithm determined K clusters that achieved the highest similarity among objects within the same cluster and the lowest similarity between objects in different clusters [77]. The clustering results obtained by the K-means clustering analysis have a clear structure and are widely used in geographical research [78][79][80].

Model Results Validation Methods
Reliable ES assessment results obtained from model simulations are crucial for subsequent studies. Hence, this study conducted a verification of the simulation results of the four models, including water resource supply, soil conservation, nitrogen export, and climate regulation (Table 3). Due to limited hydrological stations in the study area, only the measured data from some basins could be obtained. Therefore, runoff data were obtained from the Changmabao, Yingluoxia, and Dangchengwan hydrological stations to validate the water yield results. Sediment data from the Dangchengwan and Jiutiaoling hydrological stations were used for verifying the soil conservation results. River water quality data from the Yingluoxia, Jiutiaoling, and Zamusi hydrological stations were selected to verify the water quality purification results. In this study, grassland sample data were collected for different types of grassland in the study area, with a total of 115 grassland sample data collected. Based on the empirical formula, the NPP calculated from the measured data was compared with the NPP results simulated by the CASA model to verify the results of climate regulation services [63][64][65][66][67]. The watersheds and sampling sites selected for this study are shown in Figure 2.

Verification of Model Simulation Results
Verification of model simulation showed the following: (1) Regarding water resource supply (Figure 3a), when the Zhang coefficients of the InVEST model in 2000, 2010, and 2018 were 2.5, 2.5, and 3.5, the overall accuracy of the model was 86.0%, 89.2%, and 90.2%, respectively. Further, the water yield simulation results showed good credibility, and the obtained results were realistic and accurate. (2) Regarding soil conservation (Figure 3b), the accuracy of the simulation results of soil conservation in 2010 and 2018 was 60.0% and 95.0%, and the overall accuracy was 77.5%, which proved that the InVEST model met the requirements for evaluating soil conservation services in the study area. (3) Regarding water purification (Figure 3c), during the study period, the nitrogen export results simulated by the nutrient delivery ratio module of the InVEST model were all within the concentration range of water quality observations at the hydrological station, proving that the model meets the research requirements. (4) Regarding climate regulation, the sampling yield of 2021 grassland quadrats was compared with the simulation results of the 2021 CASA model (Figure 3d), and the overall simulation accuracy of the CASA model was 63.6%. The R2 value of the linear fitting of grass yield and NPP was 0.7458, proving that the model-simulated NPP results and the actual sampling data have a good degree of fit. Therefore, the CASA model is applicable to the study of climate regulation services. These verification results provide confidence in the reliability and accuracy of the simulation outcomes obtained from the four models used in this study.

Changes in Ecosystem Areas
This study constructed a transfer matrix for the different ecosystems from 2000 to 2018 (Table 4) and analyzed the area served by each ecosystem (Table 5). In 2018, the dominant ecosystem type in the Qilian Mountains was grassland, covering 44.45% of the total land area. Desert, shrubland, and forest ecosystems accounted for 36.89%, 5.33%, and 3.31% of the total land area, respectively. The human settlement ecosystem was the smallest, accounting for only 0.2-0.3% of the total land area. Over the study period, significant changes in ecosystem areas were observed in the Qilian Mountains. Notably, there was a decrease in desert and glacier areas. The areas of grassland, forest, farmland, water, wetland, and human settlement ecosystems increased, while the area of the shrubland ecosystem remained relatively constant.
During the study period, desert ecosystems had the largest net transfer outflows, with a total transfer of 4099.35 km 2 , primarily to the grassland ecosystems (Table 5). Between 2000 and 2010, there was a significant transfer of the desert ecosystem to the grassland ecosystem, totaling 3858.96 km 2 . The glacier ecosystem experienced a net outflow of 215.84 km 2 , most of which transferred to the desert and grassland ecosystems, with limited transfers to the shrubland and wetland ecosystems. The grassland ecosystem showed the largest net transfer of 3568.51 km 2 from the desert ecosystem. While the area of the grassland ecosystem increased significantly, a small part of the grassland ecosystem turned into human settlements, farmland, forest, and water ecosystems. The main source of the increase in the forest ecosystem was the transfer of the desert and grassland ecosystems, which was 68.92 km 2 and 18.05 km 2 , respectively. The area of the water ecosystem increased from the desert, grassland, and wetland ecosystems, with an inflow of 85.21 km 2 , 55.27 km 2 , and 37.80 km 2 , respectively. Except for the glacier and water ecosystems, all other ecosystems in the Qilian Mountains underwent a transformation into human settlement ecosystems. Among them, the grassland ecosystem had the largest transformation area (56.26 km 2 ). The grassland ecosystem was the main source of the increase in the area of the farmland ecosystem, with a transferred area of 229.78 km 2 . The area of the water ecosystem increased slightly, and 267 km 2 of the desert ecosystem was transferred to it.

Identification of Ecosystem Services and Hotspots in the Qilian Mountains
During the study period, the Qilian Mountains exhibited clear increases in the food supply, raw material supply, water resource supply, water conservation, climate regulation, and entertainment tourism services (Table 6). Food supply and raw material supply in the area increased by 1038.83 and 448.21 Yuan·ha −1 , respectively. Water resource supply and water conservation increased by 55.45 mm and 7.80 mm in 2018 compared to 2000. Climate regulation improved by 0.60 tc·ha −1 , indicating an enhancement in the ecosystem's capacity to regulate climate. Entertainment tourism increased to 26.76. Nitrogen export and habitat quality remained stable, and soil conservation showed a fluctuating growth trend, increasing by 40.01 t·ha −1 over the study period.
During the study period, various ecosystem services in the Qilian Mountains exhibited distinct characteristics ( Figure 4).
During the study period, the farmland ecosystem had the largest food supply capacity, with a mean capacity of 20,396.62 yuan·ha −1 in 2018, followed by the grassland ecosystem, with a capacity of 2178.36 yuan·ha −1 (Figure 4a).
In 2018, the forest ecosystem had the strongest raw material supply capacity, with a mean annual raw material supply of 3773.70 yuan·ha −1 , followed by the wetland ecosystem with 3178.51 yuan·ha −1 , and the grassland ecosystem with 2707.44 yuan·ha −1 ( Figure  4b). The raw material supply of the forest ecosystem increased significantly to 82.48% from 2000 to 2018. The increase in the raw material supply of the wetland and grassland ecosystems was 10.03% and 30.01%, respectively.
The mean annual water yield of the study area was 146.67 mm, of which the glacier ecosystem had the largest annual mean water yield of 244.29 mm·ha −1 , and the desert ecosystem had the second largest annual mean water yield of 228.47 mm·ha −1 (Figure 4c). The mean water conservation capacity of the forest ecosystem was the highest at 20.71 mm·ha −1 in 2018 (Figure 4d).
The mean climate regulation service in 2018 was 3.05 tc·ha −1 , an increase of 0.59 tc·ha −1 compared with that in 2000 ( Figure 4e). The forest ecosystem had the largest annual mean climate regulation service (4.94 tc·ha −1 ).
The mean soil conservation capacity in the study area was 1026.04 t·ha −1 , and the three ecosystems with the largest soil conservation capacities were the shrubland, forest, and desert ecosystems. Perennial mean soil conservation was 1668. 93  During the study period, the glacier ecosystem had the lowest nitrogen export, and the farmland ecosystem had the largest nitrogen export. The mean nitrogen export in 2018 was 30.20 t·ha −1 , and the nitrogen export capacity of each ecosystem remained unchanged (Figure 4g).
The habitat quality of glacial ecosystems was the highest in the study area, followed by forest ecosystems. The desert ecosystem was the natural ecological system with the lowest-quality habitat. During the period of the study, the mean habitat quality increased from 0.457 in 2000 to 0.466 in 2018 (Figure 4h).
Forest ecosystems had the highest entertainment tourism services, followed by shrubland and water ecosystems. During the study period, the entertainment tourism services of the glacier ecosystem improved significantly (Figure 4i).
These findings provide valuable insights into the diverse ecosystem services and their dynamics in the Qilian Mountains, highlighting the importance of different eco-systems in supporting various aspects of human well-being.   Table 7 shows that the global Moran I values of ES in the Qilian Mountains from 2000 to 2018 were all positive, and the p values were all less than 0.001, indicating that ES in the Qilian Mountains interacted with one another on the grid scale and showed the characteristics of aggregation distribution. Specifically, the value of the global Moran's I index for food supply service was the largest, which indicated that food supply was the most concentrated in spatial distribution and had the strongest spatial aggregation. The global Moran's I value of nitrogen export was the smallest, with the smallest spatial autocorrelation and the weakest spatial aggregation, showing the characteristics of agglomeration in a small range. During the study period, the global Moran's I value of raw material supply, climate regulation, soil conservation, and habitat quality services increased gradually, indicating an increase in spatial aggregation and a progressively stronger spatial autocorrelation. The global Moran's I value of food supply and nitrogen export gradually decreased, and the spatial autocorrelation gradually weakened. To further determine the spatial distribution of each ES aggregation area and the trade-off and synergistic relations between different ES, local Moran's I was employed, providing insights on the spatial distribution characteristics. Local Moran's I divided the spatial relations into high-high, high-low, low-high, and low-low aggregation.
The spatial distribution of ES hotspots in the Qilian Mountains remained relatively stable throughout the study period, except for notable increases in hotspots for water conservation and entertainment tourism ( Figure 5). The high-value regions of climate regulation were concentrated in the eastern Qilian Mountains. Furthermore, there was a significant spatial overlap of high-high aggregations between food supply, raw material supply, and entertainment tourism services, indicating a positive correlation among them ( Figure  5a1-a3,b1-b3,e1-e3,i1-i3). In terms of water-related ES, the high-value areas for water resource supply and water conservation exhibited spatial aggregations in the eastern region of the study area, suggesting a cooperative relationship between these two services (Figure 5c1-c3,d1-d3). However, there were distinct spatial patterns observed. For instance, the Qinghai Lake basin showed a high-high aggregation area for habitat quality but a low aggregation area for water resource supply. In contrast, the western region displayed a low aggregation area for habitat quality, while the high-altitude mountainous areas in the west showed high-value areas for water resource supply, indicating a negative correlation in their spatial distribution (Figure 5c1-c3,h1-h3).
The application of Local Moran's I analysis revealed the intricate spatial patterns and interrelationships between different ES in the Qilian Mountains. Understanding these spatial dynamics is crucial for informed decision-making and sustainable management of the region's ecosystems.

Trade-Offs and Synergies between ES
Spearman correlation analysis was conducted to assess the correlation among the nine ES in the Qilian Mountains ( Figure 6). During the study period, the strongest synergy was observed between water resource supply and water conservation services, with a correlation coefficient of 0.97. Food supply, raw material supply, climate regulation, habitat quality, and entertainment tourism showed synergistic relationships with any two of these five services. In addition, soil conservation, water resource supply, and water conservation showed synergistic relationships with one another, and the correlation coefficient showed an obvious trend of an initial increase followed by a decrease during the study period.
Nitrogen export demonstrated a synergistic relationship with six services: food supply, raw material supply, water resource supply, water conservation, climate regulation, and soil conservation. Thus, there was a trade-off relationship between these six services and the water purification service. There was a synergistic relationship between the water purification service and habitat quality, and the correlation coefficient first increased and then decreased. Soil conservation had a trade-off relationship with food and raw material supply, and the correlation coefficient first increased and then decreased. The relationship between habitat quality and three services, namely water resource supply, water conservation, and soil conservation, is a trade-off relationship. Water resource supply had a trade-off relationship with food supply, habitat quality, and entertainment tourism, and the correlation coefficient between water resource supply and entertainment tourism first increased and then decreased. dicates that the test for significance has been passed (p-value less than 0.05). FS, RMS, WS, WC, RC, SC, NL, HQ, and ET represent food supply, raw material supply, water resource supply, water conservation, climate regulation, soil conservation, nitrogen export, habitat quality, and entertainment tourism services, respectively. Red and blue indicate positive and negative correlations, respectively.

k-Means Clustering Analysis
By using k-means clustering analysis, four types of ecosystem service clusters were identified in the Qilian Mountains (Figure 7). The results showed that the same ESB has obvious aggregation relationships in space, and different bundles show significant differences. According to the characteristics of the ESB, they were named food supply, windbreak and sand fixation, water conservation, and ecological coordination bundles.
The food supply bundle was mainly located at the southeastern edge of the Qilian Mountains, and the ecosystem type was predominantly farmland (Figure 7a-c). During the study period, the area proportion of these ESB increased slightly, with an average proportion of approximately 2.5% (Figure 7g). In this service bundle, human activity was intense, agricultural production activity was frequent, and ecosystems were fragile. In addition, food supply, water purification, and climate regulation services were high, but raw material supply, habitat quality, and entertainment tourism services were relatively low (Figure 7d-f). With population growth in the Qilian Mountains, the demand for food supply services will increase, affecting the regional supply and demand structure. This may lead to the continuous expansion of the food production bundle, which would have a negative impact on the ecological balance of the surrounding areas.
The windbreak and sand fixation bundle were mainly distributed in the western Qilian Mountains and had a single ES function. The main vegetation types were desert and alpine steppe, and most glaciers in the Qilian Mountains were also distributed in this region (Figure 7a-c). During the study period, the area proportion of these ESB showed a trend of fluctuating increase, accounting for approximately 44.2% in 2018 (Figure 7g). These ESB had higher soil conservation and water resource supply services, but the food supply, water conservation, habitat quality, and entertainment tourism services were all lower than the mean level of the study area (Figure 7d-f). It is worth noting that although the vegetation coverage in the windbreak and sand fixation bundles is low, the actual wind erosion in the region is still far lower than the potential wind erosion, which is of great significance for weakening the wind erosion intensity, resisting wind sand, preventing desert expansion, and maintaining water and soil.
The water conservation bundle was mainly distributed in the middle of the Qilian Mountains and was characterized by strong water conservation ability (Figure 7a-c). In addition to the low food supply service, this ESB had an extremely high-water resource supply, water containment and soil conservation services, and relatively high raw material supply, water purification, climate regulation, habitat quality, and entertainment tourism services. This bundle was in good ecological condition, belonging to the ecological surplus area (Figure 7d-f). From 2000 to 2018, the area of the water conservation bundle increased significantly from 0.2% to 8.8% (Figure 7g). The increasing area was mainly concentrated in the eastern forest and grassland ecosystems, which had better water and thermal conditions, high vegetation coverage, and a good water conservation function.
The ecological coordination bundles were mainly concentrated in the eastern Qilian Mountains (Figure 7a-c). This bundle exhibited elevated levels of raw material supply, climate regulation, soil conservation, water purification, habitat quality, entertainment tourism services, and low food supply services, which were mainly grassland, forest, and shrubland ecosystems (Figure 7d-f). During the study period, this service cluster area showed a decreasing trend, accounting for approximately 44.5% in 2018 (Figure 7g).

Temporal and Spatial Evolution of ES
The ecosystem processes in the Qilian Mountains, including the water and carbon cycles, vegetation growth, and material and energy exchange, involve complex interactions and mutual influences. The distribution structure of ES exhibits heterogeneity due to unevenly distributed factors such as geology, topography, and hydrology, as well as social factors such as road networks, population distribution, and economic development levels [74].
Atmospheric water vapor and surface precipitation in the Qilian Mountains are affected by westerly wind belts, southerly monsoons, and the East Asian monsoon [81]. In recent decades, under the influence of global warming, the Qilian Mountains have shown a trend of warming and humidification. The characteristics of warming and humidification were obvious, and rainfall increased significantly, showing a higher spatial distribution in the southeast and lower spatial distribution in the northwest [82,83]. As a result, vegetation coverage in the eastern region was significantly higher than that in the western region. In the past 20 years, the glaciers in the Qilian Mountains have melted, and the snow line has moved up, which has a significant impact on the runoff out of the mountain pass [84]. In addition, the actual evaporation in the alpine region is small, the soil water holding capacity is weak, and the local runoff is easily formed, which leads to a strong water-producing capacity. Water conservation capacity and water resource supply were closely related in the study area, showing an obvious synergistic relationship. In the western Qilian Mountains, vegetation is sparse, and the ability to intercept precipitation is poor, resulting in a low water conservation capacity compared to areas with high vegetation coverage. Forests have the strongest water conservation capacity, which mainly includes the canopy and litter interception of precipitation and soil water storage to keep water in the soil for a long time [85]. On the one hand, farmland ecosystems are the most important source of food supply; on the other hand, farmland ecosystems have the largest nitrogen export but have a poor water purification ability, which is also the main factor leading to water eutrophication.
Climate regulation services in the Qilian Mountains increased annually during the study period, which is in accordance with the results of previous studies [86,87], and there is a good coupling relationship between the size of climate regulation services and the spatial variation trend of vegetation coverage [88]. Areas with high soil conservation value are predominantly found in the forest, shrubland, and grassland ecosystems with high vegetation cover in the eastern Qilian Mountains, as well as fragile desert and alpine grasslands in the west, which require priority protection measures. The grassland ecosystem accounted for approximately 44% of the total study area and had the highest amount of soil conservation of all ecosystem areas. Although the northwest region had less rainfall, dominant wind erosion, and low vegetation coverage, the desert vegetation in this region had a considerable amount of soil conservation.
During the study period, habitat quality in the Qilian Mountains was relatively high and stable, although a slightly increasing trend was observed. Owing to climate change, glacier retreat [84], urban expansion, and human activities, ecological degradation has occurred in some areas of the Qilian Mountains. Therefore, managing and developing land resources more scientifically and sensibly, maintaining biodiversity, and promoting sustainable development are still critical issues that need continuous attention from the academic community. Simultaneously, it is of great significance to protect the sparse grassland ecosystem in the high-altitude area of the Qilian Mountains through the further prohibition of grazing by either closing or implementing seasonal grazing of grasslands with different coverage, which will help maintain the stability of the ecosystem.
The continuous improvement of transportation facilities further improves the accessibility of tourism and expands the tourism attraction of the Qilian Mountains and surrounding areas, especially the entertainment tourism service of glaciers. Given the uniqueness and significance of the region's glaciers, they have become hotspots for scientific research and climate change investigation. Glacier resources are scarce and possess high tourism value [89]. Thus, this ecosystem plays a vital role in promoting local economic development, providing employment opportunities, increasing the income of residents, and maximizing social benefits. It is crucial to protect, plan, and develop glacier resources based on the principles of protection priority and sustainable utilization.

Trade-Offs and Synergies of ES
The Qilian Mountains exhibit a strong synergistic relationship between water conservation and water resource supply, which is in line with previous research [90]. Precipitation and evapotranspiration jointly affect the spatial distribution of regional water production, and the vegetation coverage and soil conservation of precipitation determine the capacity for water conservation. Therefore, regions with sufficient precipitation and high vegetation coverage are conducive to regional water production and conservation, leading to a synergistic relation between water resource supply and conservation. The positive correlation between water resource supply and soil conservation aligns with previous studies [91,92]. This is mainly because vegetation coverage is directly affected by the supply of water resources. In areas with good hydrothermal conditions, increased precipitation will promote the growth and development of plants. The subsequent increase in vegetation coverage will improve resistance to soil erosion and reduce water and wind erosion, which is conducive to soil conservation [65]. Water resource supply and habitat quality had a trade-off relationship, which was due to the apparent vegetation recovery in the Qilian Mountains in recent years. The areas of the forest and grassland ecosystems increased significantly, and habitat quality improved. However, with an increase in vegetation coverage, evapotranspiration increased, surface runoff decreased, and water resource supply decreased. This is consistent with previous research findings [64].
There was a strong positive correlation between food supply and nitrogen export. Therefore, compared to the other ecosystems, the farmland ecosystem had the highest food supply service and nitrogen export capacity and the lowest water purification capacity. The reason for this is that the use of chemical fertilizers increases with an increase in the regional food supply service, which leads to an increase in nitrogen export [93]. Hence, a synergistic relationship exists between food supply and nitrogen export, while a tradeoff relationship is observed between food supply and water purification services. Additionally, a weak trade-off relationship exists among soil conservation, food supply, and raw material supply, consistent with previous findings [28,94,95]. The amount of chemical fertilizer that is applied is a crucial factor in determining the trade-off relationship between food supply and soil conservation [96].
There is a strong trade-off relationship between water resources and food supply [93,97]. The main reason for this is that agricultural production consumes a significant amount of water resources, which decreases the amount of surface runoff. This reflects the inherent conflict between food production and water resource use and environmental protection in a specific region [21,26].

k-Means Clustering and Suggestions on Zoning Management
There is clear spatial heterogeneity in ES and its trade-offs [98]. Ecosystem management aims to enhance one or more ES. According to the trade-offs and synergies between ES, ecological protection and management measures are adopted according to local conditions. This approach is crucial for improving key regional ES, enhancing overall ecological benefits, and providing decision-making support for regional development and ecological construction [25]. The analysis of trade-offs and synergies between ES is the basis for the effective management of the ecosystem of the Qilian Mountains. The division of ES is an effective aid in the formulation of regional management policies [99].
ESB form different spatial relationships with the comprehensive action of the societyecological environment, and this relationship changes over time [36]. From the perspective of the evolution of ESB, the reason for this change is a combination of climate change and human activities [25]. Recent decades have witnessed a significant increase in precipitation and a warming trend in the Qilian Mountains [82,83]. Simultaneously, with the increase in regional population, the demand for food supply services will increase, and land with good water, soil, light, and heat conditions will be reclaimed for use as farmlands, which will lead to trade-offs between different ES in some areas [25]. The food supply bundle of the Qilian Mountains was dominated by the farmland ecosystem, and there is a spatial trade-off between the food supply and water resource supply services, which is associated with the large number of water resources required for food production [97]. The windbreak and sand fixation bundle were located northwest of the Qilian Mountains. Owing to the poor hydrothermal conditions in this area, a landscape with low vegetation coverage, desert, and sparse grassland was formed. In addition to soil conservation, all ES in this area were low. The ecological functions of windbreak and sand fixation guarantee ecological security in the fragile desert area of western China. The area covered by the water conservation bundle in the central region increases annually. Compared with other bundles, this bundle shows a strong water conservation ability, which is related to high vegetation coverage and the strong ability of this region to intercept, infiltrate, and store precipitation. The variety of ES in the ecological coordination bundle was greater than the average level of the study area, which was largely related to the wide distribution of grassland ecosystems in the region. Pasture, as the main food source of livestock, contributes to animal husbandry development. Grasslands also play a significant role in carbon dioxide absorption through photosynthesis, contributing to climate regulation. Hence, this bundle shows high supply and regulation services.
The purpose of ecosystem management is to enhance one or several ES and divide the study area into different ESB according to the spatial distribution characteristics of the ecosystem. The goal of this is to clarify the objectives of regional management and take protection measures according to local conditions, which will improve the value of key ES in a region and enhance the overall benefits of regional ES. Thus, it is important to maintain the sustainable development of ecosystems.
Based on our results from the four types of ESB, the following suggestions are made: (1) Food supply bundle: strict adherence to ecological protection guidelines and safeguarding food security are imperative. Enforcing laws and regulations on land use, prohibiting unauthorized cultivation, strengthening water and soil conservation, and improving ecosystem stability are recommended. (2) Windbreak and sand fixation bundle: in the future, we should optimize the scope of the red line of ecological protection, focus on protection, reduce the intensity of human activities, adopt ecological improvement, increase the vegetation coverage under appropriate natural conditions, strengthen the ecological control mode of ecological restoration, set up ecological restoration demonstration areas to improve the quality of the eco-environment and service functions, and give full effect to the role of ecological security barriers in western China. (3) Water conservation bundle: all activities that disrupt ecological balance and cause water pollution should be prohibited. The utilization of water resources should be sensibly controlled, soil and water conservation forests should be widely planted, closed management should be adopted, and grazing and other production activities should be strictly prohibited to reduce the interference of human activities as much as possible and reduce ecological risks. (4) Adhering to ecological protection principles as a priority and allowing economic development. We will continue implementing the policy of herbage and livestock balance, strengthening ecological protection and restoration, and improving regional ecological functions and service capacity.

Strengths, Limitations, and Future Research
One strength of the study is its use of correlation coefficients to analyze trade-off and synergy relationships among different ecosystem services (ES). Although correlation coefficients provide insights into these relationships, they may not fully capture the underlying mechanisms driving trade-offs and synergies [100,101]. Future research should explore these mechanisms to gain a deeper understanding of the causes and dynamics of trade-offs and synergy relationships among ES.
Another limitation of the study is its focus on analyzing trade-offs and spatial bundles of ES without considering the specific driving factors that influence the spatial and temporal distribution of different ES. It is important to investigate the specific factors that contribute to these patterns, as they can provide valuable insights for effective ecosystem management. Future research should aim to identify and analyze the driving factors that influence ES distribution, considering both natural and human-induced factors.
The use of k-means clustering in this study is a valuable approach to ecological partitioning. It overcomes some of the subjectivity associated with previous methods. However, it is important to note that the algorithm itself is a tool, and a more scientific and persuasive ecological partitioning can be achieved by combining the algorithm with expert suggestions from relevant fields. Future research should explore ways to integrate objective algorithms with expert knowledge to enhance the accuracy and robustness of ecological partitioning.
Future research should also consider the scale effect of trade-offs and synergy relationships among ES. Understanding how these relationships vary at different scales can provide valuable insights for ecosystem management and decision-making. Additionally, exploring the threshold of driving factors can help identify critical points at which tradeoffs and synergies become more pronounced or shift. This knowledge can inform targeted interventions and management strategies.
In summary, future research should delve into the underlying mechanisms driving trade-offs and synergies among ES, investigate the specific driving factors influencing ES distribution, explore the integration of objective algorithms with expert suggestions, consider the scale effect of trade-offs and synergies, and identify thresholds in driving factors.
By addressing these aspects, a more comprehensive understanding of ES dynamics and effective ecosystem management strategies can be achieved.

Conclusions
The spatial and temporal patterns and changes in the characteristics of different ES exhibit distinct patterns and changes influenced by global climate change and human activities. Regarding spatial distribution, the high-value area of food supply services was mainly concentrated in the farmland ecosystem in the southeast, which was the low-value area of water purification services. The high-value areas for water resource supply services were mainly located in the high-altitude mountainous areas. Forest ecosystems were the high-value areas for water conservation, climate regulation, and habitat quality services. The high-value areas for raw material supply services were mainly concentrated in the eastern grassland and forest ecosystems. The high-value areas for soil conservation services were concentrated in the western and eastern parts of the Qilian Mountains, where the former was a high-value area for preventing wind erosion, with the desert as the main ecosystem type, and the latter was a high-value area for preventing hydraulic erosion, with forest and shrubland as the main ecosystems. The high-value areas for entertainment tourism services were mainly distributed around Qinghai Lake and in the eastern forest and shrubland ecosystems. During the study period, in addition to habitat quality and water purification, which remained unchanged, and soil conservation, which showed a fluctuating growth trend, six ES (food supply, raw material supply, water resource supply, water conservation, climate regulation, and entertainment tourism) showed a significant increase, and the overall ecology of the Qilian Mountains continued to improve.
Trade-offs and synergistic relationships among ES affected the spatial distribution of ESB. The food supply bundle was dominated by farmland ecosystems, which exhibited high food supply, high climate regulation, and low water purification services, and the division of this bundle was related to the synergistic relationship between food supply and climate regulation services and the trade-off relationship between food supply and water purification services. The windbreak and sand fixation bundle were dominated by grassland ecosystems, which showed high soil conservation, higher water resource supply, and low habitat quality services; the classification of this bundle was related to the synergistic relationship between soil conservation and water resource supply and habitat quality had a trade-off relationship with soil conservation and water resource supply. The water conservation bundle was dominated by forest and grassland ecosystems, which exhibited high water resource supply, water conservation, and soil conservation services, and the division of this bundle was related to the synergistic relationship among water resource supply, water conservation, and soil conservation. The ecological coordinated bundle was dominated by grassland, forest, and shrubland ecosystems, which showed high raw material supply, climate regulation, habitat quality, and entertainment tourism services, and the division of this bundle was related to the synergistic relationship among raw material supply, climate regulation, habitat quality, and entertainment tourism services. Further research should focus on the drivers and thresholds of ES to address the impacts of future climate change on ecosystem management.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A Appendix A.1. Changes in Spatial and Temporal Patterns of Ecosystem Services in the Qilian Mountains
The spatial and temporal distributions of different ES in the Qilian Mountains exhibit clear heterogeneity, as depicted in Figure A1. Throughout the study period, the spatial distribution pattern of food supply services remained consistent. The high-value areas were primarily concentrated in the farmland ecosystem around the southeastern region, while the low-value areas were mainly focused on the desert ecosystem in the northwestern region ( Figure A1(a1-a3)).
Similarly, the spatial pattern of the water purification service remained unchanged. Regions with limited water purification ability were predominantly observed in the southeastern farmland ecosystem, where the output of nitrogen was extremely high. Conversely, areas with high water purification capability were mainly concentrated in the eastern forest ecosystem ( Figure A1(g1-g3)).
Regarding raw material supply, the high-value areas were distributed in the grassland ecosystem surrounding Qinghai Lake and the eastern forest ecosystem. In contrast, the alpine meadow area at high altitudes represented the low-value area ( Figure A1(b1-b3)).
The spatial distribution of the water resource supply service in the Qilian Mountains exhibited variations. High-value areas were primarily concentrated in high-altitude mountainous areas, characterized by abundant precipitation, sparse vegetation, limited soil infiltration capacity, and easy surface runoff formation. On the other hand, the lowvalue areas were mainly located in the northwestern region, characterized by low precipitation, strong evaporation, and mainly desert ecosystems ( Figure A1(c1-c3)).
Regarding the water conservation service, the high-value areas of the water conservation service were mainly concentrated in forest and grassland ecosystems in the eastern region. Throughout the study period, there was a significant increase in the high-value water conservation areas, primarily expanding in the eastern forest ecosystem and the high-coverage grassland with higher vegetation coverage and favorable water and heat conditions ( Figure A1(d1-d3)).
The spatial distribution of the climate regulation service demonstrated distinct characteristics, with low values in the western region and high values in the eastern region. The high-value regions were concentrated in the Qinghai Lake basin and the eastern forest ecosystem, which benefit from favorable hydrothermal conditions promoting vegetation growth ( Figure A1(e1-e3)).
The soil conservation service in the Qilian Mountains exhibited distinct spatial patterns. The high-value regions were concentrated in the western and eastern parts of the Qilian Mountains, while low-value regions were mainly found in glaciers and surrounding desert ecosystems ( Figure A1(f1-f3)). The western region, characterized by high sunshine intensity, a large diurnal temperature range, little precipitation, strong evaporation, and wind erosion, played a dominant role. As a result, it was identified as a high-value area for preventing wind erosion in the desert ecosystem, contributing to significant soil conservation in that area. In contrast, the eastern region received more precipitation, with hydraulic erosion being the dominant erosional process. The forest and shrubland ecosystems in the eastern region were identified as high-value areas for soil conservation.
Habitat quality showed notable variations between the eastern and western regions, with the eastern region exhibiting significantly higher habitat quality compared to the western region ( Figure A1(h1-h3)). However, it is worth noting that the continuous melting of glaciers has led to the transformation of previously glacier-covered areas into bare lands without vegetation cover, leading to a slight reduction in the quality of mountain habitats in some areas.
The areas with high values of the entertainment tourism service were mainly distributed around Qinghai Lake and the regions characterized by concentrated forest and shrubland ecosystems in the eastern part of the Qilian Mountain. These areas possess a picturesque ecological environment and convenient transportation, making them suitable for tourism activities ( Figure A1(i1-i3)). The water yield module is based on the Budyko water-heat coupling equilibrium hypothesis [53], and the calculation formula is as follows: where represents the annual actual evapotranspiration of the grid cells (mm); represents the annual precipitation of grid cells (mm); represents potential evapotranspiration; represents the non-physical parameters of natural climate-soil properties; represents the evapotranspiration coefficient of plants (vegetation) for specific land use types in the grid; the minimum value of is 1.25 [102]; and is the effective water content of the soil (mm), which is determined by the soil texture and effective depth. Z is an empirical constant that represents the regional precipitation distribution and other hydrogeological characteristics.
represents the grid's crop reference evapotranspiration; is net radiation (MJ/(m 2 ·d)); is soil heat flux (MJ/(m 2 ·d)); is a constant on the wet and dry scale; is the average daily temperature (°C); is 2 m wind speed (m/s); and are saturated and actual vapor pressure (kPa), respectively; δ is the slope of the saturated vapor pressure curve (kPa/°C); and PAWC is the available water content of plants. The calculation formula for PAWC is as follows [ where SAN is the sand content (%), SIL is the silt content (%), CLA is the clay content (%), and C is the soil organic carbon content (%). The specific parameters of the water production model adopted in this study are listed in Table A1 [52,103]. This study used the water balance equation to calculate the amount of water conservation. The computation formula is as follows [54][55][56]: where is the precipitation of the class ecosystem (mm), is the surface runoff of the class ecosystem (mm), is the evapotranspiration of the class ecosystem (mm), is the area of the class ecosystem, and is the number of ecosystem types in the study area.
The InVEST water production model calculates the water resource supply of each grid cell, i.e., the amount of water, by subtracting the actual evapotranspiration from the precipitation of each grid cell. Based on the water balance equation, the formula for calculating water conservation is as follows: where (mm/a) is the surface runoff of pixel in a class ecosystem. Surface runoff was calculated according to the runoff coefficient and rainfall. The runoff coefficient is the ratio of rainfall to runoff volume. Influenced by several factors, the larger the runoff coefficient, the less easily rainfall is absorbed by soil [104]. By referring to the domestic literature on the precipitation and surface runoff of runoff districts in different ecosystems [56,[105][106][107][108], combined with the climatic characteristics of the study area, vegetation coverage of different ecosystems, colony species, and other factors, the runoff coefficient was assigned, and the results are shown in Table A2. In this study, the soil conservation module of the InVEST model was used to calculate soil conservation in the water erosion area. Based on the soil loss equation (RUSLE), this module can evaluate the soil conservation function more intuitively and rationally [109]. The calculation formula is as follows: where SC represents soil conservation (t·ha −1 ·a −1 ); R is rainfall erosivity; K is the amount of soil loss per unit area caused by unit rainfall erodibility (t·h·MJ −1 ·mm −1 ), which represents the erodibility of soil affected by soil texture, organic matter content, soil structure, permeability, and other factors, and reflects the anti-erosion ability of soil. LS is the slope length and slope factor, C is the vegetation cover factor, and P is the soil and water conservation factor. The calculation formula of soil erodibility factor K is as follows [110] where SAN is the sand content (%), SIL is the silt content (%), CLA is the clay content (%), C is the soil organic carbon content (%), and SN = 1 − SAN/100. The mean annual rainfall was used to estimate the erosive force (R). The formula is as follows [111]: where P is the average annual precipitation (mm), and R is the average annual rainfall erosivity (MJ·mm·ha −1 ·h −1 ·a −1 ). In this study, α and β are model parameters, which are 0.0668 and 1.6266, respectively [112]. According to relevant research results [52,113,114], the values of the vegetation cover factor (C) and soil and water conservation factor (P) were assigned based on the situation in the Qilian Mountains. The results are shown in Table A3 below. Table A3. Vegetation cover factor (C) and soil and water conservation factor (P). The causes and erosion process of wind erosion are complex, and the parameters of most models are difficult to obtain. In addition, the models themselves have limitations. In this study, the drainage basin wind erosion model established by Dong et al. (1998) [57] was adopted, and the model formula was as follows:

Land Use Types
where θ is the slope (°), V is the wind speed (m/s), VCR is the vegetation coverage (%), and SDR is the damage rate of the artificial surface structure (%) that is obtained by referring to the hypothesis of the equation in the wind tunnel simulation, H is the relative humidity of air (%), d is the average particle size of soil particles (mm), F is the soil hardness (N/cm 2 ), x and y represent the distance from the reference point, and t represents time (s).
The formula for calculating the vegetation coverage is as follows: The NDVI probability distribution values of 95% and 5% were used to determine the threshold NDVI values representing pixels completely covered by vegetation and the areas covered by bare soil or no vegetation, respectively. Soil wind erosion is a form of soil erosion that occurs without considering vegetation coverage. The vegetation's contribution to soil conservation under wind erosion can be calculated as the difference between the potential soil wind erosion and the actual soil erosion, taking into account the protective effect of vegetation.
Appendix A.2.4. Nutrient Delivery Ratio (NDR) Module The principle of the NDR module of the InVEST model is to calculate the long-term steady-state flow of nitrogen in the entire watershed based on the ability of different ecosystem types to transfer nitrogen into runoff and to reflect the water quality in the watershed according to the content of nitrogen in the water body. The higher the nitrogen content, the more serious the water pollution or eutrophication in the basin and the weaker the water purification capacity of the ecosystem in the basin. The nitrogen output of each grid was calculated based on the product of the nutrient load and NDR. The calculation formula is as follows: The nitrogen content of each grid was correlated with the nutrient load. Nutrient loads for each land use type are empirical coefficients based on nutrient output measurements [115].
where is the runoff potential index of grid ; is the nutrient runoff proxy for runoff on pixel ; is the average runoff potential over the entire area; _ is the ratio parameter between the two nutrient sources; and , and , are the nutrient transport ratios of the surface and subsurface, respectively.

Surface
is the product of a delivery factor, representing the ability of downstream pixels to transport nutrients without retention, and a topographic index that represents the position on the landscape. For a pixel : and are calibration parameters, is a topographic index, and , is the proportion of nutrients that are not retained by the downstream pixels (irrespective of the position of the pixel on the landscape). Below, we provide details of the computation of each factor. , 1 (A22) where is the effective downstream retention on the pixel directly downstream from , is the maximum retention efficiency that Land use types can reach, and is the step factor defined as: where is the length of the flow path from pixel to its downstream neighbor. This is the Euclidean distance between the centroids of two pixels.
is the LULC retention length of the land cover type on pixel .
where is the connectivity index, is the average slope gradient of the upslope contributing area (m/m), is the upslope contributing area (m 2 ), is the length of the flow path along pixel according to the steepest downslope direction (m), and is the slope gradient of pixel .

Subsurface NDR
The expression for the subsurface NDR is a simple exponential decay with distance to stream, plateauing at the value corresponding to the user-defined maximum subsurface nutrient retention: where is the maximum nutrient retention efficiency that can be reached through the subsurface flow, is the distance from the pixel to the stream, and is the subsurface flow retention length.
Based on field investigations, background data collected in the study area, and relevant local and international research [52,[116][117][118][119], we built an NDR model parameter table, as shown in Table A4. Habitat quality simulated by the habitat quality module of the InVEST model can accurately reflect the status quo of habitat services [58]. The habitat quality index reflects the overall status of regional biodiversity. Roads, building land, farmland, and bare land are stress factors that threaten biodiversity. The calculation formula is as follows [120]: represents the habitat quality index of grid x for land use type j in the study area, and the habitat quality index ranges from 0 to 1.
represents the habitat fitness of land use type , represents the threat degree of land use type , constant is a halfsaturation constant, and z is a normalized constant, usually 2.5 [121]. is the number of threat sources, is the grid in the threat element, is the weight of different threat sources, is the stress value of grid , and is the anti-interference level of the habitat. is the relative sensitivity of different habitats to various threat factors, is the influence of threat source in grid y on grid , is the distance between grids and , and is the maximum influence distance of the threat source . For the assignment of related parameters, this study referred to the user manual of the invest model and the related literature [99,[122][123][124]. For specific parameter settings, see Tables A5 and A6.

Food Supply Service
The food supply service in the study area was quantified by the output values of agriculture, animal husbandry, and fisheries, and the corresponding ecosystem types are shown in Table A7. There is a significant linear relationship between crops and surface NDVI [59]. Therefore, the output value of agriculture and animal husbandry can be spatialized based on the NDVI [59,60,98]. The fishery output value was spatially equally distributed according to the area of the aquatic ecosystem. The calculation method is as follows: is the output value of item industry in the pixel (yuan/pixel), is the NDVI of item j industry in the pixel, is the sum of the NDVI values of item j industry, and is the total output value of item industry in each district and county.

. Raw Material Supply Service
The raw material supply service included wood and non-wood products provided by forests, fresh grass, and refined feed provided by grasslands. Considering the availability and visualization of data, this study quantified the raw material supply service provided by forest and grassland ecosystems in the study area according to the forestry output value and grass yield, respectively. Considering the forestry industry's main output, the forestry output value per unit area was calculated based on the forestry output value and area of the administrative regions involved in the study area. The raw material supply service of the forest ecosystem in the study area was obtained based on the NDVI. The calculation method is as follows: * (A34) where is the output value of item industry in the -th pixel (yuan/pixel), is the NDVI of item industry in the -th pixel, is the sum of the NDVI values of item industry, and is the total output value of item industry in each district and county.
The measurement of grass yield mainly included the cutting, yield simulation model, and remote sensing model methods [61]. The remote sensing model measurement method is the most commonly used method for estimating the grass yield and livestock carrying capacity. The CASA model is usually used to calculate the NPP of vegetation [62], and grassland yield is calculated according to the functional relationship between the dry weight of organic matter produced by the grassland ecosystem and the grass yield per unit area. The formula used is as follows: where is the annual total hay yield per unit area (g·m −2 ·a −1 ); NPP is the annual total NPP (g C·m −2 ·a −1 ) of grassland (see climate regulation service for the calculation method); is the conversion coefficient of grassland biomass (g) converted to NPP, which is 0.45 [61]. was the ratio coefficient of subsurface to above-ground biomass, which was 8.99 for alpine meadow, 3.81 for alpine steppe [125], 4.25 for warm steppe, 7.89 for warm desert steppe, alpine desert and warm desert, 4.42 for warm shrubland, and 15.68 for swamp [63]. The value is calculated based on the hay price using the following formula: where is the value of grassland raw materials, is the total hay yield, and is the hay price. The research team distributed 146 questionnaires to herders in the Qilian Mountains. This was followed by sorting, screening, and elimination of invalid questionnaires. A total of 133 valid questionnaires were obtained, with an effective rate of 91.10%. The prices of corn, oats, straw, and other forage ranged from 2400 to 4900 yuan/t. During this study period, the hay price was 3012.5 yuan/t. In ArcGIS, spatial analysis tools were used to conduct Euclidean distance analysis on roads, protected areas, national parks, and cultural attractions. The maximum impact distance was set to 10 km [130], and the results were normalized as 1-100 according to the principle that the farther the distance, the smaller the value. According to the actual situation of the study area, the activities that could be conducted for each land use type were assigned a value. The activities that could be conducted were assigned a value of 1, and those that could not be conducted were assigned a value of 0 [46], as shown in Table A9. The obtained raster layers were standardized to 1-100.