Using the InVEST Model to Assess the Impacts of Climate and Land Use Changes on Water Yield in the Upstream Regions of the Shule River Basin

Water yield is a key ecosystem function index, directly impacting the sustainable development of the basin economy and ecosystem. Climate and land use/land cover (LULC) changes are the main driving factors affecting water yield. In the context of global climate change, assessing the impacts of climate and LULC changes on water yield in the alpine regions of the Qinghai–Tibet Plateau (QTP) is essential for formulating rational management and development strategies for water resources. On the basis of the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model, we simulated and analyzed the spatiotemporal variations and the impacts of LULC and climate changes on water yield from 2001 to 2019 in the upstream regions of the Shule River Basin (USRB) on the northeastern margin of the QTP. Three scenarios were designed in the InVEST model to clearly analyze the contributions of climate and LULC changes on the variation of water yield. The first scenario integrated climate and LULC change into the model according to the actual conditions. The second scenario was simulation without LULC change, and the third scenario was without climate change. The results showed that (1) the InVEST model had a good performance in estimating water yield (coefficient of determination (R2) = 0.986; root mean square error (RMSE) = 3.012, p < 0.05); (2) the water yield significantly increased in the temporal scale from 2001 to 2019, especially in the high altitude of the marginal regions (accounting for 32.01%), while the northwest regions significantly decreased and accounted for only 8.39% (p < 0.05); (3) the spatial distribution of water yield increased from the middle low-altitude regions to the marginal high-altitude regions; and (4) through the analysis of the three scenarios, the impact of climate change on water yield was 90.56%, while that of LULC change was only 9.44%. This study reveals that climate warming has a positive impact on water yield, which will provide valuable references for the integrated assessment and management of water resources in the Shule River Basin.


Introduction
Ecosystem service function refers to the natural environmental conditions and utilities that the ecosystem forms and maintains to promote human survival and development [1,2], which include supply, regulation, support, and culture according better manage the water resources in the Shule River Basin have become urgent issues to be studied.
This study, using the USRB as the study area, evaluated and analyzed the spatiotemporal variations of water yield in 2001-2019 and the impacts of climate and LULC changes on variations of water yield. Specifically, the hypotheses of this study were: (1) the InVEST model has strong applicability in the alpine regions; (2) the water yield significantly increased in the temporal scale from 2001 to 2019, and the spatial distributions of water yield are basically consistent; and (3) the climate change has a positive impact on water yield, while LULC change has a negative impact on water yield for the analysis period. The results may provide scientific reference for quantitative assessment, effective management, and sustainable development of water resources in the Shule River Basin.

Study Area
The USRB (96 • 37 12"-98 • 59 24" E, 38 • 13 12"-39 • 52 12" N) is the formation region of the main stream of the Shule River Basin and covers an area of approximately 10,973.9 km 2 . It generates about 61.8% of the river flows of the entire basin, which makes significant sense on the social development of the midstream and maintains the eco-environment balance of the downstream (Figure 1) [39]. The elevation of the study area ranges from 2031 to 5763 m above sea level and increases gradually from the middle to the marginal region, with a mean elevation of 3945 m. The water outlet of the USRB is monitored by the Changmapu Hydrological Station (96 • 51 E, 39 • 49 N) [40]. The study area has a typical continental climate [41], which experiences a warm-wet climate in summer and a cold-dry climate in winter. The annual precipitation is 325.93 mm, and the mean annual temperature and potential evapotranspiration are −6.06 • C and 1249.42 mm, respectively. According to the elevation from high to low, the landscape follows a distinct vertical variation and comprises snow and ice, meadow, shrub, steppe, and desert [42]. Moreover, settlements in the study area are few and scattered, and grazing is the main industry.
Water 2021, 13,1250 3 of 20 farmers and herdsmen in the Hexi Corridor region of Northwest China [37,38]. Under the influence of the climate and LULC changes, how the water yield changes spatiotemporally and how to better manage the water resources in the Shule River Basin have become urgent issues to be studied. This study, using the USRB as the study area, evaluated and analyzed the spatiotemporal variations of water yield in 2001-2019 and the impacts of climate and LULC changes on variations of water yield. Specifically, the hypotheses of this study were: (1) the InVEST model has strong applicability in the alpine regions; (2) the water yield significantly increased in the temporal scale from 2001 to 2019, and the spatial distributions of water yield are basically consistent; and (3) the climate change has a positive impact on water yield, while LULC change has a negative impact on water yield for the analysis period. The results may provide scientific reference for quantitative assessment, effective management, and sustainable development of water resources in the Shule River Basin.

