Improvements to a Crucial Budyko-Fu Parameter and Evapotranspiration Estimates via Vegetation Optical Depth over the Yellow River Basin

Abstract

Against the backdrop of global warming and vegetation restoration, research on the evapotranspiration mechanism of the Yellow River basin has become a hot topic. The Budyko-Fu model is widely used to estimate basin-scale evapotranspiration, and its crucial parameter $\mathsf{\omega }$ is used to characterize the underlying surface and climate characteristics of different basins. However, most studies only use factors such as the normalized difference vegetation index (NDVI), which represents the greenness of vegetation, to quantify the relationship between $\mathsf{\omega }$ and the underlying surface, thereby neglecting richer vegetation information. In this study, we used long time-series multi-source remote sensing data from 1988 to 2015 and stepwise regression to establish dynamic estimation models of parameter $\mathsf{\omega }$ for three subwatersheds of the upper Yellow River and quantify the contribution of underlying surface factors and climate factors to this parameter. In particular, vegetation optical depth (VOD) was introduced to represent plant biomass to improve the applicability of the model. The results showed that the dynamic estimation models of parameter $\mathsf{\omega }$ established for the three subwatersheds were reasonable (R² = 0.60, 0.80, and 0.40), and parameter $\mathsf{\omega }$ was significantly correlated with the VOD and standardized precipitation evapotranspiration index (SPEI) in all watersheds. The dominant factors affecting the parameter in the different subwatersheds also differed, with underlying surface factors mainly affecting the parameter in the watershed before Longyang Gorge (BLG) (contributing 64% to 76%) and the watershed from Lanzhou to Hekou Town (LHT) (contributing 63% to 83%) and climate factors mainly affecting the parameter in the watershed from Longyang Gorge to Lanzhou (LGL) (contributing 75% to 93%). The results of this study reveal the changing mechanism of evapotranspiration in the Yellow River watershed and provide an important scientific basis for regional water balance assessment, global change response, and sustainable development.

1. Introduction

The Yellow River is the second longest river in China, and it is of great significance for the balance of the water cycle and supply of water resources in China [1,2]. Under the environmental background of global warming, the climate and vegetation in the Yellow River basin have changed significantly in recent years [3]. Global warming has had a significant effect on the Yellow River basin, and the temperature in the basin has risen significantly in recent decades [4], which has led to frequent meteorological and hydrological droughts [5,6]. Moreover, the overall trend of vegetation cover in the basin has increased in recent years [7], which indicates that the effect of various restoration projects, such as returning farmland to forests in the Yellow River basin, is obvious and that the ecosystem has been restored to a certain extent [8]. However, most of the Yellow River basin is located in arid and semi-arid areas, and water resources are extremely limited. With increases in the temperature and vegetation coverage in the basin, evaporation also continues to rise, thus causing potential conflicts between the ecosystem and human demand for water [9]. Therefore, under the dual impacts of climate change and vegetation restoration, in-depth research and interpretations of the mechanism underlying evaporation changes in the Yellow River basin have become hot topics of academic attention in recent years [10].

Evapotranspiration is a key link in the water cycle that connects water balance and energy balance [11]. Budyko’s theory of water heat coupling equilibrium is one of the most widely used methods in research estimating evapotranspiration in river basins [12]. The core idea is that the average annual evapotranspiration of a watershed depends on the precipitation and evaporation capacity of the watershed [13]. This theoretical framework has the advantages of a simple structure and few parameters [14]. However, the Budyko curve may deviate from the measured data of many watersheds when simulating climate-hydrology relationships [15]. This suggests that in addition to precipitation and potential evapotranspiration, other factors also affect basin evapotranspiration. For this reason, Chinese scholar Fu Baopu introduced control parameter $\mathsf{\omega }$ into the Budyko model [16], and it is used to characterize the role played by other factors in the basin water balance. The revised model has now been widely adopted.

In recent years, to improve the accuracy of the Budyko-Fu model for estimating evapotranspiration, researchers have developed many estimation formulas for parameter $\mathsf{\omega }$ [17]. Earlier studies usually considered the parameter $\mathsf{\omega }$ to be associated with underlying surface factors such as vegetation. For example, Li et al. [18] revealed linear correlations between parameter $\mathsf{\omega }$ and long-term average annual vegetation cover by analyzing data from 26 major river basins around the world. Shao et al. [19] used a multiple adaptive regression spline model (MARS) to assess the effect of forest cover on parameter $\mathsf{\omega }$. As such research progressed, other explanatory variables, such as climate, soil, and topography, were also incorporated into the estimation of parameter $\mathsf{\omega }$. Cheng et al. [20] pointed out that the drying index (DI), climate seasonality, and asynchrony index (SAI) have a significant impact on the spatiotemporal variation of parameter $\mathsf{\omega }$. Wang et al. [7] then analyzed the parameter in the framework of Budyko correlations with 24 environmental variables, such as meteorology, soil, and land use, and further explored the mechanism of environmental changes on evapotranspiration in the watershed. However, universal models for the estimation of parameter $\mathsf{\omega }$ are not yet available, and researchers generally focus on exploring the relationship between watershed underlying surface and parameter $\mathsf{\omega }$ through vegetation coverage and normalized difference vegetation index [21], without considering the impact of overall aboveground biomass of vegetation on parameter $\mathsf{\omega }$. Compared with vegetation cover, aboveground biomass is not only closely related to transpiration of vegetation, but also directly linked to the underground root function of vegetation. Therefore, the consideration of aboveground biomass is particularly important in the estimation of parameter $\mathsf{\omega }$.

