Projected future changes in the cryosphere and hydrology of a mountainous catchment in the upper Heihe River, China

Climate warming exacerbates the degradation of the mountain cryosphere, including glacier retreat, permafrost degradation, and snow cover reduction. These changes dramatically alter the local and downstream hydrological regime, posing significant threats to basin-scale water resource management and sustainable development. However, this issue is still not adequately addressed, particularly in mountainous catchments. We developed an integrated cryospheric–hydrologic model, the FLEX-Cryo model, to comprehensively consider glaciers, snow cover, and frozen soil and their dynamic impacts on hydrological processes. Taking the mountainous Hulu catchment located in the upper Heihe River of China as a case study, we utilized the state-of-the-art climate change projection data under two scenarios (SSP2-4.5 and SSP5-8.5) from the sixth phase of the Coupled Model Intercomparison Project (CMIP6) to simulate the future changes in the mountainous cryosphere and their impacts on hydrology. Our findings showed that under the medium-emission scenario (SSP2-4.5) and high-emission scenario (SSP5-8.5), by the end of the 21st century, the glacier will completely melt out around the years 2051 and 2045, respectively. The annual maximum snow water equivalent is projected to decrease by 41.4 % and 46.0 %, while the duration of snow cover will be reduced by approximately 45 and 70 d. The freeze onset of seasonally frozen soil is expected to be delayed by 10 and 22 d, while the thaw onset of permafrost is likely to advance by 19 and 32 d. Moreover, the maximum freeze depth of seasonally frozen soil is projected to decrease by 5.2 and 10.9 cm per decade, and the depth of the active layer will increase by 8.2 and 15.5 cm per decade. Regarding hydrology, catchment total runoff exhibits a decreasing trend, and the tipping point of glacier runoff occurs approximately between 2019 and 2021. Permafrost degradation will likely reduce the duration of low runoff in the early thawing season; the discontinuous baseflow recession gradually transitions into linear recessions, and the baseflow increases. Our results highlight the significant changes expected in the mountainous cryosphere and hydrology in the future. These findings enhance our understanding of cold-region hydrological processes and have the potential to assist local and downstream water resource management in addressing the challenges posed by climate change.

1 Introduction

“How will cold-region runoff and groundwater change in a warmer climate?” This question was identified by the International Association of Hydrological Sciences (IAHS) as one of the 23 unsolved scientific problems (Blöschl et al., 2019). The mountain cryosphere, which includes glaciers, snow cover, and frozen soil in high-altitude regions, has a significant impact on water resources (Adler et al., 2019; Arendt et al., 2020; Rasul et al., 2020; Zhang et al., 2022). The mountain cryosphere is considered a crucial “water tower” and a climate change indicator due to its sensitivity to climate change (Tang et al., 2023). However, the cryosphere is rapidly retreating in many parts of the world, including glacier retreat, expansion of glacier lakes, northward movement of the permafrost southern limit, and shrinkage of the snow cover area (Moreno et al., 2022; S. Wang et al., 2022; Ding et al., 2019; Wang et al., 2023). These changes have disrupted the water tower region and pose significant challenges to sustainable water resources management (Ragettli et al., 2016; Yao et al., 2022).

The degradation of the mountain cryosphere varies from region to region (Andrianaki et al., 2019; Wang et al., 2019). Lower altitudes experience a decreasing trend in snow cover days, snow depth, snow water equivalent, and snowmelt due to climate warming, while higher altitudes present a more complex picture (Connon et al., 2021; Nury et al., 2022; Yang et al., 2022). Global continental glacier mass balance from 2006 to 2015 was approximately −123 ± 24 Gt yr⁻¹, with significant losses observed in the Southern Andes, Caucasus Mountains, and central Europe, while the Karakoram and Pamir regions exhibited less loss (Intergovernmental Panel on Climate Change (IPCC), 2022; Van Der Geest and Van Den Berg, 2021). Future projections suggest a 40 % decrease in global permafrost by the end of the century, potentially transitioning into seasonally frozen soil (Chadburn et al., 2017; Martin et al., 2023). The mountain cryosphere serves as a significant freshwater reservoir, impacting water resources and the hydrological cycle (Ding et al., 2020).

In a warming climate, glacier runoff exhibits a “tipping point” characterized by an initial increase followed by a subsequent decline (Rosier et al., 2021; Zhang et al., 2012; Zhao et al., 2019). While small glaciers have already experienced this tipping point, its occurrence in large glaciers remains uncertain (Brovkin et al., 2021; Huss and Hock, 2018). Permafrost degradation leads to an increase in active layer thickness, resulting in the melting of subsurface ice and an augmentation of soil water storage capacity (Abdelhamed et al., 2022). Additionally, the degradation of the cryosphere significantly impacts the atmosphere, biosphere, surface energy balance, ecological water use, and ecosystems (Gilg et al., 2012; Miner et al., 2022; Pothula and Adams, 2022). Understanding the complex interactions between cryosphere degradation and ecosystems is crucial, but quantitatively observing the degradation process in high-altitude regions is challenging. Hydrological models provide an effective approach to analyse degradation patterns and assess the impact on future water resources (Han and Menzel, 2022).