Study Area
The USRB (96°37′12″-98°59′24″ E, 38°13′12″-39°52′12″ N) is the formation region of the main stream of the Shule River Basin and covers an area of approximately 10,973.9 km 2 . It generates about 61.8% of the river flows of the entire basin, which makes significant sense on the social development of the midstream and maintains the eco-environment balance of the downstream (Figure 1) [39]. The elevation of the study area ranges from 2031 to 5763 m above sea level and increases gradually from the middle to the marginal region, with a mean elevation of 3945 m. The water outlet of the USRB is monitored by the Changmapu Hydrological Station (96°51′ E, 39°49′ N) [40]. The study area has a typical continental climate [41], which experiences a warm-wet climate in summer and a colddry climate in winter. The annual precipitation is 325.93 mm, and the mean annual temperature and potential evapotranspiration are −6.06 °C and 1249.42 mm, respectively. According to the elevation from high to low, the landscape follows a distinct vertical variation and comprises snow and ice, meadow, shrub, steppe, and desert [42]. Moreover, settlements in the study area are few and scattered, and grazing is the main industry.

Methods
The Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) LULC Type (MCD12Q1) data product provides global LULC data at yearly intervals (2001-2019), with a spatial resolution of 500 m. Hence, the analysis period was 2001-2019 in this study. Meanwhile, the InVEST model requires the input data to have the same spatial resolution, so meteorological data were resampled and reprojected according to the resolution and projected coordinate system of LULC. In addition, the streamflow data of the Changmapu Hydrological Station were converted to the water yield with the unit of mm based on the size of the study area and compared with the values simulated by the InVEST model. Finally, the impacts of climate and LULC changes on water yield were assessed using the InVEST model for three scenarios. A workflow of the technique used in this study is given in Figure 2.

Methods
The Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) LULC Type (MCD12Q1) data product provides global LULC data at yearly intervals (2001-2019), with a spatial resolution of 500 m. Hence, the analysis period was 2001-2019 in this study. Meanwhile, the InVEST model requires the input data to have the same spatial resolution, so meteorological data were resampled and reprojected according to the resolution and projected coordinate system of LULC. In addition, the streamflow data of the Changmapu Hydrological Station were converted to the water yield with the unit of mm based on the size of the study area and compared with the values simulated by the InVEST model. Finally, the impacts of climate and LULC changes on water yield were assessed using the InVEST model for three scenarios. A workflow of the technique used in this study is given in Figure 2.

Water Yield Module
Water yield in the InVEST is defined as the amount of water that runs off the landscape to a sub-watershed in a certain period [2]. The water yield module in the InVEST model is built on the annual average precipitation and the Budyko curve [43]. Based on the principle of water balance, the model estimates the annual water yield for each pixel of the study catchment as the annual precipitation minus the annual actual evapotranspiration (AET) [28]. The model algorithm is generated with Equation (1): where is the water yield for land use type j on pixel , is the annual actual evapotranspiration for pixel , and is the annual precipitation on pixel . The InVEST approach relates AET to potential evapotranspiration (PET) using Equation (2), which was developed by Budyko and later adapted by Fu and Zhang et al. [44][45][46].

Water Yield Module
Water yield in the InVEST is defined as the amount of water that runs off the landscape to a sub-watershed in a certain period [2]. The water yield module in the InVEST model is built on the annual average precipitation and the Budyko curve [43]. Based on the principle of water balance, the model estimates the annual water yield for each pixel of the study catchment as the annual precipitation minus the annual actual evapotranspiration (AET) [28]. The model algorithm is generated with Equation (1): where Y xj is the water yield for land use type j on pixel x, AET xj is the annual actual evapotranspiration for pixel x, and P x is the annual precipitation on pixel x. The InVEST approach relates AET to potential evapotranspiration (PET) using Equation (2), which was developed by Budyko and later adapted by Fu and Zhang et al. [44][45][46]. where PET xj is the potential evapotranspiration for land use type j on pixel x, ω is related to the plant available water content, precipitation, and the Zhang parameter (Z,) which can be calculated by Equation (3) as given below: where Z is an empirical constant, which captures the local precipitation pattern and additional hydrogeological characteristics, and AWC x is the average annual values of available water capacity. A more detailed description of the water yield module is referred to in the InVEST model's user guide [47]. Both AET and PET can reflect the evapotranspiration capacity of the surface ecosystem, but there is a difference between them (Table 1).

