Quantifying the Impacts of Climate Change and Vegetation Variation on Actual Evapotranspiration Based on the Budyko Hypothesis in North and South Panjiang Basin, China

Actual evapotranspiration (Ea) plays a key role in the global water and energy cycles. The accurate quantification of the impacts of different factors on Ea change can help us better understand the driving mechanisms of Ea in a changing environment. Climate change and vegetation variations are well known as two main factors that have significant impacts on Ea change. Our study used three differential Budyko-type equations to quantify the contributions of climate change and vegetation variations to Ea change. First, in order to establish the relationship between the parameter n, which usually presents the land surface characteristics in the Budyko-type equations, with basic climatic variables and the Normalized Difference Vegetation Index (NDVI), the stepwise linear regression method has been used. Then, elasticity and contribution analyses were performed to quantify the contributions of different numbers of climatic factors and NDVI to Ea change. The North and South Panjiang basin in China was selected to investigate the efficiency of the modified Budyko-type equations and quantify the impacts of climate change and vegetation variations on Ea change. The empirical formal of the parameter n established in this study can be used to simulate the parameter n and Ea for the study area. The results of the elasticity and contribution analyses suggest that climate change contributed (whose average contribution is 149.6%) more to Ea change than vegetation variation (whose average contribution is −49.4%). Precipitation, NDVI and the maximum temperature are the major drivers of Ea change, while the minimum temperature and wind speed contribute the least to Ea change.