Glacio-hydrology is influenced by both glacier melt and glacier dynamics. Glacier melting models can be categorized into three types: energy balance, temperature index, and hybrid models (He et al., 2021; Gao et al., 2021; Negi et al., 2022; Zekollari et al., 2022). While energy balance models analyse glacier accumulation and melt processes based on solid physical mechanisms, they require extensive forcing data that may not be readily available in mountainous regions (Huss et al., 2010). In contrast, temperature index models are simpler and more effective, requiring fewer parameters (including degree-day factor and threshold temperature) and forcing data (temperature and precipitation) (Bolibar et al., 2022; Vincent and Thibert, 2023). These models perform well at both daily and monthly scales. Glaciers move slowly due to the combined effects of gravity and high viscosity of ice. Due to climate change, ice becomes thinner, and glacier loses its mass balance, which will cause the glacier morphology to evolve to a new balance status. Glacier dynamic models, with the full-Stokes approach as the most complete form, and many other simplifications, such as the shallow-ice approximation and the shallow-shelf approximation, are still computationally expensive, hindering their implications in large-scale studies. Three conceptual models are commonly used for glacier evolution: volume–area scaling (V-A) method, accumulation area ratio (AAR) method, and Δh parameterization (Michel et al., 2022; Wiersma et al., 2022). The first two approaches do not consider the detailed changes in different elevation bands, while the Δh-parameterization approaches only require glacier mass balance as forcing data to analyse changes in ice thickness at different elevation bands based on the relationship between glacier mass balance and glacier area (Huss et al., 2010). The temperature index method coupled with the Δh-parameterization approach serves as an effective module to simulate glacier evolution and its impacts on hydrology.

Permafrost hydrology models can be classified into one-dimensional models and distributed watershed models (Elshamy et al., 2020). One-dimensional hydrological models, such as the Stefan equation, the temperature at the top of permafrost (TTOP) model, CoupModel, and the SHAW model, are effective in simulating freeze depth, hydrothermal transport, and carbon or nitrogen transport, but they are unable to capture the broader impact of permafrost on hydrology at catchment scale (Kaplan Pastíriková et al., 2023; Li et al., 2022; Liu et al., 2023). On the other hand, distributed watershed models, such as the Cold Regions Hydrological Model (CRHM), the HydroGeoSphere (HGS), and the distributed water-heat coupled model (DWHC), consider the spatial variability of permafrost properties and simulate the interactions between permafrost, surface water, and groundwater (Chen et al., 2008; He et al., 2023; Pomeroy et al., 2022). These models operate on a small-scale basis and require extensive prior knowledge, following a “bottom-up” approach that relies on small-scale field observations and situational models to comprehend the effects of permafrost on hydrology (Peng et al., 2016). However, the freeze–thaw cycle is influenced by multiple interconnected factors, including climate, topography, slope orientation, snowpack, and vegetation (Chang et al., 2022). The process of upscaling would lead to the neglect of some variables and the amplification of others (Fenicia and McDonnell, 2022). In contrast, the FLEX-Cryo model is based on the FLEX-Topo-FS model, which employs a “top-down” modelling procedure that involves observed data analysis, qualitative perceptual modelling, quantitative conceptual modelling, and the testing of model realism. This model exhibits the ability to accurately and expeditiously identify key elements in permafrost hydrological processes and then simulate hydrology at the catchment scale (Beven, 2012; Gao et al., 2022).

The aim of this study is to integrate the FLEX-Topo-FS model and a glacier evolution model (Δh parameterization) to develop a landscape-based model of the mountain cryosphere, referred to as FLEX-Cryo. This model will be utilized to simulate changes in various components of the mountain cryosphere and evaluate their impacts on hydrological processes, thereby enhancing our understanding of the hydrological cycle. The model will be driven by eight bias-corrected global climate models (GCMs) under SSP2-4.5 and SSP5-8.5 scenarios obtained from the Coupled Model Intercomparison Project Phase 6 (CMIP6), which will be used to predict future changes in glaciers, snow, and frozen soil, as well as their effects on hydrology.

2 Study area and data

2.1 Study area

The Hulu catchment is located in the upper reaches of Heihe River basin (38°12^′–38°17^′ N, 99°50^′–99°53^′ E) and about 23.1 km². The elevation ranges from 2960–4820 m. The Hulu catchment belongs to a continental monsoon climate. Rainfall is the major phase of precipitation, and there is also snowfall in the winter. Four landscapes are identified, i.e. glacier (5.6 %), alpine desert (53.5 %), vegetation hillslope (37.5 %), and riparian zone (3.4 %; Fig. 1c). The landscape pattern in the Hulu catchment has typical altitude zonality. Most of the vegetation and riparian zone is distributed in the lower-elevation bands. Alpine desert and glacier landscapes are in the high-elevation bands. There is almost no human activity in the catchment (Liu and Chen, 2016; Li et al., 2014). There are two major glaciers, i.e. Glacier 1 and Glacier 2 (Figs. 1b and 2) in the catchment. Glacier 1 was also named the Shiyi Glacier in the glacier catalogue of China. Seasonally frozen soil and permafrost both exist in the catchment. The lower limit of permafrost is around 3650–3700 m. The permafrost region accounts for 64 % of the total catchment, and the rest is seasonally frozen soil. The soil generally starts to freeze in October (Gao et al., 2019). Thus, 1 October was set as the start of the hydrological year. All the interannual variations in this study were based on the hydrology year.

Figure 1 (a) The distribution of permafrost and seasonally frozen soil on the upper Heihe River basin and the location of the Hulu catchment. (b) The digital elevation model and the thickness of the two major glaciers. (c) Spatial distribution of four landscapes (glacier, alpine desert, vegetation hillslope, and riparian zone).

Figure 2 Landscape classification at different elevation bands.

2.2 Data