Data Sources and Processing
The water yield module requires some biophysical parameters as basic data to compute water yield, including LULC, reference evapotranspiration, annual precipitation, plant available water content (PAWC), root restricting layer depth, biophysical table, seasonality factor (Zhang parameter), and sub-watersheds. All of the data were resampled at a spatial resolution of 500 m and projected using the World Geodetic System 1984. The sources of these basic data are summarized in Table 2. LULC data were provided by MCD12Q1, which describes LULC based on one-year terra and aqua observations, with a spatial resolution of 500 m [49]. MCD12Q1 adopts five different land cover classification schemes, and the main technology of information extraction is the classification of a supervised decision tree. The International Geosphere Biosphere Program (IGBP) classification method was adopted and mainly divided into 17 categories in this study. The LULC in the study area consists of woodland, cropland, grassland, built-up lands, crops-natural vegetation transition, permanent snow and ice, and barren land.

Precipitation and Reference Evapotranspiration
Original precipitation and temperature data from meteorological stations around the China were interpolated into grid data by the Australian National University Splines (ANUSPLIN, Australian National University, Canberra, Australia.) package, which refers to the research of Lian et al. [11]. Then, based on the vector boundary of the study area, the needed data were cut from the interpolated meteorological grid. Yang et al. [28] has shown that the Hargreaves equation generates improved results over the Penman-Monteith model, so the annual ET 0 was calculated in this study using the Hargreaves equation. The Hargreaves equation is shown in Equation (4): where ET 0 is reference evapotranspiration, T max and T min are the daily maximum and minimum temperature, respectively, R a is extraterrestrial radiation that is estimated from the literature, and C, E, and T are the empirical parameters that refer to Hu's correction of the QTP [50].

PAWC and Soil Depth
The soil depth gridded map was generated based on the second soil survey database downloaded from the Cold and Arid Region Science Data Center at Lanzhou ( Table 2). The PAWC refers to the difference between the field capacity and the wilting point [13], which can be estimated according to the physical and chemical properties of the soil [51]. The formula is expressed in Equation (5): where Sand, Silt, Clay, and OM are the proportion of clay, sand, silt, and organic matter in the soil.

Sub-Watershed and Biophysical Table
The hydrological analysis tool in ArcGIS 10.4 (Environmental Systems Research Institute (ESRI), Redlands, CA, USA) was used to generate the sub-watershed based on the digital elevation model. The biophysical table was mainly used to reflect the attributes of soil cover and LULC, including the LULC code, plant evapotranspiration coefficient (Kc), and root depth. The biophysical coefficients of each LULC type used in the model can be found in the literature and the InVEST model's user guide [28,48].

The Zhang Parameter
The Zhang parameter is a climate seasonality factor, which indicates the local precipitation pattern and the hydrogeological characteristics, with values varying from 1 to 30. Previous studies have explored the Zhang parameter of the QTP and found that the simulated water yields of the InVEST model were closer to the observed values when the Zhang parameter was 3.33 [52,53]. In this study, the error between the water yields simulated by the InVEST model and the water flows of the Changmapu Hydrological Station are smaller (1.11%) when the Zhang parameter is equal to 3.33.

Climate Change and LULC Change Scenarios
We estimated the water yield in the USRB from 2001 to 2019 based on the InVEST model and analyzed the impacts of climate and LULC changes on the variations of water yield. To clearly analyze these impacts, three scenarios were designed: actual conditions, actual conditions without LULC change, and actual conditions without climate change. Under the actual scenario, LULC change and climate change were input in the model in accordance with the actual conditions (Table 3). Under the scenario without climate change, the precipitation and ET 0 in 2019 were controlled to be unchanged as 2001 and the water yields for 2001 and 2019 were calculated to assess the effect of LULC change. Under the scenario without LULC change, the LULC conditions in 2019 were assumed to be the same as those in 2001, and water yields were calculated to evaluate the influence of climate change. Finally, the three scenarios were compared to reveal the effects of climate and LULC change on water yield. Table 3. Simulations of water yield variation by scenario.

Conditions without
Climate Change

Conditions without Land Use Change
According to the change of water yield under different scenarios, the contribution of climate and LULC changes to the variability of water yield can be quantified by Equations (6) and (7) [2]: where G c is the contribution of climate rate to change in water yield under the scenario without LULC change, G L is the contribution rate of LULC to change water yield under the scenario without climate change, C represents the difference in mean annual water yield between 2001 and 2019 under the scenario without LULC change, and L represents the difference in mean annual water yield between 2001 and 2019 under the scenario without climate change.

