Spatial Variation and Terrain Gradient Effect of Ecosystem Services in Heihe River Basin over the Past 20 Years

: With the advent of large-scale development, extreme imbalance in the ecology of the Heihe River Basin (HRB) has caused a series of ecological problems. In order to explore the spatiotemporal variation of ecosystem services (ESs) and to assess the characteristics of ESs under the terrain gradient effect (TGE), the three key ESs were quantified based on the InVEST model using five series of land-use data obtained from remote sensing images from 2000 to 2020 in this study. The terrain index was used to analyze the influence of terrain on ESs. The results show that most of the ESs were in high numbers in the south and low numbers in the north, as well as high numbers in the middle and upper reaches and low numbers at downstream locations. It was found that high-quality habitats degrade to general-quality habitats, and poor-quality habitats evolve into general-quality habitats. It was also found that the water production volume continues to decline and soil conservation becomes relatively stable with little change. This study illustrates different ESs showing obvious TGE with changes in elevation and slope. These results indicate that the effect of land-use change is remarkable and TGE is highly important to ESs in inland watersheds. This research study can provide a scientific basis for the optimization of regional ecosystem patterns. The results are of great significance in terms of rational planning land use, constructing ecological civilizations, and maintaining the physical conditions of land cover at inland river basins.


Introduction
Ecosystem services (ESs), which are the benefits that humans receive directly or indirectly from ecosystems, have played an important role in the development of human civilization [1]. Generally, ESs include food supply, climate regulation, and cultural and supporting services. In recent years, the conflicts between the rapid socio-economic development and the consumption of natural resources and environmental change have resulted in the degradation of ESs and frequent natural disasters in some areas [2][3][4]. ESs have been considered as indispensable features in land-use planning and resource adjustment, becoming one of the research hotspots over recent decades. Currently, ES assessments mostly focus on different land-use scenarios, such as croplands, forests, oceans, lakes, wetlands, and others [5][6][7][8]. In order to ensure the sustainable provision of ESs for human beings, it is necessary to evaluate ESs in order to support human wellbeing. ES evaluations consist of a value assessment and material quality assessment [9][10][11]. The former measures the economic value of ESs by using ecological economics methods [12][13][14] that are based on the value equivalent factor of a unit area or the price of a unit service function [15][16][17][18][19][20][21]. The value measurement method is primarily used to evaluate supply services, and the results are of high economic significance [22][23][24]. With the refinement of the ES classification, material quality methods are gradually used as the mainstream assessment methods of ESs. The material quality of various ESs can be quantitatively eval-uated based on ecological principles in this method. By contrast, the credibility and accuracy of the quality assessment method are higher than that of the value assessment method, which is influenced subjectively to some extent [25]. It is necessary to conduct research on the sustainability of an ecosystem and to provide a more scientific and reliable basis for decision making [11]. Material quality assessment methods can be roughly divided into emergy methods and model methods [26]. The emergy method is based on the emergy theory for evaluating different types of emergy and material flows in the system [27]. It combines the total amount of effective energy invested directly and indirectly in the ecosystem with the energy conversion rate to calculate the final energy value. The energy flow and utilization rate of ESs can be expressed better by emergy evaluation. This method is often used to describe the regional differences with the utilization of ESs in large-scale areas and to evaluate urban ecosystems [28]. However, with increasing studies on the coupling of multiple ESs and the relationship between ESs and human wellbeing during these years, the emergence and development of evaluation models represent a major breakthrough in the field of ES assessment. The results can be presented in the form of a map incorporating more spatial and intuitive information by using model assessments. It can simulate and predict the changing scenario of future ESs.
According to the formation mechanism of ESs, studies on ES evaluation have increased due to the quality assessments of ecosystems generated through comprehensive models. The commonly used models include Integrated Valuation Ecosystem Services and Tradeoffs (InVEST), Artificial Intelligence for Ecosystem Services (ARIES), Social Value for Ecosystem Services (SolVES), and Multi-scale Integrated Models of Ecosystem Services (MIMES) [29,30]. Some models, such as SoLVES, ARIES, and MIMES, output good evaluations for specific regions and have promising application prospects, but they have not been fully developed yet. The InVEST model is the most widely used tool for the quality assessment of ESs. The advantages of the InVEST model have been highlighted by integrating various service production functions or simulating service dynamic changes. Previous assessment studies using the InVEST model are mostly concentrated in smallscale areas, such as administrative regions, urban economic zones, or lakes and rivers [31][32][33][34][35][36]. ES evaluations in medium and large watershed scales are relatively rare. In addition, most of the studies analyzed spatiotemporal pattern changes or showed cold and hot spots in the service space by spatial mapping for individual or several services, such as biodiversity conservation, water conservation, soil conservation, and carbon storage [33][34][35][36][37][38]. Relatively few studies discuss the connection and divergence between different ES supplies and specific natural environmental factors [39,40]. As an important factor in the natural environment, the terrain factor is of great significance for fully understanding the spatiotemporal difference of ESs [41,42]. It is well known that ESs are closely related to their terrain [43,44]. Some terrain-based studies analyzed the spatial distribution of landuse patterns [45,46] and some were focused on mountain areas [47,48]. However, there are only a few studies assessing the impact of topographic factors on the spatial heterogeneity of ESs. The evaluation of the characteristics of ESs with terrain gradients has not been explored in depth. Understanding the influence of topographic elements on the spatial difference of ESs comprehensively is of great significance for deepening the understanding of ES differentiation and for ensuring effective ecosystem management.
The Heihe River Basin (HRB) is not only a key area with respect to the ecological safety barrier, but it is also an important node in the layout of the national "Belt and Road" initiative [49,50]. Typical natural conditions and long-term, frequent human activities have resulted in ecological environmental deterioration in the basin, resulting in the degradation of biodiversity, water production, soil, and water conservation. Therefore, the problem of how to maintain the stability of regional ecosystems as an economy develops is an important issue for all inland river basins in arid zones. To summarize, combined with the geographical conditions and characteristics of the HRB, we discuss the relationship between ESs and land-use changes in typical arid inland river basins from a topographical perspective on the basis of previous research studies. In this paper, we selected three key ESs (i.e., habitat quality, water yield, and soil conservation) of the HRB in Northwestern China to study. Based on the proposed framework (Figure 1), our analysis focuses on addressing three questions at fine scales: (1) How have ESs changed over the past 20 years and where are the areas of high and low distribution of individual ESs? (2) How are different ESs affected by terrain factors with respect to inland watershed scales? (3) How can variations of ESs guide land-use management in a river basin? Figure 1. Framework of the methodological process in this study.

