Quantifying the Impact of Climate Change and Human Activities on Streamflow in a Semi-Arid Watershed with the Budyko Equation Incorporating Dynamic Vegetation Information

Understanding hydrological responses to climate change and land use and land cover change (LULCC) is important for water resource planning and management, especially for water-limited areas. The annual streamflow of the Wuding River Watershed (WRW), the largest sediment source of the Yellow River in China, has decreased significantly over the past 50 years at a rate of 5.2 mm/decade. Using the Budyko equation, this study investigated this decrease with the contributions from climate change and LULCC caused by human activities, which have intensified since 1999 due to China’s Grain for Green Project (GFGP). The Budyko parameter that represents watershed characteristics was more reasonably configured and derived to improve the performance of the Budyko equation. Vegetation changes were included in the Budyko equation to further improve its simulations, and these changes showed a significant upward trend due to the GFGP based on satellite data. An improved decomposition method based on the Budyko equation was used to quantitatively separate the impact of climate change from that of LULCC on the streamflow in the WRW. Our results show that climate change generated a dominant effect on the streamflow and decreased it by 72.4% in the WRW. This climatic effect can be further explained with the drying trend of the Palmer Severity Drought Index, which was calculated based only on climate change information for the WRW. In the meantime, although human activities in this watershed have been very intense, especially since 1999, vegetation cover increase contributed a 27.6% decline to the streamflow, which played a secondary role in affecting hydrological processes in the WRW.


Introduction
Climate change and land use and land cover change (LULCC) have had profound influences on global and regional hydrological processes [1][2][3].Understanding the hydrological responses in watersheds to climate change and LULCC is important for water resource planning and management throughout the world, especially in arid and semi-arid areas where water is the primary limiting factor for environmental services and social development [4][5][6].Climate change causes changes in different components of hydrological processes [2,7].These components include evapotranspiration, infiltration, streamflow, soil moisture, etc. Global evapotranspiration has shown a significant upward trend over the past three decades, caused partly by the increasing atmospheric moisture demand [8].In particular, hydrological processes are very sensitive to climate change in arid and semi-arid areas.In the Middle East, acceleration of hydrological processes induced by climate change has caused more severe droughts and flooding events, affecting the region's well-being [9].
In addition to climate change, LULCC also alters hydrological processes.Reforestation/afforestation or deforestation changes surface evapotranspiration, canopy water interception, and soil water infiltration capacity, changing the hydrological processes within watersheds.Many previous studies have shown that reforestation results in a decrease in streamflow due to greater infiltration into the soil and higher precipitation interception by vegetation [10,11].Deforestation can reduce root density and depth, and lower leaf mass, resulting in decreased vegetation water consumption, weaker evapotranspiration, and higher streamflow [12,13].These changes within a watershed lead to a redistribution among the components of hydrological processes [14].
As mentioned, climate change and LULCC are two important factors that significantly affect hydrological processes at different temporospatial scales.Streamflow observations around the world have indicated varying levels of climate change and LULCC impact, particularly in basins located in arid and semi-arid climate zones [15,16].Modeling techniques have been adopted to evaluate the contributions of climate change and LULCC to streamflow changes.The Budyko equation is a commonly used and effective tool to address such contributions due to its simplicity and physical background [17,18].The Budyko equation, based on the water and energy balance at a watershed scale, demonstrates the physical distribution among precipitation, evapotranspiration, and streamflow at a long-term temporal scale [19].Since it was established, the Budyko equation has been widely used to answer water and energy balance questions throughout the world [20][21][22].
However, limitations still exist for applications of the Budyko equation, which assumes non-changing water storage in a watershed over an application period.This assumption is often very difficult to satisfy due to the lack of sufficient observations.Yang et al. [23] used the Budyko equation to derive the elasticity of streamflow in relation to climatic variables in China at an annual timescale.Jiang et al. [24] used a time length of 11 years to satisfy the non-changing water storage assumption without observed evidence.Donohue et al. [25] asserted that 30 years of data were required to meet the criterion of the Budyko non-changing water storage for their study watersheds.In addition, many studies assume that the physical properties of a watershed do not exhibit significant changes by setting the Budyko parameter that represents such properties to a constant [26,27].Nevertheless, vegetation as a key component in the watershed often changes significantly under climate change and/or through human activities.In this study, variable vegetation was introduced to the Budyko equation to improve understanding of the influence of LULCC on hydrological processes.Thus, we applied the Budyko equation to a watershed in the Loess Plateau, China, where vegetation cover has been significantly altered by climate and human activities.In Section 2, the study methods are described, Section 3 introduces the study area and data, Section 4 describes the results, and conclusions are given in Section 5.