To further improve the accuracy of the Budyko framework for estimating evapotranspiration, this study takes the upper Yellow River basin as the study area, uses long time-series multi-source remote sensing data to calculate the crucial parameter $\mathsf{\omega }$ of each subwatershed of the upper Yellow River for the period of 1988–2015, and selects different vegetation and climatic indexes to establish a dynamic estimation model of parameter $\mathsf{\omega }$. Specifically, we focused on VOD (Vegetation Optical Depth) as an indicator to characterize the aboveground biomass of plants to explore the influence of deeper vegetation information on the changes of parameter $\mathsf{\omega }$. This study will contribute to a deeper understanding of the spatiotemporal variation characteristics and causes of evapotranspiration in the Yellow River basin, which is of great significance for regional water balance assessment, global climate change response, and sustainable development.

2. Materials and Methods

2.1. Study Area

The study area includes three subwatersheds of the upper Yellow River basin (Figure 1): the watershed before Longyang Gorge (BLG), the watershed from Longyang Gorge to Lanzhou (LGL), and the watershed from Lanzhou to Hekou Town (LHT). The Yellow River secondary basin classification data used in the study were obtained from OpenGMS. Most of the watershed before Lanzhou is located in a semi-arid and semi-humid region, and most of the watershed below Lanzhou is located in a semi-arid and arid region. Most of the BLG and the LGL are located in the Qinghai-Tibetan Plateau, where the terrain is relatively flat and less affected by human activities and the natural environment is less damaged. Most of the watersheds from Lanzhou to Hekou Town are located on the Loess Plateau, which is more affected by human activities [22].

2.2. Evapotranspiration, Potential Evapotranspiration and Temperature Dataset