Introduction
Driven mainly by solar radiation, terrestrial actual evapotranspiration (Ea) plays a key role in the water and energy exchange between land surface and atmosphere. In the context of global change, the accurate quantification of Ea and identification of different factors (such as climatic factors and vegetation conditions) that influence Ea change are crucial for evaluating the impacts of changes in climate and vegetation on water and energy cycles, improving the accuracy of hydrological forecasting and developing a suitable strategy for water resource management.
Recently, many studies have indicated significant change in Ea due to climate change based on statistical analyses [1][2][3][4][5]. These findings have increased the interest in investigating the potential climatic driving factors that control the changes in Ea. Some previous studies indicated that the impacts of changes in different climatic factors on Ea change are quite different. For example, a change in temperature can affect Ea primarily by changing the capacity of air to hold water vapor. Increased cloudiness and decreased solar irradiance can induce reductions in Ea [6]. The decrease in solar irradiance was found to be the main cause of the decline in Ea across many regions of the world [7,8].
As an important factor related to aerodynamic turbulence in Ea calculations, wind speed affects Ea by altering the process of vapor removal, which transfers large quantities of air over the evapotranspiration surface [6]. In addition to these statistical analyses, many mathematical models of Ea estimation, such as the Penman-Monteith model [9,10], the Budyko hypothesis model [11][12][13] and the complementary relationship framework [14][15][16], have been used to analyze the effects of changes in various climatic factors on Ea change.
However, an increasing number of studies have found that climate change is not the only driver of Ea change [5,17,18]. Among the potential influential factors other than climate change, land surface vegetation conditions, which are the largest contributor to total evapotranspiration [19,20] and partly control the maximum soil water available for evapotranspiration, have received considerable attention in analyses of Ea change in recent years [12,[21][22][23]. In studies of the impacts of vegetation variations on Ea change, Shao et al. [24] and Wang and Dickinson [9] found that afforestation contributes to Ea increase at the mean annual scale and plants' growth increases Ea at the seasonal scale. Jaramillo et al. [25] indicated that a main increase of Ea is mainly due to the increase of cultivated area and/or crop production in Sweden. Zhang et al. [17] found that the significant upward trend of global Ea from 1982 to 2013 was mainly driven by vegetation greening. Previous studies have suggested that it is necessary to explicitly incorporate key indices of vegetation conditions into the Budyko framework for many catchments [22,[25][26][27][28]. Consequently, it is meaningful and important to investigate how these climate and vegetation characteristics collectively influence the spatial and temporal change in Ea.
Generally, the Budyko framework [11] has been widely used to investigate the impacts of climate change and vegetation variations on Ea change in various catchments with different climate and land surface conditions because of its conceptual appeal and the fact that it requires only routinely recorded weather data [27,25,[27][28][29][30][31][32][33]. In addition to the energy and water supplies in the Budyko framework, other factors can influence the process by which water escapes from the land surface into the atmosphere [11,22,28,34], and these factors can be numerically represented by the free parameter, n, in Budyko-type equations [35][36][37]. Moreover, many studies have found that the value of n is basinspecific [12,22,30] and suggested that it is necessary to better characterize the variables that contribute to the change in n. Based on previous studies, land surface characteristics, including climatic factors, vegetation conditions, topography and soil properties, have been identified as the potential descriptors of a catchment that control the change in n [12,21,22,28,30,33,38]. The influence of these physical characteristics on n is known to be important, and many researchers have explored how these physical characteristics collectively influence the spatial and temporal change in n [12,21,30,33]. In recent years, a few studies found that the time-varying parameter in the Budyko framework can better describe the changing processes of Ea and separate the impacts of climate change and other land-related factors on the Ea change in the changing environment, such as Jiang et al. [11], Liu et al. [39] and Tian et al. [40]. At the same time, combing the elasticity method and Budyko framework (such as the equation proposed by Turc-Pike et al. [41], Fu [35], Zhang et al., [36], Yang et al. [42] and Wang-Tang [37]), some studies specifically explored the impacts and contributions of climatic factors and land-related factors on the Ea change (e.g., Yang and Yang [43], Jiang et al. [13], Xu et al. [44] and Wang et al. [45]). However, with the time-varying parameter in the Budyko framework, limited studies have specifically investigated the impacts and contributions of different numbers of climatic factors and vegetation indices on Ea change. Therefore, we try to establish the relationship between the time-varying parameter, n, and temporally changing variables and further quantify each related climatic factor and vegetation indices on Ea change in this study.
In this study, based on the time window method, the temporal change in parameter n and the associated drivers are investigated to develop an empirical model of the time-varying parameter n by using a stepwise linear regression method. The motivation for advancing the empirical model of n is to use it in the context of three Budyko-type equations described by Fu [35], Zhang et al. [36] and Wang-Tang [37] to estimate Ea and analyze the impacts of changes in climate and vegetation coverage on change in Ea. Specifically, the objectives of this paper are: (1) To explore the temporal trends of major meteorological factors and land surface vegetation coverage from 1982 to 2013 in the study area. (2) To improve Ea estimations using the Budyko-type equations described by Fu [35], Zhang et al. [36] and Wang-Tang [37] by establishing the relationships between the time-varying parameter, n, and the indices of climate and vegetation coverage. The remainder of this paper is organized as follows. Descriptions of the datasets used in this study and the study area are introduced in Section 2. The methodology, including descriptions of the Budyko framework, moving window method, stepwise linear regression analysis, and contribution analysis, is given in Section 3. The results and discussions of the model assessment in relation to the datasets are given in Section 4. Finally, the conclusions are presented in Section 5.

Study Area and Data
The North and South Panjiang basin (NSPB) (102°12′ E-107°30′ E, 23°6′ N-26°48′ N), which is one of the Pearl River sub-basins (see Figure 1), was selected as the study catchment in this study. The drainage area of the NSPB is approximately 102,199 km 2 , which occupies about 22.5% of the Pearl River basin. The climate of this catchment is a humid monsoon climate with an annual average precipitation of 1200 mm and an annual average temperature of 16.5 ℃. The flood season begins in April and ends in September, and the precipitation during this period accounts for approximately 75% of the total annual precipitation. Table 1 presents the percentages of six major land use/cover types from 1980 to 2010 in the NSPB. It is found that the changes in each type of land use/cover from 1980 to 2010 are not obvious, and the rate of change of all types of land cover is less than 1%, thus, it can be assumed that the land use/cover in the NSPB are almost constant from 1980 to 2010.  Tian-e hydrologic station is the outlet of the NSPB, with an average runoff of about 450 mm during 1982 to 2013, and the annual runoff data from 1982 to 2013 at this outlet were obtained from the Ministry of Water Resources of China in this study. The meteorological data obtained from 1982 to 2013 at the 15 meteorological stations, which include daily precipitation (P), average temperature (T), maximum temperature (Tmax), minimum temperature (Tmin), sunshine duration (S), wind speed at a height of 10 m (u) and relative humidity (RH), were provided by the China Administration of Meteorology (http://www.data.cma.cn). Additionally, the value of potential evapotranspiration (Eo) was estimated using the FAO Penman-Monteith equation [46], which can estimate Eo accurately under different climatic conditions [47,48]. The wind speed at a height of 2 m (noted as u2) was calculated based on u using the profile relationship proposed by Allen et al. [46]. The net radiation (Rn) was calculated from S using the method recommended by Allen et al. [46], and the ground heat flux (G) was neglected because the daily and long-term average values of G are close to zero [46]. All the meteorological datasets noted above were spatially averaged by the Thiessen polygons method using ArcGIS software. To quantify the impacts of changes in vegetation coverage on Ea, the 8 km Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index (NDVI) dataset from 1982 to 2013 was derived from the Advanced Very High-Resolution Radiometer (AVHRR) sensor. Similar to climatic variables, the monthly NDVI dataset was used to represent the land surface vegetation conditions [49] of the catchment in this study. To analyze the process of land surface use/cover change in the NSPB from 1980 to 2010, a 1 km resolution land use/cover map was provided by the Resource and Environment Data Cloud Platform, Chinese Academy of Sciences (CAS) (http://www.resdc.cn/). In this study, the types of land use/cover are classified as forestland, grass land, agricultural land, built-up land, water area and unused land [21].

Budyko Framework
The Budyko hypothesis indicates that terrestrial Ea is subject to the common limitations of the water supply (usually measured by P) and energy supply (usually measured by Eo) [11]. According to this hypothesis, the water supply is the limiting factor of Ea in arid catchment (Eo/P > 1), whereas energy availability is the limiting factor of Ea in humid catchments (Eo/P < 1). Based on this approach, a coupling relationship between Ea/P and Eo/P was proposed, and the general form of this relationship can be expressed as follows [11]: where, the units of the variables in Equation (1) are mm. Thereafter, many similar equations have been developed [35][36][37]. In this study, the three Budyko-type equations proposed by Fu [35], Zhang et al. [36] and Wang-Tang [37] (see Table 2) were selected to investigate the impacts of climate change and vegetation variation on Ea change. Compared with the original Budyko-type equations, which have no basin-specific parameter, the Budyko-type equations with basin-specific parameters are more flexible to account for the effects of different catchment characteristics on Ea change [50]. In this study, the basin-specific parameters in the three Budyko-type equations proposed by Fu [35], Zhang et al. [36] and Wang-Tang. [37] are expressed as n Fu , n Z and n W-T , respectively.

Moving Window Method
In general, the basic water balance in closed catchments can be expressed as follows: where P is annual precipitation, R is annual runoff measured as streamflow and ΔS is the annual change of total water storage in the closed catchment. Because the Budyko-type equations can only be applied at the long-term average scale when ΔS is approximately zero (i.e., ΔS ≈ 0), the moving window method is applied to remove the influence of annual water storage changes in this study. Specifically, in the application of the moving window method, the width of the window size is related to the length of the variable series, and the longer window can better reduce the scatter of the variable series [13,51]; however, some specific evolution processes may be missed. According to previous studies [13,38,40], the changes of total water storage can be neglected during a long time period, and different sizes of time windows, such as 11 years and 13 years, were used in studies proposed by Jiang et al. [13] and Tian et al. [40], respectively. Here, as the P-R did not change much for time window size larger than 11 years, we selected the 11-year time window to minimize the influence of changes in ΔS series. In this study, the 11-year time window is applied to annual R, Ea, P and ΔS series, and their values in the time window centered at year t are expressed as Rt, Eat, Pt and ΔSt, respectively. Thus, in the time window centered at year t, ΔSt, which mainly includes changes in ground water and the plant moisture content of the catchment, can be neglected (i.e., ΔSt ≈ 0), and the water balance equation (i.e., Equation (2)) can be simplified as: Similar to R, Ea and P, the 11-year time window is also used for other variables, including n, T, RH, u2, Tmin, Tmax, S and NDVI, the values of which, in the time window centered at year t, are expressed as nt, Tt, RHt, u2t, Tmint, Tmaxt, St and NDVIt in the following sections.

Modeling the Parameter nt in the Budyko-type Equations
For near-surface climate conditions, the basic meteorological factors, including Pt, Tt, Tmaxt, Tmint, RHt, u2t and St, were selected as indicators to describe the variability in nt in this study. For other factors related to land surface conditions, because the change in each type of land use/cover from 1980 to 2010 is not significant, and it can be assumed that the soil conditions and topography of a single catchment remain unchanged at the multiyear average scale, the natural variation in vegetation coverage described by NDVI in this study was selected as the only land surface variable other than climate to describe the variability in nt in the study area. In summary, eight candidate variables, namely, Pt, Tt, Tmaxt, Tmint, RHt, u2t, St and NDVIt, were selected as the indicators of the catchment characteristics to investigate their respective relation to the variability in nt in each Budyko-type equation.
The stepwise linear regression method was used to determine the potential relationship between the parameter nt and the variables that were selected as the indicators of the catchment characteristics.
Moreover, an empirical model is considered to explain and formulize the potential relationships, and it can be expressed as follows: where f is a function to be determined. The parameter nt obtained in this way is called the "modeled nt", which is denoted as ' t n , and the Eat value calculated by the Budyko-type equations with ' t n is called "modeled Eat", which is denoted as ' at E . The performance of this empirical model can be evaluated by comparing the values of nt calculated by the Budyko-type equations using the observations to the values of ' t n . The resulting statistic, namely, the coefficient of determination (r 2 ) and root mean squared error (RMSE), were used to evaluate the performance of this empirical formula.

Contribution Analysis
According to Wang et al. [45], assuming that each single variable is independent in the Budykotype equations (see Table 2), the elasticity method is used to analyze the contributions of climate and vegetation changes to Eat change, and the first-order approximation of the changes in Eat can be expressed as:  Substituting Equation (7) into Equation (6) yields the following expression: and the relative contributions of single variable to Eat change can be calculated as follows: where vi is a variable (such as Pt, Tmaxt, Tmint, RHt, u2t, St), and vi  is the relative contribution of changes in vi to Eat changes.

Temporal Variations in Climatic Variables and the NDVI
To investigate the temporal variations in meteorological factors and the NDVI from 1982 to 2013 in the NSPB, the nonparametric trend-free pre-whitening Mann-Kendall test developed by Mann [52], Kendall [53] and Yue et al. [54,55] was used in this section. As shown in Figure 2, significant warming can be observed in the study area, as the annual T, Tmax and Tmin values displayed an evident increasing trend at the 0.05 significance level, which is consistent with the general trends of the global warming [56] and warming in China [57]. Wind speed has significantly decreased, which supports the results given by McVicar et al. [58]. Similar to the wind speed, annual P and RH also exhibit significant decreasing trends. In addition, S presents a slight decreasing trend, and the NDVI in the study area displays a slight increasing trend (p-value > 0.05).

Modeling '
t n A stepwise regression method was used to establish the potential relationship between the parameter nt and the eight candidate variables (discussed in Section 3.3) selected as the indicators of the catchment characteristics. The results of stepwise regression show that the parameter nt simulated with four candidate variables, including St, NDVIt, Tmaxt and Pt, yields a satisfactory determination coefficient, r 2 , between ' t n and nt calculated by the Budyko-type equations using the observations. As shown in Table 3, the first candidate variable entered into the stepwise regression equation was St, which resulted in r 2 values greater than 0.5 among the three parameters, ' t n . The other three selected candidate variables, namely, NDVIt, Tmaxt and Pt, increased r 2 by more than 0.02. The use of remaining candidate variables (such as RHt, u2t, Tt and Tmint) did not effectively increase the r 2 value. Therefore, St, NDVIt, Tmaxt and Pt were identified as the decision variables and used to predict the parameter nt in this study, and the final fitted model of ' t n in each Budyko-type equation was obtained as: Fu    Additionally, the empirical formula for nt is used together with the Budyko-type equations to calibrate the Eat estimation in the NSPB. In this study, two methods of Eat estimation were considered for comparison: (i) using the time-varying ' t n calculated by Equations (11a)-(11c) with time-varying variables and (ii) using the constant parameter nt calculated with the long-term average variables. Figure 4 shows the comparison between the values of Eat obtained from the water balance equation (i.e., Equation (3)) and these two cases. The r 2 values with a constant nt vary from 0.55 to 0.61, and the r 2 values with a time-varying ' t n are all larger than 0.90. This result indicates that the performance of each Budyko-type equation with a time-varying ' t n is better than that for a constant nt, which is consistent with the findings of Yang et al. [21], Jiang et al. [13] and Tian et al. [40]. Therefore, compared with the temporally constant parameter n which was used in most previous studies [25,32,33,44,59], the time-varying parameter n improves the accuracy of Ea simulation based on the Budyko-type equations and specifically describes the evolution process of Ea in the NSPB.   [35], Zhang et al. [36] and Wang-Tang [37], respectively.

Quantitative Contributions to Eat Variation
Based on the analysis in Section 4.2, both climate change and vegetation variations will affect Eat variations in the NSPB. Moreover, combining the empirical equation for nt estimation (i.e., Equations (11a)-(11c)) and the FAO Penman-Monteith equation, which was used to estimate Eo, the Budykotype equations used in this study can be written as at Therefore, based on Equations (8) and (9) In this section, the period of 1982-1992 is treated as the benchmark (pre-change period), and the following time widows are individually treated as post-change periods. Equation (12) is used to analyze the contribution of climate change and vegetation variations to Eat change ( at E  ) in each post-change period, which is estimated by subtracting the value of Eat in the first time window (the period of 1982-1992). The results based on the equation proposed by Zhang et al. [36] were selected for analysis in this section, because the results of all three Budyko-type equations are very similar. First, the elasticity of Eat to all variables in each time window is shown in Table 4. Notably, the elasticity of Eat to Pt, Tmaxt, Tmint and u2t is positive, which indicates that an increase (decrease) in these variables will increase (decrease) Eat in the study area. To the contrary, the elasticity of Eat to RHt, St and NDVIt is negative, which indicates that an increase (decrease) in these variables will decrease (increase) Eat. In addition, the magnitude of the elasticity also reflects the sensitivity level of Eat to these variables, the greater the absolute value of the elasticity is, the more impact the variable has on Eat. Thus, Eat is more sensitive to NDVIt and Pt than other variables in the NSPB, while the absolute values of elasticity are higher for NDVIt and Pt than for other variables.
Then, considering the changes of each variable in each time window relative to the values in the prechange period (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992) and the elasticity of Eat to all variables, the contributions of Pt, Tmaxt, Tmint, RHt, St, u2t and NDVIt to at E  can be obtained, and the results are shown in Table 5. From this table, it is found that the changes in Pt, Tmaxt and NDVIt contributed most to   The contributions to Eat variation obtained in this study matched those estimated in previous studies in China, such as by Yang and Yang [43] and Wang et al. [45], who demonstrated that precipitation changes were mainly responsible for the variations of evapotranspiration and runoff among the basic climatic variables. Additionally, for the temperature index, the elasticity of Eat to Tmaxt is significantly greater than that to Tmint in the NSPB, mainly because Tmax was used to calculate several important terms in the Eo and Ea estimation process, such as net radiation at the crop surface (Rn), the saturation vapor pressure (es), and the slops of the saturation vapor pressure versus air temperature curve (Δ). Additionally, the values of Tmax are always greater than Tmint at any time of the day and in any season [60]. Furthermore, the variation in the NDVI, which was selected as the indicator of vegetation conditions, has a large influence on the change in Ea. In fact, many studies have explored how vegetation coverage influences basin-scale water and energy balances based on the Budyko framework. For instance, Zhang et al. [36] suggested that a vegetation type-specific parameter can better express the water and energy balances in different catchments. Yang et al. [21] found that the time varying NDVI can be used to capture the variability of parameter n in the small catchments (<50,000 km 2 ) of northern China. Li et al. [22] established the relationship between the NDVI and the parameter in Fu's equation in 26 global river basins (>300,000 km 2 ). In addition, Abatzoglou and Ficklin [30] and Xing et al. [33] also used the NDVI to model n in the Budyko-type equations. The relationship between nt and NDVIt obtained in this study displays good agreement with the relationships obtained in previous studies, and it can be used to quantitatively assess the contribution of NDVIt variations to changes in Eat.

Uncertainty in the Contribution Estimation
It is worth noting that adding more variables to the elasticity and contribution estimations may generate more errors and increase uncertainty [45]. Figure 5 describes the comparison between the simulated at E  estimated by Equation (12)  was associated with the empirical model for parameter nt estimation, because the empirical models of nt established in this study cannot fully simulate the observed values and therefore introduce error into the at E  simulation. Second, the approximation of the first-order Taylor expansion in Equation (12), which ignores the second-order and higher terms of the Taylor expansion [61], may be one of the reasons for the error between the modeled at E  and the observed at E  . Additionally, there are other causes of the error between the observed Eat change and the modeled Eat change. For example, the interactions between input elements could introduce uncertainty into the attribution analysis, while some of them are not fully independent [24,33]. The temporal and spatial variability in the empirical parameters used for Eo estimation may increase the uncertainty in the attribution analysis [62]. Therefore, the errors and uncertainties that exist in the elasticity and contribution analyses must be further studied in the future.

Conclusions
In this study, an empirical model of the time-varying parameter nt was established to quantify the impacts of climate change and vegetation coverage on actual evapotranspiration. Eat in the NSPB of China was simulated, and the contribution of each climatic variable and the NDVI to Eat change was assessed based on the three Budyko-type equations proposed by Fu [35], Zhang et al. [36] and Wang-Tang [37]. The conclusions are drawn as follows.
(1) Significant warming was observed in the NSPB, and P, u2 and RH significantly decreased. For land surface conditions, the land use/cover in this catchment changed slightly from 1980 to 2010, and the NDVI, which was selected as the only land surface variable to describe the variability in parameter nt in the study area, displayed a slight increasing trend. (2) The stepwise linear regression method showed that combining three basic climatic factors and NDVIt can well explain the changes in the basin-specific nt in the Budyko-type equations, with r 2 values between the modeled ' t n and nt calculated based on the Budyko-type equations using observations all being greater than 0.78. This finding indicated that the time-varying nt proposed in this study yielded better performance of the modeling Eat than a constant nt in the study area. (3) Based on the elasticity and contribution analyses of Eat, it was found that Pt, NDVIt and Tmaxt are the major driving forces to the Eat variations in the study area, and they contributed the most to the Eat variations. During the study period, climate change contributed (whose average contribution is 149.6%) more to at E  than vegetation change (whose contribution rate is -49.4%). Furthermore, the results also indicate that in addition to climatic factors, considering the impacts of vegetation coverage variation on actual evapotranspiration and the water cycle is very important in the study area.