Analysis of Spatiotemporal Variations of Research Elements
Linear regression analysis and the nonparametric Mann-Kendall test were applied to evaluate the spatial variation of water yield in this study and were conducted using MATLAB R2016a (MathWorks, Natick, MA, USA). The slope of the linear regression was the superiority index that quantified the variation of water yield in the study period. The slope is calculated by Equation (8): where n is the number of study years (19 in this study), i is the serial number of the year, wyild i is the water yield value in the year i in the grid, and slope is the trend of the water yield-when slope > 0, the trend is positive over the 19 years, and slope < 0 indicates that the trend is negative. In order to determine the variation trends of water yield and quantify the statistical significance of each pixel, the nonparametric Mann-Kendall test was also applied in this study. Generally, the values of Kendall inclination (β) and Z statistics were used to estimate the trend of variation [54,55]. β is an unbiased estimate for the trend: when β > 0, the water yield shows a significant upward trend; when β = 0, the water yield shows no significant trend; and when β < 0, the water yield shows a significant downward trend [55]. Based on the Z statistics, when Z > 1.96, the results significantly increase; when Z < −1.96, the results significantly decrease; and when −1.96 ≤ Z ≤ 1.96, the results have nonsignificant changes [54]. In addition, the linear fitting in OriginPro 9.1 (OriginLab, Northampton, MA, USA) was used to analyze the temporal variation of research elements (including air temperature, precipitation, actual evapotranspiration, and water yield) from 2001 to 2019.

Performance Assessment
The performance of the InVEST model can be evaluated by comparing the simulated water yield with the observed streamflow data [56], and the performance can be determined by the RMSE and R 2 . The R 2 (ranging from 0 to 1) represents the proportion of the variance in measured data explained by the model [57]. A higher value of R 2 indicates less error variance, and a value greater than 0.5 is considered acceptable [58]. RMSE is one of the most commonly used error index statistics, and the lower the RMSE the better the model performance [57]. In this study, linear regression analysis was used to compared estimated water yield against the observed data. The P, R 2 (all R 2 in the text represent the adjusted R 2 ), and RMSE were calculated in this process. Besides that, the slope differences between the linear regression and the 1:1 line were analyzed by the General Linear Model in IBM SPSS Statistics 22.0 (IBM, Beijing, China.).

USRB Climatic Change
The USRB is characterized by a dry and cold climate with the mean annual air temperature and annual precipitation of −6.06 • C and 325.96 mm, respectively. In terms of temporal variation, both mean annual air temperature and annual precipitation increased annually by averages of 0.02 • C·year −1 (y = 0.02x − 41.10, p > 0.05) and 5.92 mm·year −1 (y = 5.92x -11,570.50, p < 0.05), respectively, over the course of our study period, and the climate conditions of the study area were characterized by a trend of warming and wetting ( Figure 3). Under the influence of air temperature and precipitation, annual actual evapotranspiration also increased annually by an average of 3.00 mm·year −1 (y = 3.00x − 5857.00, p > 0.05). In terms of spatial variation, mean annual temperature varied spatially between −17.1 and 5.4 • C, increasing along a gradient from southeast-tonorthwest and decreasing with the rise of elevation ( Figure 4A), while annual precipitation varied between 76.60 and 628.40 mm, decreasing from southeast-to-northwest and increasing with the rise of elevation ( Figure 4B). The annual actual evapotranspiration varied spatially from 0.00 to 515.30 mm; the higher values were found in the east of the northwest and low-altitude regions of the southeast, while lower values were observed in the west of the northwest and marginal regions of the southeast of the whole USRB ( Figure 4C).