Temperature and precipitation are observed at 2920 m, near the outlet of the catchment, from 2011 to 2014. Farinotti et al. (2019) used five models which used the ice flow dynamics to invert ice thickness from surface features to estimate the ice thickness distribution of about 21 500 glaciers outside the Greenland and Antarctic ice sheets. We used the estimated data for the initial thickness distribution of Glacier 1 and Glacier 2 (https://doi.org/10.3929/ethz-b-000315707, Farinotti, 2019; Farinotti et al., 2019).

The Couple Model Intercomparison Project phase 6 (CMIP6) is widely used to predict future climate. Eight general circulation models (GCMs) (Table 1) under two climate scenarios (SSP2-4.5 and SSP5-8.5) are used for predicting future climate. The selected models have been validated well in nearby catchments (Xing et al., 2023; Yin et al., 2021; Ma et al., 2022; Zhu and Yang, 2020; R. Chen et al., 2022). The SSP2-4.5 scenario represents a medium part of the future pathways, which is usually shown by a referenced experiment comparing others (CMIP6-Endorsed MIPs), and it produces a radiative forcing of 4.5 W m⁻² in 2100. The SSP5-8.5 scenario is a high-emission scenario, and it produces a radiative forcing of 8.5 W m⁻² in 2100.

Table 1 Details of data from eight GCMs used in this study.

Although the reliability of GCMs has been verified in the previous studies, there is certain bias in the output that needs to be corrected. Firstly, outputs from eight GCMs under two climate scenarios are interpolated to 0.5°×0.5°, and then the bias correction is carried out by CMhyd software (download from https://swat.tamu.edu/software/cmhyd/, last access: October 2021) in which four methods were used, including distribution mapping of precipitation and temperature, linear scaling of precipitation and temperature, variance scaling of temperature, and local intensity scaling (LOCI) of precipitation (Teutschbein and Seibert, 2012). The bias-corrected precipitation and temperature were calculated using the equal weighted average method to obtain the multi-model ensemble average values under the SSP2-4.5 and SSP5-8.5 scenarios, which reduce the uncertainty caused by a single bias correction method and a single GCM. The method is described as follows:

$$\begin{array}{}\text{(1)}& P_{\text{ave}}={\displaystyle \frac{\mathrm{1}}{N_{\text{GCM}}}}\left(\sum _{j=\mathrm{1}}^{N_{\text{GCM}}}\left({\displaystyle \frac{\mathrm{1}}{N_{\text{bias}}}}\left(\sum _{i=\mathrm{1}}^{N_{\text{bias}}}\left(P_i\right)\right)\right)\right),\end{array}$$

where P_ave is the average value of the multi-model and multi-method, P_i is the projected climate data of an GCM, N_bias is the number of correction methods (N_bias is 3 in this research), and N_GCM is the number of GCM (N_GCM is 8 in this research).

3 Methodology

3.1 FLEX-Cryo

The FLEX-Cryo model is a landscape-based cryospheric–hydrological model that considers multiple elements of cryosphere and their impacts on hydrology, including glaciers, snow, and frozen soil. Figure 3 shows the structure of the FLEX-Cryo model.

Figure 3 Structure of the FLEX-Cryo model. The abbreviations in red indicate parameters, and the abbreviations in black indicate storage components and fluxes.

The model parameters used in this research were obtained the optimal parameter set from a previous study conducted in the same catchment (Gao et al., 2022). The selected parameters are listed in Table 2, and the other variables used in calculations (Fig. 3) are listed in Table 4.

Table 2 Model parameters and their values in this study.

Table 3 The FLEX-Cryo model equations.

Table 4 The variables in Table 3 and Fig. 3 and their meaning.

3.1.1 Glacier melt and snowmelt

The threshold temperature (T_t) determines the phase of precipitation, i.e. snowfall or rainfall. The snow reservoir (S_w) accounts for the snow accumulation, melt (M_w), and water balance (Table 3, Eq. 9). The number of days when S_w is non-zero represents the snow cover days, and the maximum S_w is the maximum snow water equivalent of a year (Giovando and Niemann, 2022). Both glacier melt and snowmelt were calculated by the temperature index method, which is based on the degree-day factor (F_dd). If there is no snow cover, the glacier starts to melt. Due to the lower albedo, the degree-day factor of ice is greater than that of snow cover, and it is multiplied by the coefficient C_g to calculate glacier melt (Table 3, Eq. 14).

3.1.2 Rainfall–runoff module

The rainfall and snowmelt enter the root zone reservoir S_u, and then runoff (R_U) is generated based on the input water and the relative root zone soil moisture ($S_{\mathrm{u}}/S_{\text{umax}}$) and the shape of root zone storage capacity distribution determined by parameter β (Table 3, Eq. 16). Actual evaporation E_a is also estimated based on the soil moisture $S_{\mathrm{u}}/S_{\text{umax}}$ and the potential evaporation by the Hamon equation (Hamon, 1961). The generated runoff (R_U) is separated, by parameter D, into two linear reservoirs, i.e. the fast-response reservoir (S_f) and the slow-response reservoir (S_s) (Table 3, Eqs. 18 and 19). The two reservoirs are, respectively, controlled by the fast-recession parameter K_f and the slow-recession parameter K_s to simulate the subsurface storm flow Q_f and groundwater runoff Q_s (Table 3, Eqs. 7 and 8).

Different landscapes, for example, alpine desert, vegetation hillslope, and riparian zone, have different sizes of root zone storage capacity (S_umax) (Aubry-Wake and Pomeroy, 2023). In the vegetation hillslope, plants have well-developed root systems, and the root zone has a larger storage capacity. Therefore, S_{umax_V} was set to a larger value. For the alpine desert and riparian zone, S_{umax_D} and S_{umax_R} were both limited due to the less developed root system and storage capacity.

3.2 Frozen soil module

The Stefan equation was employed to estimate freeze (thaw) depth. This equation is calculated by the freeze (thaw) index (F), which neglects the sensible heat. The equation is as follows:

$$\begin{array}{}\text{(2)}& \mathit{\epsilon }={\left({\displaystyle \frac{\mathrm{2}\cdot \mathrm{86}\mathrm{400}\cdot k\cdot F}{L\cdot \mathit{\omega }\cdot \mathit{\rho }}}\right)}^{\mathrm{0.5}},\end{array}$$

where ε is the freeze–thaw depth (m), k is the thermal conductivity (2 W (m K)⁻¹), F is the freeze–thaw index (°C), L is the latent heat of the fusion of ice (3.35×10⁵ J kg⁻¹), ω is the water content as a decimal fraction of the dry soil weight (0.12), and ρ is the bulk density of the soil (1000 kg m⁻³).

The Stefan equation requires ground surface temperature, which is difficult to measure and often lacks data. During freezing, the air temperature was translated into ground temperature by a multiplier of 0.6, and the ground temperature was the same as the air temperature during thawing (Gisnås et al., 2016). In this research, the freeze–thaw process was simulated at each Hydrologic Response Unit (HRU) by the Stefan equation, driven by distributed air temperature. The lower limit of permafrost was also estimated by the distributed soil freeze index and thaw index, where the freeze index is equal to the thaw index in the mountain region.

In the freezing and frozen season, there is no runoff generated due to precipitation in the form of snowfall and the soil being frozen. During this period, runoff only comes from the groundwater of the supra-permafrost, and no runoff (R_U) is generated from the root zone reservoir to the fast-response reservoir (S_f) and slow-response reservoir (Ss). Therefore, we set the R_U to zero in this season. In the freezing season, when the freezing depth is less than 3 m, the groundwater discharge in the supra-permafrost is still connected, which can be simulated with a linear groundwater reservoir (S_s) and the slow-recession coefficient (K_s). When the freezing depth is greater than 3 m in a Hydrologic Response Unit, the groundwater is frozen, and there is little runoff generated from the groundwater discharge in the Hydrologic Response Unit. So, in the FLEX-Cryo model, the groundwater reservoir (S_s) was reduced to 10 % of its storage to represent the groundwater being frozen (Eqs. 3 and 4). The other 90 % of the storage water was frozen in the groundwater system (Eq. 4). In the model, the soil begins to freeze from the high elevation to the lower elevation, which affects the groundwater. The groundwater reservoir freezes along the elevation, stopping the function of a series of cascade groundwater buckets, which is the key reason for the discontinuous recession.

In the thawing season, the freeze–thaw conditions in the lowest elevation zone control the hydraulic connectivity between groundwater system and soil. If the freeze depth is larger than thaw depth calculated by the Stefan equation, the soil is still frozen, and the connectivity between groundwater system and soil is still closed. There is no runoff generated (Ru), but the root soil moisture (Su) accumulates, and evaporation is the only outflow from the root zone. Once the thaw depth is larger than the freeze depth, the frozen groundwater reservoir is released to the groundwater discharge (Eq. 4). Complete thawing at the lowest elevation represents the end of the thawing season and the start of the completely thawed season. In the completely thawed season, the groundwater and soil are connected and not affected by the frozen soil.

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{s}}}{\mathrm{d}t}}=R_{\mathrm{s}}-Q_{\mathrm{s}}-F_{\mathrm{s}}\text{(4)}& {\displaystyle }{F}_{\mathrm{s}}=\left\{\begin{array}{ll}\mathrm{0.9}\cdot S_{\mathrm{s}}& \text{freeze depth }\left(\mathit{\epsilon }\right)\ge \mathrm{3}\mathrm{m}\\ -\mathrm{0.9}\cdot S_{\mathrm{s}}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\text{once thaw depth reaches yearly }max\\ {\displaystyle }& {\displaystyle }\text{or thaw depth}\ge \text{thaw depth}\end{array}\end{array}\right.\end{array}$$

3.3 Δh parameterization

The Δh parameterization is a mass conservation method to assess the change of ice-covered, glacier length and glacier thickness in response to global warming. The glacier mass balance (GMB) calculated by FLEX-Cryo was redistributed to glacier elevation bands. It is an observed truth that the lower-elevation bands lose more ice than higher-elevation bands. The lost ice volume, calculated by a mass balance model, is converted into a distributed ice thickness change according to the Δh parameterization (Gao et al., 2021; Huss et al., 2010).

The Δh-parameterization method was employed, which relies on empirical curves that are dependent on the size of the glacier. The study categorized glaciers into three size classes: large valley glacier (area>20 km²), medium valley glaciers ($\mathrm{5}{\mathrm{km}}^{\mathrm{2}}<\text{area}<\mathrm{20}{\mathrm{km}}^{\mathrm{2}}$), and small glaciers (area<5 km²). Both Glacier 1 and Glacier 2 had areas of less than 5 km² and were categorized as small glaciers. The small glacier equation in this study is as follows:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }h=(h_{\mathrm{r}}-\mathrm{0.30}{)}^{\mathrm{2}}+\mathrm{0.60}(h_{\mathrm{r}}-\mathrm{0.30})+\mathrm{0.09},\end{array}$$

