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

Abstract

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 (R²) = 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.

1. 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 to the Millennium Ecosystem Assessment [3]. Many ecosystem service functions are of great importance to human wellbeing, especially those related to water [4,5]. Water yield is one of the important indices of regulation function of the ecosystem in watersheds and plays a critical role in agriculture, industry, fisheries, domestic activities, etc. [6]. The water yield plays a crucial role in the sustainable development of these abovementioned sectors and directly impacts the regional economy and society [7,8]. Meanwhile, the spatiotemporal variation of water yield is also crucial for the allocation of water resources between regions [9]. Besides that, climate and LULC changes are considered to be the main factors affecting the spatiotemporal changes of water yield [10]. Hence, evaluation and analysis of the spatiotemporal variations of water yield and the driving factors are helpful to realize the effective management and protection of water resources, especially for alpine areas.

With the rapid development of remote sensing technology, the InVEST model plays the dominant role in the evaluation of water yield [11]. Since the 1970s, on the basis of the coupling of biological and geochemical processes with hydrological processes, models of distributed physical and conceptual hydrology have been applied to assess and forecast the water resource dynamics at the basin scale [12,13], such as the Soil and Water Assessment Tool (SWAT) [14], Artificial Intelligence for Ecosystem Services (ARIES) [15], Topography-Based Hydrological Model (TOPMODEL) [16], Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) [13], Variable Infiltration Capacity (VIC) [17,18], Conceptual Hydrological Model [19], etc. These models perform well with hydrological studies of various geographical and climatic characteristics. However, simpler models are potentially more suitable when available datasets are lacking [11]. Compared with other models, the InVEST model, jointly developed by the Stanford University, the World Wildlife Fund (WWF), the Nature Conservation Society (TNC), and other related institutions, is based on the Budyko curve and annual precipitation and operates with relatively low data requirements, providing an effective approach for estimating water yield at a high spatial resolution and at different scales [13]. In addition, the InVEST model can easily be used to explore the long-term impacts of climate and LULC changes on water yield with its simple interface and multiple-scenario setting function [9].

Water yield can be influenced by the combined impacts of climate and LULC changes [20]. Climate change can alter the water yield by changing the precipitation and air temperature in a basin [11,12,21]. LULC change can have a significant effect on hydrological fluxes due to variations in the physical characteristics of the land surface, soil, and vegetation, such as roughness, albedo, infiltration capacity, root depth, architectural resistance, leaf area index (LAI), and stomatal conductance [17,22]. Previous studies have successfully evaluated the water yield of different basin scales in different regions of the world, including Iran, India, Europe, and China, using the InVEST model [20,21,22,23,24,25,26] and explored the effects of climate and LULC changes on spatial and temporal variations in water yield [6,13,27,28]. However, the InVEST model has only been used in a handful of studies in alpine areas of Northwest China. Moreover, due to the limitations of terrain, transport, and economy, hydrological stations are relatively scarce in alpine areas, resulting in high uncertainty in hydrological simulations and analyses. Thus, the performance of the InVEST model on water yield in the alpine areas of Northwest China still need to be explored. The Fifth Assessment Report (AR5) was published by the Intergovernmental Panel on Climate Change (IPCC) in 2012, and it was confirmed that global warming will persist [29]. The QTP is one of most sensitive regions to global climate change; the air temperature has risen significantly by 0.3 °C per decade over 50 years and the precipitation in some regions has fluctuated greatly [9,30,31,32,33]. Previous studies have elucidated that climate change has profoundly affected the hydrogeological process of the QTP [34,35,36]. The Shule River Basin, located in the northeastern margin of the QTP and the western part of the Qilian Mountain, is the “Lifeline” and “Natural Water Tower” of 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.

2. Materials and Methods

2.1. 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². 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.

2.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.

2.3. 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):

$$Y_{xj}=\left[1-\frac{AE{T}_{xj}}{P_x}\right]\times P_x$$

(1)

where $Y_{xj}$ is the water yield for land use type j on pixel $x$, $AE{T}_{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].

