Quantitatively Calculating the Contribution of Vegetation Variation to Runoff in the Middle Reaches of Yellow River Using an Adjusted Budyko Formula

: The middle reaches of the Yellow River (MRYR) are a key area for carrying out China’s vegetation restoration project. However, the impact of vegetation variation on runoff in the MRYR is still unclear. For quantitatively evaluating the contribution rate of vegetation variation to runoff in the MRYR, this paper quantiﬁed the relationship between Normalized Difference Vegetation Index (NDVI) and Budyko parameters ( w ). Then, we used multiple linear regression to quantitatively calculate the contribution rate of different factors on vegetation variation. Finally, an adjusted Budyko formula was constructed to quantitatively calculate the inﬂuence of vegetation variation on runoff. The results showed that there is a linear relationship between NDVI and Budyko parameters ( w ) ( p < 0.05); the ﬁtting parameter and constant term were 12.327 and − 0.992, respectively. Vegetation change accounted for 33.37% in the MRYR. The contribution of climatic and non-climatic factors on vegetation change is about 1:99. The contribution of precipitation, potential evaporation, anthropogenic activities on the runoff variation in the MRYR are 23.07%, 13.85% and 29.71%, respectively.


Introduction
In the past few decades, many ecological protection projects carried out in China have led to a significant increase in vegetation [1,2]. The middle reaches of the Yellow River (MRYR) are the key area for the implementation of the vegetation restoration project in China [3,4]. Vegetation index is a simple, effective and empirical measure of surface vegetation, so the Normalized Difference Vegetation Index (NDVI) was used to characterize the vegetation change in the MRYR. Some studies show that NDVI in the MRYR has increased significantly, and its vegetation coverage has been restored [5][6][7].
Vegetation change can significantly change runoff by changing hydrological processes, such as vegetation transpiration and interception evaporation, and then affect the available amount of watershed runoff [8][9][10]. The debate on the relationship between "vegetation and water" dates back to at least the mid-19th century [11]. Although most studies show that the increase in vegetation will lead to a decrease in water yield (called negative impact here), only a few studies show that the increase in vegetation has little impact on water yield [12,13] (called no impact here) and even leads to an increase in water yield [14,15] (called positive impact here). Interestingly, the areas with positive and no impact are mostly in humid areas and large watersheds with complex terrain. In the past few decades, the spatial variability and periodicity of precipitation and runoff in the MRYR have changed [16][17][18][19]. Therefore, quantitatively calculating the influence of vegetation variation on runoff variation is helpful to deeply understand the response mechanism of the water cycle process to the underlying surface condition changes, and has theoretical and practical significance for improving watershed water resources management measures. It is useful for the utilization and development of water resources in the Yellow River, the protection of the human living environment and the sustainable development of the economy.
Many scholars have carried out studies for quantitatively calculating the influence of vegetation variation on runoff variation in the MRYR [20][21][22]. The runoff decreased with the improvement of vegetation coverage, seen by analyzing the statistical relationship between vegetation coverage and runoff in most basins in the MRYR [23]. Liang et al. [24] quantitatively analyzed the influence of ecological restoration measures for the streamflow variation of most watersheds on the Loess Plateau from 1961 to 2009, based on Budyko theory, and found that ecological restoration measures are the leading factor for the decline in runoff in most watersheds. Li et al. [25] found that the Normalized Difference Vegetation Index (NDVI) in most basins in the Loess Plateau increased significantly after 2000, while the runoff decreased significantly. A higher evapotranspiration rate from restored vegetation is the primary reason for the reduced runoff coefficient. Zhang et al. [26] quantitatively evaluated the influence of vegetation cover variation on streamflow variation in the MRYR over the past 37 years. It was shown that vegetation variation is the leading driving force of streamflow decline from 2000 to 2010. Wang et al. [27] discovered that the vegetation coverage in the MRYR increased by 29.72% from 1982 to 2018. The average runoff at each station showed a downward trend from 1961 to 2018, and ecological restoration played an important role in the decline in runoff in the MRYR. However, the contribution of vegetation variation on streamflow variation in the MRYR has not been calculated quantitatively.
Therefore, in order to quantitatively evaluate the contribution rate of vegetation variation on streamflow variation in the middle reaches of the Yellow River (MRYR), this paper analyzes its hydro-meteorological characteristics from 1982 to 2015, identifies the mutation year of streamflow, and quantifies the relationship between the NDVI and a Budyko parameter (w). Then, multiple linear regression is adopted to distinguish the contribution rate of different factors on vegetation variation. Finally, an adjusted Budyko formula is constructed to quantitatively calculate the influence of vegetation variation on runoff in the MRYR. This work is useful to understand the influence of vegetation variation on the evolution of water resources in the MRYR, and also supply important insights for water resources management.