The actual evapotranspiration (E) and potential evapotranspiration (Ep) data used in this study were obtained from the Global Land Evaporation Amsterdam Model (GLEAM), with a spatial resolution of 0.25° × 0.25° (http://www.gleam.eu) (accessed on 11 April 2023). The model is a set of algorithms [23,24] that separately estimate the different components of land evaporation (often referred to as “evapotranspiration”): transpiration, bare-soil evaporation, interception loss, open-water evaporation, and sublimation. Additionally, GLEAM provides surface and root-zone soil moisture, potential evaporation, and evaporative stress conditions. The rationale of the method is to maximize the recovery of information on evaporation contained in current satellite observations of climatic and environmental variables. This study selected the annual E and EP data (unit: mm/a) of GLEAM v3.7a version from 1988 to 2015 for the proposed model.

The temperature (TMP) data used in this study were obtained from version 4.07 of the CRU TS (Climate Research Unit Time Series) dataset [25] (accessed on 14 April 2023), which is a widely used global climate data resource provided by the Climate Research Unit (CRU) at the University of East Anglia, UK. The dataset provides global monthly data at 0.5° resolution from 1901 to 2022, covering all land regions of the world except Antarctica. In this study, data from 1988–2015 in the CRU TS v4.07 dataset were selected to correspond to the SPEIBase v2.9 dataset and bicubic interpolated to a spatial resolution of 0.25° × 0.25°.

2.3. Precipitation Dataset

The precipitation data used in this study were obtained from the China Regional Surface Meteorological Factor Driven Dataset (http://data.tpdc.ac.cn) of the National Tibetan Plateau Science Data Center [26,27] (accessed on 24 July 2023). The data are in NETCDF format and have a temporal resolution of 3 h and a horizontal spatial resolution of 0.1°. The dataset is produced by integrating the internationally available Princeton reanalysis data, GLDAS (Global Land Data Assimilation System) data, GEWEX-SRB (Global Energy and Water Exchanges, Surface Radiation Budget) radiation data, and TRMM (Tropical Rainfall Measuring Mission) precipitation data as the background field and integrating routine meteorological observation data from the China Meteorological Administration (CMA). In this study, precipitation rate (PREC) data (in mm/h) from 1988–2015 were selected, and to ensure consistency in the spatial resolution of the evapotranspiration data, the PREC data were bicubically interpolated to a spatial resolution of 0.25° × 0.25°.

2.4. Vegetation Optical Depth Dataset

Vegetation optical depth (VOD) data used in this study were obtained from the Global Long Term Microwave Vegetation Optical Depth Climate Archive (VODCA) [28,29,30].

VOD describes the attenuation of radiation by plants and is a function of frequency, vegetation water content, and biomass (by extension). VOD has many possible applications in studies of the biosphere, such as biomass monitoring, drought monitoring, phenology analyses, or fire risk management.

The dataset merged VOD observations from various spaceborne sensors (SSM/I, TMI, AMSR-E, AMSR2, and WindSat) to create global long-term VOD time series. Prior to aggregation, the data were rescaled to AMSR-E to remove systematic differences between the observations.

A product is also available for the C-band (~6.9 GHz, 2002–2018), X-band (10.7 GHz, 1997–2018), and Ku-band (~19 GHz, 1987–2017), and the data were globally sampled on a regular 0.25° grid. Each product is available as daily global netcdf4 files. In this study, products in the Ku-band from 1988–2015 were selected for the dataset.

2.5. Normalized Difference Vegetation Index Dataset

Normalized difference vegetation index (NDVI) data used in this study were obtained from the Global Inventory Monitoring and Modeling System (GIMMS) NDVI 3rd generation (3 g) dataset [31,32]. Vegetation indices are radiometric measures of photosynthetically active radiation absorbed by chlorophyll in the green leaves of vegetation canopies and are therefore good surrogate measures of the physiologically functioning surface greenness level of a region. Over 30 years, Compton J. Tucker created the NDVI time series within the framework of the GIMMS project and carefully assembled it from different AVHRR sensors, accounting for various deleterious effects, such as calibration loss, orbital drift, and volcanic eruptions [30]. The latest version of the GIMMS NDVI data set spans the period July 1981 to December 2015 and is termed NDVI3g (from AVHRR sensors), and its spatial resolution is 8 km.

The USGS Remote sensing phenology states that “NDVI values range from +1.0 to −1.0. Areas of barren rock, sand, or snow usually show very low NDVI values (for example, 0.1 or less). Sparse vegetation such as shrubs and grasslands or senescing crops may result in moderate NDVI values (approximately 0.2 to 0.5). High NDVI values (approximately 0.6 to 0.9) correspond to dense vegetation such as that found in temperate and tropical forests or crops at their peak growth stage”.

In this study, NDVI data from this dataset for 1988–2015 were selected and subjected to bicubic interpolated to a spatial resolution of 0.25° × 0.25°.

2.6. Standardized Precipitation Evapotranspiration Index Dataset

The standardized precipitation evapotranspiration index (SPEI) is a commonly used drought index. The SPEI data used in this study were obtained from the SPEI Global Database (http://spei.csic.es/database.html) (accessed on 14 September 2023) [33,34,35]. The global 0.5° gridded SPEI dataset is available under the Open Database License. Any rights to individual contents of the database are licensed under the Database Contents License. Users of the dataset are free to share, create, and adapt the data under the conditions of attribution and share-alike. The Global SPEI database, SPEIbase, offers long-term robust information on the drought conditions at the global scale, and it has a 0.5° spatial resolution and a monthly temporal resolution. Moreover, it has a multi-scale character and provides SPEI timescales between 1 and 48 months. The SPEI expresses deviations in the current climatic balance (precipitation minus evapotranspiration potential) as a standardized variate (zero mean and unit variance) with respect to the long-term balance. The reference period for the calculation in SPEIbase corresponds to the whole study period. Being a standardized variate means that the SPEI condition can be compared across space and time. Calculation of the evapotranspiration potential in SPEIbase is based on the FAO-56 Penman-Monteith method. Data are the float type, and the unit is the z-value (standard deviation) [36,37]. No land pixels were assigned a value of 1.0 × 10³⁰. In some rare cases, a good fit to the log-logistic distribution could not be achieved, resulting in a NAN (not a number) value in the database. Dimensions of the dataset: lon = 720; lat = 360; time = 1356. Resolution of the dataset: lon = 0.5°; lat = 0.5°; time = 1 month.

In this study, we selected SPEIBase v2.9. version (based on the CRU TS v4.07 dataset) of SPEI data on a 1-month scale from 1988–2015 and bicubic interpolated it to a spatial resolution of 0.25° × 0.25°.

2.7. Statistical Analysis Strategy

Budyko’s theoretical framework was proposed by Budyko, a prominent Soviet climatologist, based on the principle of coupled water-energy coupling balance, which suggests that evapotranspiration (E) within a watershed is limited by both energy (solar radiation, usually replaced by potential evaporation, Ep) and water (precipitation, P) over a long time span [38,39,40], as follows:

$${\displaystyle \frac{\mathrm{E}}{\mathrm{P}}}=\mathrm{f}\left({\displaystyle \frac{{\mathrm{E}}_{\mathrm{P}}}{\mathrm{P}}}\right)$$

(1)

The Budyko-Fu model is an extension and improvement of the Budyko hypothesis. It incorporates factors other than precipitation and potential evapotranspiration into the analysis of water and heat coupling balance in watersheds by introducing parameter $\mathsf{\omega }$, thereby improving the applicability and accuracy of the model in different watersheds. The Budyko-Fu model is as follows:

$${\displaystyle \frac{\mathrm{E}}{\mathrm{P}}}=1+{\displaystyle \frac{{\mathrm{E}}_{\mathrm{P}}}{\mathrm{P}}}-{\left[1+{\left({\displaystyle \frac{{\mathrm{E}}_{\mathrm{P}}}{\mathrm{P}}}\right)}^{\mathsf{\omega }}\right]}^{{\displaystyle \frac{1}{\mathsf{\omega }}}}$$

(2)

where E is actual evapotranspiration, Ep is potential evapotranspiration, P is precipitation, and $\mathsf{\omega }$ is the control parameter.

As the Budyko framework applies to a longer time span, to fully utilize the research data in calculating the true values of parameter $\mathsf{\omega }$, the evapotranspiration, potential evapotranspiration, and precipitation data from 1988–2015 were processed in this study by determining a sliding average every 5 years. These data were then introduced into Equation (2) to obtain the 5-year sliding averages of parameter $\mathsf{\omega }$ for the period 1988–2011.

After obtaining the long-term time series of parameter $\mathsf{\omega }$, this study further considered the factors that should be selected to establish a dynamic estimation model for parameter $\mathsf{\omega }$. We divided the influencing factors of parameter $\mathsf{\omega }$ into two categories: underlying surface influencing factors and climate influencing factors. Among the underlying surface factors, we considered the following factors: VOD and NDVI. As a traditional optical remote sensing indicator, the NDVI is usually limited to monitoring vegetation dynamics in the green canopy and thus mainly reflects the photosynthetic capacity of vegetation. VOD is different from traditional optical remote sensing indicators, which are directly proportional to the total vegetation water content of the aboveground biomass, including leafy and woody components. Rather, VOD is proportional to the aboveground biomass, including foliar and woody components, and is therefore sensitive to both photosynthetic and non-photosynthetic aboveground biomass [41,42,43]. In terms of climate, the factors we considered included SPEI and temperature (TMP). SPEI characterizes the extent to which dry and wet conditions in an area deviate from the norm through the difference between standardized potential evapotranspiration and precipitation [44]. When the SPEI is negative, it indicates that at a given time scale, the moisture conditions in an area are drier than the historical average. A larger negative value indicates a more severe drought, while a positive value indicates that the moisture condition of an area is wetter than the historical average.

To better clarify the relationship between different non-normally distributed factors and parameter $\mathsf{\omega }$, Spearman’s correlation coefficient is used to characterize the correlations between factors. Spearman’s correlations can be regarded as a nonparametric version of Pearson’s correlations, which is calculated using the data sample values themselves. However, Spearman’s correlation coefficient is calculated using the data sample rank order values, and the specific calculation formula is shown as follows:

$$\mathsf{\rho }=1-{\displaystyle \frac{6\sum {\mathrm{d}}_{\mathrm{i}}^2}{\mathrm{n}({\mathrm{n}}^2-1)}}$$

(3)

where ${\mathrm{d}}_{\mathrm{i}}$ denotes the difference between the rank values of the i-th data pair, and $\mathrm{n}$ denotes the total number of observed samples. The value of the coefficient is between −1 and 1, with values closer to 1 or −1 indicating a stronger monotonic relationship between variables. Positive values indicate positive correlations, and negative values indicate negative correlations.

In the dynamic modeling process of parameter $\mathsf{\omega }$, this study used stepwise regression equations, which introduce all the explanatory variables considered into the regression equation one by one. Variables that have been introduced into the regression equation that lose their significance after the introduction of a new variable need to be eliminated from the equation. The introduction of a variable or the exclusion of a variable from the regression equation was subjected to an F-test to ensure that the regression equation contained only those variables that had a significant effect on the explanatory variables before the introduction of the new variable and to ensure that non-significant variables had been excluded. After performing the stepwise regression, the study verified the accuracy of the model using the K-fold cross-test and quantified the contribution of the significant variables to parameter $\mathsf{\omega }$ in each watershed using the standardized coefficient method and R² decomposition method. The datasets used in the research have been described in Section 2.2, Section 2.3, Section 2.4, Section 2.5, Section 2.6 and Section 2.7, and the data processing process in the study is implemented using Python. The overall flow of the study is shown in Figure 2.

3. Results

3.1. Time Series Analysis of Parameter $\mathsf{\omega }$

In the linear analysis, evapotranspiration, potential evapotranspiration, and precipitation in the three subwatersheds of the upper Yellow River showed an increasing trend year by year from 1988 to 2015 (Figure 3a–c), which is consistent with the results of Theil-Sen’s slope estimation and the Mann–Kendall test (

Table 1. Theil-Sen’s slope estimation and Mann–Kendall test for evapotranspiration, potential evapotranspiration, precipitation, and parameter $\mathsf{\omega }$ in three sub-watersheds.

Subwatershed	Variable	Trends	Sen’s Slope	Significance	Z
BLG	Ep	increase	1.175	highly significant	2.825
	E	increase	1.275	highly significant	3.378
	PREC	increase	6.015	highly significant	3.497
	$\mathsf{\omega }$	decrease	−0.008	statistically significant	–2.257
LGL	Ep	increase	1.430	highly significant	3.141
	E	increase	2.173	highly significant	4.366
	PREC	increase	2.860	statistically significant	2.035
	$\mathsf{\omega }$	increase	0.052	highly significant	3.696
LHT	Ep	increase	0.265	insignificant	0.968
	E	increase	0.883	insignificant	1.126
	PREC	increase	0.770	insignificant	0.652
	$\mathsf{\omega }$	increase	0.038	highly significant	2.604

). The potential evapotranspiration capacity of subwatersheds closer to the source of the Yellow River was weaker, while the precipitation was higher. Most of the watersheds before Lanzhou are located in semi-arid and semi-humid regions, where the regional mean annual precipitation is greater than the mean annual potential evapotranspiration and the actual evapotranspiration is close to the potential evapotranspiration. However, most of the LHT is located in semi-arid and arid regions, where the regional mean annual potential evapotranspiration is much greater than the precipitation and the actual evapotranspiration.

The evapotranspiration, potential evapotranspiration, and precipitation data of the three subwatersheds from 1988–2015 were calculated as 5-year sliding averages and then introduced into the Budyko-Fu model to obtain the values for parameter $\mathsf{\omega }$ for the three subwatersheds after determining the sliding averages from 1988–2011 (Figure 3d).

Time series data of parameter $\mathsf{\omega }$ with the Budyko-Fu formula were used to generate a Budyko plot (Figure 4), which shows that the values of parameter $\mathsf{\omega }$ calculated using the sliding average for each subwatershed are distributed within reasonable intervals. Both the BLG and the LGL are in energy-limited regions, and the evapotranspiration in the watershed is mainly limited by energy factors, such as solar radiation. Moreover, the values of parameter $\mathsf{\omega }$ in the BLG ranged from 3.5 to 4.5 and showed a downward trend, while those in the LGL had the highest values between 4.5 and 6.5 and showed an increasing trend. The LHT is in a water-limited area, and evapotranspiration is mainly limited by precipitation. Thus, the values of parameter $\mathsf{\omega }$ were relatively low and showed a trend of initially decreasing and then fluctuating.

3.2. Dynamic Estimation Model for Parameter $\mathsf{\omega }$

In this study, stepwise regression was used to model the dynamic estimation of parameters $\mathsf{\omega }$. To simultaneously consider the effects of subsurface and climate on parameter $\mathsf{\omega }$, the VOD and NDVI were selected as independent variables to characterize the underlying surface condition of the watershed, while the SPEI and TMP were selected as independent variables to characterize the climate characteristics of the watershed.

The time series changes of independent variable factors (Figure 5) showed that in the underlying surface, the VOD of the BLG presented a slight linear downward trend and had value of 0.69–0.72, and the NDVI showed a linear upward trend and had a value of 0.32–0.33. In the LGL, VOD and NDVI showed upward trends, in which the VOD values ranged from 0.81 to 0.85 and NDVI values ranged from 0.34 to 0.36. The LGL had the highest vegetation biomass and highest vegetation cover among the three subwatersheds. The LHT showed an increasing trend in both VOD and NDVI values, with VOD values ranging roughly from 0.42 to 0.48 and NDVI values ranging roughly from 0.19 to 0.17. In general, this watershed had the lowest biomass and vegetation cover (Figure 5a,b).

In terms of watershed climate characteristics (Figure 5c,d), most of the SPEI values of the three subwatersheds were between −1 and 1 and showed a general decreasing trend, indicating that the drought conditions in the upper Yellow River increased from 1988 to 2015 relative to the previous period. In terms of temperature, the lowest temperature was found in the BLG, with the mean annual temperature roughly ranging from −1 °C to −2.5 °C. Moreover, the mean annual temperature in the LGL roughly ranged from 1.5 °C to 2.5 °C. The highest temperature was found from LHT, with the mean annual temperature roughly ranging from 7 °C to 8 °C. Against the background of global warming, the temperatures in all three subwatersheds showed a gradual upward trend.

Theil-Sen’s slope estimation and Mann–Kendall test of the subwatershed’s underlying surface and climate factors showed that the linear trends and Sen’s slopes of the factors were consistent except for the VOD values of the BLG, which showed a slight downward trend in linearity and a non-significant upward trend in Sen’s slope estimation (

Table 2. The Theil-Sen’s Slope Estimation and MK Test of Underlying Surface and Climate Factors in Three Sub watersheds.

Subwatershed	Variable	Trends	Sen’s Slope	Significance	Z
BLG	VOD	increase	0.00021	insignificant	0.929
	NDVI	increase	0.00026	insignificant	1.620
	SPEI	decrease	−0.001	insignificant	−0.020
	TMP	increase	0.015	marginally significant	1.798
LGL	VOD	increase	0.00047	insignificant	1.245
	NDVI	increase	0.00049	statistically significant	2.963
	SPEI	decrease	−0.004	insignificant	−0.533
	TMP	increase	0.015	marginally significant	1.798
LHT	VOD	increase	0.002	highly significant	5.235
	NDVI	increase	0.000	highly significant	3.477
	SPEI	decrease	−0.018	statistically significant	−2.509
	TMP	increase	0.034	highly significant	2.904

). In general, the three subwatersheds of the upper Yellow River showed a general trend of increasing vegetation biomass, an increasing trend of vegetation cover, increasing trend in temperature, and a tendency towards meteorological drought.

To further clarify the relationship between underlying surface factors, climate factors, and parameter $\mathsf{\omega }$, we performed Spearman’s correlation analysis between independent variable factors and parameter $\mathsf{\omega }$ in the different watersheds (Figure 6). Parameter ω was positively correlated with VOD to different degrees in the three subwatersheds. Parameter ω in the BLG and the LGL is negatively correlated with SPEI, while ω in the LHT is positively correlated with SPEI. Parameter $\mathsf{\omega }$ in the BLG showed negative correlations with NDVI and TMP, and parameter $\mathsf{\omega }$ in the LGL and the LHT showed positive correlations with NDVI and TMP. Among the independent variables, the SPEI in the three subwatersheds always showed negative correlations with TMP and NDVI while TMP always showed positive correlations with NDVI.

With VOD, NDVI, SPEI, and TMP as independent variables, parameter $\mathsf{\omega }$ was used as the dependent variable to build a stepwise regression model. Each of the three subwatersheds was screened to identify different significant variables (

Table 3. Variables retained as independent variables with significant effects on parameter $\mathsf{\omega }$ in the stepwise regression models for the three subwatersheds.

Subwatershed	Independent Variable	p	VIF
BLG	VOD	0.005 *	1.224
	NDVI	0.000 *	1.643
	SPEI	0.001 *	1.883
LGL	VOD	0.002 *	1.301
	SPEI	0.000 *	1.353
	TMP	0.000 *	1.100
LHT	VOD	0.001 *	1.587
	SPEI	0.039 *	1.587

p is p-vlaue, VIF is Variance Influence Factor, and * represents the independent variable passing the significance test.

).

BLG: the significant retained variables were VOD, NDVI, and SPEI, and the stepwise regression model was as follows:

$$\mathsf{\omega }=9.548+9.904·\mathrm{V}\mathrm{O}\mathrm{D}-36.029·\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-1.018·\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}$$

(4)

LGL: the significant retained variables were VOD, NDVI, and TMP, and the stepwise regression model was as follows:

$$\mathsf{\omega }=-17.741+24.683·\mathrm{V}\mathrm{O}\mathrm{D}-3.354·\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}+1.414·\mathrm{T}\mathrm{M}\mathrm{P}$$

(5)

LHT: the significant retained variables were VOD and SPEI, and the stepwise regression model was as follows:

$$\mathsf{\omega }=-5.779+20.749·\mathrm{V}\mathrm{O}\mathrm{D}+1.113·\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}$$

