Comparison of the effectiveness of four Budyko-based methods in attributing long-term changes in actual evapotranspiration

The responses of hydrological processes to climate change and anthropogenic influence have received significant attention over the past few decades. Several Budyko-based methods have been widely used to attribute hydrological variations and identify the extent of variation due to climate change and human activities. However, the accuracy of various methods has rarely been compared. This study employed four Budyko-based methods, namely the total differential method, complementary method, extrapolation method, and decomposition method, to attribute the changes in actual evapotranspiration in 13 basins in China’s Loess Plateau. We compared their performances and analysed factors that contribute to the differences in attribution results yielded by the various methods. The results showed that the total differential, complementary, and decomposition methods presented similar estimates of the contributions of climate change and human activities. However, the extrapolation method showed a large deviation in the contribution of human activities. The error of the extrapolation method was the largest, followed by that of the two-stage total differential method. The complementary method and decomposition method exhibited negligible errors.

As one of the major components of the hydrologic cycle, actual evapotranspiration (ET) drives energy and water exchanges between the hydrosphere, atmosphere and biosphere 1 . However, because of complexity of the land-plant-atmosphere system, direct measurement of ET is time consuming and expensive, and is hard to launch at regional scales 2,3 . Several methods have been proposed for evaluating regional ET, mainly including meteorology-driven diagnostic models, satellite data-driven approaches, energy-water balance methods, vegetation index-ET empirical relationship methods and statistical methods 4 . Among them, water balance method has been widely used because they are simple and sound in theory 5 . The long-term average water balance is in a steady state, and the water storage change in a catchment can be negligible. However, when it comes to shorter time scales, water storage change can be great 6 . To minimize the potential errors introduced by neglecting water storage variation, the hydrological year was introduced to estimate the annual ET [7][8][9] . In this way, the water input occurs mainly at the beginning of the year and the water is consumed within that year 9 .
Global hydrological research over the past decades has focussed on variations in the catchment water balance under changing environments where changes in processes such as runoff and ET have been assessed [10][11][12] . Climate change and human activities are considered as the major drivers that influence water resources and the catchment water cycle 13,14 . Quantification of the contributions of climate change and human activities on ET is important for efficient water resource management and ecosystem sustainability. Three kinds of approaches have been developed to evaluate the effects of climate change and human activities on hydrology, namely, the paired catchments approach, hydrological modelling, and statistical methods 12,15 . However, the accuracy of these methods has varied across applications [16][17][18][19] . For example, Jiang, et al. 20 used three different methods namely, multi-regression, total differential method, and the hydrologic model to examine runoff variation in the Laohahe basin in north SCIENtIFIC RepoRTs | (2018) 8:12665 | DOI: 10.1038/s41598-018-31036-x China 11 . They found that the contribution of human activities to runoff reduction estimated using the total differential method was slightly lower than values estimated by the other methods. Li, et al. 21 applied two approaches including a sensitivity-based approach (comprising of a non-parametric model and six Budyko-based models) and hydrological modelling approach (employing the Xinanjiang and SIMHYD models) to three catchments in Australia. They concluded that the non-parametric method and Xinanjiang model underestimated the contribution of vegetation changes and overestimated the contribution of climate variability to the runoff reduction. Zhao, et al. 22 attempted to quantify the impacts of human activities and climate change on runoff change in the middle reaches of the Yellow River. They noted that the results of the multi-regression method and total differential method presented different estimations for all the catchments assessed. Therefore, the issue of inconsistency in results derived using different methods has to be investigated. The accuracy of various methods in attributing runoff and ET variation should be compared and reasons and mechanisms contributing to the differences should be investigated as it can help us to understand the effects of climate change and human activities on the changes in hydrological processes more accurately.
The Budyko framework 23 proposes that the long-term average water balance is controlled by the supply and demand of water from the atmosphere (represented by precipitation and potential evapotranspiration, respectively). The framework has recently been extensively applied in hydrological research due to its comprehensiveness and effectiveness in studying the effects of climate and surface condition on water resources 5,24,25 . Moreover, the framework does not depend on extensive historical climate and runoff data 26 . In past half-century, a number of Budyko equations have been proposed for the Budyko curve 5 . Among them, the Fu and Choudhury-Yang equations have been used widely. Recently, Zhou et al. 27 found a generating function that can be used to produce any number of valid Budyko equations, which reflects the existence of functional relationship between these Budyko equations. Therefore, the similar performance of water balance variation estimation and its attribution could be obtained using different Budyko equations.
The Budyko-based attribution methods for catchment ET/runoff variation mainly include the total differential method 11 , extrapolation method 28 , decomposition method 29 , and complementary method 30 . However, most previous researches have only employed one of the Budyko-based methods to quantitatively identify the contributions of climate change and human activities on runoff/ET variation. Few researches have compared the similarities and differences of the attribution results estimated using different Budyko-based methods. Therefore, the objective of this study was to compare the effectiveness and applicability of these four different methods in estimating the relative contributions of climate change and human activities on ET variations in 13 basins in China's Loess Plateau.