Budyko Parameter Estimation
With the Budyko equation's assumption that changes in water storage in a watershed are negligible over a sufficiently long time, precipitation (P) is partitioned into evapotranspiration (E) and streamflow (R) for a watershed [19].The ratio of actual evapotranspiration to precipitation (θ = E/P, the evapotranspiration ratio) is controlled principally by the ratio of potential evapotranspiration to precipitation (ε = E p /P, the climatic dryness index) on a long-term timescale.For humid watersheds Water 2018, 10, 1781 3 of 17 (P > E p ), the actual evapotranspiration is controlled predominantly by the energy supply (E p ), while for non-humid watersheds (P < E p ), it is controlled mainly by the water supply (P), as shown in Figure 1.Different functional forms of the Budyko equation have been developed [28]; one of the most widely used forms, the Choudhury-Yang (CY hereafter) equation, was selected for this study [29].E p was estimated using the Penman-Monteith method [30], and P, E p , and R were used as inputs for the CY equation: where η is the Budyko parameter that represents the average state of watershed characteristics such as vegetation cover, soil properties, topography, etc. energy supply ( ), while for non-humid watersheds ( <  ), it is controlled mainly by the water supply (), as shown in Figure 1.Different functional forms of the Budyko equation have been developed [28]; one of the most widely used forms, the Choudhury-Yang (CY hereafter) equation, was selected for this study [29]. was estimated using the Penman-Monteith method [30], and ,  , and  were used as inputs for the CY equation: where  is the Budyko parameter that represents the average state of watershed characteristics such as vegetation cover, soil properties, topography, etc. Traditionally,  can be derived from climate and streamflow data [27], but  cannot be calculated for ungauged watersheds using such a method (e.g., lack of streamflow data).Thus, determining  for ungauged watersheds is a challenge.For this study, we propose a polynomial equation to calculate  using climate, soil, topography, vegetation, and other available data (e.g., remote sensing data) for ungauged watersheds (without streamflow measurements) as follows: where  represents explanatory variables defining LULCC caused by human activities,  represents explanatory variables defining climate change;  is a constant term, and  and  are the corresponding regression coefficients.Through the maximum likelihood estimation method,  and  are estimated, and  is then estimated.

Quantifying the Contributions of Different Factors to Streamflow Changes
The Budyko parameter  might change for a watershed, implying a change in the watershed's characteristics.To quantify the contribution of each factor to a change in a watershed's water-energy balance, we adopted the decomposition method [14,24], described in Figure 1.There are two assumed paths to change a watershed from Point A to Point B: (1) a move from A to C along the dashed line, and (2) a vertical move from C to B. The first (A to C) shows that the  value for the watershed does not change, implying that the watershed ecosystem automatically adapts itself to climate change.The second (C to B) indicates a change in , implying that external forcing alters the watershed's physical features such as vegetation.Such external forcing could stem from human influences and/or climate change, but in the original decomposition method this external forcing is wholly attributed to human Traditionally, η can be derived from climate and streamflow data [27], but η cannot be calculated for ungauged watersheds using such a method (e.g., lack of streamflow data).Thus, determining η for ungauged watersheds is a challenge.For this study, we propose a polynomial equation to calculate η using climate, soil, topography, vegetation, and other available data (e.g., remote sensing data) for ungauged watersheds (without streamflow measurements) as follows: where H represents explanatory variables defining LULCC caused by human activities, C represents explanatory variables defining climate change; β 0 is a constant term, and β h and β c are the corresponding regression coefficients.Through the maximum likelihood estimation method, β h and β c are estimated, and η is then estimated.

Quantifying the Contributions of Different Factors to Streamflow Changes
The Budyko parameter η might change for a watershed, implying a change in the watershed's characteristics.To quantify the contribution of each factor to a change in a watershed's water-energy balance, we adopted the decomposition method [14,24], described in Figure 1.There are two assumed paths to change a watershed from Point A to Point B: (1) a move from A to C along the dashed line, and (2) a vertical move from C to B. The first (A to C) shows that the η value for the watershed does not change, implying that the watershed ecosystem automatically adapts itself to climate change.The second (C to B) indicates a change in η, implying that external forcing alters the watershed's physical features such as vegetation.Such external forcing could stem from human influences and/or climate change, but in the original decomposition method this external forcing is wholly attributed to human activities, assuming that all the factors that influence η originate from human activities.In our study, the contribution represented by the second path is decomposed based on the polynomial equation (Equation ( 2)) in Section 2.1.

Calculation of Vegetation Fraction and Relative Infiltration Capacity
The accuracy of the Budyko equation can be improved if vegetation changes are included [31][32][33].
To study the effect of vegetation on the hydrological processes, the green vegetation fraction (F g ) was introduced in the Budyko equation.F g can be derived from the normalized difference vegetation index (NDVI) based on satellite data.In this study, a quadratic equation was adopted to calculate F g using NDVI [34]: where NDV I i is the NDVI value on a remote sensing map, NDV I s is for bare soil, and NDV I ∞ is for dense green vegetation.For this study, NDV I ∞ and NDV I s were set to 0.05 and 0.68, respectively, based on remotely sensed data and land use types [35,36].Besides vegetation, water infiltration into the soil also affects the production of streamflow.The infiltration rate is controlled by rainfall intensity and soil infiltration capacity.In this study, the relative infiltration capacity was used to describe the relationship between the soil and the parameter η.The relative infiltration capacity is defined as the ratio of the saturation hydraulic conductivity, K s , to the average rainfall intensity, i r , within a period of 24 h [37]; i r is the average value for rainy days, and K s is obtained from the soil type database for the Wuding River Watershed (WRW) [38].