Study Area
The HRB is the second largest inland river basin in arid region of Northwestern China, within the range of 96°42′-102°00′ E, 34°41′-42°42′ N ( Figure 2). The area of core drainage is approximately 128,900 km², with a mainstream length of 821 km. Different reaches of the river are distributed in the three provinces or autonomous region, belonging to Qinghai, Gansu, and Inner Mongolia, respectively. Climatic characteristics of the basin are a typical continental arid climate, little but concentrated precipitation, ample sunshine, and greater diurnal temperature range. With an average altitude over 1200 m, the topography varies significantly from south to north [51]. The regional climate is significantly different in the basin. Therefore, it is an important for ecological function in terms of water production, soil conservation, biodiversity protection, windbreak, and sand fixation in Northwestern China [52]. This basin can be divided into three parts according to the locations of the gorge stations. The upstream region is located in the Qilian Mountain, with low population density, high vegetation coverage, and good ecological environment. This segment, the main contributing area, is bounded by the Yingluo Gorge. The midstream area is a main field for developing human activities, extending from the Yingluo Gorge to the Zhengyi Gorge. This region is rich in light and heat resources. The downstream region, below the Zhengyi Gorge, constitutes areas where runoff disappears. It is mostly made up of desert and bare land [53]. More than two thousand years of different human activities in HRB have caused collision and exchange of various cultures, making natural and human processes meet together. This is an ideal region for studying the variation of different ESs.

