Attributing Evapotranspiration Changes with an Extended Budyko Framework Considering Glacier Changes in a Cryospheric-Dominated Watershed

Abstract

The retreat of glaciers has altered hydrological processes in cryospheric regions and affects water resources at the basin scale. It is necessary to elucidate the contributions of environmental changes to evapotranspiration (ET) variation in cryospheric-dominated regions. Considering the upper reach of the Shule River Basin as a typical cryospheric-dominated watershed, an extended Budyko framework addressing glacier change was constructed and applied to investigate the sensitivity and contribution of changes in environmental variables to ET variation. The annual ET showed a significant upward trend of 1.158 mm yr⁻¹ during 1982–2015 in the study area. ET was found to be the most sensitive to precipitation (P), followed by the controlling parameter (w), which reflects the integrated effects of landscape alterations, potential evapotranspiration (ET₀), and glacier change (∆W). The increase in P was the dominant factor influencing the increase in ET, with a contribution of 112.64%, while the decrease in w largely offset its effect. The contributions of P and ET₀ to ET change decreased, whereas that of w increased when considering glaciers using the extended Budyko framework. The change in glaciers played a clear role in ET change and hydrological processes, which cannot be ignored in cryospheric watersheds. These findings are helpful for better understanding changes in water resources in cryospheric regions.

1. Introduction

Catchment hydrology is fundamental to the study of interactions among climate, soil, vegetation, and terrain [1,2,3]. A full understanding of hydrological processes at the catchment scale is essential for water resource management. Global hydrological processes are undergoing significant change in response to climate change. Evapotranspiration (ET), as a significant part of hydrological processes, plays a crucial role in the water and energy exchanges in the atmosphere, hydrology, cryosphere, and biosphere [4,5]. Therefore, it is important to investigate the potential mechanisms behind ET variation to understand hydrological processes. However, quantifying the mechanisms behind ET variation remains a challenge, owing to limited long-term observations and the complexity of the processes between climate, soil, and vegetation.

There are several methods to quantify the impact of environmental changes on hydrological processes: (a) hydrological models, including physical-based and conceptual models [6,7], (b) the paired-basin method [8], and (c) statistical methods [9]. Physical-based hydrological models are powerful tools for evaluating the contributions of climate change and landscape alterations on hydrological regimes and water partitioning [6]. However, they require significant prior knowledge of soil and vegetation parameters [10], which are difficult to obtain in data-scarce regions, especially in cryospheric regions with limited observations and harsh environments. The paired-basin method, based on the similarity principle, is suitable for small catchments but not for medium or large catchments [11]. Statistical methods have been extensively applied due to their low computing requirements and high adaptability [12], among which the Budyko hypothesis [13,14] has been popular in recent years [15,16]. In addition, the Budyko hypothesis coupled with water-energy balance information provides an effective way to elucidate the contribution of each factor to hydrological processes [9,17,18,19,20,21], making the approach more feasible and robust.

The Budyko equation describes water balance by dividing precipitation (P) into runoff (R) and ET at a catchment scale [14]. Budyko [14] hypothesized that the annual ET is controlled by atmospheric demand (potential evapotranspiration, ET0) and water supply (P). Owing to the mismatch between observations and the Budyko curve, a controlling parameter was considered to reflect the impact of catchment characteristics on water partitioning [22,23,24,25]. The controlling parameter represents the combined effects of all factors other than the aridity index (ET₀/P), which is related to soil, vegetation, land use, topography, and human activity [26]. Thus, the controlling parameter was used to differentiate the contributions of climatic and anthropogenic factors to water partitioning.

Previous studies have quantified the influences of climate change and landscape factors on ET change using the Budyko framework [5,21,27,28,29]. For example, Yang et al. [21] proposed that the positive contributions of P and the normalized difference vegetation index (NDVI) offset the negative contribution of ET₀ to changes in ET in the Qilian Mountains. Ning et al. [5] suggested that ET variation in arid alpine basins was mainly controlled by rainfall variation during the growing season. Li et al. [30] found that ET variation was controlled by P in water-limited regions and dominated by the vapor pressure deficit in energy-limited regions in China, based on a modified Budyko framework. Furthermore, Li et al. [27] evaluated the relative contributions of environmental factors to global ET. These studies analyzed the relative contribution of each factor to ET in many regions using a series of Budyko equations. However, the application of the Budyko framework was limited to cryospheric regions, where glacier shrinking and permafrost degradation seriously affect water balance and hydrological processes.

