Spatial and Temporal Evolution Characteristics of Water Conservation in the Three-Rivers Headwater Region and the Driving Factors over the Past 30 Years

: The Three-Rivers Headwater Region (TRHR), located in the hinterland of the Tibetan Plat-eau, serves as the “Water Tower of China”, providing vital water conservation (WC) services. Understanding the variations in WC is crucial for locally tailored eﬀorts to adapt to climate change. This study improves the Integrated Valuation of Ecosystem Services and Trade-oﬀs (InVEST) water yield model by integrating long-term time series of vegetation data, emphasizing the role of inter-annual vegetation variation. This study also analyzes the inﬂuences of various factors on WC variations. The results show a signiﬁcant increase in WC from 1991 to 2020 (1.4 mm/yr, p < 0.05), with 78.17% of the TRHR showing improvement. Precipitation is the primary factor driving the interan-nual variations in WC. Moreover, distinct interactions play dominant roles in WC across diﬀerent eco-geographical regions. In the north-central and western areas, the interaction between annual precipitation and potential evapotranspiration has the highest inﬂuence. Conversely, the interaction between annual precipitation and vegetation has the greatest impact in the eastern and central-southern areas. This study provides valuable insights into the complex interactions between the land and atmosphere of the TRHR, which are crucial for enhancing the stability of the ecosystem.


Introduction
The Three-Rivers Headwater Region (TRHR), located in the hinterland of the Tibetan Plateau, plays a key role in ensuring ecosystem security [1].Renowned as the "Water Tower of China", the TRHR is the birthplace of the Yellow River, Yang e River, and Lancang River and is essential for water storage and conservation [2].Water conservation (WC) is a crucial ecosystem regulation service, encompassing the interception and maintenance of precipitation by various components of the ecosystem, including vegetation canopy, li er, soil, lakes, water reservoirs, and more [3].The WC services of the TRHR satisfy the water needs of the local ecosystem and provide water resources for external and downstream areas [4,5].However, the TRHR experiences an uneven and relatively low precipitation distribution.Rising temperatures also contribute to the degradation of frozen soil and increased evapotranspiration, which have a negative impact on WC [6,7].Consequently, accurately assessing the WC in the TRHR and determining its driving factors is crucial to understanding of the ecosystem's ability to maintain the stability and sustainability of water resources.
Previous studies have revealed an increasing trend in WC within the TRHR since the 1990s [8][9][10].The increase in WC results from the complex interplay between ecosystem and hydrological processes, influenced by various factors such as ecosystem type, surface characteristics, climate conditions, and human activities [11].However, differences exist in understanding the a ribution of WC variations in this region.Studies suggest the rise in WC can be a ributed to factors such as precipitation [3], potential evapotranspiration [12], the combined effect of increased precipitation and vegetation growth [13,14], or the combined effect of increased precipitation and decreased potential evapotranspiration [10].Climate, vegetation, and other factors, such as slope and altitude, also play a significant role in shaping the spatial distribution of WC [15].However, comprehensive investigations into the spatial responses of WC services in the TRHR under the synergistic influence of environmental elements are lacking.Assessing the changes in WC services and understanding the driving factors behind temporal and spatial variations are crucial.Moreover, further exploration is necessary for quantitatively examining of the interplay among distinct influencing factors [16].
Currently, the main methods used to evaluate WC include the water storage capacity [5] and the water balance methods [17].The water storage capacity method is primarily used to assess the scale of water storage capacity in sample plots, considering the contributions of vegetation, dead branches, and soil in intercepting precipitation [18].By contrast, the assessment of large-scale WC is typically accomplished through remote sensing inversion and ecological model simulation based on water balance [19].Remote sensing products provide information on WC at larger time and spatial scales but may not accurately capture local parameters such as soil water storage capacity [20].The Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model, based on the water balance, is widely used for long-term WC simulation due to its relatively straightforward data preparation and its ability to capture the physical mechanisms of soil-vegetationatmosphere interactions in terrestrial ecosystems [21].Previous research has highlighted that neglecting the influence of interannual changes in vegetation may result in inaccurate simulations of hydrological models [22,23].However, most studies that employ the In-VEST model to assess WC have overlooked the role of interannual variation in vegetation, potentially leading to inaccuracies in estimating actual evapotranspiration [24].Therefore, adequately considering vegetation dynamics is crucial to accurately representing the WC response.
This study aims to simulate and analyze the temporal and spatial characteristics of WC in the TRHR over the past 30 years.Correlation analysis and geographical detectors are used to investigate the combined impact of climate, vegetation, and other driving factors on WC variations.A key aspect of this study is improving the InVEST water yield model by integrating long-term vegetation dynamics and ground-based meteorological observations, thus highlighting the impact of interannual vegetation variation.Specifically, considering its importance as a parameter controlling land-atmosphere interactions and water cycling, a regional correction of actual evapotranspiration on alpine grasslands was conducted [25].The differences in the interaction between influencing factors on WC in different eco-geographical regions were also emphasized.Overall, this research deepens the understanding of land-atmosphere interactions in the TRHR, provides valuable insights for maintaining ecosystem stability in alpine areas and offers scientific support for ecosystem management and decision making processes.