(6)

3.3. Model Testing and Analysis of the Contribution of Independent Variables

The residual plot of the stepwise regression model (Figure 7) reveals that the R² of the stepwise regression model for the LGL was 0.8 while that of the model for the BLG was close to 0.6. The interpretability of the selected underlying surface and climate variables of the two watersheds to the parameters $\mathsf{\omega }$ was higher. The R² of the stepwise regression model in the LHT was approximately 0.4, and the model fitting effect was poorer than that in the other two subwatersheds. The Yellow River passes through the Loess Plateau in this watershed, and, based on the early agricultural development and abundant human activities, the natural environment has been greatly influenced by humans. Therefore, the hydrological cycle process in this watershed cannot fully satisfy the coupling of water and energy, and the measured data do not reflect the actual condition of this watershed. Therefore, when considering the influencing factors of parameter $\mathsf{\omega }$ in this subwatershed, further indicators related to human activities can be added to improve the explanatory power of the model for parameter $\mathsf{\omega }$.

The predictive ability of the Budyko-Fu model, which includes the dynamic parameter $\mathsf{\omega }$ for watershed evapotranspiration was determined (Figure 7). The R² between the predicted and real values of evapotranspiration in the BLG and the LGL both reached values greater than 0.9, and the R² between the predicted and real values of evapotranspiration in the LHT reached values greater than 0.8. This indicates that the estimation model of parameter $\mathsf{\omega }$ established in this study can be used for dynamic predictions of the actual evapotranspiration in these watersheds.