where Δh is normalized surface elevation change, and h_r is the normalized elevation range. Based on this equation, the glacier elevation and surface area were evolved every 5 years to avoid the circumstance of glacier advancing. The corresponding glacier melting HRU was transformed into an alpine desert landscape (Wei et al., 2023).

3.4 Spatial discretization of the catchment

The catchment area was divided into 37 elevation bands ranging from 2960 to 4820 m, with an interval of 50 m. These elevation bands were classified based on four landscapes: glacier, alpine desert, vegetation hillslope, and riparian zone (Figs. 2 and 3). As a result, there were a total of 148 Hydrologic Response Units (HRUs) in the catchment. The landscape of alpine desert was the most widespread, covering an elevation range of 3425 to 4727 m. The glacier was found in higher-altitude areas, specifically between the elevation bands of 3725 and 4727 m.

Figure 4 Taylor diagram of monthly temperature and precipitation simulation. The hollow points are the uncorrected projection, the solid points are the corrected projection, and the solid red circles are the reference values (observation).

3.5 Model evaluation metrics

The Kling–Gupta efficiency (KGE), Nash–Sutcliffe efficiency (NSE), coefficient of correlation (R), and root mean square error (RMSE) were used to comprehensively assess the model performance and the reliability for the model. The KGE, NSE, R, and RMSE are all less than 1. For KGE, NSE, and R, values closer to 1 indicate better performance. A lower RMSE value indicates less error and better model performance. These metrics can be calculated as follows:

