Impacts of climate and vegetation on actual evapotranspiration in typical arid mountainous regions using a Budyko-based framework

Evapotranspiration (ET) is a crucial component in hydrological processes, which could account for up to over 60% of the precipitation on a global scale (Oki & Kanae 2006) and more than 90% in arid regions (Glenn et al. 2007). Given that the differences between precipitation and evapotranspiration (ET) have been commonly used to evaluate the water availability at catchment scale (Falkenmark et al. 1989), it is of vital importance to quantify the actual evapotranspiration (ET_a) and its influencing factors to assist decision making in water resources management and ecological protection.

Although it is often difficult to measure the ET_a directly either at field scale or at catchment scale, several ET_a quantification methods have been proposed including meteorological (Penman 1948; Bouchet 1963; Monteith 1965; Priestley & Taylor 1972; Allen et al. 1998), water balance (Liu 2018; Rodell et al. 2004, 2011; Zeng et al. 2013), energy balance and remote sensing techniques (Menenti & Choudhury 1993; Bastiaanssen et al. 1998; Su 2002). Particularly, the Food and Agriculture Organization (FAO) Penman–Monteith method has been globally recommended to estimate the crop reference ET (ET₀) (Allen et al. 1998). Based on the results of ET₀, the ET_a could therefore be estimated by multiplying the ET₀ by an empirical crop coefficient K_c (Allen et al. 1998; Xu et al. 2006; Drexler et al. 2008; Zhan et al. 2012; Yang et al. 2017). When there are a lack of crop coefficients in determining the ET_a, especially in mountainous regions where K_c is likely hard to be obtained (Drexler et al. 2008; Reddy 2015; Yang et al. 2017; Wang et al. 2019), development of the remote sensing-based evapotranspiration datasets provides a proxy in interpreting the spatiotemporal characteristics of the ET_a at catchment scale (Mauser & Schädlich 1998; Golubev et al. 2001; Xu et al. 2006; Bai & Liu 2018; Lu et al. 2019). These datasets are generated based on the implementation of remote sensing techniques in simulating the variation of ET, e.g. the MOD16 global evapotranspiration (ET) datasets (hereafter: the MODIS ET) (Cleugh et al. 2007; Mu et al. 2007, 2009, 2011), the Global Land Evaporation Amsterdam Model (GLEAM) datasets (hereafter: the GLEAM ET) (Miralles et al. 2011; Martens et al. 2017), the CSIRO PML ET (Commonwealth Scientific and Industrial Research Organisation Penman–Monteith–Leuning ET) data collection (Zhang et al. 2016) and the TerraClimate ET (Abatzoglou et al. 2018). For example, Tang & Zhang (2011) evaluated the performances of the MODIS ET in estimating the daily variation of ET in the Yellow River Basin. Li et al. (2017) applied the MODIS ET and another two ET products in capturing interannual variation of ET_a in the middle Yellow River Basin. Yang et al. (2018) quantified the spatiotemporal characteristics using the GLEAM ET over 10 basins in China.

It has been widely accepted that evapotranspiration (ET) changes with climate variability (Brümmer et al. 2012; Feng et al. 2013, 2020; Christoffersen et al. 2014; Beaulieu et al. 2016) and underlying surface conditions, particularly vegetation conditions (Donohue et al. 2007; Huo et al. 2013). As the only link between water and energy, ET is significantly affected by climate, in drought periods, temperature and sunshine hours had positive effects on the variation of ET, while other key climatic factors (i.e. effective precipitation, vapour pressure, wind speed, relative humidity) had negative effects on ET (Feng et al. 2013). As the global climate becomes warmer and warmer, it is generally expected that the evaporation will increase (Roderick & Farquhar 2002). Goyal (2004) suggested an increase of 14.8% of total ET demand with an increase in temperature by 20%. Xu et al. (2018) suggested that air temperature was one of the main climatic factors affecting ET_a trends, on the annual scale, as the temperature increases, ET_a decreases in the upper reaches of the Yellow River Basin. However, even as increasing air temperature has been observed over the last few decades, the evaporation could show a downward trend, known as the ‘evaporation paradox’ phenomenon (Peterson et al. 1995; Brutsaert & Parlange 1998; Roderick & Farquhar 2002). Feng et al. (2020) indicated that the dominant factor of ET variability varied with timescales and dry-wet conditions, and the precipitation plays a dominant role in ET variability in most arid regions. Brümmer et al. (2012) found that as the precipitation increased, ET tended to be stabilized. Yang et al. (2003) found that in arid and semi-arid regions of China, both precipitation and evaporation are decreasing, but the precipitation decreases faster than evaporation. ET_a was influenced by the vegetation greening degree, which is manifested in the increase in evaporation caused by vegetation coverage, and vice versa (Beaulieu et al. 2016). In arid areas, Bai et al. (2019) found that the inter-annual change of ET_a was driven by the NDVI (i.e. Normalized Difference Vegetation Index) with a contribution rate of 93%. He et al. (2019) found that the average impacts of vegetation restoration on ET change was 26.9% in the Loess Plateau, China. However, impacts of climate and vegetation on the variation of the ET_a have not reached a consensus.