The mean absolute error (MAE) and mean square error (MSE) of the model were calculated, and the results showed that the MSE and MAE values of the stepwise regression models for the three subwatersheds were small, indicating that the models were well fit (

Table 4. Model accuracy and generalization ability test table.

Subwatershed	Accuracy Type	MSE	MAE
BLG	Training Accuracy	0.007	0.054
	k-fold inspection accuracy	0.014	0.268
LGL	Training Accuracy	0.055	0.166
	k-fold inspection accuracy	0.099	0.268
LHT	Training Accuracy	0.132	0.281
	k-fold inspection accuracy	0.368	0.481

). The model was tested for overfitting using the 6-fold cross-test, and the mean values of MSE and MAE were slightly higher than those of the original stepwise regression model but within the acceptable range. This indicates that the model has a strong generalization ability.

To further understand the influence of each underlying and climate factor in the model on parameter $\mathsf{\omega }$, we used the standardized coefficient method (Figure 8a) and R² decomposition method (Figure 8b) to quantify the contribution degree of each factor in the model to parameter $\mathsf{\omega }$. The total contribution of the two underlying surface actors VOD and NDVI to parameter $\mathsf{\omega }$ in the BLG reached 64% and 76% using the two methods, respectively, whereas the climate factor SPEI contributed 36% and 24% to parameter $\mathsf{\omega }$, which indicates that the underlying vegetation factors in the BLG play a dominant role in parameter $\mathsf{\omega }$. In the LGL, the contribution degree of VOD to parameter $\mathsf{\omega }$ was 25% and 7% in the two methods, while the total contribution degree of the two climate factors SPEI and TMP to parameter $\mathsf{\omega }$ reached 75% and 93%, respectively, indicating that climate factors play a dominant role in parameter $\mathsf{\omega }$ in this watershed. In the LHT, only VOD and SPEI had a significant effect on parameter $\mathsf{\omega }$, and the contribution degree of VOD to parameter $\mathsf{\omega }$ was higher than that of SPEI in the two methods. Therefore, the vegetation biomass of this watershed played a dominant role in parameter $\mathsf{\omega }$, followed by drought conditions.