$$\begin{array}{}\text{(22)}& {\displaystyle }\text{KGE}=\mathrm{1}-\sqrt{(r-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}},\text{(23)}& {\displaystyle }\text{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n(Q_{\mathrm{0}}-Q_{\mathrm{m}}{)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^n(Q_{\mathrm{0}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{0}}}{)}^{\mathrm{2}}}},\text{(24)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^n(Q_{\mathrm{0}}-Q_{\mathrm{m}}{)}^{\mathrm{2}}},\text{(25)}& {\displaystyle }R={\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n(Q_{\mathrm{0}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{0}}})(Q_{\mathrm{m}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{m}}})}{\sqrt{{\sum }_{t=\mathrm{1}}^n(Q_{\mathrm{0}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{0}}}{)}^{\mathrm{2}}}\sqrt{{\sum }_{t=\mathrm{1}}^n(Q_{\mathrm{m}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{m}}}{)}^{\mathrm{2}}}}},\end{array}$$

where r is the linear correction coefficient between simulation and observation, α is the ratio of the stand deviation of simulated variables and observed variables, β is the ratio of the average value of simulated and observed variables, Q₀ is the observation runoff, $\stackrel{\mathrm{‾}}{Q_{\mathrm{0}}}$ is the average observed runoff, and Q_m is the simulation runoff.

4 Results

4.1 Performance of bias correction and runoff depth simulation

4.1.1 Bias correction performance

The accuracy of climate projection varied with the multiple bias correction method (Fig. 4). The distance between the observation and the projection is inversely proportional to the accuracy. Before the bias correction, the distance is relatively far, especially for precipitation, indicating that there is a large error between observations and GCM projection. After the bias correction, the distance diminishes, indicating that the bias correction improves the accuracy, particularly for precipitation.

Figure 5 (a) Comparison of the modelled and observed glacier mass balance (GMB) of Glacier 1 from January 2011 to December 2014. (b) Comparison of the simulated freeze–thaw depth using the Stefan equation and observations. (c) Comparison of simulated runoff using the FLEX-Cryo model and observations from 2011 to 2014.

4.1.2 Performance of cryospheric elements and runoff simulation

We assessed the performance of the FLEX-Cryo model for glacier mass balance change, freeze–thaw depth, and runoff simulation based on historical observations. The model demonstrated strong capabilities across all evaluated aspects. For the glacier mass balance change, the model showed good accuracy throughout the entire assessment period. Monthly simulations yielded a KGE value of 0.45, an NSE of 0.83, a correlation coefficient R of 0.95, and an RMSE of 130.13 mm per month (Fig. 5a). Regarding the freeze–thaw dynamics, the model accurately captured both timing and duration. The simulated freeze onset consistently aligned with observations, typically occurring in late October and early November. Moreover, the simulated freeze–thaw cycle duration closely matched observations, with both spanning approximately 217 d and varying by no more than 15 d. Notably, the model exhibited exceptional accuracy in predicting maximum freezing depth, with a mere 2 mm error recorded in April 2013 (Fig. 5b). For the runoff simulation, the model showed very good performance over the assessment period, with a KGE value of 0.83, NSE of 0.73, R of 0.74, and RMSE of 0.77 mm d⁻¹ (Fig. 5c). These results demonstrated that the FLEX-Cryo model can effectively reproduce hydrographs and capture changes in the cryospheric elements. The good model performance in terms of various metrics demonstrates the robustness of the FLEX-Cryo model, providing confidence in its ability to accurately estimate future hydrological changes.

4.2 Future climate change

Figure 6 shows the prediction of future climate in 2015–2100 under SSP2-4.5 and SSP5-8.5 based on the average values of eight climate models (adjusted for bias). According to SSP2-4.5, the temperature will increase by 2.07 °C relatively steadily by 2100. Under SSP5-8.5, temperatures are projected to continue to rise by 5.04 °C over the course of the century. Precipitation changes are more variable than temperature, especially after the 2080s under SSP5-8.5. Overall, the precipitation under SSP2-4.5 increased by 14.25 %, and the precipitation increased by 33.50 % under SSP5-8.5. Before the 2080s, the increase in precipitation was almost the same under different scenarios, about 8.9 and 8.5 mm per 10 years, respectively. Although there are some uncertainties associated with temperature and precipitation, the increasing trends of temperature and precipitation are still distinguished, especially for SSP5-8.5.