Research Area
The TRHR (31°39′-36°12′ N, 89°45′-102°23′ E) is located in the northwestern part of Qinghai Province and is renowned as the birthplace of the Yellow River, Yang e River and Lancang River (Figure 1a).This region has an average elevation of over 4400 m (Figure 1a) and is characterized by predominantly low temperatures.A west-to-east elevation trend is exhibited in this region, with altitudes ranging from 2300 to 6600 m.Slopes in the area range from 0 to 26 degrees, with steeper slopes (>10 degrees) concentrated in the southeastern and central-southern parts of the region (Figure 1e).The annual average temperature fluctuates between -8 °C and 17 °C, leading to widespread permafrost in the TRHR.Precipitation in the area is typically concentrated between June and September, with a gradient of decreasing rainfall from southeast to northwest, ranging from over 840 mm to less than 160 mm (Figure 1b).The region also experiences high potential evapotranspiration, with a multi-year average of 600 mm and a range between 400 mm and 1000 mm (Figure 1c).The ecosystem in the TRHR is relatively simple in structure and fragile in terms of its environmental conditions, making it highly vulnerable to the direct impacts of climate change.These unique climatic characteristics in the region have resulted in the formation of distinct eco-geographical regions (Figure 1f, where I and II represent subcold and temperate, and B, C, and D represent semi-humid, semi-arid, and arid zones, respectively) [26].The dominant ecosystem type in the TRHR is alpine grassland, which covers 71.18% of the total area (Figure 1a).Thus, the region has extensive coverage of alpine shrub meadows, alpine grasslands, alpine meadow steppes, and alpine meadows (Figure 1f), collectively forming the alpine grassland.Additionally, the land use distribution in the TRHR includes farmland, woodland, water bodies, built-up land, and unused land, accounting for 0.66%, 4.26%, 5.63%, 0.1%, and 18.17% of the total land area in 2020, respectively.Unused land is located primarily in the western part of the study area, while forest land is in the central-southern and eastern regions.The leaf area index (LAI) in the TRHR exhibits a spatial gradient decreasing from southeast to northwest.The annual average LAI value ranges from 0 to 2 m 2 m −2 .The eastern region exhibits an average LAI value exceeding 1 m 2 m −2 , while the western region has an average LAI value below 0.2 m 2 m −2 (Figure 1d).