4. Discussion

According to the theory of the Budyko-Fu model, in areas with poor permeability of the underlying surface, limited vegetation, large topographic slopes, and thus strong surface runoff, the values of parameter $\mathsf{\omega }$ are small. On the contrary, in areas where the underlying surface has good permeability, more vegetation, a flat topography, and thus low surface runoff, the values of parameter $\mathsf{\omega }$ are large [16].

In this study, most of the BLG and the LGL are located on the Qinghai Tibet Plateau, which has a relatively flat terrain and is less affected by human activities. The vegetation coverage in the BLG is higher than that in the LGL, with more precipitation and less evapotranspiration, making it relatively wetter. The majority of the LHT is located on the Loess Plateau, with significant terrain fluctuations and lower vegetation coverage than the other two watersheds. It has the least precipitation and the highest potential evapotranspiration, making it relatively drier, and it is more affected by human activities. An analysis of our research results revealed that parameter $\mathsf{\omega }$ in the watershed before Lanzhou was greater than that in the LTH, which is consistent with the prediction of Fu’s theory and verifies the rationality of the research results.

This study analyzed parameter $\mathsf{\omega }$ with underlying surface and climate factors, and the results indicated that VOD exhibits varying degrees of positive correlations with parameter $\mathsf{\omega }$ in the three watersheds. As a representative of aboveground biomass of vegetation, VOD usually represents the development degree of underground roots of vegetation. The more developed the root system, the stronger the ability of the vegetation to absorb water and the stronger the transpiration; thus, the corresponding parameter $\mathsf{\omega }$ will be greater. The lower the frequency, the greater the penetration ability of the VOD product, with higher frequency VOD being most sensitive to leaves and branches (X and Ku bands) and lower frequency VOD showing increased sensitivity to the trunk (L and C bands). The VOD product used in this study was VODCA-Ku. Estimating total aboveground biomass using the lower frequency L-VOD is relatively more reliable than Ku-VOD, but L-VOD has been used in China for a relatively short period of time and is subject to relatively severe RF interference. We plan to test and apply the L-band data in future studies when they have accumulated a long enough time series to improve the incompleteness of the study. In contrast, SPEI showed a negative correlation with parameter $\mathsf{\omega }$ of the three watersheds. The SPEI values of the three watersheds generally decreased, which indicates that under the background of global warming, although precipitation increases, potential evapotranspiration and actual evapotranspiration also increase. Thus, the rate of water demand may exceed the rate of precipitation increases, resulting in a decrease in SPEI values. In addition, global warming may also lead to an increase in the spatial and temporal heterogeneity of precipitation distribution, which may trigger droughts in specific time periods or regions, thus leading to low SPEI values [45]. Accordingly, the lower the SPEI, the drier the region is relative to the historical period, which indicates limited transpiration of vegetation and decreases in evapotranspiration. Thus, SPEI showed a negative correlation with parameter $\mathsf{\omega }$.