Results
Comparison of different total differential methods. The total differential method is one of the most popular attribution methods for water balance variation under the Budyko framework. ET change induced by a certain variable was indicated by multiplying the change value of the variable by its partial derivative (Eq.(6)). The partial derivatives can be estimated using the mean annual data for the entire study period [31][32][33] or for the pre-change period 11,28 . In this study, we referred to the two methods as the whole-period method and forward-approximation method respectively. The average of the partial derivatives for the pre-change period and postchange period is also popularly used; we called it as the two-stage method. All three forms of the total differential method were employed in this study to identify the best-suited one.
The partial derivative values of ω, ET p , and P estimated using the three total differential methods were compared graphically, as shown in Figure 1. As evident, the three groups of partial derivative values were different; specifically, the values for ET p and P varied significantly between the methods. Moreover, the partial derivative values estimated using the forward-approximation method were compared with values calculated using the two-stage method; ∂ET/∂ET p was larger while ∂ET/∂P was smaller using the forward-approximation method for all the catchments with an absolute mean difference of 17.8% and 4.3%, respectively. Further, ∂ET/∂ET p and ∂ET/∂P estimated using the whole-period method was smaller than the values derived using the two-stage method for most catchments with an absolute mean difference of 8.5% and 1.6%, respectively.
The contributions of the three variables (P, ET p and ω) to ET variation in the 13 basins were estimated using Eq. (6) ( Table 1). The errors were evaluated by comparing the sum of contributions of the three variables (C_sum) and the actual change in ET (∆ET). Overall, the errors were the largest for the forward-approximation method, intermediate for the whole-period method, and smallest for the two-stage method at average values of 5.2 mm, 1.7 mm, and 0.5 mm across the 13 catchments, accounting for 16.4%, 5.7% and 1.4% of the observation respectively. The two-stage method appeared to be superior to the other methods in attributing ET variation. Thus, we selected it for comparison with the other three Budyko-based methods in section 3.2.
Comparing the impacts of climate change. The estimated contributions of climate change on the variation in ET (ΔET c ) in the 13 basins calculated using the four Budyko-based attribution methods were compared. The results are listed in Table 2. Generally, the ΔET c estimated using the four methods were close. Moreover, the results of the total differential and complementary methods were the same, suggesting that all four methods were appropriate to calculate the contributions of climate change on ET variation. However, differences existed between these methods. The absolute mean relative difference between ΔET c calculated using the decomposition method and total differential method was 2.8%. Further, their values were higher than 5% in basins 1, 10, and 13. ΔET c calculated by the extrapolation method was smaller than the value calculated by the total differential method and complementary method in most of the basins with a mean relative difference of 5.4%; the value was up to 13% in basin 13.
During the assessment of the contribution of climate change on catchment water balance change, researchers attempt to identify whether the water condition (represented by P) or heat condition (represented by ET p ) SCIENtIFIC RepoRTs | (2018) 8:12665 | DOI:10.1038/s41598-018-31036-x dominated the water balance change [34][35][36] . However, in this case, the contributions of P and ET p cannot be separated further using the decomposition method and extrapolation method.
Comparing the impacts of human activities. Similarly, ∆ET h values were also obtained by the total differential, complementary, and decomposition methods, and gaps were noted (Table 3). Unlike the values of ∆ET c , ∆ET h values calculated by the total differential method were smaller than values calculated by the complementary method in most basins, with an absolute mean relative difference of 1.4%. The largest gap appeared in basins 5 and 6 with a difference of −4.4% and 4.1%. However, ∆ET h values calculated by the extrapolation method exhibited large difference compared to values calculated by the other three methods. For example, the absolute mean relative difference of ∆ET h between the extrapolation method and complementary method was 31.7% and the relative difference in basins 8 and 15 were as high as 90.2% and 116.9%, respectively, suggesting that the extrapolation method was not suitable for estimating the contributions of human activities on ET variation in this study.   Comparing errors of the four methods. The changes in ET calculated by the four methods, i.e. the sum of ΔET c and ΔET h were compared to the observed ET (∆ET). The results showed that the errors of the extrapolation method were the largest with mean relative value of 13.3%; the error exceeded 40% in basins 7, 13, and 14. The two-stage total differential method ranked next, with a mean relative value of 1.4%. The complementary and decomposition methods exhibited negligible errors. However, it should be noted that the decomposition method directly used the difference between ΔET and ΔET h as ΔET c , and its error was hidden in ΔET c .