Data Source and Preprocessing
Monthly observational data were gathered from 20 weather stations located in the TRHR.These data cover the period from 1991 to 2020 and include parameters such as temperature, precipitation, wind speed, and relative humidity.A thin plate spline function was employed to spatially interpolate these site-based observational data, utilizing elevation as a covariate, to obtain regional-scale data across the entire region.Potential evapotranspiration data were derived from the 1 km monthly potential evapotranspiration dataset in China provided by the National Tibetan Plateau Data Center [27].Previous studies have demonstrated the effectiveness of this dataset in simulating the WC in the Tibetan Plateau [16].
The vegetation dataset used in the current study included the LAI and the homogenized vegetation index (NDVI).The LAI datasets were utilized to calculate the annual vegetation evapotranspiration coefficient (Kc) of the region.In contrast, the NDVI datasets were employed to assess the explanatory power of vegetation factors regarding the spatial distribution of WC.The Global Mapping (GLOBMAP) LAI version 3 dataset was acquired to obtain the required LAI data [28].This dataset provided the necessary LAI data with a temporal resolution of 16 days before 2001 and 8 days after.The 1 km spatial resolution NDVI data were downloaded from Zenodo (h ps://zenodo.org/record/6295928,accessed on 1 September 2022).Monthly vegetation index data from 1991 to 2020 were generated using the maximum synthesis method.
Spatial data for the velocity coefficient were obtained based on empirical values associated with different land cover types to calculate water yield using the InVEST model [29].Root depths for different vegetation types were determined by referring to relevant data from the Food and Agriculture Organization (FAO).Woodland areas were assigned a root depth of 3 m, grassland areas were assigned a root depth of 1 m, and non-vegetated areas were not assigned a root depth.
The annual total water yields  in the Qinghai Water Resources Bulletin were obtained to assess the accuracy of the model simulation.All the data used in this study are listed in Table 1.The WC for each grid was calculated using the water yield obtained from the InVEST model.This calculation considers various factors, including soil permeability, topography, and flow velocity coefficients, which were specific to different land-use types [29]: where WC (mm) is the WC for each land use type (j) in each grid (i).The influence of different land use types on surface runoff was considered by using the flow velocity coefficient (Velocity ).The topographic index (TI ), expressed as a dimensionless parameter, was used to assess the topographic characteristics of each grid (i).The Cosby model [30] was utilized to calculate the soil-saturated hydraulic conductivity (K ) (cm/d) in centimeters per day.Additionally, the water yield (Y ) was determined using the InVEST model and measured in millimeters (mm).
where D is the number of grids in the catchment area (dimensionless), Soil is the soil depth (mm), and P is the percentage slope (%) calculated using the suite of Hydrology Tools (Arc Hydro) in ArcGIS 10.5.
where sand and clay represent the content of sand and clay in the soil (%), respectively.Water yield in InVEST is defined as the amount of water that runs off the landscape.It calculates water yield by applying the water balance principle at the sub-watershed level.The conceptual diagram of the water balance model is illustrated in Figure 2. The water yield model is based on the Budyko curve and annual average precipitation.The annual water yield for each pixel was calculated as follows: where AET(x) is the annual actual evapotranspiration of grid x and P(x) is the annual precipitation of grid x.
In the water balance formula, the vegetation evapotranspiration is calculated based on the Budyko hypothesis proposed by Zhang et al. [32]: where ET (x) is the potential evapotranspiration of grid x.K (l ) represents the evapotranspiration coefficient of vegetation for land use type j and grid x.An estimation of Kc values for farmland and woodland was obtained using the relationship Kc = LAI/3 (Kc equal to 1 for LAI > 3).The Kc values for water bodies, built-up land, and unused land were set to 1.1, 0.3, and 0.5, respectively [31].Kc for grassland was estimated using the method in Section 2.3.2.w(x) is an empirical parameter proposed by Donohue et al. [33]: where Z is the Zhang coefficient, which corresponds to the number of precipitation occurrences every year.When Z = 3.5, the error between the simulated water yield and actual runoff is the smallest (7.9%) [34].AWC(x) represents the available water content of soil (mm), which is used to determine the total water stored and provided by soil for plant growth.AWC(x) is calculated as follows: AWC(x) = Min(root restricting layer depth, root depth) × PAWC (7) where AWC(x) is estimated as the product of the plant available water capacity (PAWC) and the minimum of root-restricting layer depth and vegetation rooting depth.The calculation of AWC is performed using the following formula for PAWC: PAWC = 54.509-0.132sand%− 0.003(sand%) − 0.055silt% − 0.006(silt%) − 0.738clay% + 0.007 (clay%) − 2.688OM% + 0.501(OM%) Actual evapotranspiration (AET) for other land use types (water bodies, built-up land, unused land) is calculated directly with ET(x): where ET represents potential evapotranspiration and K (l ) represents the evapotranspiration coefficient for the land use type without vegetation cover.The value of K (l ) is taken in accordance with a fixed value (x) [35]: the value of the water body is 1.1, the value of built-up land is 0.3, and the value of unused land is 0.5.

Calibration Method of Evapotranspiration Coefficient in Grassland Based on Interannual Variation in Vegetation
Alpine grassland covers more than 70% of the total area of the TRHR, making it the dominant terrestrial ecosystem type.Therefore, the main objective of this study is to calibrate the Kc coefficient specifically for alpine grasslands.Capturing the dynamic evapotranspiration pa erns of alpine grasslands within the InVEST model accurately is crucial.The growth process of alpine grassland can be categorized into different stages, which exhibit interannual variability influenced by climatic changes.
Monthly grid data of LAI were employed to estimate the Kc coefficient for each grid.Subsequently, the monthly average Kc coefficient was calculated for the entire study area by averaging the grid-based Kc coefficients.Finally, the 12 monthly regional average Kc coefficients for each year were combined to generate the annual value of Kc coefficients for grassland, following the approach recommended by the InVEST model.The specific procedure is as follows.
First, the monthly average Kc coefficients of grassland were computed based on the different growth stages of vegetation.Previous research and actual meteorological conditions identify three distinct stages for the growth season of alpine grasslands: the initial growth season (May), the middle growth season (June, July, and August), and the late growth season (September and October) [36].Following the FAO-56 guidelines, the recommended values for Kc during the initial, middle, late, and non-growing stages of grassland are 0.4, 1.05, 0.85, and 0.4, respectively [35].However, these coefficients can be adjusted to reflect the actual climatic conditions of the study area.The Kc coefficients for the middle and late growth stages were modified.Notably, the Kc coefficients for the early growth and non-growth stages remained unadjusted.The adjustment equation employed in this study is as follows: where RH is the mean value for the daily minimum relative humidity (%), u is the wind speed at 2 m height (m/s), and K is the recommended Kc value in FAO-56.h is the mean plant height (m), fi ed by wind speed, relative humidity, and LAI [37]: Second, the annual vegetation evapotranspiration coefficients (K _ ) was determined using the calculation formula recommended by the InVEST model.Each month's regional average Kc coefficients was calculated by averaging the grid-based monthly Kc coefficients across the entire study area.Then, each month's regional average potential evapotranspiration was computed.Finally, the annual Kc coefficient was derived by using the following formula: where K represents an average Kc coefficient of month m (1-12) and ET is the corresponding ETo.
The differences between the observed yearly total water yield (TWY) and the simulated TWY were evaluated to investigate the accuracy of the simulations by the improved InVEST model using the regionally calibrated Kc coefficient, including the coefficient of determination (R 2 ) and the root mean square error (RMSE).The differences between the observed and simulated TWY obtained through three different methods of calculating Kc are further compared to illustrate the effectiveness of the improved model (Figure 3).The approach proposed in Section 2.3.2 was utilized in the first method to estimate the Kc coefficients (Figure 3a).In the second method, the Kc coefficient was determined using the relationship Kc = LAI/3 (where Kc equals 1 when LAI > 3) (Figure 3b).This method is mentioned in the InVEST model manual.The calculation process involves synthesizing the monthly values of LAI using the maximum value synthesis method, followed by averaging the values for each month across each year to derive the annual average LAI.The Kc value for each grid is computed based on this relationship Kc = LAI/3, and the overall Kc value required by the model operation is obtained by averaging across the entire region.The third method uses a fixed Kc coefficient of 0.65, corresponding to the recommended Kc value for grassland in the InVEST model (Figure 3c).The observed TWY data are obtained from the Qinghai Water Resources Bulletin, covering 1994 to 2020.
Comparatively, the simulation results of TWY based on the improved model display better performance, with an increase in the R 2 ranging from 0.51% to 14.6% and a reduction in RMSE ranging from 0.478 × 10 10 m 3 to 2.193 × 10 10 m 3 .The improved model enhances the performance of the InVEST water yield model by incorporating information on interannual vegetation changes.

Geographical Detector
The geographical detector method is a powerful tool for analyzing the driving forces by detecting the heterogeneity of the spatial stratification of elements [38].Factor detectors and interaction detectors were utilized to examine the driving factors of WC.This study focused on two important aspects: identifying the factors that significantly influence WC and investigating the interactions between these factors.
The factor detector is employed to detect the spatial differentiation of dependent variables and assess the explanatory power of independent variables for WC.The q-value quantifies the strength of this explanatory power.The specific calculation formula is as follows: where h is the stratification of the dependent or independent variable, with values ranging from 1 to L. N and N are the number of cells within stratum h and the whole region, respectively.σ and σ are the variances of the values of the dependent variable in the stratum and the whole region, respectively.SSW is the sum of the variances within the stratum, and SST is the total variance of the whole region.q indicates the explanatory power of independent variables for dependent variables, with a value range between 0 and 1.A higher value of q indicates a stronger explanatory power of independent variables for dependent variables.The main influencing factors of WC can be determined by examining the value of q.The q value can be tested for significance using the geographic detector method because it follows a non-central distribution after a simple transformation.This study assesses whether the two factors have an interaction.It determines the strength and characteristics of the interaction by comparing the values obtained from individual factors and the interaction values obtained from the geographic detector.For example, if q(X1 ∩ X2) < Min(q(X1), q(X2)), it indicates that the interaction is characterized by nonlinear weakening.If Min(q(X1), q(X2)) < q(X1 ∩ X2) < Max(q(X1), q(X2)), the interaction shows the characteristics of one-factor nonlinearity weakening.If q(X1 ∩ X2) > Max(q(X1), q(X2)), the interaction displays the characteristics of two-factor enhancement.If q(X1 ∩ X2) = q(X1) + q(X2), it signifies that the interaction is characterized by two-factor independence.Finally, if q(X1 ∩ X2) > q(X1) + q(X2) , it indicates the interaction is characterized by nonlinear enhancement.
The objective of this study is to examine the factors that influence the spatial heterogeneity of WC.These factors include climatic factors such as annual average precipitation and annual average potential evapotranspiration, topographic factors such as slope and digital elevation model, and vegetation factors such as annual average growing season LAI.Continuous variables need to be converted into classified data to meet the requirements of the geographic detectors.Therefore, the natural breakpoint method in ArcGIS 10.5 was utilized to divide altitude, slope, annual average precipitation, annual average potential evapotranspiration, and LAI into nine categories.The analysis process of the geographic detector using NDVI as the vegetation factor was repeated to ensure the reliability of the results.

Trend Changes and Correlation Analysis
Simple linear regression analysis was initially employed to examine the interannual variation trends within the TRHR at regional and grid scales.Subsequently, Pearson correlation analysis was used to investigate the impact of individual environmental factors on the interannual variability in WC.

Total Research Approach
The specific research schemes are outlined as follows: (1) first, the monthly Kc coefficients are calculated using the monthly LAI grid data to correct for the grassland's middle and late growth stages.The method recommended in the InVEST model manual is used to synthesize the monthly Kc coefficients into an annual value.Other parameter values necessary for model execution can be found in Sections 2.3.1 and 2.3.2.(2) Next, the InVEST model is utilized to simulate the water yield of the TRHR from 1991 to 2020.The TWYs obtained from the model output are compared with the recorded TWYs in the water resources bulletin to verify the accuracy of the simulation results.Once the verification is complete, we estimate the WC based on the grid-scale water yield results. (3) The influence of vegetation and climate on the temporal changes in WC is analyzed using correlation analysis.Additionally, the geographic detector method is employed to explore the effects of climate, vegetation, and topography on the spatial differentiation of WC and the interactions among these factors (Figure 4).Implementation flow chart of this study (Formula a is derived from Zhao, et al. [37].Formula b is derived from Allen, et al. [35].Formula c is derived from Sharp, et al. [31]).

Characteristics of Interannual Variation in Water Conservation
After regional correction, the multi-year average Kc coefficient for the grassland is 0.72, which is higher than the recommended value of 0.65 suggested by the InVEST model (Figure 5).Overall, the Kc coefficient showed a slight downward trend over the years, which was not very pronounced.The coefficient of variation of Kc from 1990 to 2000 was 0.018, indicating significant fluctuation of the Kc coefficient of grassland during that period.By contrast, the overall coefficient of variation from 1990 to 2020 decreased to 0.013, indicating decreased fluctuation in the 2000s and 2010s.Among the years analyzed, the evapotranspiration coefficient of grassland reached its maximum value of 0.75 in 1997 and its minimum value of 0.71 in 1993.WC in the TRHR exhibited a significant increase (p < 0.05) from 1991 to 2020, with an annual growth rate of 1.4 mm (Figure 6).A series of ecological projects have been implemented, and the ecological environment in the TRHR has improved since 2000.Accordingly, WC has changed from a weak increase to a strong increase since 2000.Specifically, WC demonstrated a nonsignificant upward trend (1.16 mm/yr) from 1991 to 1999.However, a significant increase was observed in WC, with a growth rate of 1.

Spatial and Temporal Distributions of Water Conservation
The spatial distribution of WC in the TRHR shows significant heterogeneity (Figure 7a).The areas with high WC were found mainly in the southeast of the Yellow River and the source area of the Lancang River, specifically in Banma County, Nangqian County, and Jiuzhi County.The average value of WC over the years was 77.82 mm.Among the sub-regions, the mean value of WC was highest in the HIB1 region (125.32mm), followed by the HIIC2 (112.97 mm) and the HIC1 regions (71.97 mm).The HIB1 region is characterized by lower altitude and a more humid climate favorable conditions for WC.Moreover, this area boasts relatively high vegetation coverage, pivotal in regulating runoff and increasing water yield, making it a promising region for WC.WC exhibited an increasing trend in 96.12% of the region during the study period, and among these trends, 78.17% were statistically significant (p < 0.05) (Figure 7b).These findings highlight substantial WC improvements across most TRHR over the past three decades.Notably, the southeastern region of the TRHR exhibited an even more pronounced trend, demonstrating a rate exceeding 2 mm/yr.However, approximately 3.87% of the region still faces challenges maintaining WC efforts.Moreover, approximately 1.72% of the total area of the TRHR experienced a significant decrease in WC throughout the study period, with a rate of -2.77 mm/yr (p < 0.05).

Drivers of Interannual Change
The TRHR experienced an increase in average annual precipitation during the period from 1991 to 2020, with a rate of 3.6 mm/yr (p < 0.5) (Figure 8a).Before 2000, the increase in precipitation was not statistically significant, with an annual increase of only 3.71 mm/yr.However, the rate of precipitation increase has noticeably accelerated since 2000, with an annual increase of 3.80 mm/yr.By contrast, potential evapotranspiration demonstrated no significant change throughout the study period, with a rate of 0.48 mm/yr (Figure 8b).The rate of increase in potential evapotranspiration slowed after 2000.Additionally, the average annual LAI from 1991 to 2020 exhibited an upward trend, indicating improved vegetation growth (Figure 8c).The TRHR generally had a low baseline LAI, but a fluctuating upward trend was observed from 1991 to 2020, with a growth rate of 0.0008 m 2 m −2 /yr (p < 0.05).In particular, the annual fluctuation in LAI was more pronounced during 2000-2020 than in the 1991-2000 period.Table 2 shows a significant positive correlation between precipitation and WC.The increase in precipitation in the study area over the past 30 years has played an important role in improving WC.However, no significant correlation between potential evapotranspiration and WC was observed throughout the study.Before 2000, potential evapotranspiration and WC had a positive correlation, but this relationship was not statistically significant.However, this correlation changed to a significant negative correlation after 2000.The increase in potential evapotranspiration indicates that meteorological factors such as temperature, radiation, and wind speed have a stronger ability to demand atmospheric water, ultimately leading to an increase in actual evapotranspiration and a decrease in WC.The increase in LAI also promotes greater WC.After 2000, the influence of LAI on WC increased.

Drivers of Spatial Heterogeneity
The geographical detector method was employed to identify the primary factors influencing the spatial distribution of the annual average WC in the region, given the significant spatial variability of WC within the TRHR.Analysis of the q values of each factor indicated that annual average precipitation had the highest explanatory power for WC (0.73), followed by LAI (0.582), annual average potential evapotranspiration (0.207), digital elevation model (DEM) (0.221), and slope (0.109), as determined.These factors were statistically significant (p< 0.01) (Table 3).The annual average precipitation accounted 73% of the spatial distribution of WC, which was the primary driving factor in the region.LAI demonstrated moderate explanatory power for WC, while potential evapotranspiration and altitude had weak influences.The slope had the weakest influence on the spatial heterogeneity of WC, contributing only 10.9%.A multi-year average NDVI from 1991 to 2020 was used as a vegetation factor, yielding consistent findings with those obtained using LAI to ensure the reliability of the results.An interaction detector was used to analyze the interaction between two factors in the TRHR and various ecological regions to explore further the impact of different factors on the spatial variation in WC (Figure 9).The NDVI data were applied to replicate the entire interaction detection process and validate the reliability of the interaction detection findings.The entire interaction detection process was also replicated using NDVI data to validate the reliability of the findings.The interaction between two factors provided a better explanation for the spatial differentiation of WC compared to a single factor.In particular, the interaction between precipitation and other factors played a key role in the spatial differentiation pa ern of WC.The interaction between precipitation and potential evapotranspiration had the strongest influence on WC in the north-central and western regions (HIC1) (explaining at least 63.1% of the spatial distribution of WC), where the background value of LAI was relatively low.Precipitation and vegetation together accounted for 72.8% and 60% of the distribution of WC in regions HIIC2 and HIB1, respectively.

Interannual Variation in Water Conservation
Vegetation dynamics affect the regional WC through their influence on rainfall-runoff capacity [39] and evapotranspiration [40].Thus, considering the impact of vegetation dynamics is crucial in quantifying WC services.Some hydrological models, including the Soil and Water Assessment Tool (SWAT), Water and Energy transfer Processes in Large river basins (WEP-L), and Water and Energy Budget Distributed Hydrological Model (WEB-DHM), also consider the effects of vegetation growth on precipitation distribution, soil moisture, and groundwater recharge [41][42][43].This study focused specifically on the influence of vegetation on WC and improved the InVEST model to reflect interannual variation in vegetation and enhance the accuracy of the simulation.
The grassland growth and actual evapotranspiration in the Tibetan Plateau exhibit evident stage and interannual variations [44].An empirical coefficient called the Kc coefficient is commonly used to determine vegetation water requirements.This coefficient is typically calculated by dividing actual evapotranspiration (AET) by reference evapotranspiration (ETo).Previous studies have employed this ratio to estimate the annual Kc coefficient.However, obtaining AET data in advance is necessary for the corresponding research, which limits the generalizability of this approach [45].The Kc coefficient for a specific region can be determined by establishing an empirical relationship between vegetation indexes, such as LAI and NDVI, and the corresponding Kc coefficients.However, collecting and obtaining the observed Kc coefficients beforehand to establish this relationship is crucial [46].This study first categorized grassland Kc coefficients into three growth stages to be er understand the variations in grassland evapotranspiration in this research.Then, the Kc coefficients recommended by FAO-56 for each specific stage of grassland growth were utilized.The Kc coefficients for the middle and late growth stages were also refined by incorporating local climate data.These adjusted coefficients were combined to derive annual-scale grassland Kc coefficients.The methodology developed in this study enables rapid and accurate estimation of Kc coefficients across a wide range of temporal and spatial scales.However, when estimating the Kc coefficient, the absence of extensive plant height data requires estimating plant height using LAI and an empirical formula.However, this approach is specific to alpine grasslands because of the regional nature of the empirical formula employed and may not be apply to other regions.If simulated or observed plant height data become available for a broader range of times and locations, the method proposed in the current research could potentially be applied to grassland areas on a global scale.

Factors Affecting Temporal and Spatial Differentiation of Water Conservation
In the past 30 years, the WC in the TRHR has increased significantly, which is consistent with previous studies [8,13,47].WC is a complex process influenced by various factors.Climatic factors, including precipitation, temperature, and humidity directly impact on the hydrological cycle and the ability to conserve water [48,49].The relationship between interannual variations in WC and climatic factors was examined in the TRHR.The results indicated a significant positive correlation between precipitation and WC over the past 30 years, in line with previous research [3,13,16].Precipitation emerged as the primary factor influencing the spatial variation in WC in the TRHR, followed by vegetation (measured by NDVI or LAI) consistent with previous studies in the region [14,16].The evapotranspiration capacity of the terrestrial ecosystem, indicated by potential evapotranspiration, has a great effect on WC in the TRHR, along with vegetation.The TRHR displays significant spatial variations in WC, demonstrating high WC values found in areas with abundant precipitation and well-developed vegetation.Topography regulates WC by impacting vegetation structure, soil properties, and water and heat conditions [50].Nevertheless, slope has a weaker effect on WC than other influencing factors.
The implementation of ecological projects in the TRHR also contributes to enhancing WC.Grassland degradation has been observed in this region since the 1980s, reducing effectiveness in intercepting runoff and storing water [51].Since 2000, the government has initiated several ecological projects to address this issue, including the grassland ecological reward mechanism, the grazing withdrawal project, and the first and second stages of ecological conservation and restoration projects [52,53].These projects transformed approximately 1.98 × 10 4 km 2 of unused land into grassland in the TRHR by 2020.Furthermore, a remarkable increase in vegetation coverage has been observed in the entire region over the past 30 years, stabilizing surface temperature and enhancing soil moisture conditions [54].Implementing ecological projects has effectively improved WC in the TRHR, strengthening its ecological security barrier function [55].
A single factor alone is often insufficient to comprehend the entire ecosystem's complexity fully.The complexity of geographical processes typically involves the interplay of multiple factors.While the interaction between precipitation and vegetation enhances their explanatory power for WC, the degree of interaction may vary across eco-geographical regions.In the TRHR, HIC1 experiences relatively low precipitation and has a low background value of vegetation.Therefore, the interaction between precipitation and potential evapotranspiration effectively explains the distribution of WC in this area.By contrast, HIB1 and HIIC2 are relatively moist and have abundant vegetation.Vegetation plays a crucial role in regulating the distribution of precipitation and ensuring WC through transpiration and interception [56].Thus, the interaction between vegetation and precipitation in these areas has a high explanatory capacity for WC.

Adaptation to Future Climate Change
The Tibetan Plateau is predicted to experience increased precipitation in the future, which will positively impact WC in most areas of the TRHR, such as the Yang e River and Lancang River sources [3,57].However, certain areas, particularly some counties at the source of the Yellow River, face a decline in WC because of reduced precipitation and increased potential evapotranspiration resulting from a weakening summer monsoon.This increase in rainfall is accompanied by higher rainfall erosivity, leading to soil erosion and uncertainties in WC [58].The Tibetan Plateau is recognized as one of the most sensitive regions to global climate change, with its warming rate twice exceeding the global average [59].The TRHR, located in the hinterland of the Tibetan Plateau, has observed even more significant warming changes.The increasing frequency of extreme heat events in the future will have a profound impact on the terrestrial ecosystem balance in this region [60].Rising air temperatures lead to a decrease in soil moisture [61] and an increase in the saturated vapor pressure deficit [62], thereby increasing the vulnerability of WC services in this region.Climate change is anticipated to intensify vegetation expansion (greening) in the TRHR [63].The development of grasslands is closely linked to WC services [64].This large-scale vegetation greening will cause changes in evapotranspiration and precipitation pa erns [65].Permafrost in the TRHR has also been declining since 1980, while the thickness of the active layer has been increasing [51,66].Hence, future research should monitor the impact intensity of vegetation greening and permafrost degradation on WC.Changes in WC also affect the terrestrial ecosystem in alpine regions [67].Therefore, understanding WC services and their interactions with various environmental factors is crucial to developing effective measures to adapt to climate change.

Limitations and Future Work
WC requires a comprehensive understanding of the ecohydrological processes within a basin, including factors such as precipitation interception and evaporation.It also involves exploring the interaction between multiple hydrological processes and the ecosystem, such as runoff timing, water quality, flood protection, and temperature regulation.Runoff timing can indicate the capacity for WC.If the surface and groundwater systems have a strong water storage capacity, runoff from rainfall or snowmelt will slowly enter rivers or lakes, delaying runoff formation.In such cases, the runoff time will be relatively late [68].WC processes reduce soil erosion and prevent pollutants from entering water bodies.For example, vegetation slows down the flow of rainfall, allowing more time for water to permeate the soil, which reduces surface runoff and helps decrease the risk of soil erosion and the introduction of pollutants into water bodies [69].Additionally, a basin with a good WC capacity mitigates the likelihood of natural disasters such as floods and droughts by intercepting rainfall and storing soil water [70,71].These aspects of WC function will be considered to further enhance the understanding of WC in future research.
Research methods employed in this study have several limitations.First, the WC simulation does not differentiate between the Kc coefficients at grid scales.Instead, a regional average Kc coefficient is utilized, which may introduce certain uncertainties to the results.Second, the InVEST model does not account for the influence of runoff in the calculation of regional WC, which is another limitation [31,72].the model also fails to consider the impact of freezing and thawing on WC [14].Therefore, future research should focus on developing more precise methods for simulating vegetation evapotranspiration coefficients at the grid scale and exploring ways to optimize WC models in plateau regions.
The environment in the TRHR is highly complex, and this study has not delved deeply into the various factors contributing to changes in WC.Human activities such as overgrazing, expansion of agricultural land, road construction projects, and ecological project initiatives alter the surface, directly or indirectly affecting the WC services of the region [73][74][75].Future studies should focus on quantitatively analyzing the positive and negative impacts of human activities on WC to gain a more comprehensive understanding.Differentiating between the effects of human activities and natural environmental factors on WC is crucial.

Conclusions
In this study, we improved the InVEST model by integrating the long-term series of LAI data to capture the interannual variation in vegetation.The temporal and spatial distribution characteristics of WC in the TRHR were assessed quantitatively from 1991 to 2020.Trend and correlation analyses were used to examine the temporal variation trend and driving factors of WC.We also employed the geographical detector method to explore the influencing factors and their interactions in the spatial distribution of WC.The main findings are as follows: (1) WC in the TRHR exhibits a spatial pa ern with high values in the south and low values in the north.WC significantly increased (1.4 mm/yr, p < 0.05) from 1991 to 2020.Compared with 1991-1999, the annual average value of WC in 2000-2020 increased by 28.17%.Over the past 30 years, WC improved in 78.17% of the regions (p < 0.5).(2) The increase in precipitation over the past three decades has had a significant positive impact on WC (R = 0.97, p < 0.01).However, the growth in potential evapotranspiration has had a significant inhibitory effect on WC since 2000 (R = -0.5, p < 0.05).(3) The spatial variation in WC is influenced primarily by precipitation, followed by vegetation and potential evapotranspiration.The interaction among these factors has a stronger explanatory power for WC than individual factors alone.The interaction between precipitation and other influencing factors demonstrates the greatest explanatory power.The combined influence of precipitation and vegetation accounts for approximately 79.1% of the WC distribution across the study area.However, the dominant interaction factors for WC vary in different eco-geographical regions.In the north-central and western regions (HIC1) with low vegetation, the interaction between annual precipitation and potential evapotranspiration explains 65% of the variation in WC, making it the dominant interaction factor in the region.In the eastern and central-southern areas (HIB1 and HIIC2), the interaction between annual precipitation and vegetation exhibits the strongest influence.
The method proposed provides valuable insights into simulating WC.Incorporating processes such as permafrost degradation is crucial for the accurate assessment of WC in plateau regions in future work.

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

Figure 2 .
Figure 2. Conceptual diagram of the water balance model used in the InVEST water yield model [31].Only parameters shown in color are included, and parameters shown in grey are ignored.

Figure 3 .
Figure 3. Evaluation of simulated total water yield (TWY) versus observed TWY (TWY was simulated by the InVEST model using (a) calibration Kc values with interannual variation in vegetation, (b) Kc values obtained using the relationship Kc = LAI/3, or (c) Kc coefficients estimated at Kc = 0.65.Shadow represents a 95% confidence interval).

Figure 4 .
Figure 4. Implementation flow chart of this study (Formula a is derived from Zhao, et al.[37].Formula b is derived from Allen, et al.[35].Formula c is derived from Sharp, et al.[31]).

Figure 5 .
Figure 5. Regional corrected Kc coefficient of grassland in the Three-Rivers Headwater Region from 1990 to 2020.
74 mm/yr (p < 0.05) from 2000 to 2020.Over the past 30 years, the average annual WC in the TRHR was 77.82 mm.During the 1990s, the average annual WC in this region was 65 mm.The annual average WC has increased to 83.31 mm since 2000, representing an overall growth of 28.17%.The change trend of WC is similar to that of precipitation, with a change rate of 3.71 mm/yr before 2000 and 3.8 mm/yr from 2000 to 2020.The lowest value for WC occurred in 2015 (46.31 mm), while the highest was observed in 2019 (114.76 mm).The average annual total WC over the years was 2.79 × 10 10 m 3 .The overall trend of total WC aligns with the trend of WC over time.

Figure 6 .
Figure 6.Temporal variation in water conservation in the Three-Rivers Headwater Region from 1991 to 2020 (relative to the 1991-2020 anomaly).

Figure 7 .
Figure 7.The spatial distribution (a) and trend (b) of water conservation in the Three-Rivers Headwater Region from 1991 to 2020 (the inset panels in the bo om right of (a) display the water conservation values in different ecosystem zones using a violin diagram.The inset panels in the bo om right of (b) indicate the significance level (p < 0.05).The percentages of increasing (I) and decreasing (D) trends (percentage of significant trends in parentheses) are shown at the bo om of (b)).

Figure 8 .
Figure 8. Interannual variability of climatic and vegetation elements from 1991 to 2000 ((a) Average annual precipitation.(b) Average annual potential evapotranspiration.(c) Annual average leaf area index).

Figure 9 .
Figure 9. Explanatory power of interactive influencing factors on water conservation in the Three-Rivers Headwater Region and its sub-regions ((a-d): LAI represents vegetation factor, (e-g): NDVI represents vegetation factor).

Author Contributions:
Conceptualization, data curation, formal analysis, methodology, writingoriginal draft, Y.P.; methodology, conceptualization, funding acquisition, supervision, writing-review and editing, Y.Y.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the Second Tibetan Plateau Scientific Expedition Program (2019QZKK0403).Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.

Table 2 .
Correlation between climate, vegetation, and water conservation.

Table 3 .
The q values of influencing water conservation factors in the Three-Rivers Headwater Region.