Research Region
The middle reaches of the Yellow River (MRYR) are in the area of the Toudaoguai-Huayuankou hydrological station ( Figure 1). The length of this reach is about 1234.6 km, its area accounts for about 45.7% and its runoff accounts for about 44.3% of the Basin. Guanzhong Plain and Hetao Plain are in the MRYR, which produce a large number of agricultural products every year. The wheat output ranks first, soybean output ranks second, and cotton output ranks fourth in China. The decreasing runoff in the MRYR has brought great harm to agricultural production and food security.

Data
(1) The annual runoff observation data of Toudaoguai and Huayuankou hydrometric stations from 1982 to 2015 were obtained from the hydrological statistical yearbook of China. In this study, the annual runoff data in the MRYR was equal to a value, which was the annual measured runoff of Huayuankou hydrometric station minus the annual measured runoff data of Toudaoguai hydrological station.
(2) 82 meteorological stations data  were downloaded from China Meteorological Administration. We adopted the Penman-Monteith formula for calculating the potential evaporation of 82 meteorological stations, and then the annual precipita-tion and potential evaporation were interpolated by Kriging. Finally, the average annual precipitation and potential evaporation in the MRYR can be obtained.
(3) The NDVI dataset (8 × 8 km 2 ) was sourced from the National Aeronautics and Space Administration. The time span of the dataset is from July 1981 to December 2015, its time resolution is 15 d, and its downloaded files are in NC4 format. The annual NDVI are calculated by the average value composites method, and is equal to the arithmetic mean of all grids in the study area. This method further eliminates some interference of cloud, atmosphere, solar altitude, and sensor sensitivity. Finally, the annual NDVI from 1982 to 2015 can be obtained.

Data
(1) The annual runoff observation data of Toudaoguai and Huayuankou hydrometric stations from 1982 to 2015 were obtained from the hydrological statistical yearbook of China. In this study, the annual runoff data in the MRYR was equal to a value, which was the annual measured runoff of Huayuankou hydrometric station minus the annual measured runoff data of Toudaoguai hydrological station.
(2) 82 meteorological stations data  were downloaded from China Meteorological Administration. We adopted the Penman-Monteith formula for calculating the potential evaporation of 82 meteorological stations, and then the annual precipitation and potential evaporation were interpolated by Kriging. Finally, the average annual precipitation and potential evaporation in the MRYR can be obtained.
(3) The NDVI dataset (8 × 8 km 2 ) was sourced from the National Aeronautics and Space Administration. The time span of the dataset is from July 1981 to December 2015, its time resolution is 15 d, and its downloaded files are in NC4 format. The annual NDVI are calculated by the average value composites method, and is equal to the arithmetic mean of all grids in the study area. This method further eliminates some interference of cloud, atmosphere, solar altitude, and sensor sensitivity. Finally, the annual NDVI from 1982 to 2015 can be obtained.

Trend Analysis Method (Slope)
The formula of Slope method is as follows [28,29]: If Slope is greater than 0, the variable shows ever-growing properties. If Slope < 0, it implies that the variable gets less with the day; n means the number of years; X indicates the annual precipitation, and potential evaporation runoff.