Discussion
In this study, we compared the effectiveness of four Budyko-based methods in estimating the contributions of climate change and human activities on the variation in ET. The results highlighted that the estimated changes in the ET varied between methods because of their assumptions and structure. This section discusses some possible reasons for the difference in performance of each method.
Uncertainties of the decomposition and total differential methods. The decomposition method hinges on two assumptions. Firstly, ET/P follows the same Budyko curve as the pre-change period, which implies that the controlling parameters in the Budyko equations are constant. However, Ning, et al. 9 reported that the values of the parameter ω in all the basins have exhibited an upward trend, and this trend has been more dramatic after the 1980 s. For example, in the Huangfu basin (no. 1), the mean annual ω in the pre-change and post-change periods was 2.32 and 2.56, respectively, and their difference was 0.24. The mean partial derivative coefficient of ET to ω was 46.9. Thus, if the ω values of different periods were used to calculate the contribution of human activities to ET variation, the error would increase to 11.3 mm. Therefore, the constant controlling parameter is an error source, and also applies to the total differential method and complementary method. Secondly, the decomposition method also assumed that P and ET P are stationary during the pre-change and post-change periods 29 . However, both natural climate variability and human activities can induce change in P and ET P at the interdecadal or interannual scales.
Moreover, the partial derivative coefficients in the total differential method should be calculated using the complete Taylor expansion to deliver accurate results. However, in our study, the partial derivative coefficients were mostly estimated as a first-order approximation. Yang, et al. 37 demonstrated that ignoring the higher orders of the Taylor expansion will underestimate the contribution of climate change to hydrological processes when precipitation increases or potential ET decreases. Furthermore, the forward-approximation and whole-period total differential methods actually assume that the partial derivative coefficients of ET to P, ET p , and ω did not change during the study period. However, the annual partial derivative coefficients for the period 1961-2010 for the 13 basins calculated in this study showed that the partial derivative coefficients of ET to ET p and ω in all the basins exhibited significant downward trend from 1960 to 2010 (p < 0.05). Conversely, the partial derivative coefficients of ET to P in all the basins showed significant upward trend (p < 0.05). This implies that using the constant partial derivative coefficients could result in large errors in the estimated ∆ET.