Various studies have quantified the impact of cryospheric elements on runoff [3,31,32,33,34,35]. For example, Zhang et al. [31] indicated that runoff is sensitive to snow change in northern mountainous and high-altitude areas. Wang et al. [36] analyzed the response of streamflow to permafrost degradation in the source region of the Yellow River. Saydi et al. [37] and Liu et al. [35] reported that the effect of glaciers on runoff should be considered in the Budyko framework for cryospheric catchments. However, few studies have investigated the elasticity and contribution of glaciers to ET in cryospheric regions.

Glaciers, as natural solid reservoirs, significantly contribute to water supply in many catchments worldwide [38]. Global warming causes glacier retreat [39], which further changes hydrological processes in cryospheric regions [40]. The Tibetan Plateau (TP), known as the “Asian tower”, has approximately 36,763 glaciers with an area of approximately 49,873.44 km² [41]. Owing to limited observations and complicated cryospheric processes, the responses of hydrological elements (i.e., R and ET) to climate change and glacier change are still unclear. Therefore, it is crucial to investigate the response of ET to climate change and glacier change in cryospheric regions, which will be helpful for understanding the potential driving mechanisms of ET variation.

The main purpose of this research was to investigate the potential mechanisms of climate change, glacier change, and landscape alterations in response to ET change in a typical cryospheric-dominated watershed. First, the spatiotemporal variations in ET and environmental factors during 1982–2015 were evaluated. Second, the elasticity coefficients of ET for the climatic factors, glacier change, and landscape alterations were assessed. Third, the contribution of all environmental factors to ET change were quantitatively analyzed. This study provides a foundation for better understanding the response of hydrological elements in cryospheric regions to glacier change, which is of great significance for water resource management.

2. Materials and Methods

2.1. Study Area

The upper reach of the Shule River Basin (URSRB) is located on the northeastern Tibetan Plateau (Figure 1), where the famous Hexi Corridor is also located [42,43]. Throughout, the Shule River in the basin flows from southeast to northwest; its total length is approximately 670 km [41]. Usually, the river section upstream of the Changmabao hydrological station is the URSRB, which covers an area of 1.1 × 10⁴ km². It is a typical cryospheric watershed fed by glaciers and snowmelt. There were 486 glaciers in the study area during 2006–2010 [44]. The glaciers cover an area of 403.5 km² and account for 3.7% of the total area based on data extracted from the Second Chinese Glacier Inventory [44,45]. Permafrost covers about 83% of the total area [46]. The elevation ranges from 2100 m to 5637 m (Figure 1).

2.2. Data Collection

Annual runoff data at the Changmabao Hydrologic Station were obtained from the Hydrological Bureau. Due to limited meteorological observations in the URSRB, P and ET₀ were extracted from ERA5-Land reanalysis data. The ERA5-Land dataset has a higher horizontal resolution of 9 km compared to ERA5 [47] and its performance has been verified in previous studies [48,49,50]. Glacier data were obtained from the Cold and Arid Region Science Data Center (http://www.crensed.ac.cn/portal/ (accessed on 3 September 2022)) and the National Tibetan Plateau Scientific Data Center (http://data.tpdc.ac.cn/zhhans/ (accessed on 26 August 2022). NDVI were extracted from the Global Inventory Modelling and Mapping Studies (GIMMS) dataset (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/ (accessed on 12 October 2021) with a spatial resolution of 8 km, which provides an index representing vegetation dynamics. The collection period in this study was from 1982 to 2015.

2.3. Methods

2.3.1. Water Balance and Budyko Framework

The water balance in the watershed is written as follows:

$$\mathrm{ET}=\mathrm{P}-\mathrm{R}-\mathsf{\Delta }\mathrm{S}$$

(1)

where ET, P, R, and ∆S are the evapotranspiration (mm), precipitation (mm), runoff (mm), and terrestrial water storage change (mm), respectively. For long-term periods, ∆S can be ignored because of the small anomalies [20,35,51] on an annual scale.

In the cryospheric-dominated watershed, the water change caused by glacier mass balance cannot be ignored; it represents the effect of climate change on glacier volume and the fluctuation of glacier meltwater. The average annual glacier mass balance in the URSRB ranged from −595.9 to 222.8 mm w.e. (mm water equivalent), with an average value of −116.1 mm w.e. during 1970–2015 (Figure S1). It decreased significantly by −5.64 mm w.e. yr⁻¹. Therefore, the glacier mass balance should be considered in the water balance of glaciated watersheds. Thus, the water balance for the glaciated watershed was modified as follows:

$$\begin{array}{l}\mathrm{ET}=(\mathrm{P}-\mathsf{\Delta }\mathrm{W})-\mathrm{R}\\ \mathsf{\Delta }\mathrm{W}=mb\times f\end{array}$$

(2)

where ∆W is the glacier mass balance with equivalent water height (mm), which can be calculated from the glacial mass balance (mb; mm w.e.) and the glacier fraction (f; %); mb can be calculated using the degree-day factor and positive accumulated temperature [45,52].

The Budyko framework was used for ET estimation to couple the concepts of water and energy balance. This framework divides the steady-state water balance as a function of the relative importance of atmospheric water supply (P), water demand (ET₀), and the controlling parameter (w). There are many classical equations for describing the Budyko framework, such as the Fu equation [23,53], Zhang equation [24], Choudhury–Yang equation [22,54], and the Wang–Tang equation [25]. These equations produce similar performances in ET estimation [20]. Therefore, the Fu equation [53] was employed to estimate ET.

$$\frac{\mathrm{ET}}{\mathrm{P}}=1+\frac{{\mathrm{ET}}_0}{\mathrm{P}}-{[1+{(\frac{{\mathrm{ET}}_0}{\mathrm{P}})}^w]}^{\frac{1}{w}}$$

(3)

where P and ET₀ can be extracted from ERA5-Land data and w is the controlling parameter, representing the combined effects of the landscape characteristics [55].

2.3.2. Extended Budyko Framework

By considering ∆W, the Budyko equation can be written as [36]:

$$\frac{\mathrm{ET}}{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}=1+\frac{{\mathrm{ET}}_0}{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}-{[1+{(\frac{{\mathrm{ET}}_0}{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})})}^w]}^{\frac{1}{w}}$$