When modeling parameter $\mathsf{\omega }$, the significant variables that were retained in the BLG were VOD, NDVI, and SPEI. According to Shen et al. [46], the watershed BLG is the source area of the Yellow River and is located in the hinterland of the Tibetan Plateau, where grasslands account for 80% of the surface cover and changes in vegetation cover have a very important influence on the hydrological process of the whole Yellow River watershed. The results of this study showed that there was no significant trend in the VOD values in this watershed, although the NDVI values showed an increasing trend. We hypothesized that this might be due to the fact that although the overall biomass of alpine grassland ecosystems in this watershed is not sensitive to changes in hydrothermal conditions, the start of the growing season (SOS) of grassland vegetation under the influence of warmer temperatures and more precipitation is advanced and the end of the growing season (EOS) is delayed [47], which leads to an increase in the NDVI. Parameter $\mathsf{\omega }$ showed a decreasing trend and a negative correlation with the NDVI, which is consistent with the study of Li et al. [48].

In the analysis of the watershed from LGL, the significant variables that were retained were VOD, SPEI, and TMP, and climate factors had the greatest total contribution to the change of parameter $\mathsf{\omega }$. First, this region is an energy-limited watershed, and evapotranspiration in the watershed is mainly limited by energy factors, such as solar radiation; second, this watershed flows through the Qilian Mountains, and its vegetation types are grassland and woodland and VOD and NDVI values are also the highest among the three watersheds. Du et al. [49] demonstrated that the timberline of this area is very sensitive to the climate change response, and, under the background of global warming, the vegetation type in this watershed is transiting from grassland to woodland [50], which corresponds to the increasing trend of VOD values in this study. Therefore, we hypothesize that climate factors dominate the change of parameter $\mathsf{\omega }$ in this watershed for two reasons: the lower temperatures in the watershed limit the evapotranspiration process and the climate in this region has a significant influence on the vegetation. Thus, changes in parameter $\mathsf{\omega }$ were affected by the combination of climate and vegetation.