Mann-Kendall (MK) Mutation Detection Algorithm
For X containing n data, we first construct a rank sequence. Then, we calculated S k .
E(S k ) represents the average value of S k , Var(S k ) represents the variance of S k . Finally, UF 1 = 0, the Uf k statistics can be calculated.
Constructing the reverse order of X, and making UB k = UF k (k = n, n − 1, . . ., 1), UB 1 = 0, we calculated the UB k statistics in this way.

Budyko Hypothesis
There are two hypotheses in the attribution analysis of runoff change using Budyko hypothesis: (1) assuming that human activities, climate factors and vegetation do not affect each other, and are independent factors; (2) assuming that the base period is only affected by climate factors.
The multi-year scale water balance leads to Formula (8).
where, Q represents runoff depth; P represents precipitation; E represents actual evaporation. According to the Buduko hypothesis, there is a nonlinear function relation between E/P and dryness ratio (phi = E 0 /P) of a specific watershed.
E 0 is the potential evaporation (mm). The Fu-type Budyko equation is [30]: w represents the variable representing the underlying surface information of the basin, and is used to characterize human activities. Human activities can affect the runoff change of the Yellow River Basin through many aspects, including vegetation change, water conservancy project construction, domestic water for urban residents and agricultural irrigation water, etc. The larger the value of w, the less precipitation is transformed into runoff.
Li et al. [31] found that there is a good linear relationship between NDVI and Budyko parameter (w) in the 26 major global river basins, which makes it possible to quantitatively evaluate the contribution rate of vegetation variation on runoff change.
a and b are constants for a specific watershed and can be calculated by the least squares method.
The elasticity coefficient is used to measure the relative change of another variable caused by the change of one variable. ε P , ε E0 , ε w and ε NDV I represent the elastic coefficient of P, E 0 , w and NDVI on runoff, respectively.
∆R P , ∆R E0 , ∆R w , ∆R NDV I and ∆R H represent the value of runoff change caused by precipitation, potential evapotranspiration, underlying surface change, vegetation change and human activities.
∆P, ∆E 0 , ∆w and ∆NDV I represent the changes of multi-year average P, E 0 , w and NDVI from T 1 to T 2 period.
∆NDV I is attributed to NDVI change value caused by climate factors (∆NDV I c ) and human activities (∆NDV I H ): Average NDVI is significantly correlated with P and E 0 [32][33][34], so we calculated the multiple linear regression equation between NDVI, P and E 0 in T 1 period. There are two assumptions in the attribution analysis of vegetation change using multiple linear regression method: (1) Human activities and climate factors are independent factors. (2) The base period is only affected by climate factors. Therefore, except for climate factors, other factors affecting vegetation change are classified as human activities.
where, NDV I T2,S represents average NDVI value under climate change in the T 2 period. NDV I T1 and NDV I T2 represent average NDVI value of the two periods, respectively. η NDV I H and η NDV I C represent the influence of anthropic and climatic factors on vegetation variation.

Mutation Year Identification
Because the Mann-Kendall mutation test algorithm is not disturbed by a few outliers and has the advantages of simple calculation, it is widely used to identify the mutation points of hydro-meteorological elements [35][36][37][38].
In the Mann-Kendall mutation detection algorithm, UF statistics are a sequence of statistics calculated by the order of time series X, and UB statistics are a sequence of statistics calculated by the inverse order of time series X. In practical application, it is generally considered that the result is credible at the significance level of 0.05. If the UF statistics line and UB statistics line have intersections and these intersections are within the range

Mutation Year Identification
Because the Mann-Kendall mutation test algorithm is not disturbed by a few outliers and has the advantages of simple calculation, it is widely used to identify the mutation points of hydro-meteorological elements [35][36][37][38].
In the Mann-Kendall mutation detection algorithm, UF statistics are a sequence of statistics calculated by the order of time series X, and UB statistics are a sequence of statistics calculated by the inverse order of time series X. In practical application, it is generally considered that the result is credible at the significance level of 0.05. If the UF statistics line and UB statistics line have intersections and these intersections are within the range of 0.05, this is significant. The year corresponding to the intersection is the mutation year of the time series data.
UF statistics of runoff time series data in the MRYR intersected with UB statistics in 1989 (Figure 3), and this intersection is within the range of 0.05 (significant level), indicating that 1989 is a runoff mutation year in the MRYR.