LULC Change
It can be seen from the spatial distribution of LULC that the area of grassland coverage is the highest, accounting for approximately 51.72% of the USRB, being mostly distributed in the east of the northwest and low-altitude regions of the southeast ( Figure 5). The distribution pattern of the grassland is consistent with the high-value regions of annual actual evapotranspiration. The second most common LULC is barren land, encompassing more than 31.95% of the total area and is mainly distributed in the west of the northwest and marginal regions of the southeast. The areas of permanent wetland, cropland, built-up land, crops-natural vegetation transition, and permanent snow and ice accounted for about 2.91%, 2.58%, 2.54%, 3.79%, and 4.51%, respectively. Permanent snow and ice are mainly distributed in the marginal high-altitude regions; permanent wetland, cropland, and cropsnatural vegetation transition are mainly distributed in the northwest, while urban and built-up land are scattered in the study area. The spatial distribution of LULC changed little from 2001 to 2010, and the major changes happened mainly from 2010 to 2019.
In Table 4, the rows display the areas of the seven LULC types in 2019, whereas the columns display the area in 2001. Hence, Table 4 indicates the area of the LULC types that experiences a transition from one type to another type between 2001 and 2019. The main diagonal elements indicate the area of LULC types that show persistence of LULC types in 2019, and off diagonal entries represent a transition from one LULC type in 2001 to another LULC type in 2019. Data show that the total conversion area within the USRB was 2127.53 km 2 from 2001 to 2019, encompassing 19.42% of the total area. Barren land, grassland, and permanent snow and ice were converted in the highest amounts. A further detailed analysis of LULC type changes can reveal more important information, especially in relation to barren land. There was a decline in barren land (1039.73 km 2 ), which was attributed mainly to the transitions from barren land to grassland and permanent snow and ice. The percentage of newly converted grassland derived from barren land was 37.04%, which was mainly distributed in the northwest of the USRB, and 21.37% of the newly converted permanent snow and ice was derived from barren land, which mostly occurred in the marginal high-altitude regions of the study area. There was also a small transition of grassland (15.28 km 2 ) and permanent snow and ice (10.05 km 2 ) to barren land, while the surplus transition included permanent snow and ice to grassland (134.46 km 2 ) and grassland to permanent snow and ice (39.83 km 2 ). Overall, except for barren land, the areas of all other LULC types showed increasing trends from 2001 to 2019. Statistically, the area of grasslands had the largest proportional gain (77.55%), followed by permanent snow and ice (10.33%).

LULC Change
It can be seen from the spatial distribution of LULC that the area of grassland coverage is the highest, accounting for approximately 51.72% of the USRB, being mostly distributed in the east of the northwest and low-altitude regions of the southeast ( Figure  5). The distribution pattern of the grassland is consistent with the high-value regions of annual actual evapotranspiration. The second most common LULC is barren land, encompassing more than 31.95% of the total area and is mainly distributed in the west of the northwest and marginal regions of the southeast. The areas of permanent wetland, cropland, built-up land, crops-natural vegetation transition, and permanent snow and  In Table 4, the rows display the areas of the seven LULC types in 2019, whereas the columns display the area in 2001. Hence, Table 4 indicates the area of the LULC types that experiences a transition from one type to another type between 2001 and 2019. The main diagonal elements indicate the area of LULC types that show persistence of LULC types in 2019, and off diagonal entries represent a transition from one LULC type in 2001 to another LULC type in 2019. Data show that the total conversion area within the USRB was 2127.53 km 2 from 2001 to 2019, encompassing 19.42% of the total area. Barren land, grassland, and permanent snow and ice were converted in the highest amounts. A further detailed analysis of LULC type changes can reveal more important information, especially in relation to barren land. There was a decline in barren land (1039.73 km 2 ), which was attributed mainly to the transitions from barren land to grassland and permanent snow and ice. The percentage of newly converted grassland derived from barren land was 37.04%, which was mainly distributed in the northwest of the USRB, and 21.37% of the newly converted permanent snow and ice was derived from barren land, which mostly occurred in the marginal high-altitude regions of the study area. There was also a small transition of grassland (15.28 km 2 ) and permanent snow and ice (10.05 km 2 ) to barren land, while the surplus transition included permanent snow and ice to grassland (134.46 km 2 ) and grassland to permanent snow and ice (39.83 km 2 ). Overall, except for barren land, the areas of all other LULC types showed increasing trends from 2001 to 2019. Statistically, the area of grasslands had the largest proportional gain (77.55%), followed by permanent snow and ice (10.33%).

Model Validation
To evaluate the performance of the InVEST water yield model, we analyzed the relationship between the InVEST simulated mean annual water yield and gauged mean annual water yield. The results showed that the points were symmetrically and dispersedly distributed on both sides of the 1:1 line. Besides that, there was a strong linear relationship between the simulated mean annual water yield and the observed corresponding values (y = 0.962x + 3.377, R 2 = 0.986, p < 0.05, RMSE = 3.012) ( Figure 6). Although the intercept indicated that the InVEST model underestimated by 3.38 mm per year, the slope of linear regression was not significantly different from one. annual water yield. The results showed that the points were symmetrically and dispersedly distributed on both sides of the 1:1 line. Besides that, there was a strong linear relationship between the simulated mean annual water yield and the observed corresponding values (y = 0.962x + 3.377, R 2 = 0.986, p < 0.05, RMSE = 3.012) (Figure 6). Although the intercept indicated that the InVEST model underestimated by 3.38 mm per year, the slope of linear regression was not significantly different from one. Simulated mean annual water yield (mm) Observed mean annual water yield (mm) Figure 6. InVEST modeled mean annual water yield against gauged mean annual water yield. Red line indicates a relationship with intercept = zero and slope = one.