According to the Budyko hypothesis, water and energy supply conditions (represented by precipitation and potential evapotranspiration, respectively) are the two main factors controlling the actual evapotranspiration in the basin (Xu & Singh 2005), and are widely used in the study of evapotranspiration in the basin due to their concise expression and clear physical meaning. The Budyko framework provides a simple first-order relationship to estimate water availability represented by the ratio E/P as a function of the aridity index (Ep = P, with Ep denoting potential evaporation). The Budyko method is widely used in analyzing climate and vegetation factors, including arid regions and mountainous areas. Previous studies focus on the influence of ET in arid regions with a long series of climatic changes, while others focus on the influence of climatic extremes on ET in arid regions. However, there are few studies on the influence of vegetation and seasonal climatic changes on ET. Greve et al. (2016) proposed an analytically derived modification of the Budyko framework, that explicitly accounts for conditions under which ET is also driven by other water sources than precipitation (i.e. changes in water storage). Ning et al. (2018) 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. The results 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 will be estimated with a large error. Gao et al. (2017) used the Budyko framework to simulate the actual water consumption of 161 sub-basins from 1990 to 2014. The study indicates that with the vegetation recovery, the percentage of ET in the total precipitation keeps decreasing. However, at present, in arid mountainous regions research on these coefficients and parameters is still limited. It should be noted that large differences still exist in these studies, especially in quantifying the impacts of climate factors and vegetation conditions on ET_a. Particularly, in arid mountainous regions, there is still an urgent need to conduct further studies on the application of the Budyko method in revealing the climate-vegetation-hydrology interactions.

The objectives of this study were to investigate impacts of climatic factors and vegetation conditions on the variations of ET_a. Firstly, we quantify the changes of the climatic factors and vegetation conditions using ground-based and satellite-based observation data. Variations of the ET_a were validated by applying a Budyko-based empirical formula, which was fitted by the data of climatic factors and vegetation conditions. The water-heat coupling control parameter ω was therefore quantified according to the empirical formula. Based on the above results, we made an attribution analysis on the variation of ET_a to the climatic factors and vegetation coverage by applying a total differential method. Case studies are conducted in two reservoir watersheds in the Daqing River Basin (DRB), which act as two of the most important water supply regions in this river basin. It is interesting to note that the Xiong'an New Area, which depends largely on these two watersheds in the river basin, has received widespread attention from the scientific community, leading to an emerging focus on the hydrological processes and attribution anlaysis in this region.

We selected the Wangkuai Reservoir Watershed (WRW) and Xidayang Reservoir Watershed (XRW) as the study area, which are two typical arid mountainous regions located in the DRB, North China. These two watersheds act as the main water storage in the river basin and provide more than 80% of water resources to downstream of the DRB. These two watersheds are of great significance as the two reservoirs are two of the multiple water sources of the Xiong'an New Area, which faces large challenges in guaranteeing its water security in the coming future. The geophysical location of the study area is shown in Figure 1. The two watersheds are located between 113°40′–115°E and 38°30′ ∼ 39°45′N in the mountainous region of DRB, where the vegetation coverage area accounts for more than 70% of the watersheds. Mountains and plains characterize the topography of DRB. The northwest is mainly composed of mountainous areas, the middle region is a plain area, and the eastern low-lying area is the Baiyangdian Lake area. This region has a semi humid continental monsoon climate.

