Vegetation dynamics and climate seasonality jointly control the interannual 1 catchment water balance in the Loess Plateau under the Budyko framework

Within the Budyko framework, the controlling parameter (ω in the Fu equation) is widely 11 considered to represent landscape conditions in terms of vegetation coverage (M); however, some 12 qualitative studies have concluded that climate seasonality (S) should be incorporated in ω. Here, we 13 discuss the relationship between ω, M, and S, and further develop an empirical equation so that the 14 contributions from M to actual evapotranspiration (ET) can be determined more accurately. Taking 13 15 catchments in the Loess Plateau as examples, ω was found to be well correlated with M and S. The 16 developed empirical formula for ω calculations at the annual scale performed well for estimating ET by 17 the cross-validation approach. By combining the Budyko framework with the semi-empirical formula, 18 the contributions of changes in ω to ET variations were further decomposed as those of M and S. 19 Results showed that the contributions of S to ET changes ranged from 0.1% to 65.5 % (absolute values); 20 therefore, the impacts of climate seasonality on ET cannot be ignored. Otherwise, the contribution of M 21 to ET changes will be estimated with a large error. The developed empirical formula between ω, M, and 22 S provides an effective method to separate the contributions of M and S to ET changes. 23


Introduction
The water cycle has been influenced greatly by human activities and climate change since the 1960s, and considerable variability in hydrological processes has been observed in many basins around the world; this has led to a series of problems concerning essential water resources (Stocker et al., 2014).Analyses of the mechanisms of the interactions among the water balance, climate, and catchment surface conditions are important for understanding these complex processes at different spatio-temporal scales (Zhang et al., 2008), and such work has practical significance in regard to the improvement of water resources and land management (Rodriguez-Iturbe, 2000;Xu et al., 2014).Budyko (1948Budyko ( , 1974) ) postulated that precipitation (P, represents the water supply from the atmosphere) and potential evapotranspiration (ET 0 , represents the demand by the atmosphere) are the two dominant variables that control the longterm average water balance.The Budyko framework is considered one of the most abiding frameworks linking climatic conditions to the runoff (R) and actual evapotranspiration (ET) of a catchment (Donohue et al., 2007), and it has been used successfully to investigate interactions between hydrological processes, climate variability, and landscape characteristics (e.g.Milly, 1994;Woods, 2003;Yokoo et al., 2008;Yang et al., 2009;Liu et al., 2012).A series of empirical formulas have been developed for the Budyko curve based on theoretical research and case studies of regional water bal-T.Ning et al.: Interannual catchment water balance ance over the past 50 years.Among them, the Fu (Fu, 1981;Zhang et al., 2004) and Choudhury-Yang equations (Choudhury, 1999;Yang et al., 2008) have been used widely; furthermore, the controlling parameters ω (in the Fu equation) and n (in the Choudhury-Yang equation) are related linearly (Yang et al., 2008).
Deviations from the Budyko curve have been detected in previous studies, which indicates that in addition to climate conditions, other variables can also influence the variability of regional water balances (Yang et al., 2007;Wang and Alimohammadi, 2012).Two kinds of factors have been identified to be responsible for the deviations.The first type of factor is related to land surface conditions, and includes vegetation dynamics (Donohue et al., 2007(Donohue et al., , 2010;;Yang et al., 2009;Li et al., 2013;Zhang et al., 2016), soil properties, and topography (Yang et al., 2007;Peel et al., 2010).The second type of factor includes seasonal climate variability (in addition to P and ET 0 ), such as storm depth (Donohue et al., 2012;Shao et al., 2012;Li et al., 2014), frequency of daily rainfall (Milly, 1994), and differences in the timing of P and ET 0 (Budyko, 1961;Potter et al., 2005).All of these factors can be encoded into the controlling parameter of the Budyko equations (e.g.ω in the Fu equation and n in the Choudhury-Yang equation).So far, a great deal of attention has been paid to the relationships between land surface conditions and the controlling parameter.Based on satellite products of vegetation such as the normalized difference vegetation index (NDVI), vegetation has been found to correlate well with the controlling parameter, and some empirical relationships have been successfully developed (Yang et al., 2009;Li et al., 2013).In particular, the controlling parameter can be better represented by vegetation when higher spatiotemporal resolution products are used.Therefore, the impacts of dynamic changes in vegetation on hydrology can be effectively quantified.
Many current studies attribute any effects of the controlling parameter to landscape characteristics (Roderick and Farquhar, 2011;Zhou et al., 2015;Zhang et al., 2016).However, both empirical evidence and modelling tests have demonstrated the important function of climate seasonality on catchment water yield, and thereby, evidence exists that climate seasonality also strongly affects the controlling parameter in the Budyko equations (Berghuijs and Woods, 2016).Some indices and models have thus been developed to address this issue, and several potential solutions have been discussed (Milly, 1993(Milly, , 1994;;Potter et al., 2005;Yokoo et al., 2008;Feng et al., 2012;Li, 2014).Yang et al. (2012) introduced the climate seasonality index into the Budyko framework and proposed an empirical equation to include its effect in the estimation of the long-term controlling parameters; however, by focusing on the mean annual scale, the effects of vegetation dynamics were not considered.Therefore, how the vegetation dynamics and climate seasonality jointly control the interannual variability in the controlling parameters needs further interpretation.
Therefore, the primary motivation behind this study was to detect the potential linkages between the controlling parameter and surface condition change, as well as climate seasonality at the annual scale.The specific objectives were to derive an appropriate analytic formula between parameter ω in the Fu equation and the above two factors for typical catchments in the Loess Plateau, and then quantify the impacts of vegetation change and climate seasonality variability on the catchment water balance.

Annual water balance definition
The Budyko framework assumes that the long-term average water balance is in a steady state (Wang and Alimohammadi, 2012), and the water storage change in a catchment can be negligible.The interannual variability of the water balance in individual basins can also be studied by overlooking the interannual variation of the catchment water storage (Sankarasubramanian and Vogel, 2002;Yang et al., 2007;Potter and Zhang, 2009).However, water storage change can be great when analysing the interannual variability of the water balance (Wang, 2012).To minimize the potential errors introduced by neglecting water storage variation, the hydrological year (Sivapalan et al., 2011;Carmona et al., 2014) and moving windows (Jiang et al., 2015) were introduced to the time series of annual hydrological variables.Similar to Sivapalan et al. (2011) andCarmona et al. (2014), the hydrological year rather than the calendar year is introduced to calculate the annual ET in our study, and this is called the "measured" ET in the subsequent discussion.Specifically, as the study area has a continental monsoon climate with most rainfall occurring in summer and autumn (July-September), a hydrological year is defined as July to June of the following year.In this way, the water input occurs mainly at the beginning of the year and the water is consumed within that year.

Identification of factors determining parameter ω in Fu's equation
The Fu equation is used in this study with the following expressions: where ω is the controlling parameter of the Budyko curve.ET 0 is calculated by using the equation of Priestley and Taylor (1972).
The important issue regarding the parameterization of ω in Fu's equation is to choose factors with physical meanings.According to the results from related studies, land surface conditions can be mainly represented by vegetation, which was also true in this study.Although soil properties and topography also influence the variation of catchment water balance, they are relative stable for a catchment.Therefore, vegetation dynamics (i.e.vegetation coverage) were chosen to represent the variations in surface conditions and the effects of soil and topography were not considered in our study.The vegetation coverage (M) was estimated by the following equation (Yang et al., 2009): where NDVI max and NDVI min are the NDVI values of dense forest (0.80) and bare soil (0.05) in the previous studies (Yang et al., 2009;Li et al., 2013), respectively.The maximum NDVI in our study is 0.85, thus the value of NDVI max in Eq. ( 2) was defined as 0.85.The effects of seasonal variations in coupled water and energy are important for better understanding the water cycle for different climatic regions.For example, two catchments with the same annual P and ET 0 may have different evapotranspiration partitioning or runoff if the balance between P and ET 0 within a year are different, Solar radiation was considered as the dominant factor that controls the climate seasonality and thus the seasonality of ET 0 can be expressed by sine function.Similarly, the seasonality of P also exhibited the characteristic of sine distribution (Milly, 1994;Woods, 2003;Yang et al., 2012): where δ P and δ ET 0 are the ratios of the amplitudes of the monthly harmonics to the monthly averages of precipitation and potential evapotranspiration, respectively.Larger absolute values of δ P and δ ET 0 mean larger variability of climate seasonality.ϑ is the duration of the seasonal cycle, 2π/ϑ equal to 1 year.Woods (2003) summarized the modelled climate of Eq. ( 3a) and (b) in dimensionless form and defined the climate seasonality index (S) and here it was used to reflect the non-uniformity in the intra-annual distribution of water and energy in our study: where ∅ is the dryness index, ∅ = ET 0 /P .If S = 0, there is no seasonal fluctuation of the difference between P and ET 0 .Larger values of S indicate that the larger changes in the balance between P and ET 0 during the seasonal cycle.

Evaluating the contributions of climate change and surface condition alterations to ET changes
Based on the climate elasticity method, which was introduced by Schaake (1990) and improved by Sankarasubramanian et al. (2001), the contribution of change for each factor to runoff was defined as the product of the sensitivity coefficient and the variation of the factor (Roderick and Farquhar, 2011): (5) However, due to ignoring the higher orders of the Taylor expansion in Eq. ( 5), this method will result in high errors (Yang et al., 2014b).Recently, Zhou et al. (2016) proposed a new method to partition climate and catchment effect on the mean annual runoff based on the Budyko complementary relationship, called "the complementary method".The algebraic identities in their work can ensure that the change in runoff can be decomposed into two components precisely without any residuals.Here, we extend "the complementary method" to conduct attribution analysis of ET changes for each basin by further incorporating the effects of vegetation coverage and climate seasonality: where α is a weighting factor that varies from 0 to 1, which can determine the upper and lower bounds of the climate and the controlling parameter effect.In this study, we defined α = 0.5 according to the recommendation of Zhou et al. (2016).The difference operator ( ) refers to the difference of a variable from period I (1981 to the changing point detected by Pettitt's test; Pettitt, 1979) to period II (period-I end to 2012), e.g.ET 0 = ET 0,2 − ET 0,1 .Then the contributions of P, ET 0 , and ω changes to the ET change can be expressed as follows: T. Ning et al.: Interannual catchment water balance After obtaining the contribution of parameter ω to the ET change, the contributions of vegetation coverage (M) and climate seasonality (S) changes to the ET change can be further decomposed as follows.
First, the contributions of M and S changes to the variation of parameter ω are calculated by using the total differential method similar to Eq. ( 5) based on the relationship between ω and M as well as S we built: Ignoring the higher orders of the Taylor expansion, Eq. ( 8a) can be rewritten as Furthermore, the individual relative contribution (RC) of M and S to ω can be calculated.Then, the contributions of M and S changes to ET variation can be obtained as follows: 3 Study area and data The Loess Plateau, which is located in the upper and middle reaches of the Yellow River in China, experiences a sub-humid and semiarid continental monsoon climate (Ning et al., 2016).Frequent heavy summer storms, sparse vegetation coverage, easily erodible wind-deposited loess soil, and a long agricultural history have all contributed to severe drought and soil erosion problems in this region (Li et al., 2012).To recover and preserve the ecosystem, the Chinese government has launched numerous soil and conservation measures since the 1950s, and these include biologic measures ("Grain to Green" project) and engineering measures (building terraces and sediment trapping dams) (Mu et al., 2007).As a result, the hydrological processes of this area have undergone significant changes (Huang and Zhang, 2004;Zhang et al., 2008).Thirteen catchments on the Loess Plateau were selected as our study area (Fig. 1).Monthly runoff data for the 13 catchments were supplied by the Yellow River Conservancy Commission.Detailed information about the catchment characteristics and data durations are shown in Table 1.Daily meteorological data (1960Daily meteorological data ( -2012) ) comprised of precipitation, daily maximum and minimum temperatures, atmospheric pressures, wind speeds, mean relative humidity values, and sunshine durations, which were recorded at 96 stations, were provided by the China Meteorological Administration.The new NDVI third-generation (NDVI3g) dataset was used to represent the vegetation characteristics of the study area, and detailed information about this dataset was presented earlier by Fensholt and Proud (2012).The maximum value composite (MVC) procedure (Holben, 1986) was applied to produce the annual NDVI values.

The variability of parameter ω
The Budyko framework is usually used for analyses of long-term average catchment water balance; however, in this study, it was employed for the interpretation of the interannual variability of the water balance by using the hydrological year approach described earlier.To validate the feasibility of using Fu's equation for interannual variability, the evapotranspiration ratio (ET /P ) and dryness index (ET 0 /P ) on the annual scale for 13 basins are presented in the supporting information (Fig. S1 in the Supplement), and it can be seen that almost all points are focused on Fu's curves in each basin.Therefore, Fu's equation was considered adequate for the analysis of the interannual variability of the water balance.
If the controlling parameter ω on the annual scale can reflect the combined impacts of vegetation change and climate seasonality, it should also exhibit interannual variability with the seasonal variation in vegetation and climate, especially in those basins affected significantly by climate change and human activities.Obviously, this is true for basins in Loess Plateau (Fig. 2).During 1961-2012, ω values in all 13 basins had an upward trend.Along with such a changing trend in ω, ET should increase for the same levels of P and ET 0 .Before the 1980s, the variation in ω for each basin was relatively gentle; however, since that time, it has increased dramatically.

Development of the semi-empirical formula for parameter ω
The relationships between the annual parameter ω and vegetation coverage M as well as the climate seasonality index S were first explored in each study basin during the period 1981-2012, and the results are shown in Figs.S2 and S3.We can see that the parameter ω generally had a positive correlation with M, which implies that evapotranspiration increased with improvements in the vegetation conditions.However, ω was correlated negatively with S, which means that larger seasonal variations of coupled water and energy resulted in less evapotranspiration in this area.The relationships between ω and M as well as S imply that the annual variation in parameter ω can be estimated by the changes in vegetation dynamics and climate seasonality.To expand the sample size and span a wider range of climate conditions, as well as to make the derived semiempirical formula of parameter ω more representative, relationships were then developed based on the combined dataset from the 13 basins (Fig. 3).These results also indicate a good relationship between ω and M (R 2 = 0.43, p < 0.01) as well as S (R 2 = 0.24, p < 0.01).
To develop the semi-empirical formula of parameter ω, the limiting conditions of the two variables were considered as follows.
In Eq. ( 1), when ω → 1, ET → 0. After further considering the range of M and S, we can get that if ET→ 0, which indicates that T → 0 (transpiration), and thus M → 0, or implies that R → P , the match of P and ET 0 tends to be worst, i.e. S → ∞.Considering the relationships shown in Fig. 3 and given the above limiting conditions, the general form of parameter ω can be expressed as follows: where a, b, and c are constants.Using the least linear squares regression method, the semi-empirical formula of parameter ω is derived as follows: The coefficient of determination R 2 and the statistics for the F test of the modelled ω were 0.53 and 214.54, respectively.A cross-validation approach was chosen to calibrate and test the above semi-empirical formula for parameter ω.Specifically, the dataset for the 13 basins in our study was separated into two groups.One was applied to build the semiempirical formula, and it consisted of 12 basins for each time; the other was used for testing the performance of the semi-empirical formula, and it consisted of the remaining basin.In total, the cross-validation process was conducted 13 times.After building the semi-empirical formula by using the vegetation coverage and climate seasonality index data for the 12 basins, the parameter ω for the validated basin was modelled by using this fitted formula, and the annual ET for the validated basin was evaluated with the modelled ω, which is referred to as the "modelled" ET.Then, the "modelled" ET was compared with the "measured" ET.
Table 2 shows the cross-validation results for each basin.Except for basins 4 and 12, the MAE (mean absolute error) and RMSE (square root of the mean square error) values for each cross-validation process were relatively low, with mean values of 14.6 and 18.8 mm, respectively.The NSE coefficient (Nash-Sutcliffe coefficient of efficiency) for each process was greater than 0.85, thus suggesting that vegetation changes and climate seasonality can well explain the variation in the controlling parameter of the catchment water balance on the shorter timescale.

Quantitative attribution of the variation in ET
The impacts of vegetation changes on ET have been widely studied with the Budyko framework by assuming surface conditions can be represented by the controlling parameter.However, according to the developed relationships in our study, the controlling parameter is not only related to surface condition change, but also to climate seasonality.The contributions of changes in climate (P, ET 0 , and S) and vegetation (M) to the ET change were thus estimated by using the semi-empirical formula for parameter ω in the context of Fu's framework.
Trends in hydrometeorological variables and vegetation coverage were first analysed for each basin (Table 3).ET 0 , M, and S in almost all basins exhibited an upward trend, though with different significances.Based on the sensitivity coefficients of ET (Table S1 in the Supplement) and the changes in mean annual P, ET 0 , ω, M, and S from period I to period II (Table 4), the changes in ET due to those in P, ET 0 , M, and S were estimated using the method described in Sect.2.3.The contributions of four variables to ET change for each basin were presented in Table 4.In basins 1, 3-4, and 6, the ET changes were controlled by vegetation improvement; however, in the other basins, the dominant factor was precipitation.Except for basins 6, 9, and 12, elevated vegetation in most basins positively contributed to ET changes, which is consistent with Feng et al. (2016).ET in several basins showed a downward trend even though M positively contributed to ET changes; which is due to the offsetting effect of the other factors.
It should be noted that the climate seasonality (represented by S) played an important role in the catchment ET variation.The contributions of S to ET changes ranged from 0.1 to 74.8 % (absolute values).Besides basins 6, 9, and 12, the climate seasonality had a negative effect on ET variation in most of the basins, which means that larger seasonality differences between seasonal water and energy will lead to smaller amounts of evapotranspiration.Accordingly, if ω is supposed to only represent the landscape condition, the effects of landscape condition change on ET variation will be underestimated in basins 1, 3, 6-7, 9, and 11.Except for basin 9, the area of these basins is relative smaller; while its effects will be overestimated in the other basins, and the error would be equal to the contributions of S to ET changes.

Discussion
Although the controlling parameter ω showed a good relationship with the vegetation change and climate seasonality index, two groups of deviations around the regressed curves were detected (Fig. 3).The deviation points for the relationship between ω and M were mainly located at the top of the curve, i.e. corresponding to the same M values, where ω values were greater.We checked those points and found that precipitation and vegetation coverage in those years were normal, but runoff was very low compared to normal years.Excluding abrupt climate change, possible reasons for the extremely low runoff in those years include dam and reservoir operations, as well as irrigation diversions.A study conducted by Liang et al. (2015) on the same basins that we investigated in the Loess Plateau showed that check-dams increased continuously starting from the 1960s.By the year 2006, the numbers of dams along basins 10 and 5 reached 482 and 181, respectively.Dams can intercept stormwater runoff for a short period during flood seasons and allow more time for infiltration (Polyakov et al., 2014).A total of 21 large and 136 medium-sized reservoirs were installed along the Yellow River by 2001.Such infrastructure can also influence the runoff change by controlling the flooding, regulating the water discharge, and diverting the water to other regions (Chen et al., 2005).Agricultural production is heavily dependent on irrigation throughout the entire Yellow River basin, and it has been reported that water consumption by agricultural irrigation accounted for nearly 80.0 % of the entire water consumed from 1998 to 2011 (Wang et al., 2014).Thereby, water withdrawn for irrigation also plays an important role in the changing trends in runoff.In this study, the deviation points around the relationship curve between the annual ω and S fell in the upper left, and they were likely influenced by the low runoff.However, separation of the impacts on runoff from vegetation change, climate seasonality, and engineering works will have to await future work.
The relationships of parameter ω with vegetation dynamics and climate seasonality in some single basins were not significant in this study.Similarly, Yang et al. (2014a) also found a weak relationship between parameter n and vegetation coverage in 201 basins in China.This implies that the parameter ω might represent the combined effects of some other factors.For example, strong interactions among vegetation, climate, and soil conditions will lead to specific hydrologic partitioning at the catchment scale.In dry years, with low soil water contents, plants are trying to adapt by making use of hydrological processes, e.g.ground water dynamics and plant water storage mechanisms (Soylu et al., 2011).Therefore, the relationship between the parameter ω and vegetation dynamics can be influenced by climate and soil conditions.However, it is difficult to separate the climatic and soil components from the vegetation change.Moreover, Zhang et al. (2001) reported that the impact of different vegetation types on catchment water balance can be vastly different, and the plant-available water coefficient in their function, which is similar to parameter ω in Fu's equation, is related to vegetation type.Therefore, the vegetation type may also be an important variable that influences the parameter ω.
Despite that catchment-scale water storage changes are usually assumed to be zero on a long-term scale, the interannual variability of water storage can be an important component in annual water budget during dry or wet years (Wang and Alimohammadi, 2012), and cannot be ignored.However, the Loess Plateau has a sub-humid to semiarid climate, the water storage and its annual variation are relatively small compared with humid regions (see Fig. 5 from Mo et al., 2016).For example, using GRACE (Gravity Recovery and Climate Experiment), the water storage variations in the Yangtze, Yellow, and Zhujiang from 2003 to 2008 were analysed by Zhao et al. (2011), and the values for the Yangtze and Zhujiang basins were 37.8 and 65.2 mm, while no clear annual variations are observed in the Yellow River basin (3.0 mm).Furthermore, Mo et al. (2016) found that the water storage in Yellow River kept decreasing from 2004 to 2011, whereas it was changing slowly with a rate of 1.3 mm yr −1 .Therefore, considering the small water storage change in study area, ignoring water storage change in a period of hydrologic year is reasonable.
There were still errors in the attribution analysis of ET changes.As the changes in evapotranspiration were decomposed without residual by the complementary method (Eqs.6-7), the errors were induced from the developed empirical formula for ω (Eq.11).This suggested that ω cannot be completely explained by M and S, and it might include some other factors.Therefore, further discussion of factors influencing ω remains a goal for future work.

Conclusions
This study explored the concomitant effects of vegetation dynamics and climate seasonality on the variation in interannual controlling parameter ω from Fu's equation in the Loess Plateau.First, to reduce the impact of ignoring the water storage change on annual catchment water balance, the hydrological year approach was introduced to examine the interannual variability of the controlling parameter ω for the 13 basins in the Loess Plateau from 1961 to 2012.The findings showed that parameter ω in all these basins presented an increasing trend, especially after the 1980s.Furthermore, we checked the relationship between ω and vegetation dynamics (represented by the annual vegetation coverage, M) as well as climate seasonality (represented by the climate seasonality index, S).The interannual changes of parameter ω were found to be related strongly to M and S. As such, a semi-empirical formula for the annual value of ω was developed based on these two parameters, and it has proven superior for estimating the actual evapotranspiration (ET) by a cross-validation approach.Finally, based on the proposed semi-empirical formula for parameter ω, the contributions of changes in climate (P, ET 0 , and S) and vegetation (M) to ET variations were estimated.The results showed that the improved vegetation conditions in almost all basins made a positive contribution to the ET change, but these effects were largely offset by other variables in some basins.The contribution of landscape condition changes to ET variation will be estimated with a large error if the effects of climate seasonality are ignored.

Figure 1 .
Figure 1.Locations of the study area and hydrometeorological stations.

Table 1 .
Long -term hydrometeorological characteristics and vegetation coverage.Because a few runoff data points were missing for several basins, the data length in these basins was less than 32.Each item represents the mean annual value.

Table 2 .
Cross-validation results for each basin.

Table 3 .
Trend analysis for the hydrometeorological variables and vegetation coverage.ID BasinPeriod ET, mm yr −2 ET 0 , mm yr −2 P, mm yr −2 M, yr −1 S, yr −1 * and * * indicate the trend is significant at the level of p = 0.05 and p = 0.01 by the Mann-Kendall test, respectively.

Table 4 .
Attribution analysis for ET changes for each basin.The relative contribution of a certain variable to the ET change (ϕ(x)) was calculated as follows: ϕ (x) = (C_(x)/ ET) × 100 %, where C_(x) represents the contribution of each variable.and * * indicate the significance level of p = 0.05 and p = 0.01, respectively; ns indicates that the significance level exceeds 0.05. *