Temporal Variation of Water Yield
We focused on the temporal variations of water yield in 2001-2019, and the results are given in Figure 7. Based on the slope of the water yield trend line, the mean annual water yield increased significantly in the USRB, with an average rate of 2.36 mm·year −1 (y = 2.361x − 4621.892, p < 0.05).

Temporal Variation of Water Yield
We focused on the temporal variations of water yield in 2001-2019, and the results are given in Figure 7. Based on the slope of the water yield trend line, the mean annual water yield increased significantly in the USRB, with an average rate of 2.36 mm·year −1 (y = 2.361x − 4621.892, p < 0.05).

Spatial Distribution of Water Yield and its Dynamic Change
In terms of geographical distribution (Figure 8), the spatial pattern of water yield showed strong heterogeneity across the USRB, and the spatial patterns of the high and low values of water yield distribution were consistent in the whole study period. The

Spatial Distribution of Water Yield and Its Dynamic Change
In terms of geographical distribution (Figure 8), the spatial pattern of water yield showed strong heterogeneity across the USRB, and the spatial patterns of the high and low values of water yield distribution were consistent in the whole study period. The distribution of water yield seems to be related to the elevation variation of the USRB. In detail, the water yield values were found to be higher at the marginal high-altitude regions and lower at the low-altitude regions; the lowest values were located at the low-altitude regions of the southeast and east of the northwest.

Spatial Distribution of Water Yield and its Dynamic Change
In terms of geographical distribution (Figure 8), the spatial pattern of water yield showed strong heterogeneity across the USRB, and the spatial patterns of the high and low values of water yield distribution were consistent in the whole study period. The distribution of water yield seems to be related to the elevation variation of the USRB. In detail, the water yield values were found to be higher at the marginal high-altitude regions and lower at the low-altitude regions; the lowest values were located at the lowaltitude regions of the southeast and east of the northwest. Over the course of our study period, the variation trend of water yield to all the pixels is shown in Figure 9. More than 32.01% of the water yield in the study area experienced an extremely significant (p < 0.01) or significant (p < 0.05) increasing distribution largely in the high-altitude regions in the margin. Only 8.39% of the study area experienced an extremely significant (p < 0.01) or significant (p < 0.05) decrease, most of which was in the northwest of the study area. In contrast, the areas of nonsignificant changing of water yield were approximately 59.60%, which mainly occurred in the low-altitude regions of the southeast and parts of the northwest of the USRB.

Difference in Water Yield among LULC Types
There were obvious differences in mean water yield under different LULC types ( Figure 10A). Generally, the mean water yields from permanent snow and ice, barren land, and permanent wetland were the highest, reaching more than 270 mm; these were followed by crops-natural vegetation transition, which was up to 133.7 mm. In contrast, the mean water yields on grassland, croplands, and built-up land were relatively small, approximately 50.13, 79.17, and 42.80 mm, respectively. From 2001 to 2019, the mean water yields of grassland, permanent wetland, cropland, and barren land showed an increasing trend, built-up land and crops-natural vegetation transition fluctuated upward of the mean water yield, and permanent snow and ice decreased slightly. Similar variation trends also occurred in the comparison of the total water yield of different LULC types in 2001-2019.
Over the course of our study period, the variation trend of water yield to all the pixels is shown in Figure 9. More than 32.01% of the water yield in the study area experienced an extremely significant (p < 0.01) or significant (p < 0.05) increasing distribution largely in the high-altitude regions in the margin. Only 8.39% of the study area experienced an extremely significant (p < 0.01) or significant (p < 0.05) decrease, most of which was in the northwest of the study area. In contrast, the areas of nonsignificant changing of water yield were approximately 59.60%, which mainly occurred in the low-altitude regions of the southeast and parts of the northwest of the USRB.

Difference in Water Yield among LULC Types
There were obvious differences in mean water yield under different LULC types ( Figure 10A). Generally, the mean water yields from permanent snow and ice, barren land, and permanent wetland were the highest, reaching more than 270 mm; these were followed by crops-natural vegetation transition, which was up to 133.7 mm. In contrast, the mean water yields on grassland, croplands, and built-up land were relatively small, approximately 50.13, 79.17, and 42.80 mm, respectively. From 2001 to 2019, the mean water yields of grassland, permanent wetland, cropland, and barren land showed an increasing trend, built-up land and crops-natural vegetation transition fluctuated upward of the mean water yield, and permanent snow and ice decreased slightly. Similar variation trends also occurred in the comparison of the total water yield of different LULC types in 2001-2019. In terms of the total water yield ( Figure 10B), barren land produced the largest total water yield, reaching 10.80 × 10 8 m 3 , followed by grassland, permanent wetland, and permanent snow and ice, which were up to 2.88 × 10 8 , 1.03 × 10 8 , and 1.41 × 10 8 m 3 , respectively. In contrast, the total water yields of cropland, built-up land, and cropsnatural vegetation transition were relatively small, up to just 0.24 × 10 8 , 0.12 × 10 8 , and 0.38 × 10 8 m 3 , respectively.