(4)

By eliminating the denominator, Equation (4) can be written as:

$$\mathrm{ET}=(\mathrm{P}-\mathsf{\Delta }\mathrm{W})+{\mathrm{ET}}_0-{[{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w]}^{\frac{1}{w}}$$

(5)

2.3.3. Sensitivity Analysis

The elasticity coefficient method was used to evaluate the sensitivity of ET to environmental factors. The elasticity coefficient of ET for a particular variable can be expressed as [56]

$${\mathsf{\epsilon }}_{x_i}=\frac{\partial \mathrm{ET}}{\partial x_i}\times \frac{x_i}{\mathrm{ET}}$$

(6)

where ε_xi means the elasticity coefficient of ET for x_i, where x_i represents each environmental factor (P, ET₀, ∆W, or w).

Based on Equation (5), the first-order partial derivative coefficients for P, ET₀, ∆W, and w can be written as follows:

$$\frac{\partial \mathrm{ET}}{\partial \mathrm{P}}=1-{[{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w]}^{\frac{1}{w}-1}{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^{w-1}$$

(7)

$$\frac{\partial \mathrm{ET}}{\partial {\mathrm{ET}}_0}=1-{[{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w]}^{\frac{1}{w}-1}{\mathrm{ET}}_0{}^{w-1}$$

(8)

$$\frac{\partial \mathrm{ET}}{\partial \mathsf{\Delta }\mathrm{W}}=-1+{[{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w]}^{\frac{1}{w}-1}{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^{w-1}$$

(9)

$$\begin{array}{ll}\frac{\partial \mathrm{ET}}{\partial w}=& -[\frac{{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w\ln ({(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w)+{\mathrm{ET}}_0{}^w\ln ({\mathrm{ET}}_0)}{w({(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w)}-\frac{\ln ({(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w)}{w^2}]\times \\ & {[{(\mathrm{P}-\mathsf{\Delta }\mathrm{W})}^w+{\mathrm{ET}}_0{}^w]}^{\frac{1}{w}}\end{array}$$

(10)

A positive (negative) elasticity coefficient means that ET increases (decreases) with the variable. Based on the above equations, the elasticity coefficients of ET for P, ET₀, ∆W, and w can be obtained and the changes in the elasticity coefficients can be derived.

2.3.4. Contribution Analysis of ET Change

A differential equation method was adopted to quantify the contribution to ET variation. This method assumes that all the variables are independent of each other. Thus, the contribution for each variable to ET variation can be calculated as follows:

$$\frac{d\mathrm{ET}}{dt}=\frac{\partial \mathrm{ET}}{\partial {\mathrm{ET}}_0}\frac{d{\mathrm{ET}}_0}{dt}+\frac{\partial \mathrm{ET}}{\partial \mathrm{P}}\frac{d\mathrm{P}}{dt}+\frac{\partial \mathrm{ET}}{\partial \mathsf{\Delta }\mathrm{W}}\frac{d\mathsf{\Delta }\mathrm{W}}{dt}+\frac{\partial \mathrm{ET}}{\partial w}\frac{dw}{dt}+\delta$$

(11)

The equation is simplified as:

$$\mathrm{L}(\mathrm{ET})=\mathrm{C}(\mathrm{P})+\mathrm{C}({\mathrm{ET}}_0)+\mathrm{C}(\mathsf{\Delta }\mathrm{W})+\mathrm{C}(w)+\delta$$

(12)

where L(ET) is the slope of ET change; C(P), C(ET₀), C(∆W), and C(w) are the contributions of change in P, ET₀, ∆W, and w to change in ET, respectively; and δ is the systemic error.

The relative contribution of each variable to ET change can be expressed as:

$$\mathrm{RC}(x_i)=\frac{\mathrm{C}(x_i)}{\mathrm{C}(\mathrm{P})+\mathrm{C}({\mathrm{ET}}_0)+\mathrm{C}(\mathsf{\Delta }\mathrm{W})+\mathrm{C}(w)}$$

(13)

where x_i means one of the four variables and RC(x_i) represents the relative contribution of x_i to the ET change.

3. Results

3.1. Variations of ET, P, ET₀, and ∆W

To better understand the influences of climate change, glacier change, and landscape alteration on ET change, the temporal variations of ET and environmental factors were analyzed (Figure 2). Annual ET and P showed significant increasing trends, whereas the increasing trend in ET₀ was insignificant. The average values of annual ET, P, and ET₀ were 394.2, 488.3, and 719.8 mm during 1982–2015, with increases of 1.158, 2.472, and 0.716 mm yr⁻¹, respectively. Moreover, ∆W showed a significant decreasing trend with a slope of −0.392 mm yr⁻¹, with an average value of −4.78 mm.

3.2. Sensitivity of ET to Environmental Variables

The elasticity coefficients of ET with respect to precipitation, potential evapotranspiration, glacier and landscape are shown in Figure 3. The elasticity coefficient ranged from 0.57 to 0.75 for P, from 0.25 to 0.42 for ET0, from −0.01 to 0.03 for ∆W, and from 0.31 to 0.7 for w. These values indicate that an increase of 1% in P, ET₀, ∆W, or w would result in a 0.57~0.75% increase, 0.25~0.42% increase, −0.01 ~0.03% increase, or 0.31~0.7% increase in ET, respectively.

On average, an increase of 1% in P, ET₀, ∆W, or w increased ET by 0.678, 0.325, 0.007, and 0.471%, respectively (Figure 4). The elasticity coefficient of ET for P was the largest, followed by the controlling parameter w, indicating that ET was the most sensitive to climate change.

The largest absolute value of the four elasticity coefficients was ε_P, whereas the smallest absolute value of the four elasticity coefficients was ε_∆W. The absolute values of ε_ET0 and ε_w increased significantly, whereas that of ε_P decreased significantly. In addition, the absolute value of ε_∆W increased insignificantly. This suggests that ET became more sensitive to variations in landscape and glacier changes, whereas the sensitivity of ET to climate change decreased.

3.3. Attribution Analysis of ET Change

The calculated changes in ET (1.202) were similar to that observed for the ET trend (1.158), indicating that the method used in this study was appropriate for evaluating the contribution of environmental factors to ET variation (Figure 5a). The increased P and ET₀ caused an increase in ET of 1.353 mm and 0.127 mm, respectively; the corresponding contribution rates were 112.64% and 10.55%, respectively (Figure 5b). The decrease in ∆W resulted in an increase in ET of 0.215 mm, with a contribution rate of 17.86%. The decrease in w caused a decrease in ET of −0.493 mm, with a contribution rate of −41.05%. Therefore, climate change increased ET, whereas landscape alteration decreased ET, and the contribution of climate change was greater than that of the latter. These results suggest that changes in P were the dominant factor in the increase in ET, followed by changes in landscape. Changes in ∆W had an impact on ET change, while changes in ET₀ played a limited role in the change in ET. The combined effects of the positive contribution of P and the negative contribution of w resulted in an increase in ET.

4. Discussion

4.1. Performances of the Budyko Framework Considering Glaciers

The original Budyko framework assumes that the potential water supply for ET is from P. However, in areas covered by glaciers, the water supply originates from P and from glacier melting. To better interpret the influence of glaciers on ET, the original Budyko framework, without considering glaciers, was also analyzed (Figure 6). Parameter w in the original Budyko framework ranged from 2.05 to 3.39, with an optimal value of 2.518. The estimated w in this study ranged from 2.15 to 3.32, with an optimal value of 2.548. A larger optimal w value appears when considering meltwater from snow and glaciers [34].

In addition, the contribution rates of environmental factors to ET variation were different between the original and extended Budyko frameworks (Figure 7). The contribution rates of P and ET₀ in the original Budyko framework were higher than those in the present study. The absolute value of the contribution rate of w was higher than that observed in this study, which considered the effect of glaciers, indicating that there were other factors that influenced the ET change.

4.2. Relationship between the Parameter w and NDVI

Parameter w is related to landscape conditions, including the topography, soil properties, and vegetation conditions [55,57]. The soil properties and topography of the study area have hardly changed in recent decades, whereas the vegetation has varied with time. The increasing NDVI resulted in an increase in transpiration and, thus, an increase in ET (Figure 8a). The decrease in w was negatively correlated with the increase in NDVI (Figure 8b), indicating that ET increased with decrease in w and the change in w negatively contributed to the change in ET. This further indicates that w might be influenced by factors other than NDVI. For example, Ning et al. [20] suggested that climate seasonality influences the ET. Liu et al. [35] reported that w reflected the combined results of NDVI, terrestrial slope, and glacier fraction. Future studies should focus on exploring the relationships between w and other related factors (i.e., topography, permafrost degradation, and climate seasonality).

4.3. Uncertainty Analysis

The Budyko framework assumes that all factors should be independent of each other, whereas some factors may exhibit autocorrelation [58], which leads to uncertainties in the contributions of environmental factors to ET change. Moreover, the interactions of related factors may have an indirect effect on ET, but the framework cannot quantify the indirect effects and needs to be further improved. The lack of long-term observations of ET limits the verification of model accuracy. Moreover, P, ET₀ and NDVI were extracted from the existing data products (ERA5-Land and GIMMS), and the mismatches of these products may introduce uncertainties in the estimation of ET based on the water balance equation. In addition, terrestrial water storage has been ignored over long periods, which may introduce uncertainties in water-balance estimation [59].

In cryospheric regions, permafrost degradation directly influences hydrological processes [60], potentially affecting the redistribution of water resources. Previous studies have proposed that permafrost thawing and lengthening of non-freezing days may lead to a significant increase in ET [61,62]. However, the influence of permafrost (i.e., increasing ground temperature, deepening of the active layer thickness, and thawing of ground ice) on hydrological processes is difficult to determine directly. For example, the thawing of ice consumes energy, which was neglected in the ET estimation [63]. In addition, the ground ice content is difficult to quantify and the results estimated by different algorithms vary greatly [64,65]. Moreover, the spatial distribution of ground ice content is yet to be clearly investigated, which makes assessment of the impact of ground ice content on hydrological processes challenging [66]. Therefore, the role of permafrost is not considered in this study. Further studies should consider the impacts of permafrost freezing and thawing on hydrological processes.

5. Conclusions

Hydrological processes have been strongly affected by climate change and landscape alteration, especially in cryospheric-dominated regions. This study evaluated the long-term change in ET from 1982 to 2015 in the upper reach of the Shule River Basin using an extended Budyko framework involving consideration of glaciers. The responses of ET to the direct impact of climate change, the impact of glacier change, and landscape alteration were then also determined.

Annual ET, P, and ET₀ exhibited increasing trends, whereas ∆W exhibited a significant decreasing trend. The elasticity coefficients were highest for P, followed by w, ET₀, and ∆W. In addition, the elasticity coefficients varied with time. The increase in ET was dominated by changes in P, accounting for 112.64% of the total ET variation. The contributions of w, ∆W, and ET₀ change accounted for −41.05%, 17.86%, and 10.55% of ET change, respectively. The positive contribution of P largely offset the negative contribution of w. Compared with the original Budyko framework, the contributions of P and ET₀ to ET change decreased, while that of w increased when considering glaciers in the extended Budyko framework. This indicates that glaciers had a significant impact on annual ET which cannot be ignored.

These results provide insights for understanding the response of ET to climate change and landscape alteration in cryospheric-dominated watersheds, which can also serve as a reference for other cryospheric regions. Nevertheless, further studies should consider the effects of permafrost on hydrological processes in cryospheric regions.