Figure 6 (a) The annual average temperature and (b) the annual precipitation mean of bias-adjusted multi-global climate model from 2015–2100. The blue and red areas indicate the uncertainty caused by eight climate change models of SSP2-4.5 and SSP5-8.5 scenarios.

4.3 The change of cryosphere in the future

4.3.1 Predicting glacier retreat

Figure 7 shows changes in glacier thickness at the highest-elevation band and volume for Glacier 1 and Glacier 2 under two SSPs from the years 2015–2100. Starting from the 2020s, the glacier volume showed a rapid decline, and after the 2030s, the glacier entered a phase of rapid thinning. Around 2040, the glacier degradation reached a stabilization period, during which glaciers were only present in the highest-elevation band. According to the SSP2-4.5 scenario, Glacier 1 and Glacier 2 are projected to completely melt and disappear by 2051 and 2046, respectively. Under the SSP5-8.5 scenario, the complete melt-out time is slightly earlier, occurring in 2045 and 2044 for Glacier 1 and Glacier 2, respectively. After the glaciers completely melt, the ablated glacier area will transform into alpine desert.

Figure 7 The glacier thickness (a) and the glacier volume (b) change from 2015 to 2100 for Glacier 1 and Glacier 2.

Figure 8 Changes in seasonally frozen soil and permafrost from 2015–2100 under SSP2-4.5 and SSP5-8.5 scenarios. (a, b) Freeze–thaw onset. (c, d) Freeze–thaw duration of frozen soil and permafrost. (e, f) The maximum freezing depth and active layer thickness.

Taking the glacier changes in 2025, 2035, and 2045 as examples, under the SSP2-4.5 scenario, the area of Glacier 1 is projected to decrease to 5.49×10⁵, 1.52×10⁵, and 0.26×10⁵ m², with corresponding volume reductions to 5.27×10⁶, 1.03×10⁶, and 0.26×10⁶ m³, respectively (Fig. 9). Comparatively, the retreat trend is more pronounced under the SSP5-8.5 scenario. The area of Glacier 1 is projected to be 4.00×10⁵, 0.81×10⁵, and 0.26×10⁵ m², with volumes of 4.86×10⁶, 0.71×10⁶, and 0.03×10⁶ m³, respectively. The degradation of Glacier 2 follows a similar pattern to that of Glacier 1, except that Glacier 2 experiences less ice loss. According to the SSP5-8.5 scenario, Glacier 2 is projected to completely melt by 2045. In 2025 and 2035, the area of Glacier 2 is 1.67×105 and 0.51×105 m² for both scenarios, respectively. These glaciers are only distributed within the elevation bands from 4625 to 4727 m and from 4675 to 4727 m.

4.3.2 Forecasting the degradation of frozen soil

The degradation of seasonally frozen soil and permafrost is projected by the FLEX-Cryo model (Fig. 8). Under SSP2-4.5, by the end of 21st century, the freeze onset of seasonally frozen soil will be delayed by 10 d, and the freeze–thaw cycle duration will shorten by approximately 1 month. The maximum freeze depth of seasonally frozen soil is expected to decrease by 5.17 cm per decade. The thaw onset of permafrost will be advanced by 19 d, and the freeze–thaw cycle duration will increase nearly 50 d. The active layer thickness will rise by approximately 8.24 cm per decade. Meanwhile, the degradation trend of permafrost is more severe under the SSP5-8.5 scenario. Under SSP2-4.5, the freeze onset of seasonally frozen soil will be shortened by 22 d, and the freeze–thaw cycle duration will reduce by over 2 months. The thaw onset of permafrost will occur approximately 1 month earlier, and the freeze–thaw cycle duration of permafrost will increase by nearly 3 months. Compared with SSP2-4.5, the decreasing trend of the maximum freeze depth and the increasing trend of the active layer thickness are approximately twice as pronounced under SSP5-8.5. By 2100, seasonally frozen soil will begin to freeze around mid-November and late November, while permafrost will start to thaw in mid-May and early June by the year 2100 under two SSPs.

Under SSP2-4.5 and SSP5-8.5, the lower limit of permafrost gradually expands along the altitudinal gradient, with rates of 4.30 and 8.75 m yr⁻¹, respectively (Fig. 9). In SSP2-4.5, the lower limit of permafrost is projected to reach altitudes of 3685, 3795, 3835, 3865, 3985, and 4015 m in the years 2025, 2035, 2045, 2055, 2075, and 2095, respectively. The lower limit of permafrost in 2095 under the SSP2-4.5 scenario is comparable to the lower limit of permafrost (3965 m) in 2055 under the SSP5-8.5 scenario. The lower limit is projected to increase to 4355 m by 2095 under the SSP5-8.5 scenario.

Figure 9 Changes in ice thickness and the lower limit of permafrost in 2025, 2035, 2045, 2055, 2075, and 2095 under SSP2-4.5 and SSP5-8.5.

4.3.3 Snow change in the future

The duration of snow cover is projected to decrease continuously in the future (Fig. 10). Under SSP2-4.5, snow cover days are likely to be shortened by 45 d, and snow water equivalent will decrease by 0.24 mm yr⁻¹. Compared with SSP2-4.5, snow cover reduces more under SSP5-8.5. Under SSP5-8.5, snow cover day is expected to be around 76 d, and snow water equivalent will decrease by 0.35 mm yr⁻¹.

Figure 10 Changes in snow water equivalent and snow cover days from 2015–2100 under SSP2-4.5 and SSP5-8.5.

4.4 Projected future runoff