Influence of LULC and Climate Changes on Water Yield
LULC and climate changes are important drivers to the variations of water yield in the USRB. The total water yield in the USRB increased by 58.75% between 2001 and 2019 in the actual scenario (Table 5). In the scenario without LULC change, climate change led to a 92.62% increase in regional total water yield, which was consistent with the actual situation. In the scenario without climate change, LULC conversions resulted in a 9.65% In terms of the total water yield ( Figure 10B), barren land produced the largest total water yield, reaching 10.80 × 10 8 m 3 , followed by grassland, permanent wetland, and permanent snow and ice, which were up to 2.88 × 10 8 , 1.03 × 10 8 , and 1.41 × 10 8 m 3 , respectively. In contrast, the total water yields of cropland, built-up land, and cropsnatural vegetation transition were relatively small, up to just 0.24 × 10 8 , 0.12 × 10 8 , and 0.38 × 10 8 m 3 , respectively.

Influence of LULC and Climate Changes on Water Yield
LULC and climate changes are important drivers to the variations of water yield in the USRB. The total water yield in the USRB increased by 58.75% between 2001 and 2019 in the actual scenario (Table 5). In the scenario without LULC change, climate change led to a 92.62% increase in regional total water yield, which was consistent with the actual situation. In the scenario without climate change, LULC conversions resulted in a 9.65% reduction in total water yield. Additionally, the actual conditions and two scenarios were compared, we found that the contribution of climate change to total water yield was 90.56%, while that of LULC change only accounted for 9.44%. It should be acknowledged that the impact of climate change was far bigger than that of LULC change. Meanwhile, consistent with the previous hypotheses, climate change had a positive impact on the variation of water yield, while LULC change had a negative impact on it.

Temporal Variation Characteristics of Water Yield and Its Influencing Factors
Water yield is a crucial component of the regulation function of the ecosystem in basins, which relates to the regional sustainable development of water resources and ecological security [6,11]. In recent years, drought conditions in the QTP, which is called the "Roof of the World" and the "Water Tower of Asia", have received a lot of attention due to climate warming [55]. In this context, some studies have explored the temporal variations of water yield on the QTP. Among them, the Three-River Headwaters Region is the largest ecological function region of water source supply and conservation on the QTP. Lü et al. [59] assessed the temporal variation of water yield in the Three-River Headwaters Region from 1980 to 2016 and showed an insignificant increasing trend, while Pan et al. [52] found that the water yield in the region decreased from 1980 to 2005. In addition, Qilian Mountain is an important region of water supply and a priority region for biodiversity protection in Northwest China [60]. Zhao et al. [61] focused on the water yield in the upstream regions of the Shiyang River Basin and showed a variation pattern of increase first, then a decrease, and then showed an increase in 1986-2015. Zhang et al. [62] explored the variation trend of water yield in the upstream regions of the Heihe River Basin and indicated a decreasing trend from 2001 to 2015. In this study, we found that water yield significantly increased in the USRB from 2001 to 2019, which was similar to the study of Lü et al. and was at odds with the studies of Pan, Zhao, and Zhang et al. The different outcomes often related to the difference in the climatic conditions of the study areas and the time scales of the studies. In the USRB, the annual precipitation and annual actual evapotranspiration increased at the rates of 5.92 and 3.00 mm·year −1 , respectively, from 2001 to 2019. Thus, the recharging of precipitation to water yield exceeded the consumption of actual evapotranspiration, resulting in a significant increase in water yield. Interestingly, the water yield of permanent snow and ice presented a decreasing trend in 2001-2019, which could be due to permanent snow and ice being more sensitive to climate warming. Throughout the whole study period, the annual precipitation and annual actual precipitation in the permanent snow and ice region increased at the rates of 8.77 and 12.02 mm·year −1 , respectively. So, the evapotranspiration of the water yield was greater than the recharging of precipitation to it, leading to the water yield decreasing in the permanent snow and ice region.
Through the three scenarios analysis, we found that the proportional contribution of climate change to the variation in water yield was 90.56%, while LULC accounted for 9.44%. Compared with LULC change, climate change played a more important role in affecting the regional variation of water yield, which was consistent with many previous studies [2,11,63]. There are some reasons for the smaller contribution of LULC change. Firstly, the size of LULC change was small and the process of that was complex. In addition, different LULC conversion patterns cause both positive and negative impacts on water yield [11]. In contrast, climate change could directly alter surface runoff and produce a significant impact on water yield [2].