Possible reason for large error of the extrapolation method. A comparison of the error in ΔET h
derived using the four methods indicated that the largest derivation was produced by the extrapolation method. The extrapolation method used the sensitivity analysis and the trend analysis to evaluate the effects of climate effect and human activities, respectively. And the large error of this method was mainly induced by the trend analysis. Specifically, early occurrence of the abrupt year resulted in small sample size for the relationship development in Eq. (9), and subsequently led to errors when the relationship was applied to the post-change period. Thus, the unstable relationship between ET and P in the pre-change period is the most dominating factor leads to the inaccurate result in this study. To address the problem of small sample size in the pre-change period, Zhang, et al. 28 used monthly runoff and precipitation to build the relationship. In most cases, annual ET was obtained by ignoring water storage change, which was acceptable if this was carried out for multiple hydrological years [7][8][9] . If the relationship between ET and P is built on a monthly scale, ET estimation would exhibit higher uncertainty as water storage variations on a monthly scale cannot be ignored. This argument indicates that even when the sample size is large enough and the determining coefficient of the developed relationship is high, the errors would be large. Further, climate change may make some contribution to ΔET h in this method and induce the attribution error. For instance, the P post series was used to predict the ET post ; however, the P post itself must have already been influenced by the climate change.
Other potential factors influencing attribution results. Many previous studies 25,31,32 have considered that the controlling parameters of the Budyko model represent the effects of the land use/cover change (i.e. human activities) completely. However, several researches have highlighted that climate seasonality also strongly affects the controlling parameters of the Budyko models 9,30,37,38 . In particular, our previous study found that parameter ω was well-correlated with vegetation coverage and climate seasonality, and if the impacts of climate seasonality on ET were ignored, the contribution of vegetation (i.e. land use/cover) will be estimated with a large error 9 . Thus, the effects of climate seasonality should be considered when attribution analysis of catchment runoff or ET variation is conducted. Further, although the attribution equation requires the variables to be independent, major factors such as potential evapotranspiration, precipitation, and controlling parameters can interact with each other, increasing uncertainties.

Conclusions
The water cycle in many regions of the world has been affected by climate change and human activities and has experienced drastic changes in the recent decades. Budyko-based methods are popularly used to detect hydrological variation attributed to climate change and human activities due to their simplicity and clear physical foundations. The methods mainly include the total differential, complementary, extrapolation, and decomposition methods. However, little attention has been paid to assessing and comparing the effectiveness of these methods. The study obtained the annual evapotranspiration data of 13 typical basins in China's Loess Plateau from 1961 to SCIENtIFIC RepoRTs | (2018) 8:12665 | DOI:10.1038/s41598-018-31036-x 2012 using the hydrological year approach, and compared the performance of the various Budyko-based methods in quantifying the relative contributions of related environmental/anthropogenic factors to long-term ET variation. The following main conclusions were drawn: in the total differential method, use of average values of the partial derivatives for the pre-change and post-change periods to conduct attribution analysis of ET variation is better than using long-term mean annual data for the whole study period or for the pre-change period. The contribution of climate change (∆ET c ) and human activities (∆ET h ) estimated by the total differential, complementary, and decomposition methods were similar. Similar ∆ET c could also be obtained by the extrapolation method. However, ∆ET h values calculated by this method exhibited large difference compared to the other three methods. The error of the extrapolation method was largest, followed by the two-stage total differential method. The complementary and decomposition methods recorded negligible errors.

Data and Methods
Study area. The Loess Plateau is located in the upper and middle reaches of the Yellow River in North China ( Figure 2), and extends over an area of 6.4 × 10 5 km 2 . Most areas experience semi-arid and sub-humid climate. Based on the daily meteorological data of 96 stations in the Loess Plateau from 1961 to 2012, the spatial and temporal change trend of precipitation (P), mean temperature (T) and potential evapotranspiration for whole Loess Plateau were exhibited. Spatially, P decreases along the southeast-northwest direction, ranging from 200 mm to 750 mm, 80% of which is distributed during June to September. T decreases from the southeast to northwest, ranging from 1.5-15 °C. ET P , calculated by the method of Priestley and Taylor 39 , distinctly decreases in three directions, namely, northwest-east, northwest-south, and northwest-southwest, ranging from 812-1235 mm. However, temporal variation in these variables has been noted; P is decreasing slowly by 5.2 mm/decade (p > 0.05); while T has increased significantly by 0.3 °C/decade (p < 0.05). ET P has exhibited a downward trend at −3.1 mm/decade (p > 0.05). The soil erosion in the study area is very heavy as it is influenced by scarce and unevenly distributed rainfall, steep landscape, highly erodible loessial soil, low canopy cover degree, and long history of intensive cultivation. Human activities in this region are also intensive, mainly refers to the implementation of soil and water conservation measures. The construction of energy bases and new towns have also significantly altered the land surface and consequently affected the regional water cycle.