In the LHT, the significant variables retained were VOD and SPEI, and the contribution of VOD to the change of parameter $\mathsf{\omega }$ was larger than that of SPEI. This watershed flows through the Loess Plateau and Inner Mongolia Plateau, which presents sparse vegetation and fragile ecosystems; therefore, parameter $\mathsf{\omega }$ was sensitive to the response of vegetation. In all three watersheds, the variables that had a significant effect on parameter $\mathsf{\omega }$ included VOD and SPEI; therefore, the inclusion of VOD and SPEI in the dynamic estimation model of parameter $\mathsf{\omega }$ is generalizable to the whole upper Yellow River watershed. The NDVI was also retained in the BLG, and TMP was retained in the LGL, which reflects the differentiation among the three subwatersheds.

This study has made some achievements in clarifying the dynamic changes of parameter $\mathsf{\omega }$ and its influencing factors in the upper reaches of the Yellow River Basin; however, certain limitations were observed, which provide a further direction for future research. First, the study area was limited to the upper reaches of the Yellow River Basin. However, the hydrological process of the lower reaches may be significantly different from that of the upper reaches due to the more significant impact of human activities. In subsequent studies, factors related to human activities can be included in the dynamic estimation model of parameter $\mathsf{\omega }$.

Second, this study only considered the relationship between the underlying surface and parameter $\mathsf{\omega }$ by using factors related to vegetation; thus, the effects of soil characteristics, topography, and other factors were generally ignored. Soil and terrain factors may have an important impact on water infiltration, surface runoff, and evapotranspiration. In future research, these factors should be comprehensively considered to build a more comprehensive hydrological model.

Of course, this study also clearly demonstrated the influence of vegetation biomass on the crucial parameter $\mathsf{\omega }$ in the hydrological cycle of the basin and revealed the higher applicability of VOD in explaining the variation of parameter $\mathsf{\omega }$ compared with the NDVI. Therefore, the study provides an important scientific basis for understanding the role of vegetation biomass in the hydrological cycle of the basin, which is of great practical significance in guiding water resources management and ecological protection in the Yellow River Basin and other similar basins.

5. Conclusions

In summary, this study used the Buydko-Fu model to calculate parameter $\mathsf{\omega }$ for three subwatersheds in the upper reaches of the Yellow River and established dynamic estimations model of parameter $\mathsf{\omega }$ through stepwise regression. Moreover, this study also quantified the contribution of underlying surface factors and climate factors to parameter $\mathsf{\omega }$ in the model. The results showed: (1) Both the BLG and the LGL are in energy-limited regions, the LHT is in a water-limited area. (2) Parameter $\mathsf{\omega }$ of the three subwatersheds was significantly correlated with VOD and SPEI, which indicates that the parameter is affected by both vegetation and climate and indicates that the inclusion of VOD and SPEI in the dynamic estimation models of parameter $\mathsf{\omega }$ in the upper reaches of the Yellow River are somewhat generalizable. In addition, VOD has a higher applicability in the explanation of the change of parameter $\mathsf{\omega }$ than the NDVI. (3) Furthermore, parameter $\mathsf{\omega }$ for the BLG was also affected by the NDVI while that for the LGL was affected by TMP, thus revealing the differences in the ecohydrological responses of different basins in the upper Yellow River Basin.