The Relationship between NDVI and w
In Budyko's hypothesis, the parameter (w) reflects the change in the underlying surface, so the parameter (w) should change on an annual or multi-year scale. For the Fu-type Budyko equation, if the elements (Q, P, E0) in the equation can be determined, the parameter (w) can be solved. Using the hydrological and meteorological data in the MRYR from 1982 to 2015, the parameters (w) from 1982 to 2015 are obtained. Figure 5 displayed the variation tendency of NDVI and underlying surface parameters (w) in the MRYR from 1982 to 2015. The NDVI and underlying surface parameters (w) in the middle reaches of the Yellow River increased significantly (p < 0.05), and the slopes were 0.017/10a and 0.268/10a, respectively.

The Relationship between NDVI and w
In Budyko's hypothesis, the parameter (w) reflects the change in the underlying surface, so the parameter (w) should change on an annual or multi-year scale. For the Fu-type Budyko equation, if the elements (Q, P, E0) in the equation can be determined, the parameter (w) can be solved. Using the hydrological and meteorological data in the MRYR from 1982 to 2015, the parameters (w) from 1982 to 2015 are obtained. Figure 5 displayed the variation tendency of NDVI and underlying surface parameters (w) in the MRYR from 1982 to 2015. The NDVI and underlying surface parameters (w) in the middle reaches of the Yellow River increased significantly (p < 0.05), and the slopes were 0.017/10a and

The Relationship between NDVI and w
In Budyko's hypothesis, the parameter (w) reflects the change in the underlying surface, so the parameter (w) should change on an annual or multi-year scale. For the Fu-type Budyko equation, if the elements (Q, P, E 0 ) in the equation can be determined, the parameter (w) can be solved. Using the hydrological and meteorological data in the MRYR from 1982 to 2015, the parameters (w) from 1982 to 2015 are obtained. Figure 5 displayed the variation tendency of NDVI and underlying surface parameters (w) in the MRYR from 1982 to 2015. The NDVI and underlying surface parameters (w) in the middle reaches of the Yellow River increased significantly (p < 0.05), and the slopes were 0.017/10a and 0.268/10a, respectively.

Quantifying the Influence of Vegetation Variation on Runoff
Because we want to quantitatively calculate the contribution rate of different factors on runoff change (Figure 7), 1982-2015 is separated into T1 (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) and T2 , according to the catastrophe analysis results of runoff data in the MRYR. Firstly, we analyzed the correlation between NDVI, P and E0 (Figure 8), implying that NDVI is significantly correlated with P and E0. Next, we calculated the linear function relation between NDVI, P and E0 in T1 period. Finally, we analyzed the contribution rate of different factors on vegetation change, using formulas 27-31 ( Table 1). The results showed that the contribution of climatic and non-climatic factors on vegetation variation is close to 1:99. In this paper, the values of underlying surface parameters (w) and NDVI are processed by a 10-year moving average, and then the regression coefficients (a and b) are obtained by least square fitting, a = 12.327, b = −0.992. The R 2 of equation fitting is 0.4465, and is significant (p < 0.05), indicating that the fitting result is good ( Figure 6).

Quantifying the Influence of Vegetation Variation on Runoff
Because we want to quantitatively calculate the contribution rate of different factors on runoff change (Figure 7), 1982-2015 is separated into T1 (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) and T2 , according to the catastrophe analysis results of runoff data in the MRYR. Firstly, we analyzed the correlation between NDVI, P and E0 (Figure 8), implying that NDVI is significantly correlated with P and E0. Next, we calculated the linear function relation between NDVI, P and E0 in T1 period. Finally, we analyzed the contribution rate of different factors on vegetation change, using formulas 27-31 (Table 1). The results showed that the contribution of climatic and non-climatic factors on vegetation variation is close to 1:99.