The areas of WRW and XRW are 3,662 and 1,296 km², respectively. The average annual precipitation of WRW is about 416 mm, while it is about 394 mm in XRW. The average outbound flow of WRW is 240 million m³, and that of XRW is 130 million m³.

Daily meteorological data for the period 1961–2018, including precipitation, temperature, sunshine hours and relative humidity, were obtained from the China Meteorological Administration (CMA, http://data.cma.cn/). Digital elevation data, at 30 m resolution, was provided by the Geospatial Data Cloud site, Computer Network Information Centre, Chinese Academy of Sciences (www.gscloud.cn/).

In this study, MODIS satellite remote sensing data MOD13A2 were selected to analyze the NDVI variations characteristics in the DRB (Daqing River Basin). The MOD13A2 product provides Vegetation Index (VI) values at a per pixel basis at 1 km spatial resolution. It contains the NDVI, which is referred to as the continuity index to the existing National Oceanic and Atmospheric Administration-Advanced Very High-Resolution Radiometer (NOAA-AVHRR) derived NDVI. The algorithm for this product chooses the best available pixel value from all the acquisitions from the 16-day period (Holben 1986). The criteria used is low clouds, low view angle and the highest NDVI value. The MODIS product NDVI data from 2000 to 2015 were obtained and used to estimate the dynamic of vegetation activity in the DRB. Data is obtained from geospatial data cloud sites (www.gscloud.cn/).

In this study, MODIS satellite remote sensing data MOD16A2 were obtained and used to analyze the actual evapotranspiration from 2000 to 2015 in the DRB. Provided in the MOD16A2 product are layers for composited evapotranspiration (ET), latent heat flux (LE), potential ET (PET) and potential LE (PLE) along with a quality control layer. The pixel values for the two evapotranspiration layers (ET and PET) are the sum of all 8 days within the composite period and the pixel values for the two latent heat layers (LE and PLE) are the average of all 8 days within the composite period. ET was extracted for analysis of actual evapotranspiration in this paper. The datasets were 500 m resolution and 8-day evapotranspiration products. Data was obtained from the net of EarthData (https://ladsweb.modaps.eosdis.nasa.gov).

According to Fu's formula, when → 1, ET → 0. Therefore, the values of M and S are: when ET → 0, vegetation transpiration T → 0, that is, M → 0. Or when runoff R → P, the matching between precipitation and potential evapotranspiration tends to be the worst, that is, S → ∞.

In order to verify the performances of the empirical formula, we tested the above formulas by estimating the annual ET of the two watersheds against the MODIS ET measurements combining with the Fu's function (Zhang et al. 2004).

represents the contribution of meteorological factors to ET change. represents the change of ET on the multi-year scale, mm/yr. In order to simplify the calculation process, the contribution can be approximately calculated by multiple linear regression method.

The moving average method was used to test the variations trend of meteorological data and NDVI sequence, and the Mann–Kendall method was used to test the significance of the trend (p ≤ 0.05). Analysis of variance (ANOVA) was used for the difference test of M and S. The linear correlation method was used to analyze the correlation between vegetation and climate seasonality factors to ET_a.

The annual temperature and precipitation of WRW and XRW showed a rising trend (Figure 2).

Similar increasing trends are witnessed in the two watersheds in terms of the vegetation coverage index (M), with the mean M value of the WRW and XRW being 0.49 and 0.45. The slope of M is 0.01 in both WRW and XRW, and the fluctuation is similar. (Table 1 and Figure 3(a)). The vegetation index changes from 0.41 to 0.54 in the past 16 years for the WRW, while the index varies from 0.37 to 0.51 for the XRW. Significant increasing trends are observed for these two watersheds from 2000 to 2015 (p < 0.05) (Table 1 and Figure 3(a)). It can be seen from the ANOVA analysis results that significant differences have been found between the variations of the vegetation cover trends (p < 0.05, F > F-crit, Table 1).

In terms of the climate seasonality index (S), it is shown that a contradictory trend has been found in the two watersheds. In WRW and XRW, the slope of S is 0.02 and –0.001 respectively (Figure 3(b)). For the WRW, the S tends to show an overall upward trend, ranging from 0.06 to 0.90 with a mean value of 0.49. Large annual variations are found in the WRW, the maximum differs from the minimum by 0.85 (Figure 3(b)), while a slight drop of the S in the XRW is shown, ranging from 0.02 to 0.19. Fewer annual variations exist compared with those of the WRW. It should be noted that large significant difference between the S values of the two watersheds have been found (p < 0.01, F > F-crit) (Table 1).

These formulas are verified by comparing the estimated ET by Fu's equation and measured ET from the MODIS ET data. The verification results showed that both the formulas could estimate the ET_a well at annual scale, with a low value of MAE and RMSE and good performances in terms of NSE (Table 2 and Figure 4). There are many methods to estimate ET, among which the Romanenko (Rom1, Rom2) (Romanenko 1961) and Penman method simplified by Linacre (L-P) (Linacre 1977) have been widely applied. We estimated the variations of ET by Rom1, Rom2 and L-P. The correlation coefficient R² of the three methods in WRW were 0.69, 0.65 and 0.67, and in XRW were 0.81, 0.72 and 0.80. Although the overall fitting effect of these methods was good, the accuracy of ET_a is less than that of Fu's method (Table 3). Therefore, we deduce that Fu's method could be more applicable in estimating actual evapotranspiration in this study.

The results showed that the ET_a presented a significant increasing trend in the two watersheds based on the MODIS ET data (Figure 5). The ET_a of these two watersheds rises swiftly from 289 to 415 mm within 16 years, with a mean value of 388.94 and 360.41 mm for the WRW and XRW, respectively. In this study, we attempted to determine the effects of vegetation and climatic index on the annual variation of ET_a in the two watersheds via attribution analysis.

The correlation analysis shows that the ET_a increases with the rising of vegetation coverage at annual scale (Figure 6). Significant positive correlation has been found between the variation of ET_a and the vegetation index (M) in these two watersheds (Figure 6). However, no significant linear relationships exist between ET_a and the climate seasonality index (S) (Figure 7). Particularly, the climate seasonality index is positively related with the ET_a in the WRW, while a negative relationship is found in the XRW. These results indicated that the variation of vegetation coverage may exert a consistent effect on the variation of ET_a, while the effects of climatic factor may change across watersheds (Zhang et al. 2001; Donohue et al. 2007, 2010).

We also investigate the contribution of variation of climate and vegetation factors to the changes of ET_a by applying a total difference method. Results showed that the vegetation index explained up to 42.15–58.56 mm of the variations of the ET_a in the two watersheds, which is determined as the main influencing factor to the change. However, the climate seasonality index tends to lead to a decreasing trend of ET_a in the two watersheds, with negative contribution rates of 11.48 and 13.47 mm, respectively (Table 4).

According to the results analysis of the contribution, the contribution of precipitation is second only to vegetation (Table 4). The precipitation is more sensitive than other climate factors in this area (Figure 8).

In this study we used the observed meteorological data from two reservoir watersheds in DRB. Additionally, we applied the data of climate factor and vegetation coverage to analyze the variation of ET_a in WRW and XRW of DRB. The results showed that there was a larger correlation between vegetation factors and ET_a, and the contribution of climate factors was slightly less than that of vegetation factors. We deduced that the vegetation factor is likely to be the main contributing factor to the variations of ET_a.

There have been many methods for the quantitative attribution study of evapotranspiration time change, such as linear multiple regression (Shan et al. 2015), detrending method (Xu et al. 2006) and sensitivity analysis (Wang et al. 2013). However, the research on the spatial variation attribution of evapotranspiration is mostly in the qualitative description stage (Li et al. 2012; Zhang et al. 2015). The Budyko method considers the effects of vegetation, topography and meteorological factors on evapotranspiration, which is conducive to a more detailed study. The above findings are generally consistent with previous studies on the ET_a attribution analysis across arid regions. For example, Ning et al. (2020) found that vegetation coverage averaged across all the catchments has increased obviously, and has further led to the increase in ET_a in China, especially in the dry catchments, while the variation amplitude of precipitation was small. Liu et al. (2019b) found that greening was the major reason for the ET_a increases in the Ziya-Daqing Basins, which are partially included in the mountainous regions of the Haihe River Basin. Xu et al. (2018) found that large-scale coherent vegetation greening would be the primary factor affecting inter-annual ET_a in comparison with climate factors in the upper reaches of the Yellow River Basin. Ma et al. (2019) found that vegetation greening was the dominant factor to the increases of ET_a in mixed forest and grassland in the Loess Plateau. In these studies, we found that the variation amplitudes of their vegetation were larger than those of precipitation and temperature. Moreover, compared with the correlation between precipitation and temperature on ET_a, vegetation has a greater correlation with ET_a. Particularly, in the Daqing River Basin where the two reservoir watersheds are located, Cui et al. (2019) revealed that there exists a significant increasing trend of vegetation coverage, especially in the mountainous regions. This may serve as an important reason for the increase of ET_a in the DRB. Recent studies have attempted to elaborate the mechanism of vegetation and ET_a interactions in arid regions. Findings of these studies indicate that the greening of vegetation could lead to the rising of ET_a in a number of patterns, e.g. increasing the water storage in the watersheds and reducing runoff to downstreams (Chang et al. 2017). Experimental studies like Feng et al. (2016) and Jaramillo et al. (2018) deduced that the changing of the leaf stomatal conductance, surface roughness and leaf area index and biomass of trees under the greening circumstance could substantially influence ET_a.

Some studies, however, have come to the opposite conclusion. For example, Feng et al. (2020) indicated that precipitation plays a dominant role in ET_a variability in most arid regions. This may due to the large variability of precipitation and the small variability of vegetation in these regions. Particularly, in arid regions such as Xinjiang Province in Northwest China and the Loess Plateau, the large precipitation amplitude is the main controlling factor affecting ET_a. He et al. (2019) found that the variations of precipitation contributed the largest proportion of the ET_a changes in the Loess Plateau, where significant fluctuations of annual precipitation were observed in the study area. However, a recent study by Ning et al. (2017) is basically consistent with the conclusion of our results, indicating that vegetation, rather than precipitation, is the leading factor of the variations of ET_a via a comparative study in 13 watersheds in arid areas of Northern China. Therefore, we deduce that it is of great necessity to incorporate the variations of climatic factors, e.g. precipitation, and vegetation coverage in arid regions into the attribution analysis of the ET_a.

Two reservoir watersheds in mountainous regions of the DRB are selected as case studies. Mountain regions are experiencing more pronounced climate change than other global land areas. Xu et al. (2020) indicate that vegetation greening and/or earlier growing season are imposing positive effects on evapotranspiration in mountainous areas. In this study, climate and vegetation factors affecting actual evapotranspiration in two major reservoirs of the DRB were analyzed, and attribution analysis was conducted using the Budyko framework, which has guiding significance for water resources regulation and management, agricultural production and ecological construction in the basin. Similar results are found in the mountainous regions of the North China Plain, one of the most severe water scarce regions in the world. The results indicated that the improved vegetation conditions could positively contribute to the ET_a changes. For example, Liu et al. (2019a) found that crop coefficients (K_c) and leaf area index (LAI) were positively correlated with ET_a in the North China Plain. Wu & Xing (2016) pointed out that the change of ET_a in Haihe River Basin was positively correlated with the change of vegetation cover area. Cao et al. (2018) found that vegetation change had a promoting effect on ET_a, the increase of forest area enhanced the evapotranspiration capacity of the basin in the Chaohe Basin of the North China Plain. In other mountainous areas, e.g. the Loess Plateau, climate change has been reported to have significant impacts on the hydrological processes. Nevertheless, changes of vegetation could also exert great impacts on the variations of ET_a. For instance, Xu et al. (2018) quantified the sensitivity of climate and vegetation growth changes to ET_a in the upper reaches of the Yellow River Basin, revealing that both the regional averaged NDVI and ET_a showed an increasing trend, where positive relationships were detected between these two variables. Ma et al. (2019) reported that under the circumstance of ecological restoration in the Loess Plateau, recovering vegetation has affected the ET_a in the past 16 years significantly.

There are certain uncertainties in our study, which may come from the uncertainty of methods, data and parameters. We use the total differential method for attribution analysis under the Budyko framework. The partial derivative coefficients in the total differential method should be calculated by the complete Taylor expansion method so as to obtain accurate results. However, in our study, the partial derivative coefficient is estimated as a first-order approximation. Yang et al. (2014) showed that ignoring the higher-order coefficients of the Taylor expansion would lead to an underestimation of the contribution of climate change to hydrological processes as precipitation increases. And since the R² of the semi-empirical formula is less than 0.5, the change of vegetation coverage M and climate seasonal index S cannot fully explain the interannual variation of parameter ⁠. In this study, we used MODIS datasets to analyze the actual evapotranspiration. Results showed that ET_a data of MODIS is slightly higher than the ground-based observed data. This may lead to a higher value for the study of the contribution of vegetation coverage to ET_a. Since the difference between MODIS-ET_a and ground-based observed ET_a is very small and the linear fitting result is very good, the error between MODIS-ET_a and ground-based observed ET_a can be ignored. The relationships of parameter with vegetation coverage and climate seasonality index in the two watersheds were not significant in this study (Figure 9). This implies that 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 (Ning et al. 2017).

This study focused on the impacts of climatic factors and vegetation coverage on ET_a, particularly in arid mountainous regions. Nevertheless, there may be other factors affecting ET_a, such as human activities, irrigation and water use, which may not be quantified by the Budyko-based framework in one study. Liu et al. (2019b) investigated the attribution of evapotranspiration changes in the west area of China and found that human activities (e.g. air pollution, urbanization, and water supply project) may be one factor to influence ET_a. Sun & Ren (2013) found that ET_a from irrigated farmland increased in Haihe River Basin, and land use played the dominant role in the ET_a changes. Djaman et al. (2016) pointed out that crop water use efficiency is related to ET_a and irrigation water demand, and ET_a is positively correlated with irrigation. Jung et al. (2010) found that the decrease of the annual ET_a from 1998 to 2008 resulted from the decrease of the soil moisture in many regions of the world. Further studies are warranted in quantifying the impacts of other factors, e.g. human activities, construction and regulation of water supply infrastructures, on the variations of ET_a in these regions.

This study, using a Budyko-based framework, explored the impacts of climate and vegetation on actual evapotranspiration typical arid mountainous regions. This study quantified the climatic and vegetation indices, validated the accuracy of the estimated ET_a, and researched the attribution of ET_a variations to climatic and vegetation factors. The conclusions can be summarized as follows:

• (1)

  The increase of vegetation coverage led to the increase of ET_a in arid mountainous river basins such as the upper reservoir watersheds of the Daqing River Basin.

• (2)

  In arid mountainous regions, vegetation coverage may be the main impact factor relative to climate factors, especially for regions with large change of vegetation coverage but little change of climate factors. The contributions of vegetation were 42.15 and 58.56 mm, and that of climate factor were –11.48 and –13.47 mm, in the WRW and XRW, respectively.

The actual evapotranspiration ET_a in the study area showed an increasing trend from 2000 to 2015, and the vegetation cover also showed a green trend. The increase of vegetation led to the increase of evapotranspiration, which had a negative impact on the water resources and water supply in the downstream. Therefore, the extent of vegetation restoration in the upstream needed to be limited. Moreover, a proxy strategy is to choose crops that use less water and are suitable for growing in arid regions. Findings of this study tend to be a reference in exploring the complex climate-vegetation-hydrology relationships in arid mountainous regions.

This work was jointly supported by the National Key Research and Development Program of China (Grant No. 2018YFC0406506 and 2018YFC0406501), the National Natural Science Foundation of China (No. 51909275 and 51779271) and Guangxi Graduate Education Innovation Project (YCBZ2018023).

All relevant data are included in the paper or its Supplementary Information.