Data Sources
The datasets used in this study included: (1) land-use types in 2000, 2005, 2010, 2015 and 2020 (with a spatial resolution of 1km) and NDVI (Normalized Difference Vegetation Index). These were both obtained from the Resource and Environment Science and Data Center (http://www.resdc.cn accessed on 21 March 2021). According to the model input acquirements and the reference of land-use classification system of the Resources and Environment Database, the land-use data were reclassified into six categories: cropland, forest, grassland, waters, built-up land, and unused land; (2) DEM (Digital Elevation Model) and potential evapotranspiration data (with a spatial resolution of 250m) were derived from the National Tibetan Plateau / Third Pole Environment Data Center (https://data.tpdc.ac.cn/zh-hans/ accessed on 24 March 2021); (3) basic geographic information data (including regional boundaries, roads, settlements, etc.) were acquired from the "Digital Heihe River" project of the Cold and Arid Region Scientific Data Center; (4) soil attribute data were derived from Chinese soil datasets (1995) in the Food and Agriculture Organization's Harmonized World Soil Database; and (5) meteorological data were retrieved from the China Meteorological Data Service Center (http://www.nmic.cn/ accessed on 26 March 2021) and the National Earth System Science Data Center (http://www.geodata.cn accessed on 06 April 2021), including annual, monthly, and daily value data from 2000 to 2020.

Methodology
In this study, terrain analysis and integrated maps were carried out by using the ArcGIS (Ver.10.2) software. Regression analyses to detect trend significance were performed in SPSS (Ver. 21.0), and visualized charts were generated by using Microsoft Excel (Ver. 2017). All the spatial assessments of ESs were processed by using InVEST (Ver. 3.8.0). As a strong technical support in ES research, the InVEST model is used for a quantitative assessment of ESs on the basis of ecosystem processes, and a full demonstration of their spatially distributed characteristics [54,55]. The three sub-models, habitat quality, water yield, sediment delivery and retention were selected to quantitatively evaluate and analyze the changes of three ESs (including habitat quality, water production and soil conservation) in this work.

Habitat Quality (HQ)
The HQ in this model can be described as the connection between different types of land use and threat sources. The results and scarcity of HQ can be obtained in accordance with the response degree of different habitats to various threat sources.
The formula for calculating habitat quality ( ) is: where is the habitat suitability of land-use type ; is the total threat level; k is the half-saturation value.
The input data of this sub-module include current land-use and threat data. ArcGIS was used to process the original data (i.e., mosaic, reclassify, merge, etc.) in order to obtain the land-use and threat factor layer data. The HQ correlation function in the model corresponds to the four variables: the relative impact of each threat, the relative sensitivity of each habitat to each threat, the distance between habitat and threat and impact from the threat, and the level at which a grid cell of the habitat is legally protected. In order to reduce the interference of human factor on results, the impacts of the first three variables were discussed in this study. We selected cropland, town, residential area, and road as habitat threats. The corresponding parameters of the model (Tables 1 and 2) were set by referring to the user guide and previous works with similar eco-environments in arid region [56,5657], as well as the suggestions from experts and professors in the ecological field. The half-saturation constant defaults to 0.5 according to the user guide. The WY estimation is based on the principle of water balance, which simplifies the flow process. It does not distinguish between surface runoff, runoff in the soil, and base flow. WY defines the amount of water produced within each grid range by the difference between precipitation and evapotranspiration (including plant transpiration and surface evaporation). Its calculation is based on the acquired data, such as topography, land use, meteorological factors, etc.
The water yield ( ) is calculated using the following formulas: where is the annual actual evapotranspiration for land-use type in category j on grid x (mm); is the average annual value of precipitation on grid x (mm); is the ratio of the year's available water for modified vegetation to the expected amount of water; is the ratio of potential evaporation to rainfall; Z is the Zhang coefficient, which represents the parameters of seasonal rainfall distribution and precipitation depth, ranging from 1 to 30; is the average annual value of available water capacity on grid x (mm); 0 is the potential evapotranspiration on grid x (mm); is the evapotranspiration coefficient of vegetation. The input data of the model include precipitation, reference evapotranspiration, depth to root restricting layer, plant available water fraction, land use, and watersheds. The parameters include Z parameter and evapotranspiration coefficient of vegetation. The precipitation data were interpolated by using ArcGIS. The reference evapotranspiration was calculated by the modified formula. The soil depth data were derived by rasterizing the soil spatial attribute data. The plant available water fraction was set according to the reference results [58]. The watershed boundary was extracted by hydrological analysis with DEM in ArcGIS. The Zhang coefficient was set as 2.2 referring to the instruction manual and previous reference [59]. The plant evapotranspiration coefficient was defined according to the calculation method proposed by FAO and the reference value was given by the model.

Soil Conservation (SC)
The SC can be calculated with soil erosion reduction and sediment retention using the sediment delivery and retention sub-module. The decrease in soil erosion was expressed by the difference between potential erosion and actual erosion. The amount of sediment retention refers to the amount of uphill sand retained by plot. On the basis of cell scale USLE (Universal Soil Loss Equation) calculation method, the calculation on the grid unit scale was simulated using the following data, including land-use data, soil attribute, DEM, rainfall data, vegetation cover factor, and soil and water protection measure factor. The rainfall erosivity was calculated using observation data from meteorological stations in the basin and then interpolated with Kriging method in ArcGIS. The soil erodibility factor was calculated from the soil attribute data by formula-10. The vegetation cover factor and soil and water protection measure factor can be defined by consulting the user guide.
The formulas for calculating SC are defined as the following: where SC is the amount of soil conservation (t· ℎ −2 · −1 ); USLE is the amount of potential soil loss in the original land-use cover (t· ℎ −2 · −1 ); RKLS is the amount of potential soil loss for bare soil (t· ℎ −2 · −1 ); SEDR is the retention from the upstream sediment (t· ℎ −2 · −1 ); R is the rainfall erosivity (MJ· mm· ℎ −2 ·ℎ −1 · −1 ); K is the soil erodibility factor (t· ℎ 2 ·h·ℎ −2 · −1 · −1 ); LS is the slope length gradient factor; P is the soil and water protection measure factor; is the average monthly rainfall (mm); is the average annual precipitation (mm); SAN, SIL and CLA are the content values of sand grains, powder grains and sticky grains (%); 0 indicates the organic carbon content value (%); c is the vegetation cover; NDVI is the normalized difference vegetation index; C is the crop management factor (value range is between 0 and 1). If C is 0, the vegetation cover of the land surface is good and almost not eroded; if C is 1, there is almost no vegetation cover on the surface.

Terrain Niche Index
Terrain is one of the important factors affecting spatial distribution and changes in land use. TGE is also manifested as the variations in ES supply that are caused by landuse changes to some extent. The terrain differences are outstanding in HRB, which has a variety of mountain, plain, basin, and hilly geomorphological types. In order to mirror the comprehensive relationship between terrain condition and ESs, the terrain niche index was selected as a geographical factor to analyze TGE between land-use patterns and terrain gradient in this study. The terrain niche index (T) can be calculated using the following formula: where E and ̅ , respectively, indicate the elevation value of any point in the region and that of the areas where the point is located; S and ̅ , respectively, represent the slope value of any point in the region and the average slope value of the areas where the point is located. If terrain level index of a point is low, the elevation and the slope of this point are low values; if that of a point is centered, the values are high and low. Both are high values if that of a point is large.

Spatial Distribution and Variation of HQ
HQ score is a comprehensive indication of the influence of each threat factor on itself and the impact of habitat on threat resistance. The higher the number is, the better the HQ is. The HQ in this study was classified as higher (0.8-1), high (0.5-0.8), general (0.3-0.5) and poor (0-0.3). Figure 3 shows the proportion of different HQ grades in the basin. More than half of the HRB was in poor and average HQ over the past 20 years. Higher and high HQ account for less than 30% of the total area, and the overall proportion composition is relatively stable from 2000 to 2020. In order to further clearly explore the changes of HQ spatial distribution over time, we calculated the difference based on the HQ data and obtained the interannual changes in four five-year periods ( Figure 4). Spatial differences and variations of HQ are obvious throughout the basin. Most of the upstream region has highand higher-grade HQ. There have been many changes over the past two decades, especially in the two phases of 2005-2010 and 2015-2020. The area of higher-grade HQ increased and the HQ showed an upward trend. The overall ecological environment became better, which might be the benefit of a sparse population and dense vegetation in the upper mountains. A different layout was presented in the middle stream. The south-central region was dominated by general-and high-grade HQ. The area of general HQ level increased, which was due to the human reclamation of wasteland, grasslands, and forest land, and the expansion of cultivated land and construction land in the active areas. At the same time, the area of high-level HQ also reduced with the reclamation and the degradation of grassland and forest land. By incorporating the changes of four stages, it can be found that mostly the HQ in the northern part of the midstream area was at the poor level and it increased over the past 20 years, due to large distribution of unused land. Downstream areas were clearly affected by land use. The HQ was higher in most of the east, with grasslands, woodlands and waters. Obvious increases can be seen through the four periods in this region, while the central and western regions were dominated by unused land with low coverage and the HQ levels therefore dropped. There were few changes in the areas during 2000-2005 and 2010-2015, while the HQ changed significantly and mainly increased in the other two periods. With the transformation of unused land into grasslands and woodlands, as well as the increase in water area over the past 20 years, the area of high-grade HQ in the downstream region increased and that of poor-grade HQ decreased gradually. This indicates that the overall HQ of the downstream area improved, which may be related to the implementation of the Heihe Ecological Water-divide Project (2000).  Totally, the quality of habitats at all levels changed by a small margin, with a slight decrease in poor and high HQ areas and a slight increase at general and high levels. The overall HQ in the basin became better and the improvement was mostly concentrated in the upstream and downstream areas in the northeast, while the quality levels deteriorated in the midstream and lower-west regions.

Spatial Distribution and Variation of WY
The annual WY in this study was classified into one, two, three and four levels from low to higher weight ( Figure 5). In 2000-2015, the proportion of low-grade WY gradually increased, while that of high-grade WY decreased bit by bit. However, during the fiveyear period from 2015 to 2020, the proportion of low-value WY decreased, and that of high-value WY increased. However, the WY in most areas was at a lower level and the depth of water production was not large. As can be seen in Figure 6, the spatial distribution of WY shows great difference. The upstream region was the principal origin for highvalue water production, and the area of high-value WY gradually decreased. In particular, some of them degraded to medium-yielding areas. However, with the exception of largescale decline in the five years from 2010 to 2015, the changes of WY still increased in the other years. The overall trend also caused an increase in the upstream water production. The mid-range showed a decreasing trend from south to north. The high-value waterproducing region in the southern gradually decreased, yet some low-value water-producing areas increased to medium-yielding areas. The central and northern regions were dominated by low-value water production and continued to decrease in most areas over the past 20 years. Compared with the upper and middle reaches, the downstream area is a region with a low-value WY. Although the WY partially increased between 2005 and 2015, there has still been a slight decline in the past 20 years. The overall change was not significant.  Combined with precipitation and land use, it can be found that the spatial distribution shows consistency for WY, precipitation and vegetation in the basin. The region with strong WY capacity has high average annual precipitation and high vegetation cover and low steaming emission. On the other hand, the areas with low annual precipitation and low vegetation cover and strong steaming had a low-value WY. In terms of land-use changes from 2000 to 2020, except rainfall, the variation of WY capacity in HRB was dominated by the area decrease in grassland and forest and the area increase in unused land. Considering the effects of grassland and forest degradation through the whole basin, the overall WY finally showed continuous reduction. Compared to land use in recent years, it can be found that grassland provided the largest proportion among the overall WY services in the basin, reaching more than half of the total WY. This was followed by forest land and cropland.

Spatial Distribution and Variation of SC
The annual SC was classified into one, two, three and four levels, from low to higher values (Figure 7). The soil-conservation composition strongly varied through the whole basin. The low-grade SC accounts for about 70% of the total area, but there was a downward trend year by year. The proportion of SC changed little at middle and high levels. As shown in Figure 8, the spatial distribution of SC in HRB showed that high-value SC had more in the east and less in the north in the southern area. From a regional perspective, the soil holds in the upstream area was relatively small. The SC decreased in 2000-2005 and 2010-2015 and increased significantly with a wide range in the five years from 2005 to 2010. The mesh distribution of high-value SC increased over the past 20 years. The high-value SC focused on the southern part of the mid-range areas. It expanded year by year and the amount of hold increased. Conversely, soil erosion in the north region was more serious and the changes mainly occurred in the mutual transformation of low SC into medium-grade SC. As for the four stages, the changes slightly increased or decreased but were not obvious. So, the soil-keeping ability in mid-range areas has been improved since 2000. High soil holds in the downstream area were few. SC in the first decade showed an upward trend and declined in the next decade. Some soil-hold areas previously of a general level have degraded to that of low level.  Incorporating the analysis of land use, the expanded utilization of cropland and forest land has brought an increase in SC in the upper and the middle regions. However, the main body to provide SC supply was grassland throughout the whole basin. The areas decrease in grassland caused the decline of SC. At the same time, the growth of cropland and the reduction in unused land compensated the loss of SC caused by the area drop of grassland to some extent. So, in terms of total volume, the soil holds in the study area showed a slight upward trend from 2000 to 2015, followed by a downward trend from 2015 to 2020.

TGE of Different ESs
In order to analyze the TGE of different ESs, 2015-HRB was selected to conduct this practice.

Terrain Niche Index of HRB
The slope and the elevation data were extracted from DEM, and the topographic index of HRB was calculated by the terrain niche index formula-14 ( Figure 9). It was divided into six levels: one (0.179-0.305), two (0.305-0.463), three (0.463-0.709), four (0.709-1.005), five (1.005 -1.282), and six (1.282-1.786). The topography was dominated by one and two levels, the area of which accounts for 70.96% of the whole basin. As seen from the spatial distribution (Figure 10), the topographic index was low in the south (east) and high in the north (west). In other words, the high-grade topographic level was concentrated in the upper and middle reaches of the basin, while low-and medium-level terrain areas were heavily distributed in the north of midstream region and throughout the downstream region.

TGE of HQ
The HQ was reclassified into four levels: one (0-0.3), two (0.3-0.5), three (0.5-0.8), and four (0.8-1), which correspond to poor, general, high, and higher grades (see Section 3.1.1), respectively. The scatter plot was obtained by the statistical re-regression analysis on the HQ and topographic index ( Figure 11). The HQ in the HRB has an upward trend with the growth of the topographic index, and there is a significant positive correlation between the two. The curve fitting (R²) is 0.441, indicating that the logarithmic function can effectively describe their relationship. Figure 11. Correlation between terrain index and HQ in HRB in 2015.
As can be seen from Figure 12, the majority of low-grade habitat mass was distributed in the low-grade topographic index in HRB in 2015. With the growth of the topographic level index, the low-quality habitat area gradually reduced with the area increase in high-quality habitats. Additionally, it can be clearly seen that the proportion of highgrade HQ significantly increased with the growth of the terrain niche index from the proportion of change matrix during 2000-2015 (Table 3). The dominant interval was the third HQ level for the fourth and fifth levels of the terrain index, while the other levels of terrain were the first level of HQ. Over the past 15 years, the HQ of second-and fourth-grade topography positions declined. By incorporating land-use changes, this appearance was caused by the decrease in grassland and the increase in built-up land and unused land. Conversely, the increase in HQ at other topographical levels was caused by the increase in cropland and the decrease in unused land. It is sufficient to show that the transformation of nature has a considerable impact on HQ. Therefore, under the impact of rapid urbanization on ecological environment in this basin, relevant departments should take some ecological measures to deal with this.

TGE of WY
The water depth in the basin was divided into four levels based on the raster data obtained from model (see Section 3.1.2). The scatter plot was obtained by the statistical reregression analysis of WY and the topographic index ( Figure 13). The WY in the HRB has an upward trend with the growth of the topographic index. There is a significant positive correlation between the two. The curve fitting (R²) is 0.566, which indicates that the logarithmic function can successfully describe their relationship. The overlaid spatial distribution of the terrain niche index and WY shows the distinct discrepancy in the distribution of water depth levels on the local shape level index in HRB in 2015 ( Figure 14). Similar to the distribution of HQ with the local level index, the WY distributed in the low-terrainlevel index areas significantly decreased. With the growth of the topographic level index, the depth of water production also gradually increased. There was little low-grade water depth in the areas of the sixth topographic level index.   (Table 4). WY in the low terrain level was dominated by the low depth of first and second levels. The water production depth increased to three or four levels in this basin when it raised to fifth and sixth topographic indexes. This is enough to show that WY is closely related to the topographic gradient. The reason for this is that land use in low-terrain areas is majorly unused land with low grassland cover and low water conservation. High-terrain areas mostly have with medium-and high-density grassland cover, and the climate is humid at higher altitudes, resulting in higher water production. Comparing the water depth of 2015 with that of 2000, a decline in the water depth of the low-grade topographic index and the increase in water-producing depth of medium-and high-level topographic index areas was found. In addition to the differential rainfall, this difference may be due to the increasing area of grassland and forest land cover in high-altitude and high-terrain areas. The ecological environment in the unused land of the desert-dominated region continues to deteriorate in the north due to the continuous and extensive human activities. This needs to arouse people's vigilance with regard to, e.g., reducing human extrusion of original animal and plant living space, improving bad ecological environments, and maintaining good surroundings together.

TGE on SC
The scatter plot was obtained by the statistical re-regression analysis of SC and the topographic index ( Figure 15). The SC in the HRB has an upward trend with the increase in the topographic index, and there is a significant positive correlation between the two. The curve fitting (R²) is 0.687, indicating that the logarithmic function can properly describe their relationship.

Figure 15. Correlation between terrain index and SC in 2015.
As seen in Figure 16, it was in a poor condition for the whole basin in 2015. However, it can also be seen that the TGE of SC is similar to that of HQ and WY, although the SC in different locations is primarily of a low grade. The SC gradually became better with the increase in the terrain niche index. Furthermore, in order to analyze the proportions of each part, the SC was divided into four levels using the natural break point method (Table  5) and was then overlaid with the terrain niche index. Most of the low-grade areas transformed into high-grade areas and the proportion decreased gradually with the increase in the terrain niche index. This is mainly due to the concentration of high-terrain zones in the upper reaches of the basin, which have large areas of forest land, grassland and cropland, as well as dense vegetation and thick and fertile soil. The downstream area was dominated by desert ecosystems, with sparse vegetation and poor soil. The area proportion of low soil levels in various shapes decreased slightly, while that of high-grade soil levels increased. However, soil conditions in the basin remained poor and large amounts of soil were at risk of erosion, particularly in the southern mountainous areas and unused land in the north. Administrators should develop reasonable management measures. Only by strengthening the protection of woodland and grassland, planting trees and expanding vegetation cover can we reduce the threat of soil erosion and decrease the harm.

Discussion
In this study, we evaluated the changes in ES patterns in the HRB from 2000 to 2020 by using the three sub-modules in the InVEST model. The spatial variations of HQ, WY and SC were revealed for the past 20 years. We also analyzed the distribution characteristics of three ESs from the perspective of topography. This can provide the scientific basis for optimizing regional ecosystem patterns and land-use planning and management. The results are of great significance for rational planning of land use and constructing ecological civilizations.

Validation of Results with Previous Works
In general, the overall trends of the spatiotemporal changes of ESs are the same compared with the similar studies in the region [52,57]. Therefore, the reliable results can objectively mirror the spatial differentiation characteristics of ESs in the HRB. Otherwise, elevation and slope were used to analyze the changes in ESs on different terrain gradients and the results illustrate different ESs showing obvious TGE in the HRB. The distribution pattern of ESs is generally high in the mountainous areas and low in the plain, which is consistent with the relevant research [60,61]. In order to further verify the accuracy and reliability of the results, some specific cross-validations were conducted between this paper and previous studies [62,63], including the values of different ESs and spatial distributions ( Table 6). The range of HQ in this result is the same with the referenced work. It shows complete consistency with the distribution of high values and low values. The WY value in our study is slightly higher than that in the referenced work, but the difference is within a reasonable range. The major spatial differentiation is that the high value in our work is continuously distributed in the south of middle reaches, while it is dispersed with fragments in the referenced study. This may be due to different methods used to calculate the potential evapotranspiration in the two studies. Additionally, owing to lacking SC in previous studies in HRB, research from 2011 was selected to compare with our results in 2010. The result of this study is numerically low. The spatial difference lies in the shortage of high-value areas for SC in the northern basin and the distribution is not obvious along the riverside. This may be the factor parameter difference in the biophysical coefficient table. In addition, different formulas used to calculate the K factor in the two studies may also cause a discrepancy in the final outcomes. Table 6. Comparison between the results in this paper and those in the referenced work.

ESs Supply Results in this Paper
Year Area Range Spatial Distribution

Driving Forces of ESs Variation
The changes of ESs are affected by natural and human factors. Climate change is the main driving issue for natural factors. Furthermore, human factors, such as activity intensity on different topographic gradient will also have an important impact on ESs [64]. We qualitatively analyze the natural and human factors of ES changes in the HRB in terms of temperature, precipitation and human activities.
Research indicated that the temperature in the HRB shows a significant upward trend in the recent decades [65]. The increase rate of temperature in the study area is much higher than that of the average in the northern hemisphere during the same period. In particular, a large increase in temperature around 2010 was rare [66]. Rising temperature will promote the conversion of some unused land covered by glaciers and snow into grassland and grassland into woodland in the basin. The general trend of precipitation is increasing but the magnitude is negligible. The change is not obvious. The precipitation in this basin is not enough to balance the increase in evapotranspiration caused by the rising temperature. This gave rise to a further reduce in surface water, while the glacier melt water alleviated the decrease in runoff [66]. So, runoff changes little as temperature rises. The water supply gradually decreased, and some rivers and lakes dried up. The water area shrunk and desertification increased significantly [67].
Human activity is another essential driving force for land-use change in the watersheds. As be seen from Table 7, urban and rural industrial and other residential land continue to expand with the gradual deepening of human influence over the past 20 years. A large amount of unused land and grassland were reclaimed as arable land due to more intensive human activities in the middle reaches, which causes the area growth of cropland significantly. On the other hand, the manual allocation of water resources has given rise to a large amount of water diversion to the downstream region of the HRB. As a result, some unused land has been restored to grasslands. The ecological environment has been remarkably improved. A series of ecological management measures, such as natural forest and grass fence enclosures, have provoked the area increase in forest land in the upstream and the midstream areas. However, the growth of the population has brought about an increase in water demand. Then, adding up the decrease in rainfall, forest land and grassland shows a degradation trend in some areas. Therefore, the land-use change in the watersheds is a result of the joint power of human activities and climate change [37]. Furthermore, the land-use change will have an impact on the variation of different ESs.

Suggestions for Sustainable Improvement of ESs
The implementation of ecological protection policies and effective management in the Qilian Mountain reserve have enhanced the grassland coverage and water conservation capacity in the upstream area over the years. The Forestry Bureau should continue to strengthen the protection and surveillance of forest and grassland [68]. In the upper mountain region, climate change is always a dominant element for water resource allocation. Furthermore, except for water conservation, the increase in forest and grassland will also play an important role in adjusting runoff in the form of transpiration interception due to the warming and drying trend in recent years [69,70]. Grass planting and forestation measures should also be taken under the guidance of the Forestry Bureau within a scientific and reasonable range to stabilize runoff according to local conditions. In addition, thoughtful reservoir construction by the Water Resources Bureau is also useful in the regulation of seasonal runoff and in the efficient implementation of the ecological water diversion project.
The optimization and adjustment of the agricultural planting structure have been implemented in the midstream area in recent years, which yielded preliminary results for the construction of a water-saving society. However, the rapid development of the agricultural economy has brought the expansion of the irrigation area, which has caused the increase in water consumption, offsetting the consequences of saving water [70]. According to the analysis of long-term climate change in the basin, runoff may turn into the dry season in the future [71]. The implementation of the water transfer project may not be enough to support the rapid growth of population and fast development of the agricultural economy in the middle reaches [69][70][71][72]. In order to adapt to the changes in water resources in various stages, it is necessary to control the cropland expansion in the midstream region as soon as possible and the Bureau of Land Management needs to adjust the industrial structure. Moreover, during the implementation of the Heihe River "97" water diversion scheme, the continuous expansion of the oasis has begun to consume a large amount of groundwater and the ecological environment has a tendency of degradation [70]. It is essential to strictly limit the exploitation of groundwater, make zoning planning and continue to implement ecological protection projects in fragile ecological environment areas by the Natural Resources Planning Bureau. This basin needs to develop urbanization without destroying natural environment and follow the principles of sustainable development.
Since the implementation of the Heihe River "97" water diversion scheme in 2000, the amount of water flowing into the downstream region has increased significantly. The water area gradually recovered and steadily expanded, and the growth of vegetation improved. This alleviated the deterioration of the ecological environment in the downstream area, remarkably restored the ecological environment and initially realized the control goal [69,73]. However, the continuous expansion of cropland and construction land still brings great pressure on the environment. So, to improve ecological conditions, the Environmental Protection Agency should strictly control overgrazing and deforestation to allow the self-healing function of the ecosystem to take place. According to the local conditions, planting trees and grass can alleviate the ecological deterioration of desert. The government should further optimize and adjust the amount of water allocation and conduct a more reasonable allocation of resources for the sustainable development of the basin.

Limitations and Future Study
Although some reliable results have been obtained in this study, we recognize some limitations of this research. First of all, more attention should be paid in future studies. For example, only three key services provided by the ecosystems in the HRB were considered in our study due to the urgency to solve eco-environmental problems, such as the deterioration of the ecological environment, shortage of water resources and soil erosion. At the same time, the eco-environmental characteristics and data availability of the research area were also taken into consideration. However, the temporal and spatial changes of other ESs in the basin cannot be ignored, and they can be further discussed in future research. Then, the terrain niche index was used to comprehensively mirror the effect of topography on ESs. However, the respective weights of elevation and slope in the impact is not clear, and this needs to be further explored so as to formulate more practical ecological protection policies. Additionally, distinct operation mechanisms of different models may result in dissimilar outputs. It is necessary to focus on the comparative verification of the evaluated results between various models in future work.

Conclusions
In this study, three key ESs (e.g., HQ, WY and SC) were evaluated using the InVEST model in HRB from 2000 to 2020, and the TGE of different ESs was analyzed based on the terrain niche index. First of all, HQ and SC in the basin improved overall, while WY was in continuous decline. The values of various ESs changed little and the spatial distribution of them always presented similar discrepancy laws. The distribution of ESs in the basin shows lower values in the north and higher values in the south, high values at upstream areas and low values at downstream area. Furthermore, the correlations between different ESs and topography are strong, showing obvious TGE. HQ, WY and SC increased with the growth of the topographic level index. A high-quality habitat has high rates of water production. A low-quality habitat has a low yield of water. High SC is mostly distributed in the mountainous hilly areas with high grassland coverage, while a low SC is related to the accumulation of built-up land and unused land.
We analyze the variation of different ESs on the two scales: watersheds and substream. Relevant suggestions should be adopted for the distinct changes of ESs on different scales. In addition, topography and land use were combined to analyze the ES variation; in particular, the terrain niche index was specifically used to analyze the influence of TGE on different ESs. The results provide useful information for policy proposals to comprehensively consider the effect of terrain factors on the spatial differences of Ess and its drivers. The method provides a new perspective for studying the spatial differentiation of various ecological problems in the inland watersheds and can be applied to other similar areas, and then be incorporated into the formulation of watershed ecological zoning management policies.