The runoff in the catchment was predicted by the FLEX-Cryo model under SSP2-4.5 and SSP5-8.5. The tipping point of the glacier melting has already occurred (around 2020). After the turning point, glacier runoff and runoff of the total basin decrease dramatically until the glacier completely melts. Then, the runoff of the total basin will increase moderately. After the glacier completely melts, runoff of the total basin will decrease by 15.56 % and 18.05 %, respectively. The runoff coefficient, which represents the proportion of precipitation that becomes runoff, follows a similar pattern to the glacier runoff changes. It initially increases until the turning point of glacier melting occurs, then it decreases, and it eventually reaches a relatively stable state after the glacier completely melts (Fig. 11c). Before the turning point, the runoff coefficient is almost equal to or even greater than 1. The maximum values of the runoff coefficient occur in 2021 and 2019, coinciding with the tipping points of the glacier runoff. By the end of the 21st century, the runoff coefficient is projected to be dramatically reduced to approximately 0.42. These results indicate that glaciers play a key role in water resource supply.

Figure 11 (a) The predicted runoff depth of the total basin. (b) Runoff in the glacier and in the non-glacier from 2015–2100. (c) Project runoff coefficient under SSP2-4.5 and SSP5-8.5 scenarios.

Two hydrological phenomena observed in permafrost mountainous catchments, namely the low runoff in the early thawing season (LRET) and discontinuous baseflow recession (DBR) (Gao et al., 2022), are expected to persist in the future (Fig. 12). Meanwhile, baseflow, which represents the sustained flow of water from groundwater, shows an increasing trend. The duration of the early thawing season is projected to be further reduced. The first recession coefficient remains unchanged, while the second recession coefficient progressively increases. Under the SSP2-4.5 scenario, the second recession coefficient is equal to 74 d, which is consistent with the recession coefficient in 2060 under the SSP5-8.5 scenario. This suggests that the permafrost area undergoes less significant change under the SSP2-4.5 scenario than the SSP2-8.5 scenario according to Fig. 9. The baseflow gradually increases, especially in the SSP5-8.5 scenario, as indicated by the runoff depth on a logarithmic scale (Fig. 12).

Figure 12 Temperature, precipitation, runoff depth, and freeze–thaw cycle in 2020, 2040, 2060, 2080, and 2100 under SSP2-4.5 (a) and SSP5-8.5 scenarios (b).

5 Discussion

5.1 The effect of the mountain cryosphere degradation on runoff

Glaciers and snow are sensitive to climate change and cover play a crucial role in water retention, with meltwater contributing significantly to downstream water resources and the ecological environment (Stecher et al., 2023; Nan and Tian, 2024; Liu et al., 2022). The turning point of glacier runoff represents a critical tipping point that signifies not only the rapid thinning of glaciers but also the irreversible stage of water resources in the basin (Brovkin et al., 2021; X. Chen et al., 2022). After the turning point, the glacier thickness and glacier volume rapidly decrease (Fig. 7). But the glacier thickness shown in this paper represents the change in the highest-elevation band, which means the turning point would lag in this band for a change in glacier thickness. In the Hulu catchment, the proportion of glacier runoff reached 51 % to 55 % between 2019 and 2021, indicating that it is in the turning point period (Fig. 11). Subsequently, the contribution of glacier runoff gradually decreases until complete melting occurs. Temperature is the primary factor influencing glacier runoff, while precipitation and temperature together determine the proportion of glacier runoff in relation to total runoff. Although the highest contribution of glacier runoff and the tipping point of glacier runoff may not align precisely, after the tipping point, the capacity of glacier runoff to contribute to overall runoff continuously diminishes. From 2015 to 2021, there has been a decreasing trend in precipitation, leading to a corresponding decline in non-glacier runoff (Figs. 6 and 11). Thus, while glacier runoff has increased, the total runoff has decreased. However, between 2032 and 2038, even though rainfall continues to decline, the contribution of glacier runoff to overall runoff becomes negligible due to the limited volume of ice remaining ($\text{glacier volume}<\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ m³), resulting in minimal glacier melting runoff (Figs. 7 and 11). On the other hand, once the glaciers have completely melted, the total runoff in the Hulu catchment is reduced by 16 % to 18 %, and the runoff coefficient is halved (Figs. 9 and 10). This highlights the critical role of glaciers as solid freshwater reservoirs in regulating water sources and mitigating droughts (McCarthy et al., 2022; Mohammadi et al., 2023).

The freeze–thaw cycle has a significant impact on runoff yield and hydrological response routines in the Hulu catchment (Sun et al., 2022; Wang et al., 2020). Precipitation in the Hulu catchment is primarily concentrated in the summer when soil moisture is high and even close to saturation, making saturation excess flow the main mechanism for runoff generation (Li et al., 2016). During the freeze–thaw cycle, the weak permeability of frozen soil affects both surface runoff and infiltration. Soil runoff primarily occurs through underground in hillslope areas and surface water flow in riparian areas, resulting in a faster response to rainfall and snowmelt and contributing to a higher runoff coefficient (Hu et al., 2022; Jones et al., 2023). However, it is important to note that shallow frozen soil does not completely block the interaction between deeper soil layers and the surface. Frost heave in the soil creates large pores, allowing snowmelt water and precipitation to bypass the matrix layer and reach the deeper soils (Jiang et al., 2021; Zhang et al., 2023). This phenomenon is considered one of the significant reasons for low runoff in the early thawing season (Mohammed et al., 2021). Low runoff is observed between the frozen season and complete thawing season (Fig. 12). The duration of freeze–thaw cycles in seasonally frozen soils is shortening, and freeze onset is being delayed due to the warming climate, resulting in a decreasing duration of low runoff. However, the temperature during the freezing season remains lower than the initial frost heave temperature of the soil, and there is still a deficit of soil water in the early thaw, indicating that the prevalence of low runoff will persist in the future (Teng et al., 2022; Wen et al., 2024).