$$\frac{AE{T}_{xj}}{P_x}=1+\frac{PE{T}_{xj}}{P_x}-{\left[1+{\left(\frac{PE{T}_{xj}}{P_x}\right)}^{\omega }\right]}^{1/\omega }$$

(2)

where $PE{T}_{xj}$ is the potential evapotranspiration for land use type $j$ on pixel $x$, $\omega$ is related to the plant available water content, precipitation, and the Zhang parameter (Z,) which can be calculated by Equation (3) as given below:

$${\omega }_x=Z\times \frac{AW{C}_x}{P_x}$$

(3)

where $Z$ is an empirical constant, which captures the local precipitation pattern and additional hydrogeological characteristics, and $AW{C}_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. Difference between AET and PET.

Item	Difference [48]
AET	actual amount of water that evaporates from the Earth back to the atmosphere
PET	combined evaporation from the soil surface and transpiration from plants, which is the maximum value of evaporation that would occur under ideal conditions

).

2.4. 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. Dataset sources and statistics for model inputs for water yield.

Data	Data Description	Data Source
Climate data	Daily precipitation data	China Meteorological Science Data Center (http://data.cma.cn (accessed on 7 October 2019)))
	Daily mean temperature data
	Daily maximum temperature data
	Daily minimum temperature data
Soil data	Soil texture and root depth	National Cryosphere Desert Data Center (http://data.casnw.net/portal/ (accessed on 5 December 2020))
Digital elevation model	Shuttle Radar Topography Mission (SRTM) with a resolution of 90 m	Resource and Environmental Science and Data Center (http://www.resdc.cn/ (accessed on 18 November 2020))
LULC	500 m spatial resolution	NASA Earth Science Data Systems (https://search.earthdata.nasa.gov/search (accessed on 16 October 2020))
Streamflow data	Annual water yield in 2001–2019	Measured data of Changmapu Hydrological Station

.

2.4.1. LULC

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.

2.4.2. 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₀ was calculated in this study using the Hargreaves equation. The Hargreaves equation is shown in Equation (4):

$$E{T}_0=C{R}_a{\left(T_{max}-T_{min}\right)}^E\left[\frac{T_{max}-T_{min}}{2}+T\right]$$

(4)

where $E{T}_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].

2.4.3. 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):

$$\begin{gathered} PAWC=54.509-0.132\times Sand-0.003\times {\left(Sand\right)}^2-0.055\times Silt-0.006\times {\left(Silt\right)}^2-0.738\times Clay \\ +0.007\times {\left(Clay\right)}^2-2.668\times OM+0.501\times {\left(OM\right)}^2 \end{gathered}$$

(5)

where $Sand$, $Silt$, $Clay$, and $OM$ are the proportion of clay, sand, silt, and organic matter in the soil.

2.4.4. 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].

2.4.5. 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.

2.5. 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. Simulations of water yield variation by scenario.

	Actual Conditions	Conditions without Climate Change	Conditions without Land Use Change
2001	2001 precipitation	2001 precipitation	2001 precipitation
	2001 ET0	2001 ET0	2001 ET0
	2001 land use	2001 land use	2001 land use
2019	2019 precipitation	2001 precipitation	2019 precipitation
	2019 ET0	2001 ET0	2019 ET0
	2019 land use	2019 land use	2001 land use

). Under the scenario without climate change, the precipitation and ET₀ 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.

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]:

$$G_c=\frac{C}{C+L}\times 100\%$$

(6)

$$G_L=\frac{L}{C+L}\times 100\%$$

(7)

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.

2.6. 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):

$$slope=\frac{n\times {\sum }_{i=1}^{n}i\times wyil{d}_i-\left({\sum }_{i=1}^{n}i\right)\left({\sum }_{i=1}^{n}wyil{d}_i\right)}{n\times {\sum }_{i=1}^n{i}^2-{\left({\sum }_{i=1}^{n}i\right)}^2}$$

(8)

where $n$ is the number of study years (19 in this study), $i$ is the serial number of the year, $wyil{d}_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.