Data.
Monthly runoff records for the 13 catchments from 1961 to 2012 were obtained from the Yellow River Conservancy Commission. Detailed characteristics of each catchment are presented in Table S1 in the Supporting information. Daily meteorological data for 96 stations, including precipitation, temperature, sunshine hours, and relative humidity for the period 1961-2012, was acquired from the China Meteorological Administration.

Methods
Abrupt point test. The non-parametric Pettitt method was used to detect the abrupt year of ET; detailed procedure can be found in Mallakpour and Villarini 40 . Once an abrupt year is identified, the entire study period is divided into two sub-periods: 1961 to the abrupt year, referring to as the pre-change period, represented the initial state of ET when no significant human activities occurred; and the residual data period, referring to as the post-change period, represented changed ET and is associated with significant human activities. Attributing changes in ET. The catchment water balance can be expressed as: where R and ΔS are the runoff and change in water storage respectively. ΔS can be ignored for long-term assessments but can be significant for assessments conducted for shorter time scales. Thus, the hydrological year concept (July to June of the following year) 9 was introduced to assess the annual P and R to minimize the errors that may occur in the annual ET calculation due to the exclusion of ΔS To estimate ET, the Budyko formula expressed by the Fu equation was used 41 : where ω is the controlling parameter that determines the shape of the Budyko curve, which is used to represent the effects of land use/cover change induced by human activities in the long-term. ET P is the potential evapotranspiration calculated using the equation suggested by Priestley and Taylor 39 . The change in ET from pre-change period to post-change period can be described as: where ET post and ET pre are the mean annual ET values during the post-change and pre-change periods respectively. The total contributions of climate change and human activities can be expressed as: where ΔET c and ΔET h are climate-induced and human-induced ET changes, respectively. The absolute and relative errors between C_sum and ΔET are calculated as: Budyko-based attribution methods. Four methods, namely, the total differential method, complementary method, extrapolation method, and decomposition method were used to calculate ΔET c and ΔET h .
(1) Total differential method Assuming that ET P , P, and parameter ω in Eq. (1) were independent variables, the total differential of ET can be used to evaluate the contributions of major climate factors and human activities on the long-term scale 2 : Ignoring the higher orders of the Taylor expansion in Eq. (6a), the equation can be rewritten as: where ΔP, ΔET p , and Δω are the changes in the annual P, ET p, and ω from the pre-change period to post-change period, respectively.  30 to reduce the errors of the total differential method caused by the exclusion of the higher orders of Taylor expansion 21 . It can be expressed as: where α is a weighting coefficient that varies from 0 to 1, which can determine the bounds of P, ET p, and ω effect. In this study, we defined α = 0.5 according to the recommendation of Zhou, et al. 30

Extrapolation method
The relationship between P and ET during the pre-change period is considered as: Specifically, the function f for each catchment was fitted by the annual ET and P series using the linear regression method.
Next, the function f was applied to the post-change period as: where ET post and ′ ET post are the mean annual values of actual and predicted ET during the post-change period respectively. Further, ΔET c was calculated using the forward approximation method.

(4) Decomposition method
The decomposition method is presented in Fig. S1 in the Supporting information. Ignoring the intra-annual variability of P, the contribution of human activities to runoff change was defined as: Similarly, the contribution of human activities to variation in ET can be expressed as: where ′ ET post is the mean annual ET calculated by Eq. (2) using the mean annual P and ET p during the post-change period and parameter ω during pre-change period.
ΔET c can be calculated as: c h