The freezing state has a significant impact on the recession process of baseflow, and permafrost plays a crucial role in discontinuous baseflow (Cooper et al., 2023; J. Wang et al., 2022). During the freezing season, baseflow follows a linear recession process (K_s=80 d), with contributions from both permafrost and seasonally frozen soil regions (Fig. 12). In the frozen season, the groundwater under the supra-permafrost layer becomes inactive, and baseflow is solely derived from the seasonally frozen soil regions, causing a discontinuous recession. With climate warming, the lower limit of permafrost gradually moves upward along the elevation, resulting in the shrinking of the permafrost region. This suggests that in the future, an increased proportion of baseflow will originate from the expanding area of seasonally frozen soil, leading to a gradual decrease in the influence of permafrost on baseflow. Consequently, the discontinuous recession of baseflow will gradually transition into a linear recession. Furthermore, an increase in the thickness of the active layer enhances the soil water storage capacity, contributing to a gradual rise in baseflow (Yao et al., 2021).

5.2 Comparison with other studies

The cryosphere, including glaciers (ice sheets), seasonal snow cover, frozen soil, and permafrost, plays a vital role in storing approximately 75 % of the world's freshwater resources (Qin et al., 2021). Although there are some differences in the driving data and models, the trends of the cryospheric elements and runoff changes are still comparable and consistent. In this study, the small glaciers are projected to completely melt in the mid-21st century, which is also reported in other areas (Mukhopadhyay and Khan, 2015; Baraer et al., 2012; Schwank et al., 2014). The projected maximum freeze depth of seasonally frozen soil calculated in this research is 5.2 cm per decade, similar to the 5.4 cm per decade predicted by Wang et al. (2018). Ni et al. (2021) showed that Qinghai–Tibet Plateau permafrost is at risk of disappearance based on statistical and machine learning (ML) modelling approaches. This shift in regions with permafrost impacts hydrological connectivity, fostering improved hydrothermal conditions that enhance vegetation growth (Han and Menzel, 2022; Jin et al., 2022). Few studies have focused on future change in snow cover days and snow water equivalent in the Heihe River basin, but many research studies have indicated that the snow-free period increases and the snow water equivalent decreases due to climatic warming in the Tibetan Plateau (Zhang and Ma, 2018). The reduction of the snow cover period may result in an earlier peak in spring snowmelt floods, thereby increasing the risk of flooding (Chai et al., 2022). Simultaneously, the decrease in snow water equivalent may impact plant water supply, placing pressure on ecosystems (Guan et al., 2022). Although cryospheric elements have a trend of degradation in different regions, the impact on runoff may differ. However, on a longer timescale, the degradation of the cryosphere will lead to a decrease in runoff (Xu et al., 2024). This study confirmed that runoff from cryospheric melting is one of the main factors controlling runoff, and degradation of the cryosphere may exacerbate the risk of future droughts.

5.3 Uncertainty and limitations

Uncertainty in this study comes from the GCMs, the downscaling and bias correction methods, and the structure and parameters of the FLEX-Cryo model. The temperature and precipitation projections from different GCMs at the basin scale introduce uncertainty. Moreover, four bias correction methods were used to correct the GCMs based on the observation, which may ensure consistent relative trends but not improve the accuracy of precipitation and temperature frequency distribution and seasonal variations. This may cause some uncertainty in the simulation results (Jia et al., 2023).

In this research, the time-variant albedo information and the aspect are worth being taken into account for improving glacier melting simulations, which require further observation and quantitative studies (Arnold et al., 2006; Feng et al., 2024). The change in elements is sensitive to energy. The snow cover and the effect of topographic shading may also have an effect on the degradation and thus hydrologic response, which warrants further investigation (Zhang, 2005). On a long timescale, the degradation of frozen soil and glaciers may result in thaw lake generation and other landscapes changes, which may have an effect on the runoff yield and baseflow recession (Serban et al., 2021).

6 Conclusions

In this study, we employed the FLEX-Cryo model and data from eight global climate models (GCMs) under the SSP2-4.5 and SSP5-8.5 scenarios to predict the potential impacts of climate change on the mountain cryosphere and hydrology. Results from the projected change of mountain cryosphere elements, glacier, snow, and frozen soil are expected to undergo degradation. The glacier will completely melt by the middle of the 21st century. Snow cover days will decrease by 45 and 76 d, and snow water equivalent will decrease by 0.24 and 0.35 mm yr⁻¹. The thaw onset is expected to advance 19 and 32 d. The active layer thickness will increase by 8.24 cm per 10-year period.

The degradation of the mountain cryosphere has significant implications for water resources. The tipping point of glacier runoff is projected to occur in the 2020s. Once the glaciers have completely melted, the runoff is projected to decrease by approximately 16 % and 18 % under the SSP2-4.5 and SSP5-8.5 scenarios, respectively. Importantly, the duration of low runoff will shorten, baseflow will increase, and the discontinuous recession of baseflow will gradually transform into a more linear pattern.

This study provides insights into the potential impacts of climate change on the mountain cryosphere and hydrology. The projected changes in glacier retreat, snow cover, and frozen soil dynamics highlight the urgent need for proactive water resource management strategies in the face of a changing climate. Further modelling research and monitoring efforts are necessary to refine these projections and guide effective adaptation measures to sustainably manage water resources in mountainous regions.