Quantifying the Influence of Vegetation Variation on Runoff
Because we want to quantitatively calculate the contribution rate of different factors on runoff change (Figure 7), 1982-2015 is separated into T 1 (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) and T 2 (1990-2015), according to the catastrophe analysis results of runoff data in the MRYR. Firstly, we analyzed the correlation between NDVI, P and E 0 (Figure 8), implying that NDVI is significantly correlated with P and E 0 . Next, we calculated the linear function relation between NDVI, P and E 0 in T 1 period. Finally, we analyzed the contribution rate of different factors on vegetation change, using formulas 27-31 ( Table 1). The results showed that the contribution of climatic and non-climatic factors on vegetation variation is close to 1:99.   Shown in Table 2, compared with the T1 period, the potential evaporation in the T2 period increased by 41.70 mm, precipitation and runoff depth decreased by 25.89 mm and 24.49 mm, respectively, underlying surface parameters (w) and NDVI increased by 0.52 and 0.023, respectively.     Shown in Table 2, compared with the T1 period, the potential evaporation in the T2 period increased by 41.70 mm, precipitation and runoff depth decreased by 25.89 mm and 24.49 mm, respectively, underlying surface parameters (w) and NDVI increased by 0.52 and 0.023, respectively.   Shown in Table 2, compared with the T 1 period, the potential evaporation in the T 2 period increased by 41.70 mm, precipitation and runoff depth decreased by 25.89 mm and 24.49 mm, respectively, underlying surface parameters (w) and NDVI increased by 0.52 and 0.023, respectively.  Table 3 displays the influence of climatic and anthropogenic factors on the annual runoff of the post-change period. Contrasted to T 1 , runoff variation amount caused by precipitation, potential evapotranspiration, vegetation variation and anthropogenic activities in the T 2 is −5.42 mm, −3.25 mm, −7.84 mm, and −6.98 mm. Their contributions were 23.07%, 13.85%, 33.37%, and 29.71%, respectively. Table 3. Attribution analysis of runoff variation in MRYR.

Conclusions and Discussion
Using NDVI, meteorological data and measured runoff data, based on an adjusted Fu-type Budyko equation, this paper quantitatively analyzed the contribution of vegetation variation on runoff in the MRYR from 1982 to 2015. The results showed that the influence of climatic and anthropogenic factors on vegetation variation is about 1:99. Vegetation change played a major role on runoff in the MRYR. Precipitation, potential evaporation and human activities account for 23.07%, 13.85% and 29.71% in runoff variation, respectively.
Although this paper carries out strict quality control on the data and model, some non-determinacy remains. Firstly, precipitation and potential evaporation data are interpolated from the data of meteorological stations in and around the MRYR, which will bring some non-determinacy. Secondly, the precondition of the study is that all variables are independent and do not affect each other. For example, when calculating the partial derivative of multi-year evapotranspiration to runoff, it is required that the factors, such as vegetation change, precipitation and human activities, are independent. However, in practice, the underlying surface and climate system are a whole. These factors interact and affect each other, which will have some non-determinacy on the research results.
This paper quantitatively calculated the influence of vegetation variation on runoff in the MRYR, but in order to reveal the influence mechanism of vegetation on water cycle, it is necessary to analyze how vegetation affects the hydrological process of "precipitationinterception-infiltration-evapotranspiration-runoff". In the follow-up study, we will combine the distributed hydrological model HLMS with the PML model, coupling leaf area index information to build a distributed Hydrothermal Coupling Model [39][40][41], for clarifying the response of water heat balance to vegetation variation. In addition, in the process of quantitatively calculating the contribution rate of different factors on runoff change, in addition to climate factors and vegetation change, water conservancy project construction, urban residents' domestic water and agricultural irrigation water are classified as human activity factors. We will quantitatively calculate the contribution rate of urban water, industrial water, agricultural water and water conservancy project construction on runoff change in the follow-up study.