Spatial Variation Characteristics of Water Yield and Its Influencing Factors
The spatial distribution pattern of water yield in the USRB was consistent and had little inter-annual change in 2001-2019, which was similar to the study of the Qinghai Lake Watershed [11]. The spatial pattern of water yield is directly related to the distribution characteristics of the regional climate elements and LULC, i.e., the region has a higher water yield with more precipitation and less actual evapotranspiration [64]. Under the comprehensive influence of climate factors and LULC, the water yield of the USRB presented that relatively higher values were distributed in the marginal high-altitude regions than in the low-altitude regions, with the southeast and east of the northwest having the lowest values. To be more specific, the high-altitude regions in the margin of the USRB have high precipitation and low evapotranspiration. Furthermore, the main LULC types are permanent ice and snow and barren land in high-altitude areas; the permanent ice and snow belongs to the water area and is most prone to forming runoff, while the barren land can only infiltrate a small part of the precipitation, and more precipitation directly forms the water yield [2,65]. However, the low-altitude regions have a relatively low water yield due to low precipitation and high evapotranspiration. In particular, the main LULC type in the low-altitude regions of the southeast and east of the northwest is grassland, which has the highest value of vegetation evapotranspiration and water infiltration, resulting in the lowest water yield [65].
In 2001-2019, the precipitation in the marginal high-altitude regions increased significantly, while the actual evapotranspiration increased insignificantly, with both leading to significant increases in the water yield in these regions. In addition, barren land in the northwest region has been largely replaced by grassland and cropland. Grassland and cropland have strong evapotranspiration compared to barren land and consume more water for plant growth than barren land. Hence, the water yield represented a significant decreasing trend in some northwest regions.

Uncertainties in Model-Based Assessment
Some uncertainties exist in the assessment of the water yield due to data limitations and the complexity of the model structure and parameters. Precipitation and reference evapotranspiration constitute the pivotal dataset to water yield simulation. However, due to the complexity of climate change or the sparse climate data, it is difficult to interpolate precipitation exactly and use the Hargreaves method to calculate reference evapotranspiration accurately, which largely affects the accuracy of the water yield simulation [13]. Hence, these reasons lead to data limitations. Additionally, the water yield is closely related to the development of the social economy and human activities; however, the input data of the model are natural data and the socio-economic-related data are rarely considered [2]. At the same time, the model represents biophysical processes in a simplified way; it assumes that all the water yield of a pixel reaches one point and does not distinguish the surface and the subsurface water [11]. In particular, it should be pointed out that the USRB is a typical high cold region and the processes of glaciers and permafrost are not considered by the InVEST model, which increases the uncertainty of the hydrological simulation. Thus, the model needs to be further improved. Despite these uncertainties, the results of our study can still reflect the temporal and spatial distribution patterns and variations of water yield and reflect the relationship between water yield and the changes of climate and LULC. It can provide a scientific reference for the utilization of water resources and the protection of the ecological environment around the Shule River Basin. In the future, it will be necessary to further expand the time scale of the study and evaluate the impacts of climate and LULC changes on water yield in the whole QTP.

Conclusions
Using the InVEST model, we assessed the variation of water yield in the USRB from 2001 to 2019 and analyzed its spatiotemporal responses to the changes of the climate and LULC. In terms of the temporal scale, the water yield increased significantly with the combined effects of climate and LULC changes in the study area, and the regional total water yield increased from 49.76 × 10 8 m 3 in 2001 to 79.00 × 10 8 m 3 in 2019, an increase of 58.75%. Additionally, the spatial distribution of water yield was heterogeneous, i.e., the high-altitude areas located in the margins of the region had relatively higher values than in the low-altitude areas. By varying the conditions, we simulated the water yield in different scenarios and analyzed the relative contributions of climate and LULC to water yield, indicating that climate change plays a dominant role in affecting water yield, while LULC change has a small impact. Generally, the proportional contribution of climate change to the variation in water yield was 90.56%, while LULC change accounted for 9.44%.
In conclusion, the InVEST model has great applicability and performance in the Shule River Basin. In addition, the results of this study will help to understand the impacts of climate warming and LULC change on water provisions and provide a foundation for the effective management of water resources and scientific strategies in alpine areas.