2.7. 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². The R² (ranging from 0 to 1) represents the proportion of the variance in measured data explained by the model [57]. A higher value of R² 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² (all R² in the text represent the adjusted R²), 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.).

3. Results

3.1. Change Characteristics of Climatic Elements and Land Use

3.1.1. 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⁻¹ (y = 0.02x − 41.10, p > 0.05) and 5.92 mm·year⁻¹ (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⁻¹ (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-to-northwest 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). Furthermore, about 93.46% of the grids had a temperature rise in 2010 than 2001. The mean annual air temperature presented an increasing trend (p > 0.05) in the low-altitude areas from 2001 to 2019, which account for 76.08% of the USRB. The annual precipitation increased over the study area in 2001–2010 and 2010–2019. In particular, the annual precipitation significantly increased (p < 0.05) from 2001 to 2019 and extremely significantly increased (p < 0.01) in the southeast. In addition, the area of actual evapotranspiration in 2010 is greater by 49.86% than that in 2001, and 2019 is greater than 2010 by 61.23%. The annual actual evapotranspiration also significantly increased in 2019 than 2001 in the east of the northwest and low-altitude regions of the southeast, accounting for 33.39%, while only 4.01% of areas emerged as significantly decreasing over the study area (p < 0.05).

3.1.2. 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 crops–natural 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. Transition matrix of LULC types between 2001 and 2019 in the USRB (unit: km²).

LULC Type		2001
		Grassland	Permanent Wetland	Cropland	Built–Up Land	Crops–Natural Vegetation Transition	Permanent Snow and Ice	Barren
2019	Grassland	5141.17	134.46	108.18	96.00	80.09	77.03	432.7
	Permanent wetland	39.83	122.22	17.48	18.47	13.36	15.06	114.8
	Cropland	25.97	10.16	107.11	14.06	10.22	16.68	109.7
	Built–up land	17.45	9.19	10.35	105.13	14.36	11.29	119.4
	Crops–natural vegetation transition	12.87	4.89	5.58	8.62	98.99	17.82	142.1
	Permanent snow and ice	10.76	6.24	8.57	7.83	7.51	257.13	249.7
	Barren	15.28	10.05	14.95	19.56	23.54	45.33	2997

, 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² 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²), 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²) and permanent snow and ice (10.05 km²) to barren land, while the surplus transition included permanent snow and ice to grassland (134.46 km²) and grassland to permanent snow and ice (39.83 km²). 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%).

3.2. 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² = 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.

3.3. Temporal and Spatial Variation Characteristics of Water Yield

3.3.1. 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⁻¹ (y = 2.361x − 4621.892, p < 0.05).

3.3.2. 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.

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.

3.4. 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⁸ m³, followed by grassland, permanent wetland, and permanent snow and ice, which were up to 2.88 × 10⁸, 1.03 × 10⁸, and 1.41 × 10⁸ m³, respectively. In contrast, the total water yields of cropland, built-up land, and crops–natural vegetation transition were relatively small, up to just 0.24 × 10⁸, 0.12 × 10⁸, and 0.38 × 10⁸ m³, respectively.

3.5. 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. Total water yield under different scenarios (unit: 10⁸ m³).

LULC Type	Scenario
	Actual Conditions		Conditions without Climate Change	Conditions without Land Use Change
	2001	2019	2019	2019
Grassland	0.85	4.21	1.01	4.57
Permanent wetland	0.73	1.30	0.80	0.84
Cropland	0.06	0.44	0.07	0.55
Built–up Land	0.10	0.19	0.10	0.54
Crops–natural vegetation transition	0.34	0.46	0.37	0.58
Permanent snow and ice	1.25	1.18	1.56	1.37
Barren	9.12	11.96	7.33	15.51
Total water yield	12.45	19.74	11.24	23.96

). 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.

4. Discussion

4.1. 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⁻¹, 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⁻¹, 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].

4.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.

4.3. 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.

5. 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⁸ m³ in 2001 to 79.00 × 10⁸ m³ 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.