Study Area
To control soil erosion and improve environmental conditions, many soil conservation measures have been applied in the Loess Plateau (Figure 2) since the 1960s, one of which is the Grain for Green Project (GFGP) [39].This project has remarkably increased the vegetation cover in the Loess Plateau through afforestation/reforestation [40].Furthermore, this water-limited, environmentally fragile area is vulnerable to climate change at different temporospatial scales [41].For this study, we selected a typical watershed in the Loess Plateau, the WRW (Figure 2), to explore how afforestation/reforestation due to the GFGP affects hydrological processes under climate change.
Covering an area of approximately 30,261 km 2 , the WRW, located at 37.04 • -39.03 • N and 108.04 • -110.57• E, is in the center of the Loess Plateau.The Wuding River is a first-order tributary of the Yellow River.Streamflow data for this study were obtained from the Baijiachuan gauging station, which is located 100 km from the outlet of the WRW and has a drainage area accounting for 98% of the WRW.The WRW is in a semi-arid temperate continental climate zone, with average annual precipitation of 405 mm, a mean annual temperature of 8.0 • C, and potential evapotranspiration of 1007 mm, based on observational data over 1960-2011 (http://data.cma.cn/).Affected by the East Asian monsoon, approximately 75% of the annual rainfall occurs between June and September and is characterized by a significant number of heavy rain events.The topography is a typical loess hilly/gullied landscape with elevation ranging from 579 m to 1824 m.
have been applied in the Loess Plateau (Figure 2) since the 1960s, one of which is the Grain for Green Project (GFGP) [39].This project has remarkably increased the vegetation cover in the Loess Plateau through afforestation/reforestation [40].Furthermore, this water-limited, environmentally fragile area is vulnerable to climate change at different temporospatial scales [41].For this study, we selected a typical watershed in the Loess Plateau, the WRW (Figure 2), to explore how afforestation/reforestation due to the GFGP affects hydrological processes under climate change.

Hydrometeorological Data
Monthly streamflow data from gauge stations located in the main stream and first-order tributaries in the WRW were obtained from the Yellow River Hydrological Bureau.Only data covering at least 50 years were used in this study; data from eight stations met this criterion.Thus, all streamflow data used in this study covered the period from 1960 to 2011.Daily meteorological data from 12 stations in and around the WRW were obtained from the National Meteorological Information Center, China Meteorological Administration (http://data.cma.cn/), for the study period.These meteorological data include precipitation, maximum and minimum temperature, relative humidity, wind speed, sunshine duration, and solar radiation.We used the nonparametric Mann-Kendall (MK) test to detect the significance of temporal trends with a 95% confidence interval [42].

Digital Elevation Model (DEM) and Soil Data
A DEM dataset at 30-m resolution was provided by the Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn).A soil dataset at 1-km resolution, containing soil property data and the spatial distribution of each soil type in the WRW, was provided by the Ecological Environment Database of the Loess Plateau (http://www.loess.csdb.cn/pdmp/index.action).The saturation hydraulic conductivity was verified with site observations from the WRW.

Satellite Remote Sensing Data
As one of the most useful indices for vegetation monitoring in terrestrial ecosystems, NDVI derived from remote sensing data was used.This study selected the Global Inventory Modeling and Mapping Studies Normalized Difference Vegetation Index 3rd generation dataset (NDVI3g) for the WRW [43].The NDVI3g covers the period from 1982 to 2011 at a 0.083 • spatial resolution and a semi-monthly time step.NDVI3g data have been examined and compared with other NDVI products [44] and were found to be consistent with these data.The maximum value composite method was used to obtain the monthly and annual NDVI values [45].Therefore, this dataset was used to analyze the long-term vegetation trends and the relationship between vegetation and climate variability.

Temporal Trends of Streamflow
Figure 3a shows the changes in annual streamflow in the WRW from 1960 to 2011.The annual streamflow in the WRW experienced a significant decrease over this 52-year period.The observed downward trend of 5.2 mm/decade passes the 95% significance level using the MK test.Moreover, the annual streamflow shows two significant abrupt points in 1972 and 1998, which were detected using the nonparametric multiple change-points detection method [46].These abrupt points divide the study period into three stages, i.e., 1960-1972, 1973-1998, and 1999-2011, defined as Stages 1, 2, and 3, respectively.Figure 3a also shows that the amplitude of streamflow variation over the study period becomes weaker with time.
Water 2018, 11, 13 FOR PEER REVIEW 6 of 16 downward trend of 5.2 mm/decade passes the 95% significance level using the MK test.Moreover, the annual streamflow shows two significant abrupt points in 1972 and 1998, which were detected using the nonparametric multiple change-points detection method [46].These abrupt points divide the study period into three stages, i.e., 1960-1972, 1973-1998, and 1999-2011, defined as Stages 1, 2, and 3, respectively.Figure 3a also shows that the amplitude of streamflow variation over the study period becomes weaker with time.These three stages are consistent with water and soil conservation activities in the WRW, according to a survey of the WRW [47] (p.385), which shows that soil and water conservation activities over the WRW involve approximately three stages.The first stage spans from the 1950s to early 1970s, during which small-scale experimental field tests were performed to explore suitable ways of controlling soil erosion.The second stage lasts from the mid-1970s to the end of the 1990s, when the WRW was used as a national water and soil erosion management area.The last stage begins in 1999, when the WRW was one of the first GFGP pilot areas and more intensive conservation was WRW involve approximately three stages.The first stage spans from the 1950s to early 1970s, during which small-scale experimental field tests were performed to explore suitable ways of controlling soil erosion.The second stage lasts from the mid-1970s to the end of the 1990s, when the WRW was used as a national water and soil erosion management area.The last stage begins in 1999, when the WRW was one of the first GFGP pilot areas and more intensive conservation was performed.Watershed management records prove the validity of the abrupt statistical tests employed; thus, streamflow changes are closely related to human activities in the watershed.

Temporal Trends of Precipitation and Temperature
Climate change is one of the main factors affecting hydrological processes in the WRW.In Figure 3b, the annual precipitation shows a downward trend of 10 mm/decade, but this trend does not reach the 95% statistical significance level.A comparison of Figure 3a,b demonstrates that streamflow variability is controlled partly by precipitation changes.The correlation coefficients between precipitation and streamflow over the three stages are 0.8, 0.5, and 0.4, respectively, implying that the response by streamflow to precipitation becomes weaker.There must be other factors causing the decline in streamflow.
Figure 3c shows the time series of area-averaged annual temperature for the WRW.An upward trend of 0.27 • C/decade at the 95% significance level was detected by the MK test.This remarkable trend is five times the global average temperature change [48].A rising temperature could contribute to the reduced streamflow in this area by increasing evapotranspiration [49], as will be discussed again in Section 4.5.

Determination of Timescale in the Budyko Equation
In the Budyko equation, water storage change in a watershed is assumed to be zero or very close to zero over a long-term period.However, the length of this period is watershed dependent, and it is impossible to accurately measure water storage change in almost any watershed.In some studies, researchers have arbitrarily set water storage change to be zero over a period ranging from 1 to 30 years with no support from observed evidence [23][24][25].For this study, we made a series of sensitivity tests to determine the timescale at which the water storage change is reasonably close to zero in the WRW.In these tests, we calculated the Budyko parameter η on timescales of 1 to 52 years with increments of one year.For each of the 52 tests, the water storage change was set to zero.We found that the Budyko parameter η stabilized on timescales longer than 13 years, although there was a slight upward trend between timescales of 13 and 52 years (Figure 4).Therefore, we derived the value of η on a timescale of 13 years, at which water storage change can be reasonably assumed to be zero.the 95% statistical significance level.A comparison of Figure 3a,b demonstrates that streamflow variability is controlled partly by precipitation changes.The correlation coefficients between precipitation and streamflow over the three stages are 0.8, 0.5, and 0.4, respectively, implying that the response by streamflow to precipitation becomes weaker.There must be other factors causing the decline in streamflow.
Figure 3c shows the time series of area-averaged annual temperature for the WRW.An upward trend of 0.27 °C/decade at the 95% significance level was detected by the MK test.This remarkable trend is five times the global average temperature change [48].A rising temperature could contribute to the reduced streamflow in this area by increasing evapotranspiration [49], as will be discussed again in Section 4.5.

Determination of Timescale in the Budyko Equation
In the Budyko equation, water storage change in a watershed is assumed to be zero or very close to zero over a long-term period.However, the length of this period is watershed dependent, and it is impossible to accurately measure water storage change in almost any watershed.In some studies, researchers have arbitrarily set water storage change to be zero over a period ranging from 1 to 30 years with no support from observed evidence [23][24][25].For this study, we made a series of sensitivity tests to determine the timescale at which the water storage change is reasonably close to zero in the WRW.In these tests, we calculated the Budyko parameter  on timescales of 1 to 52 years with increments of one year.For each of the 52 tests, the water storage change was set to zero.We found that the Budyko parameter  stabilized on timescales longer than 13 years, although there was a slight upward trend between timescales of 13 and 52 years (Figure 4).Therefore, we derived the value of  on a timescale of 13 years, at which water storage change can be reasonably assumed to be zero.

Temporospatial Changes in Vegetation
NDVI is an effective parameter representing vegetation cover in the WRW. Figure 5 shows that the area-averaged annual NDVI for the WRW increased from 1982 to 2011, indicating a growth in vegetation over this period.There was a pronounced change around 2000, which divided the period into two stages.These two stages fall within Stages 2 and 3, characterized by significant water

Temporospatial Changes in Vegetation
NDVI is an effective parameter representing vegetation cover in the WRW. Figure 5 shows that the area-averaged annual NDVI for the WRW increased from 1982 to 2011, indicating a growth in vegetation over this period.There was a pronounced change around 2000, which divided the period into two stages.These two stages fall within Stages 2 and 3, characterized by significant water conservation activities in the watershed.The NDVI trends for these two stages are 3.6 × 10 -3 /yr and 11.8 × 10 -3 /yr, respectively.The significant increase in the latter stage indicates remarkable vegetation growth in the WRW associated with the GFGP since 1999 [41].
Water 2018, 11, 13 FOR PEER REVIEW 8 of 16 and their 95% significance levels (black dots) for the same two stages, where 29% of the WRW in 1982-1998 and 83% of the WRW in 1999-2011 pass the 95% significance level.Particularly in the second stage, the middle and lower reaches of the WRW have the most significant NDVI increases, where the most severe soil erosion often occurs, and thus where reforestation/afforestation has been focused.In addition, pixels that did not pass the 95% significance level are predominantly urban areas.Based on the above analysis, the WRW has experienced remarkable vegetation growth, particularly from 1999 to 2011, due to reforestation/afforestation.Such a substantial landscape change goes against the rules of the Budyko equation application, which assumes minimal landscape changes in a watershed.In this study, we made a significant effort to include landscape changes in the Budyko equation, with a focus on vegetation.The average NDVI spatial distributions over the two stages are shown in Figure 6a,b.Generally, the NDVI increases from southeast to northwest in the watershed during both stages, consistent with a change from a humid climate to a semi-arid one.The vegetation cover increases quite dramatically from the 1982-1998 to 1999-2011 periods.Figure 6c,d show the trends of NDVI in the WRW (pixels) and their 95% significance levels (black dots) for the same two stages, where 29% of the WRW in 1982-1998 and 83% of the WRW in 1999-2011 pass the 95% significance level.Particularly in the second stage, the middle and lower reaches of the WRW have the most significant NDVI increases, where the most severe soil erosion often occurs, and thus where reforestation/afforestation has been focused.In addition, pixels that did not pass the 95% significance level are predominantly urban areas.
Based on the above analysis, the WRW has experienced remarkable vegetation growth, particularly from 1999 to 2011, due to reforestation/afforestation.Such a substantial landscape change goes against the rules of the Budyko equation application, which assumes minimal landscape changes in a watershed.In this study, we made a significant effort to include landscape changes in the Budyko equation, with a focus on vegetation.
Based on the above analysis, the WRW has experienced remarkable vegetation growth, particularly from 1999 to 2011, due to reforestation/afforestation.Such a substantial landscape change goes against the rules of the Budyko equation application, which assumes minimal landscape changes in a watershed.In this study, we made a significant effort to include landscape changes in the Budyko equation, with a focus on vegetation.

Estimation of Parameter 𝜂 in the Budyko Equation
By considering the vegetation changes in the WRW, we used covariate analysis with the Akaike information criterion [50] to develop an empirical scheme to estimate  .In this scheme, we parameterized  as a function of explanatory variables including vegetation cover,  , the relative infiltration capacity, irrigation area, and terrace area.The Budyko parameter  was optimized using the above method to quantify the relationship between  and the explanatory variables.Finally,  was estimated as follows: where  reflects the vegetation conditions as one of the most important landscape factors in a watershed and is derived from NDVI through the conversion model discussed above.The relative infiltration capacity denotes the infiltration property that influences streamflow generation.The multiple R-squared of the regression equation is 0.86 and passes the 95% significance level, indicating that these factors can realistically explain .These factors represent vegetation, soil, and climate conditions, in which vegetation changes are induced mainly by human activities.The result reveals a significant positive correlation (0.76) between  and  in the WRW. Figure 7 illustrates the  estimated with Equation (4) versus the  calculated based on the Budyko equation with the observed input variables.For the WRW, the  value generated with the above regressed polynomial equation agrees very well with that derived from the Budyko equation.
By inputting this estimated  into the Budyko equation, we calculated the streamflow for the WRW.As shown in Figure 8, the root mean square error and Nash Sutcliffe efficiency coefficient are 1.22 mm and 0.91, respectively.The streamflow results calculated by a constant  are also displayed

Estimation of Parameter η in the Budyko Equation
By considering the vegetation changes in the WRW, we used covariate analysis with the Akaike information criterion [50] to develop an empirical scheme to estimate η.In this scheme, we parameterized η as a function of explanatory variables including vegetation cover, F g , the relative infiltration capacity, irrigation area, and terrace area.The Budyko parameter η was optimized using the above method to quantify the relationship between η and the explanatory variables.Finally, η was estimated as follows: where F g reflects the vegetation conditions as one of the most important landscape factors in a watershed and is derived from NDVI through the conversion model discussed above.The relative infiltration capacity denotes the infiltration property that influences streamflow generation.The multiple R-squared of the regression equation is 0.86 and passes the 95% significance level, indicating that these factors can realistically explain η.These factors represent vegetation, soil, and climate conditions, in which vegetation changes are induced mainly by human activities.The result reveals a significant positive correlation (0.76) between F g and η in the WRW. Figure 7 illustrates the η estimated with Equation (4) versus the η calculated based on the Budyko equation with the observed input variables.For the WRW, the η value generated with the above regressed polynomial equation agrees very well with that derived from the Budyko equation.
Water 2018, 11, 13 FOR PEER REVIEW 10 of 16 also substantially demonstrates the streamflow trend.This case is useful for situations where vegetation data are insufficient, especially on the large timescale of future climate scenarios.

Contributions of Climate Change and Vegetation to Streamflow
To quantify the contributions of different factors to streamflow changes, the improved decomposition method mentioned in Section 2 was applied.In view of the good performance of explanatory variables at interpreting the Budyko parameter , Equation ( 4) was used to calculate the change in mean annual streamflow in each 13-year period, together with the Budyko equation (Equation ( 1)).Therefore, the streamflow changes in each period are compared with the baseline period 1970-1982, which is the first 13-year period containing vegetation information.The baseline period is denoted as the pre-stage, and other lengths are denoted as the post-stage.The calculated By inputting this estimated η into the Budyko equation, we calculated the streamflow for the WRW.As shown in Figure 8, the root mean square error and Nash Sutcliffe efficiency coefficient are 1.22 mm and 0.91, respectively.The streamflow results calculated by a constant η are also displayed in Figure 8, and the root mean square error and Nash Sutcliffe efficiency coefficient are 2.95 mm and 0.49, respectively.The constant η was derived on a timescale of the entire period, indicating no watershed landscape change over the study period.A comparison between these two calculated results indicates that the Budyko equation is more accurate when changes in landscape factors, especially vegetation, are included.The streamflow calculated by the η that considers vegetation changes reflects not only the streamflow trend but also the magnitude.Conversely, the streamflow calculated by a constant η greatly underestimates streamflow during the first several years and overestimates streamflow in the last several years of the period.This implies that a constant η cannot reflect dynamic changes in watershed landscape characteristics.However, it is worth noting that the constant η case also substantially demonstrates the streamflow trend.This case is useful for situations where vegetation data are insufficient, especially on the large timescale of future climate scenarios.

Contributions of Climate Change and Vegetation to Streamflow
To quantify the contributions of different factors to streamflow changes, the improved decomposition method mentioned in Section 2 was applied.In view of the good performance of explanatory variables at interpreting the Budyko parameter η, Equation (4) was used to calculate the change in mean annual streamflow in each 13-year period, together with the Budyko equation (Equation ( 1)).Therefore, the streamflow changes in each period are compared with the baseline period 1970-1982, which is the first 13-year period containing vegetation information.The baseline period is denoted as the pre-stage, and other lengths are denoted as the post-stage.The calculated result of the decomposition method is shown in Figure 9, indicating that a combination of climate and human activities (mainly from vegetation changes) led to the streamflow decline in recent years.From the average contributions of climate and vegetation during different periods (Figure 9a), the conclusion can be made that climate change is the dominant factor affecting streamflow, accounting for nearly 76% of the total streamflow reduction.Vegetation changes are also important factors, accounting for about 24% of the streamflow decrease.Further, the streamflow reduction induced by climate increased substantially after 1999 (Figure 9a), which is attributed to the increasingly dry climate.The relationship between drought and streamflow change is discussed in the following.

Contributions of Climate Change and Vegetation to Streamflow
To quantify the contributions of different factors to streamflow changes, the improved decomposition method mentioned in Section 2 was applied.In view of the good performance of explanatory variables at interpreting the Budyko parameter , Equation (4) was used to calculate the change in mean annual streamflow in each 13-year period, together with the Budyko equation (Equation ( 1)).Therefore, the streamflow changes in each period are compared with the baseline period 1970-1982, which is the first 13-year period containing vegetation information.The baseline period is denoted as the pre-stage, and other lengths are denoted as the post-stage.The calculated from the change in average rainfall intensity, which influences the relative infiltration capacity.This is a crucial factor to consider in landscape characteristics, because infiltration excess overland flow is the main mechanism for streamflow generation in a typical loess soil watershed [51].In order to distinguish the impact of climate and vegetation changes on η, the streamflow reduction induced by these two factors via altering η is compared in Figure 9b.In Figure 9b, streamflow reduction caused by climate change remains steady with little variation and is smaller than that caused by vegetation changes.This indicates that changes in η induced by climate change are not negligible, which has not previously been considered [52].Figure 9b also indicates that vegetation is the primary factor affecting η; streamflow reduction induced by vegetation changes represents the majority of streamflow reduction caused by altering parameter η.This implies that vegetation is vital to the hydrology in this semi-arid watershed, and growth in vegetation cover increases the evapotranspiration ratio and reduces the streamflow ratio to precipitation.It also demonstrates the significance of introducing a vegetation factor into streamflow estimates in the Budyko equation.
From the average contributions of climate and vegetation during different periods (Figure 9a), the conclusion can be made that climate change is the dominant factor affecting streamflow, accounting for nearly 76% of the total streamflow reduction.Vegetation changes are also important factors, accounting for about 24% of the streamflow decrease.Further, the streamflow reduction induced by climate increased substantially after 1999 (Figure 9a), which is attributed to the increasingly dry climate.The relationship between drought and streamflow change is discussed in the following.In order to test the validity and rationality of the improved decomposition method, the contribution was also quantified using another mainstream method called the elasticity method [27,53].The elasticity method Water 2018, 10, 1781 13 of 17 results are not shown here, but our results with this method are similar to those with the improved decomposition method.
There are 74 dams with a storage capacity greater than one million cubic meters in the WRW with the purpose of flooding control [54] (p.705).Nevertheless, these dams mainly affect the seasonal variations of streamflow in the WRW, and they do not have a significant influence on the volume of annual streamflow.In addition, most of these dams lost normal function in the end of 1980s due to the sediment deposition caused by severe soil erosion [55] (pp.428-429).Thus, our study did not include the influence of dam regulation on the streamflow, and focused on the change in annual streamflow in WRW over the period of 1982 to 2011.
Precipitation is the only source of water input to a closed watershed and is partitioned into different parts, such as soil water storage and evapotranspiration.The results with our improved Budyko equation application indicate that hydrological processes are the result of the long-term co-evolution of a watershed's vegetation and climate [14].The contribution analysis results of the WRW demonstrate the dominant role of climate in this complex evolved system.These findings were further confirmed with the Palmer Drought Severity Index (PDSI) [56,57]

Conclusions
In this study, we diagnosed hydrometeorological changes in the WRW, with a focus on vegetation cover changes.Over recent years, streamflow has dramatically declined, regional climate change has become evident, and the watershed has experienced more severe drought.Vegetation cover changes are the main reason for underlying surface changes in the WRW.The timing of abrupt changes indicates that NDVI changes are closely tied to water and soil conservation activities in the WRW and streamflow changes.Intense variations of NDVI in such a short time reveal that human activities are the main driving force of vegetation cover changes.

Conclusions
In this study, we diagnosed hydrometeorological changes in the WRW, with a focus on vegetation cover changes.Over recent years, streamflow has dramatically declined, regional climate change has become evident, and the watershed has experienced more severe drought.Vegetation cover changes are the main reason for underlying surface changes in the WRW.The timing of abrupt changes indicates that NDVI changes are closely tied to water and soil conservation activities in the WRW and streamflow changes.Intense variations of NDVI in such a short time reveal that human activities are the main driving force of vegetation cover changes.
Using the moving average method with a timescale of 13 years, an optimized model was established, incorporating the Budyko parameter, vegetation cover, and relative infiltration capacity.The main factors that influence watershed landscape characteristics were then determined, i.e., climate change and vegetation changes.The good performance of the estimated streamflow implies that the Budyko parameter can be explained by these variables.Based on this optimized model, an improved decomposition method was used to separate the impact of climate change and vegetation cover changes on streamflow.It should be noted that we considered the climatic impact on the Budyko parameter η, which has previously been ignored.Furthermore, introducing the main factors that affect the Budyko parameter improved the performance of the Budyko equation by incorporating physical mechanisms.

Figure 1 .
Figure 1.Schematic of water-energy balance changes in a watershed as indicated by the Choudhury-Yang (CY) Budyko-type equation, and the decomposition method.

Figure 1 .
Figure 1.Schematic of water-energy balance changes in a watershed as indicated by the Choudhury-Yang (CY) Budyko-type equation, and the decomposition method.

Figure 2 .
Figure 2. Locations of the study area and hydrometeorological stations (asterisk shows the river outlet).

Figure 3 .
Figure 3. Annual changes in hydrometeorological variables in the Wuding River Watershed (WRW) from 1960 to 2011: (a) streamflow; (b) precipitation; and (c) temperature.These three stages are consistent with water and soil conservation activities in the WRW, according to a survey of the WRW[47] (p.385), which shows that soil and water conservation activities over the

Figure 4 .
Figure 4. Sensitivity tests determining the timescale used to calculate the Budyko parameter .

Figure 4 .
Figure 4. Sensitivity tests determining the timescale used to calculate the Budyko parameter η.

Figure 5 .
Figure 5. Variation in area-averaged annual normalized difference vegetation index (NDVI) from 1982 to 2011 in the WRW (dashed lines are the trends of periods before and after the Grain for Green Project (GFGP)).

Figure 5 .
Figure 5. Variation in area-averaged annual normalized difference vegetation index (NDVI) from 1982 to 2011 in the WRW (dashed lines are the trends of periods before and after the Grain for Green Project (GFGP)).

Water 2018, 11 , 16 Figure 6 .
Figure 6.Spatial distribution of average annual NDVI in the WRW: (a) average annual NDVI from 1982 to 1998; (b) average annual NDVI from 1999 to 2011; (c) linear regression slope of NDVI from 1982 to 1998; and (d) linear regression slope of NDVI from 1999 to 2011 (dots denote the slope at the 95% significance level).

Figure 6 .
Figure 6.Spatial distribution of average annual NDVI in the WRW: (a) average annual NDVI from 1982 to 1998; (b) average annual NDVI from 1999 to 2011; (c) linear regression slope of NDVI from 1982 to 1998; and (d) linear regression slope of NDVI from 1999 to 2011 (dots denote the slope at the 95% significance level).

Figure 7 .
Figure 7.Comparison of  calculated by observed input variables with  estimated by the regression equation (dashed line is the 1:1 line).

Figure 8 .
Figure 8.Comparison of two types of calculated streamflow with the observed streamflow.

Figure 7 .
Figure 7.Comparison of η calculated by observed input variables with η estimated by the regression equation (dashed line is the 1:1 line).

Figure 7 .
Figure 7.Comparison of  calculated by observed input variables with  estimated by the regression equation (dashed line is the 1:1 line).

Figure 8 .
Figure 8.Comparison of two types of calculated streamflow with the observed streamflow.

Figure 8 .
Figure 8.Comparison of two types of calculated streamflow with the observed streamflow.Climate and landscape changes in a watershed have an important effect on hydrological processes; this effect can be reflected in the Budyko equation by altering parameter η.The contribution from climate can be divided into two parts: the first part is caused by change in the meteorological input of precipitation and potential evapotranspiration to the WRW (direct climate change); the other part is caused by climate change through modification of the watershed landscape characteristics (indirect climate change).In this study, the impact of climate change on η (indirect climate change) originates from the change in average rainfall intensity, which influences the relative infiltration capacity.This is a crucial factor to consider in landscape characteristics, because infiltration excess overland flow is the main mechanism for streamflow generation in a typical loess soil watershed[51].In order to distinguish the impact of climate and vegetation changes on η, the streamflow reduction induced by these two factors via altering η is compared in Figure9b.In Figure9b, streamflow reduction caused by climate change remains steady with little variation and is smaller than that caused by vegetation changes.This indicates that changes in η induced by climate change are not negligible, which has not previously been considered[52].Figure9balso indicates that vegetation is the primary factor affecting η; streamflow reduction induced by vegetation changes represents the majority of streamflow reduction caused by altering parameter η.This implies that vegetation is vital to the hydrology in this semi-arid watershed, and growth in vegetation cover increases the evapotranspiration ratio and reduces the streamflow ratio to precipitation.It also demonstrates the significance of introducing a vegetation factor into streamflow estimates in the Budyko equation.

Figure 9 .
Figure 9. Contributions separation results of the decomposition method: (a) comparison of the contributions of climate change and vegetation; (b) comparison of the contributions of climate change and vegetation by altering the Budyko parameter; and (c) comparison of the contributions of two types of climate change.

Figure 9 .
Figure 9. Contributions separation results of the decomposition method: (a) comparison of the contributions of climate change and vegetation; (b) comparison of the contributions of climate change and vegetation by altering the Budyko parameter; and (c) comparison of the contributions of two types of climate change.The two different components of climate change contribution to the streamflow are shown in Figure 9c.One represents the contribution directly induced by climate change (direct climate change), and the other is the contribution induced by altering η by climate change (indirect climate change).The direct climate change contribution accounts for the majority (88%) of the total climate change contribution, and the indirect climate change contribution accounts for 12%.In order to test the validity and rationality of the improved decomposition method, the contribution was also quantified using another mainstream method called the elasticity method[27,53].The elasticity method , a physically based hydrometeorological index.The calculation of PDSI does not consider interference from human activities in this watershed, and thus this index explains hydrological drought patterns regardless of human influences [58].The annual changes in streamflow and PDSI in the WRW from 1960 to 2011 are shown in Figure 10.These two variables derived from independent datasets exhibit similar trends and variations.The downward trend of streamflow is −0.048 mm/yr and that of PDSI is −0.047.The MK test results indicate that both show a significant downward trend at the 95% significance level.The decreasing PDSI indicates that the WRW has experienced increasingly serious droughts since the early 1980s.Moreover, this similarity shows that PDSI captures the trend of streamflow change and the dominant role of climate in streamflow reduction.However, the performance of PDSI deteriorates in Stage 2 and Stage 3 compared to Stage 1, which demonstrates that human activities play a non-negligible role in streamflow reduction in the WRW.Water 2018, 11, 13 FOR PEER REVIEW 13 of 16 in Stage 2 and Stage 3 compared to Stage 1, which demonstrates that human activities play a nonnegligible role in streamflow reduction in the WRW.

Figure 10 .
Figure 10.Annual streamflow and Palmer Drought Severity Index (PDSI) changes in the WRW from 1960 to 2011.

Figure 10 .
Figure 10.Annual streamflow and Palmer Drought Severity Index (PDSI) changes in the WRW from 1960 to 2011.