Future snow changes and their impact on the upstream runoff in Salween

Understanding the hydrological processes related to snow in global mountainous regions under climate change is necessary for achieving regional water and food security (e.g., the United Nation's Sustainable Development Goals 2 and 6). However, the impacts of future snow changes on the hydrological processes in the high mountains of the “Third Pole” are still largely unclear. In this study, we aimed to project future snow changes and their impacts on hydrology in the upstream region of the Salween River (USR) under two shared socioeconomic pathway (SSP) scenarios (SSP126 and SSP585) using a physically based cryosphere–hydrology model. We found that the climate would become warmer (0.2 ^∘C per decade under SSP126 and 0.7 ^∘C per decade under SSP585) and wetter (5 mm per decade under SPP126 and 27.8 mm per decade under SSP585) in the USR in the future under these two SSPs. In this context, the snowfall, snow cover, snow water equivalent, and snowmelt runoff are projected to exhibit significant decreasing trends during 1995–2100, and the decreases are projected to be most prominent in summer and autumn. The future (2021–2100) snowmelt runoff is projected to significantly increase in spring compared with the reference period (1995–2014), which would benefit the availability of water resources in the growing season. The annual total runoff would significantly increase in all of the future periods due to increased rainfall, which would increase the availability of water resources within the basin, but the high peak flow that occurs in summer may cause rain flooding with short duration and high intensity. Compared with the reference period (the contribution of snowmelt runoff to the total runoff was determined to be 17.5 %), the rain- and snow-dominated pattern of runoff would shift to a rain-dominated pattern after the near term (2021–2040) under SSP585, whereas it would remain largely unchanged under SSP126. Climate change would mainly change the pattern of the snowmelt runoff, but it would not change the annual hydrograph pattern (dominated by increased rainfall). These findings improve our understanding of the responses of cryosphere–hydrological processes under climate change, providing valuable information for integrated water resource management, natural disaster prevention, and ecological environmental protection at the Third Pole.

1 Introduction

Snow, a key component of the cryosphere, is widely distributed in high mountainous regions around the world, which are particularly sensitive to climate change; thus, it is an important indicator of regional climate change (Nepal et al., 2021; Pulliainen et al., 2020). The snowpack can store a large amount of solid precipitation in the cold season and can melt in the warm season. Accordingly, it not only has a strong effect on the regional hydrological cycle (Huning and AghaKouchak, 2020; Musselman et al., 2021) but also provides abundant water resources for industry, agriculture, and residents in basins as well as supporting about one-sixth of the world's population (Yan et al., 2022; Barnett et al., 2005). In addition, snow cover (SC) affects the radiation balance and the thermal regime of the underlying ground due to the high albedo of snow, thereby changing the regional energy balance and, in turn, affecting regional climate (Jia et al., 2021; You et al., 2020; Henderson et al., 2018; Xiao et al., 2017). Changes in snowfall, snow storage, and snowmelt under climate change would not only change the total annual and seasonal runoff at mountain outlets but would also affect and change mountain glaciers and the availability of water resources in downstream regions (Immerzeel et al., 2020; Li et al., 2019). These factors may increase the frequency, intensity, and range of natural disasters and may further threaten the security of the water supply, flood/drought control, and ecological security in downstream areas (Qin et al., 2020; Biemans et al., 2019). Therefore, an in-depth understanding of the impact of climate change on the snow-related hydrological regime in mountainous areas is crucial to regional ecological protection, water resource management, disaster prevention, and sustainable socioeconomic development (Nepal et al., 2014; Biemans et al., 2019; Yao et al., 2019; Tang et al., 2019; Viviroli et al., 2020; Qi et al., 2020).

Over the last few decades, the Tibetan Plateau (TP) has experienced intense warming, which has exceeded the global warming rate (IPCC, 2019). In addition, precipitation has also exhibited an increasing trend in the central TP and a decreasing trend in the southern and eastern TP under climate warming (Chen et al., 2015). Moreover, several studies have predicted that the TP will continue to warm in the future, accompanied by increased precipitation, especially in the monsoon-controlled regions, with an increased frequency and intensity of extreme events (Panday et al., 2015; Sanjay et al., 2017). These influences not only change the spatiotemporal distribution and magnitude of precipitation but also alter several of the key variables (e.g., solid precipitation) that drive snow occurrence and development as well as the variables (e.g., radiation and temperature) that control snow ablation (Hock et al., 2019). Yao et al. (2019) and Bibi et al. (2018) reported that the melting of the snowpack, especially in low- and mid-elevation areas (because more of the precipitation occurs as rainfall in these regions), has accelerated in recent decades on the TP, which has increased the instability and uncertainty of the inter-annual and seasonal runoff and will cause future changes in the spatiotemporal pattern and availability of water resources in the TP river basins (Tang et al., 2019; Immerzeel et al., 2010; Kraaijenbrink et al., 2021). Previous studies have assessed several snow variables (e.g., snowfall, snow storage, SC, and snowmelt) and the hydrological processes related to snow under climate change on the TP based on in situ observations and land surface snow–hydrology models (Bian et al., 2020; Xu et al., 2017; Barnhart et al., 2016; Su et al., 2016). However, the impact of future snow changes on runoff is still unclear due to a lack of reliable data. Moreover, most of these models did not fully consider the physical processes of snow accumulation and ablation (Wang et al., 2017; Liu et al., 2018). Therefore, a comprehensive cryosphere–hydrological model and high-quality forcing datasets are urgently needed to improve our understanding of snow-related hydrological processes on the TP in order to better support the sustainable development of this region.

To better understand the effect of snow changes on runoff on the TP under climate change, the upstream region of Salween River (USR) was selected as the study area for the present work. The USR is located in the alpine region of the TP and Hengduan Mountains, which has a complex underlying surface and is very sensitive to climate change (Liu et al., 2017; Chen et al., 2020). Moreover, as a transboundary river with great influence, the reasonable utilization and coordination of water resources in the Salween Basin have received a great deal of attention from relevant countries and organizations worldwide (Yao et al., 2012). Previous studies have shown that the complex underlying surface, ecosystem, and hydrological processes of the basin have undergone intense changes in recent years, which have gradually changed the total water resources and runoff within the basin, further affecting the rationing and allocation of water resources and hydropower development (Liu et al., 2017; Fan and He, 2012; Luo et al., 2017; Hong and He, 2019). However, the working conditions in the USR are very difficult due to the complexity of the topography and the harsh environment, which has resulted in an irregular and relatively sparse distribution of national meteorological stations in the basin. In addition, due to a lack of continuous snow and hydrological observations in the basin resulting from strict national data policies, it is difficult to quantify and reveal the mechanisms of snow-related hydrological processes and verify the results of numerical simulations (Wang et al., 2021). Moreover, the precipitation and snow products used in previous studies have low accuracy and poor applicability in this basin and are often insufficient for the analysis of trends due to the short periods of the records. These factors have seriously restricted our understanding of snow-related hydrological processes and the development and utilization of water resources in this basin (Mao et al., 2019; Liu et al., 2016; Ding et al., 2015).

Previous studies have used simple statistical methods to study the variations in the historical runoff in the USR based on observation data from individual sites; the dry season and annual runoff have been found to exhibit an increasing trend as a result of the increase in precipitation and meltwater (glacier and snow), wherein the precipitation contributes the most (Zhang et al., 2007; You et al., 2008; Yao et al., 2012; Cuo et al., 2014; Luo et al., 2016; Liu et al., 2017; Zhang et al., 2019). However, these studies did not separate glacier melt and snowmelt, consider the corresponding physical processes, nor predict future changes in runoff and its components. In addition, Su et al. (2016), Lutz et al. (2014), Zhao et al. (2019), and Khanal et al. (2021) used hydrological models to predict the future runoff on the TP, and they all found that the future runoff would exhibit an increasing trend; however, there are great uncertainties in the meltwater contribution and seasonal variations in these models. This may be due to the different descriptions of the cryospheric processes used in these models and the large differences in the driving datasets used. Moreover, these studies did not consider the intermediate snow change processes, which may lead to a partial understanding of snow and hydrological processes.

The main objective of this study was to simulate the changes in the snow-related hydrological processes in the USR on the southeastern TP, China. First, based on meteorological observation data, we evaluated the performance of four reanalysis precipitation products and selected the most reliable products. Second, we constructed a distributed cryosphere–hydrological model of the USR basin, namely the Water and Energy Budget-based Distributed Hydrological Model with improved Snow and Frozen ground processes (WEB-DHM-sf) (Wang et al., 2009a, b, 2010, 2016, 2017), and evaluated its performance using observed discharge and remote sensing data (land surface temperature, LST, Moderate Resolution Imaging Spectroradiometer, MODIS, SC). Third, we used WEB-DHM-sf driven by a Coupled Model Intercomparison Project Phase 6 (CMIP6) global climate model (GCM) to predict the changes in the snowfall, snow cover, snow water equivalent (SWE), total snowmelt, snowmelt runoff, and total runoff during different periods under the following different shared socioeconomic pathway (SSP) scenarios:

• SSP126, which is a combination of the RCP2.6 low-emission scenario and the SSP1 sustainable socioeconomic pathway, representing low vulnerability, low mitigation pressure, and low radiative forcing (2.6 W m⁻²);

• SSP585, which is a combination of the RCP8.5 high-emission scenario and the SSP5 fossil fuel development pathway and is the only SSP that can achieve an anthropogenic radiative forcing of 8.5 W m⁻² by 2100, representing the worst development pathway.

We also further analyzed the impact of the snow changes on runoff.

Figure 1 (a) The topography and distribution of meteorological and hydrological stations, (b) soil type, and (c) land use types in the upstream region of the Salween River (USR) basin.

2 Study area

The Salween River, also known as the Nujiang River in China, is a free-flowing international river, and it is also an important strategic hydropower and water resource reserve area in Southeast Asia (Lu et al., 2021). It originates from the southern foot of the Tanggula Mountains on the Tibetan Plateau in China, flows across Myanmar and Thailand from north to south, and finally flows into the Andaman Sea (He, 2004). In China, its length and basin areas are 2013 km and 137 800 km², respectively (Guo, 1985). The basin above Jiayuqiao (JYQ) station is called the USR, with an average elevation and area of 4800 m and 72 975 km², respectively (Fig. 1a). This region is a typical alpine mountainous zone and is controlled by the westerlies and the South Asian monsoon. The precipitation in the area is mainly concentrated from May to October (wet season). The temperature and precipitation are highly dependent on elevation and exhibit significant spatial heterogeneity (Mao et al., 2019). Therefore, this region is also sensitive to climate change (Liu et al., 2017). The snow depth exhibited a significant decreasing trend in the USR during 2000–2018. The maximum monthly snow depth occurred in December–January, and the lowest monthly snow depth occurred in July–August (Yan et al., 2021). The annual mean fraction of snow cover was about 50 % during 2001–2014, and it was mainly distributed at 4000–6000 m (Li et al., 2017). The area covered by glaciers was greater than 800 km², accounting for 1.2 % of the entire basin. The soil was mainly Lithosol soil, followed by frozen soil, and a small amount of Planosol soil (Fig. 1b). The vegetation types were mainly agriculture/C₃ grassland, with some shrubs (Fig. 1c).

3 Data and methods

3.1 Distributed cryosphere–hydrological model

WEB-DHM-sf has been further improved compared with the first version of the Water and Energy Budget-based Distributed biosphere Hydrological Model (WEB-DHM; Wang et al., 2009a, b), which coupled and improved the three-layer energy balance snow module of the Simplified Simple Biosphere 3 (SSib3) model (Shrestha et al., 2010) and the empirical frozen soil parameterization scheme (Wang et al., 2010). The three-layer snow module divides the snow pack (at each model grid) into layers, with snow depths of greater than 5 cm being divided into three layers and shallower depths being regarded as single-layer snow. It considers the energy exchange between the snow layers and the influences of the incidence angle of the solar radiation and the age of the snow on the snow albedo. Therefore, this model can describe the physical processes of snow in detail, including the phase transition, compaction, albedo, temperature, and melt runoff of each layer. Thus, we can obtain the SC, snowfall, SWE, snowmelt, and other variables from the model outputs (Shrestha et al., 2014). The frozen ground processes are characterized by the frozen soil hydrothermal transfer parametrization scheme, in which the thermal conductivity scheme used is the Johansen scheme (Wang et al., 2017). In both snow parameterization schemes, enthalpy is used instead of temperature to establish the energy equation or as a new predictive variable, which reduces the uncertainty when calculating the latent heat released by the changes in the water phase and also enhances the stability of the model (Wang et al., 2017; Song et al., 2020). In this version of the model, glaciers are considered to be snow with a thickness of 100 m (the details of this model have been described by Wang et al., 2016, and Shrestha et al., 2010, 2014). Due to the continuous improvements in recent years, this model can better describe the cryospheric–hydrological processes in the alpine region, and it has been verified to have a good performance in several basins in the TP (Liu et al., 2018; Qi et al., 2019, 2022; Zhong et al., 2020; Wang et al., 2021; Zhou et al., 2022), especially in areas where observations are very scarce. In addition, the outputs of the model can be verified using multisource data, such as in situ observations and satellite remote sensing data. Therefore, this model was used in this study to analyze the impact of the snow changes on the runoff with a temporal resolution of 1 h and a spatial resolution of 5 km.

3.2 Data

3.2.1 Meteorological forcing data

The meteorological forcing data included the near-surface air temperature, total precipitation, downward shortwave and longwave radiation, wind speed, surface pressure, and specific humidity. For the historical period (1995–2014), except for the precipitation data, the data for the meteorological variables were obtained from the China Meteorological Forcing Dataset (CMFD), which has a high spatiotemporal resolution (3 h and 0.1^∘) from 1979 to 2018 and has been widely used in hydrometeorological research because of its good applicability (He et al., 2020). To better evaluate the applicability of four reanalysis precipitation products – the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth-generation reanalysis data product (ERA5), the Global Land Data Assimilation System (GLDAS) product, version 2 of the Modern-Era Retrospective Analysis for Research and Applications (MERRA-2) product, and the Multi-Source Weighted-Ensemble Precipitation (MSWEP) product (Beck et al., 2017) – in the USR, we mainly used several statistical indices to evaluate the products based on the meteorological observation data in this work, including correlation coefficients (CCs), mean bias (MB), and the root-mean-square error (RMSE). Moreover, the probability of detection (POD), the false alarm ratio (FAR), the missing alarm rate (MAR), and the critical success index (CSI) were calculated to check the capture capability of rainfall events from reanalysis products (Ebert et al., 2007; Tian et al., 2009). It can be seen from Figs. S1 and S2 that, except for GLDAS, the statistical indicators of other precipitation products are relatively consistent; however if the international exchange station is excluded (Naqu and Dingqing), ERA5 and MSWEP have good performance (CC is about 0.5). In addition, ERA5 has a higher POD (average of about 0.97) and a lower FAR (average of about 0.44) as a whole, which is also shown in the CSI (average of about 0.56), implying that ERA5 has better characterization of rainfall events. Therefore, we selected ERA5 based on its 1 h and 0.25^∘ resolutions from 1979 to the present, its better spatiotemporal distribution, and the fact that it has been widely used in hydrological simulations of high mountainous regions (Hersbach et al., 2020; Yang et al., 2021). The projected climate forcing variables used in the hydrological model were obtained from the Inter-Sectoral Impact Model Intercomparison Project 3b (ISIMIP3b) (Lange and Büchner, 2021), of which four GCMs (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, and MRI-ESM2-0) that have a consistent experimental protocol (historical, SSP126, and SSP585) and atmospheric climate variables (spatial and temporal resolutions of 0.5^∘ and 1 d, respectively) were selected as the model input. The ISIMIP3BASD (v2.4.1) method was used for the bias adjustment and statistical downscaling of these variables to ensure their consistency with the long-term statistics of the observation reference dataset, i.e., W5E5 v2.0 (WFDE5 v0.1 over land merged with ERA5 over the ocean) (Lange, 2019, 2021a, b; Cucchi et al., 2020). To maintain the relative and absolute trends in these variables during the historical and future periods at the basin scale, we further corrected these variables on the monthly scale using the delta method (a simple linear bias correction method), which has been widely used and can effectively avoid a large model bias (Cuo et al., 2011; Su et al., 2016); the reader is referred to the bias correction methods in the Supplement for details. The historical and future meteorological forcing data were interpolated to the model's resolution (5 km) via bilinear interpolation. The time downscaling method used for the daily forcing data was consistent with that used by Song et al. (2020).

3.2.2 Digital elevation model, soil, land use, vegetation, and glacier data

A digital elevation model (DEM) was used to extract the basin's boundary and calculate the topographic parameters. The DEM was the National Aeronautics and Space Administration (NASA) Shuttle Radar Topography Mission (SRTM) DEM with a spatial resolution of 90 m. The soil type and land use maps with a spatial resolution of 1 km used to calculate the soil hydraulic parameters and land use types were obtained from the Food and Agriculture Organization (FAO, 2003) and the U.S. Geological Survey (USGS) (Fig. 1a, b), respectively. The fraction of photosynthetically active radiation (FPAR) and the leaf area index (LAI) used to calculate the vegetation dynamics were downloaded from the Global Land Surface Satellite (GLASS) datasets, both with respective temporal and spatial resolutions of 8 d and 0.05^∘ (Xiao et al., 2014). The second glacier inventory of China was provided by the National Tibetan Plateau Data Center (Guo et al., 2015). These data were resampled to the model's resolution (5 km).

3.2.3 Observation data

The daily observed discharge data at the outlet (JYQ) of the USR basin were obtained from the National Hydrology Almanac of China and were used to calibrate and validate the hydraulic parameters of the hydrological model during 1981–1987. There are six meteorological stations within the basin. Among them, Naqu and Dingqing stations are the international exchange stations (Fig. 1a). The meteorological observation data used to evaluate the precipitation products were obtained from the China Meteorological Administration (CMA). The land surface temperature (LST) data were obtained from the MODIS and MOD11A2 products, with respective spatial and temporal resolutions of 1 km and 8 d from 2001 to 2018 (Wan et al., 2014). The SC data were obtained from an improved Terra–Aqua MODIS snow cover and Randolph Glacier Inventory 6.0 combined product (MOYDGL06^∗) for high-mountain Asia with respective temporal and spatial resolutions of 8 d and 500 m from 2002 to 2018, which can be accessed from PANGAEA (Sher and Amrit, 2019).

3.3 Model evaluation criteria

The evaluation criteria used to evaluate the performance of the hydrological model outputs (e.g., discharge, LST, and snow cover) mainly included the Nash–Sutcliffe coefficient (NSE), the Kling–Gupta coefficient (KGE), the correlation coefficient (CC), mean bias (MB), root-mean-square error (RMSE), and relative bias (RB). The NSE and KGE can complement each other and make the evaluation more reasonable (Gupta et al., 2009). The equations used to calculate these evaluation criteria are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left({\text{Sim}}_i-{\text{Obs}}_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left({\text{Obs}}_i-\stackrel{\mathrm{‾}}{\mathrm{Obs}}\right)}^{\mathrm{2}}}},\text{(2)}& {\displaystyle }\mathrm{\text{KGE}}=\mathrm{1}-\sqrt{{\left(\text{CC}-\mathrm{1}\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\delta }}_{\text{sim}}}{{\mathit{\delta }}_{\text{obs}}}}\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathrm{Sim}}}{\stackrel{\mathrm{‾}}{\mathrm{Obs}}}}\right)}^{\mathrm{2}}},\text{(3)}& {\displaystyle }\text{CC}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left({\text{Sim}}_i-\stackrel{\mathrm{‾}}{\mathrm{Sim}}\right)\left({\text{Obs}}_i-\stackrel{\mathrm{‾}}{\mathrm{Obs}}\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^n{\left({\text{Sim}}_i-\stackrel{\mathrm{‾}}{\mathrm{Sim}}\right)}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^n{\left({\text{Obs}}_i-\stackrel{\mathrm{‾}}{\mathrm{Obs}}\right)}^{\mathrm{2}},}}}\end{array}$$

$$\begin{array}{}\text{(4)}& {\displaystyle }\text{MB}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left({\text{Sim}}_i-{\text{Obs}}_i\right)}{n}},\text{(5)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left({\text{Sim}}_i-{\text{Obs}}_i\right)}^{\mathrm{2}}}{n}}},\text{(6)}& {\displaystyle }\text{RB}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left({\text{Sim}}_i-{\text{Obs}}_i\right)}{{\sum }_{i=\mathrm{1}}^n{\text{Obs}}_i}}\times \mathrm{100}\mathit{\%}.\end{array}$$

Here, Sim_i is the simulated value at time i; Obs_i is the observed values at time i; n is the number of samples; $\stackrel{\mathrm{‾}}{\text{Sim}}$ and δ_sim are the mean value and standard deviation of the simulated values, respectively; and $\stackrel{\mathrm{‾}}{\text{Obs}}$ and δ_obs are the mean value and standard deviation of the observed values, respectively. The closer the values of NSE, KGE, and CC are to one, the better the simulation. The closer the values of MB, RMSE, and RB are to zero, the better the simulation.

Figure 2 Simulated and observed (a) daily and (b) monthly discharges at Jiayuqiao (JYQ) station from 1981 to 1987. The calibration and validation periods were 1981–1983 and 1984–1987, respectively.

Figure 3 Simulated (a) daytime and (b) nighttime 8 d land surface temperature (LST) compared to MODIS satellite observations, averaged over the USR basin during 2001–2018. Here, the input air temperature (T_air) has been compared with simulated LST (purple line).

4 Results

4.1 Model calibration and validation

4.1.1 Model calibration and validation using observed discharge

First, the model parameters were calibrated using the daily observed discharge data measured at JYQ station from 1981 to 1983, which mainly included the following soil hydraulic parameters: the saturated soil moisture content (θ_s=0.48), the residual soil moisture content (θ_r=0.07), the hydraulic conductivity anisotropy ratio (anik = 2.17), the van Genuchten parameter (α=0.014, n=2.385), the saturated hydraulic conductivity for soil surface (K_sat1=141.45 mm h⁻¹), the saturated hydraulic conductivity for unsaturated areas (K_sat2=21.22 mm h⁻¹), and the hydraulic conductivity for an unconfined aquifer (K_g=0.7 mm h⁻¹). Then, the daily observed discharge data for 1984–1987 (observations for 1986 were missing) were used to validate the calibrated hydrological model. The results are shown in Fig. 2. The time series of the simulated discharge matched the observations during 1981–1987 well and can capture the daily and monthly changes in the observed discharge, especially in the dry season. The values of the evaluation criteria during the calibration and validation period were determined to be as follows: 0.7 (NSE), 0.75 (R²), 0.84 (KGE), and 3.18 % (RB) for calibration, and 0.79 (NSE), 0.81 (R²), 0.88 (KGE), and 3.29 % (RB) for validation (Fig. 2a). The monthly simulated and observed discharge were in better agreement during the peak flow (Fig. 2b), accurately reproducing the inter-annual and intra-annual variations in the discharge at JYQ station, although there were overestimates and underestimates during the wet seasons in some years.

4.1.2 Model validation using LST

Due to the lack of observed discharge data from 1988 to the present, remote sensing data (MODIS LST and snow cover) were used to further validate the performance of the hydrological model at the basin scale. Figure 3 shows the comparison of the results for the basin-averaged time series between the simulated and MODIS LST values at daytime (10:30 LT, local time) and nighttime (22:30 LT) during 2001–2018. The results show that the simulated LST is in good agreement with the MODIS LST at daytime (R²=0.69) and nighttime (R²=0.91). The simulated nighttime LST has a higher R² value and lower MB value, but the RMSE is better than that of the daytime LST, which indicates that the LST is mainly affected by the incoming solar radiation and land surface absorption during the daytime and is affected very little at night.

From the perspective of the spatial distribution patterns (Fig. 4), the simulated multiyear seasonal average LST during the daytime and nighttime were similar to those of the MODIS LST from 2001 to 2018. Similarly, the simulated nighttime LST was better than the simulated daytime LST. The basin-average value of simulated seasonal LST during the nighttime was closer to the MODIS LST (December–January–February, DJF: MB $=-\mathrm{2.25}$ K; June–July–August, JJA: MB = −1.07 K; March–April–May, MAM: MB $=-\mathrm{1.5}$ K; September–October–November, SON: MB $=-\mathrm{0.66}$ K) than that at daytime, and the simulated seasonal LST during the daytime was underestimated in winter (MB $=-\mathrm{5.38}$ K), spring (MB $=-\mathrm{7.24}$ K), and autumn (MB $=-\mathrm{2.76}$ K), and it was slightly overestimated in summer (MB = 0.86 K). The reason for this difference between the daytime and nighttime LST simulations may be that the complex interactions of the surface energy balance during daytime are also one of the reasons that the simulated LST at the nighttime is closer to MODIS LST than that at daytime. The LST is controlled by solar radiation absorbed by the surface canopy and the ground during daytime, and these factors lead to greater uncertainty in the simulated daytime LST, whereas the solar radiation can be ignored at nighttime, and the LST is mainly affected by the downward longwave radiation (Xue et al., 2013; Zhou et al., 2015; Zhong et al., 2021). In addition, the complex terrain and strong convective activities (clouds) in the USR basin may also affect the accuracy of satellite observations (clouds and terrain may lead to increased missing observations) and the reanalysis products (complex terrain may lead to lapse rate errors in reanalysis products and affect temperature calculations).

4.1.3 Model validation using MOYDGL06^∗ SC

The fraction of simulated snow-covered area (FSCA) was also verified using the MOYDGL06^∗ FSCA data at the basin scale from 2003 to 2019 at 8 d and monthly temporal scales (Fig. 5). As shown in Fig. 5a, there was a good correlation between the simulated 8 d FSCA and the MOYDGL06^∗ 8 d FSCA (CC = 0.73). It can be seen from Fig. 5b that the simulated FSCA was overestimated, except for in October and November, especially in summer (by almost 100 %). Of course, this could reflect the intra-annual changes in the MOYDGL06^∗ FSCA values at the basin scale. Figure 5c illustrates the good agreement between the simulated and MOYDGL06^∗ FSCA at a monthly temporal scale during 2003–2018 (R²=0.56), which generally reflects the inter-annual changes in the MOYDGL06^∗ FSCA in the USR basin.

As can be seen from Figs. 6 and S3, the simulated SC reproduced the spatial distribution of the seasonal evolution of the SC in the USR basin, and it also captured the snow ablation (April–August) and accumulation (September–March) processes well. Moreover, the simulated SC also reflected the rapid melting of snow after 2 June and the rapid accumulation of snow after 8 October, but there was a difference on some days – that is, the SC was overestimated in the high mountains in summer and underestimated in the low valleys in winter. If the uncertain grids in MOYDGL06^∗ were considered (red grids), the difference would be smaller. The differences in some months may partly be caused by the bias in the forcing data, especially the biases in the temperature and precipitation. In addition, the resolution of MOYDGL06^∗ (500 m) is not consistent with that of the hydrological model (5 km), which may also lead to differences in the calculation of the average SC for the entire basin.

Figure 4 Comparison of the spatial distributions of the land surface temperature (LST) simulation and MODIS data for the USR basin from 2001 to 2018 at daytime (first and second rows) and nighttime (third and fourth rows). From left to right, the columns show winter (DJF), summer (JJA), spring (MAM), and autumn (SON).

Figure 5 Comparison of the simulated fraction of the snow-covered area (FSCA) and the MOYDGL06^∗ time series in the USR basin for the (a) 8 d FSCA, (b) multiyear mean monthly FSCA, and (c) variations in the monthly FSCA during 2003–2018.

Overall, the verification results show that the simulation results of the hydrological model are reasonable and can be used to simulate the historical and future snow and hydrological processes in the USR basin.

4.2 Future climate change

4.2.1 Evaluation of projected precipitation and temperature

Due to the systematic bias and coarse spatial resolution of the GCM, there was great uncertainty in using it directly to drive the hydrological model at the basin scale (Zhang et al., 2016; Liu et al., 2018). Thus, it was necessary to evaluate these data before conducting the hydrometeorological analysis. Four GCM datasets with bias adjustment conducted using the delta method were used as the meteorological forcing data to drive the WEB-DHM-sf hydrological model; the results are shown in Fig. 7. Although the simulated values of the single GCM in the wet season were overestimated (GFDL-ESM4 and IPSL-CM6A-LR) or underestimated (MPI-ESM1-2-HR and MRI-ESM2-0), the daily simulated discharge of each GCM was relatively consistent with the inter-annual variations in the observations during 1981–1987 (Fig. 7a). Moreover, the daily hydrograph curve of the multi-GCM ensemble mean (MEM) reproduced the discharge at JYQ well, with reasonable NSE (0.61), KGE (0.77), R² (0.64), and RB (−9.62 %) values. The monthly hydrograph curve of the MEM and observed data exhibit a good overall agreement (NSE = 0.78, KGE = 0.77, R²=0.79, and RB = −9.7 %) during 1981–1987 (Fig. 7b), and the intra-annual changes in the discharge at JYQ station were captured. In summary, the MEM was better than the single GCM; thus, the MEM was used in the analysis of the snow-related hydrological processes in the USR basin.

4.2.2 Projected changes in temperature and precipitation

As shown in Table S1, the temperature exhibits significant growth trends under SSP126 (0.2 ^∘C per decade) and SSP585 (0.7 ^∘C per decade) during 1995–2100. It also shows a consistently increasing trend during the different periods, except in the long term for SSP126 (during which time a slightly significant decrease occurs). In addition, precipitation also exhibits growth trends throughout the entire study period under SPP126 (5 mm per decade) and SSP585 (27.8 mm per decade) as well as an increasing trend in the near term, midterm, and long term, except during the reference period when it exhibited an obvious decreasing trend (−5.8 mm per decade for SSP126 and −9 mm per decade for SSP585). The climate in the USR generally undergoes a warming and wetting trend in the future under these scenarios. For the different periods, both the warming and wetting trend of the basin gradually slow down under SSP126, whereas they become much stronger under SSP585.

Figure 6 Comparison of the MODIS 8 d snow cover product (MOYDGL06*) and the simulated 8 d snow cover in the USR basin in 2017; the time interval used is 16 d.

Figure 8 and Table S2 show the relative changes in the annual precipitation and mean annual temperature in the USR during 1995–2100, with 1995–2014 as the baseline, under SSP126 and SSP585. Compared with the reference period, the warming amplitude increases in the near, mid- and long-term periods, but the rate of increase slows down and does not exceed 2 ^∘C at the end of this century under SSP126. However, warming exceeds 2 ^∘C in the future in the midterm and reaches about 6 ^∘C by the end of this century under SSP585 (Fig. 8a). The annual precipitation increases by 4.64 %–5.57 % under SSP126 and by 1.34 %–26.99 % under SSP585 during the different periods (Fig. 8b).

Figure 7 Comparison of simulated discharge of four GCMs and the observed discharge at JYQ station during 1981–1987. The red line is the multi-GCM ensemble mean discharge, and the gray line is the observed data on (a) daily and (b) monthly timescales.

The seasonal cycles of the projected temperature and precipitation in the USR exhibit different changes under different scenarios (Fig. 9). Intense warming occurs in all of the seasons under SSP126 and SSP585. Compared with the reference period, the changes in the seasonal warming rates for each period are relatively small under SSP126, but these changes are greater under SSP585, with a distinct gradient. The warming rates during each season do not exceed 2 ^∘C in the near term, and the smallest change occurs in spring under all of the scenarios; however, it exceeds 4 ^∘C in the long term under SSP585 (Fig. 9a, b). As can be seen from Fig. 9c and d, the precipitation in the other periods (compared to the reference period) is projected to increase from April to September and to decrease from October to March under all of the SSPs, and its growth amplitude gradually increases from the near term to the long term. The above results indicate that the USR is a typical monsoon region that will become warmer and wetter in the monsoon season and warmer and dryer in the non-monsoon season in the future.

Figure 8 Relative changes in the annual precipitation and mean annual temperature from 1995 to 2100. The black, blue, and green lines represent the precipitation and temperature during the reference period (1995–2014) and under the CMIP6 SSP585 (2015–2100) and SSP126 (2015–2100) scenarios, respectively. The rectangular shaded areas are the near-term (pink: 2021–2040), midterm (yellow: 2041–2060), and long-term (red: 2081–2100) periods. The shading around the lines represents the interquartile range of the data, and the upper and lower ranges are 95 % and 5 %, respectively.

4.3 Snow-related hydrological response to climate change

4.3.1 Changes in snowfall

There is a significant decreasing trend in the annual snowfall and the snow / precipitation ratio (Snow / Pre) in the USR under all of the SSPs during 1995–2100 (Fig. 10a, b). The rate of decrease in snowfall under SSP585 (−13.6 mm per decade) is much faster than that under SSP126 (−4.2 mm per decade). In the near term, snowfall is projected to decrease by more than 13 % under SSP126 and by 16 % under SSP585. By the end of the century, snowfall is reduced to about 118 mm under SSP585, which is less than half of the snowfall in the reference period, while snowfall decreases by approximately 16 % under SSP126. Similarly, Snow / Pre exhibits changes similar to those of the snowfall, and Snow / Pre is projected to decrease to approximately 10 % under SSP585 and to 20 % under SSP126 by the end of the century.

Figure 9 Seasonal changes in the temperature and precipitation in the reference period, the near term, the midterm, and the long term under (a, c) SSP126 and (b, d) SSP585. The shading around the lines represents the interquartile range of the data, and the upper and lower ranges are 95 % and 5 %, respectively. The dotted lines in panels (c) and (d) show that the precipitation in April and October during future periods is higher than that in the reference period.

Figure 10c and d show the changes in the monthly snowfall and rainfall at the basin scale. The rainfall mainly occurs from May to October (monsoon season), accounting for more than 70 % of the total annual precipitation, and the changes in rainfall and the pattern of rainfall in the future are similar to those of the total precipitation. However, snowfall exhibits a distinct bimodal pattern, with the first peak appearing in May (accounting for about 21 % of the annual snowfall) and the second occurring in October (accounting for about 14 % of the annual snowfall). Snowfall is projected to decrease by less than 10 % from November to April compared with the reference period, and there will likely be no snowfall in July and August after the midterm. Moreover, the projected Snow / Pre exhibits a consistent decrease in all of the months compared with the reference period (Fig. 10e, f). Further analysis of the seasonal variations in snowfall revealed that snowfall would be the heaviest in spring, and it is projected to decrease by about 30 % (1 %) under SSP585 (SSP126) by the end of the century. The largest reduction in snowfall occurs in the summer and autumn seasons under all of the SSPs, and it is likely to decrease by approximately 85 % (44 %) and 60 % (21 %) in these respective seasons under SSP585 (SSP126) by the end of the century (Table 1). The above results indicate that precipitation in the USR is less likely to occur in the form of snow in the future under climate warming.

Table 1 Changes and trends in seasonal snowfall in the USR basin during different periods under SSP126 and SSP585 compared to the reference period (1995–2014).

“^∗” represents a p value of less than 0.05.

Figure 10 Snowfall changes and changes in the snowfall / precipitation ratio (Snow / Pre): (a, b) annual changes in snowfall and the proportion of snowfall, (c, d) seasonal changes in snowfall and rainfall, and (e, f) the ratio of snowfall to total precipitation at seasonal scales in the USR basin during different periods under SSP126 and SSP585 compared to the reference period. The shading around the lines represents the interquartile range of the data, and the upper and lower ranges are 95 % and 5 %, respectively. “^∗” represents a p value of less than 0.05 and indicates that the trend is statistically significant.

4.3.2 Changes in FSCA

The trends and relative changes in the annual FSCA at the basin scale under the different SSPs during 1995–2100 are shown in Fig. S4 and Table S3. The annual FSCA exhibits a significantly decreasing trend at the basin scale during 1995–2100, and this trend is more obvious under SSP585 (−0.42 % yr⁻¹). Compared with the reference period, the annual FSCA is projected to increase by 4.9 % in the near term and to decrease by 8.59 % and 2.77 % in the midterm and long term, respectively, under SSP126. In contrast, the annual FSCA is projected to drastically decrease under SSP585, and it will be less than 40 % by the end of the century according to this scenario compared with the reference period. Moreover, it should be noted that the SC will continue to decrease from the near term to the midterm under SSP126, but it will begin to increase in the long term, which is significantly different from the pattern under SSP585. This indicates that future SC in the USR will decrease further under climate warming, and the severe warming projected under SSP585 would cause the SC to melt rapidly.

Figure 11 Changes in the monthly and seasonal FSCA in the USR basin during different periods under SSP126 and SSP585 compared to the reference period: (a, b) decadal FSCA changes and (c, d) FSCA changes in different periods, and (e, f) variations in seasonal FSCA during 1995–2100. The shading around the lines in panels (a)–(d) represents the interquartile range of the data, and the upper and lower ranges are 95 % and 5 %, respectively. The numbers in panels (e) and (f) represent the trends of the corresponding lines (the same color); “^∗” represents a p value of less than 0.05 and denotes that the trend is statistically significant; and the rectangular shaded areas indicate the near term (pink), midterm (yellow), and long term (red).

The results obtained from the analysis of the monthly-scale decadal FSCA from 1995 to the 2090s are shown in Fig. 11a and b. Figure 11 shows that the cycle of accumulation and ablation of the FSCA increases gradually from September to March and decreases rapidly from April to August. Under SSP585, the FSCA is projected to decrease in almost all months. However, the FSCA increases in June–November and decrease in December–April under SSP126. This pattern also occurs for the monthly-scale FSCA during the different periods (Fig. 11c, d). Similarly, the changes in the FSCA are more pronounced under SSP585 than under SSP126. For example, compared with the reference period, the FSCA becomes snow-free from July to September by the end of the century under SSP585.

Figure 12 Change in the monthly fraction of snow-covered area (FSCA) in the different elevation bands and periods under SSP126 and SSP585 compared to the reference period. The shading around the lines represents the interquartile range of the data, and the upper and lower ranges are 95 % and 5 %, respectively.

Table 2 Changes in the seasonal SWE in the USR basin under SSP126 and SSP585 compared to the reference period (1995–2014).

Figure 11e and f and Table S4 present the change trends in the seasonal FSCA, i.e., in spring (MAM), summer (JJA), autumn (SON), and winter (DJF). The results show that there is a decrease in all of the seasons during 1995–2100, but the decrease is more significant in winter and spring. Compared with the reference period, the FSCA is projected to decrease most (by about 61 % on average) in winter and spring under SSP585. Even under SSP126, the FSCA will decrease by about 15 % on average in winter and spring. However, there are differences between SSP126 and SSP585 in summer and autumn. From the near term to the long term, the FSCA is projected to increase in summer and autumn under SSP126, but this would not offset the large losses in spring and winter. This shows that the reduction in the SC in winter and spring would cause a reduction in the annual SC, which may further lead to a decrease in the snow storage in the basin, thereby affecting the amount of snowmelt during the ablation period.

From the perspective of different elevations, the SC is mainly distributed at middle–high elevations (≥4500 m), accounting for about 44 % of the basin area during the reference period, whereas less than 7 % of the SC is located at low elevations (<4500 m) (Table S3, Fig. 12). During 1995–2100, the annual FSCA at all elevations (3500–6000 m) significantly decreases under all of the SSPs, except in the area above 5500 m a.s.l. under SSP126, where no change occurs (Fig. S4). Further analysis revealed that the rate of decrease in the FSCA is faster at middle–high elevations (4500–5500 m) than at low elevations. However, the magnitude of the decrease in the FSCA is much greater at low elevations than at middle–high elevations in all of the periods compared with the reference period. For example, the FSCA in the low-elevation areas is projected to decrease by about 78 % under SSP585 and by 14 % under SSP126 by the end of the century. This illustrates that the amount of SC lost at middle–high elevations will be much greater than that lost at low elevations throughout the 21st century, but the response of the SC to climate warming at low elevations will more rapid and pronounced.

Figure 13 Seasonal changes in the SWE in the USR basin during different periods under SSP126 and SSP585 compared to the reference period.

Further statistical analysis of the monthly-scale FSCA at all elevations during all of the periods revealed that there is a consistent decrease in the FSCA at all elevations compared with the reference period (1995–2014), except in areas with elevations of 4500–5500 m a.s.l. in the near term under SSP126 (Table S3, Fig. 12). Under all of the SSPs, the FSCA at middle–high elevations (≥4500 m) is projected to decrease in almost all months from the midterm to the long term. At low elevations (<4500 m), the FSCA decreases by approximately 25 % from December to April under SSP126 by the end of the century, whereas there is a slight increase in the FSCA during June–November (with no snow cover in July and September). However, the FSCA is projected to decrease by approximately 80 % in almost all months under SSP585 by the end of the century. Moreover, the decrease in the FSCA is more pronounced in winter and spring at low elevations than at high elevations under all of the SSPs, and it further strengthens towards the end of the century. Remarkably, we found that the snowmelt starts earlier at high elevations than at low elevations under SSP585 by the end of the century. This was confirmed at elevations of 4500–5500 m a.s.l., where the starting time of the decrease in the FSCA shifts from May to March, and it continues until the end of August Correspondingly, the FSCA in this elevation range also decreases by more than 60 %. In addition, there is almost no snow at elevations below 5000 m from July to September, and the snow is very sparse and broken at high elevations. For all of the SSPs, the projected climate would not be conducive to the occurrence and development of SC.

Figure 14 Changes in the snowmelt in the USR basin during the different periods under SSP126 and SSP585 compared to the reference period: comparison of snowmelt runoff (SR) and total runoff at the basin outlet under (a–e) SSP126 and (e–h) SSP585, and comparison of snowmelt runoff and total melt (snowmelt runoff and snowmelt that seeps into the soil) under (i–l) SSP126 and (m–p) SSP58.

4.3.3 Changes in SWE

Figure 13 and Table 2 show the monthly- and seasonal-scale changes in the SWE in the different periods under SSP126 and SSP585. Compared with the reference period, the SWE is projected to decrease by about 39 % in the near term and midterm, and it slightly rebounds in the long term (about 31 %) under SSP126. However, there is a clear decrease in the SWE from the near term to the long term under SSP585, and the SWE decreases by more than 85 % by the end of the century. Similarly, the characteristics of the abovementioned changes are also reflected in the seasonal changes in the SWE (Fig. 13c, d). The SWE in the USR is the highest in spring and the lowest in autumn. Further analysis of the data revealed that the SWE decreases in all of the seasons and all of the periods under all of the SSPs, and the largest decrease occurs in summer and autumn, with a decrease of more than 91 % (30 %) under SSP585 (SSP126) by the end of the century. Furthermore, the peak in the SWE shifts from June to May under SSP585 by the middle of the century, and the SWE is almost zero from July to September. However, this does not occur under SSP126. Moreover, it can be seen that the peak in the SWE after the near term shifts from June to May under SSP585, whereas it remains in June under SSP126, and the peak value of the SWE weakens after the reference period.

Figure 15 Changes in the total runoff (blue line) and snowmelt runoff (orange line) during the different periods under SSP126 (a–b) and SSP585 (i–l). The green line is the percentage of snowmelt runoff in total runoff. The numbers represent the trends of the corresponding lines (the same color). “^∗” represents a p value of less than 0.05 and denotes that the trend is statistically significant. The shaded band around the line denotes the interquartile range, and the upper and lower ranges are 95 % and 5 %, respectively.

4.3.4 Changes in snowmelt and the snowmelt runoff

The trends in the annual total runoff, total snowmelt, and snowmelt runoff during 1995–2100 are shown in Fig. S5. The total runoff increases under all of the SSPs, and it is faster under SSP585 (4.4 mm yr⁻¹) than under SSP126 (0.31 mm yr⁻¹). However, the total snowmelt, including the snowmelt runoff and snowmelt that seeps into the soil, significantly decreases under SSP585 (−1.04 mm yr⁻¹) and SSP126 (−0.34 mm yr⁻¹). The snowmelt runoff changes very slightly and insignificantly under SSP126, whereas it significantly decreases (0.56 mm yr⁻¹) under SSP585. During the reference period, the snowmelt runoff (0.1 mm yr⁻¹) and its contribution to the total runoff significantly increased (0.3 % yr⁻¹), but the total runoff exhibited a nonsignificant downward trend (−0.13 mm yr⁻¹). The trends in the total runoff, snowmelt runoff, and the contribution of snowmelt runoff in the midterm and long term are the inverse of those during the reference period under all of the SSPs, except in the long term under SSP126 (Fig. 15). This indicates that the snowmelt and snowmelt runoff likely peak in the near term (2021–2041).

Table S5 presents the changes in the annual total runoff, total snowmelt, and snowmelt runoff from the near term to the long term under SSP126 and SSP585. In the near term, the total runoff slightly increases under SSP126 (1.07 %) and SSP585 (1.40 %). In the midterm, the total runoff increases by twice as much under SSP585 (22.86 %) compared with SSP126 (10.4 %). By the end of the century, the total runoff is projected to increase by 112.29 % under SSP585 compared with the reference period (341.0 mm), whereas there is only a small increase under SSP126 (4.7 %). The total snowmelt consistently decreases from the near term to the long term under all of the SSPs. The reduction in the snowmelt remains below 20 % under SSP126, whereas it decreases by 22.5 % and 43.01 % in the midterm and the long term under SSP585, respectively. The snowmelt runoff increases in the near term and long term and decreases in the midterm (by −20.3 %) under SSP126, while the snowmelt runoff decreases significantly from the near term to the long term under SSP585, and the projected decrease is about 75 % in the long term.

The intra-annual changes in the total runoff are very similar to those of the precipitation, and the hydrograph remains largely consistent in all periods under SSP126 and SSP585, with 60 % of the annual total runoff occurring in summer and the peak flow occurring in July (Fig. 14a–h). There is a very obvious change in the total snowmelt – that is, the peak snowmelt shifts from July to June after the reference period (Fig. 14i–p). The snowmelt runoff in summer is the largest compared with the other seasons, accounting for about 65 %–85 % of the annual snowmelt runoff, and it sharply increases in May, with a peak in June (Fig. 14i–p). From Table 3, compared with the reference period, it can be seen that there is a consistent decrease in the spring and winter total runoff in the future periods under all of the SSPs, except in the long term for SSP585, whereas total runoff increases in summer and autumn in all of the periods. The snowmelt consistently increases in winter and spring in the future periods under all of the SSPs, whereas it markedly decreases in summer and autumn. Similarly, the snowmelt runoff in spring is projected to increase in all of the periods under all of the SSPs, whereas it decreases in summer and autumn, except in the near term and long term under SSP126.

Table 3 Changes in the seasonal total runoff, snowmelt, and snowmelt runoff in the USR basin under SSP126 and SSP585 compared to the reference period (1995–2014).

Table  S6 and Fig. 15 show the contributions of the multiyear average and seasonal snowmelt runoff to the annual total runoff (SR / TR) in all of the periods under SSP126 and SSP585. The SR / TR was about 17.6 % in the reference period, and this value reaches about 19 % by the end of the century under SSP126. However, the SR / TR constantly decreases in all periods under SSP585, reaching about 2 % by the end of this century. Compared with the reference period, the SR / TR in spring was projected to increase in all of the periods, except in the long term under SSP585. Although the SR in summer contributes the most to the TR at the annual scale, the spring SR contributed the most to the spring TR at the seasonal scale, which is more obvious after the midterm under all of the SSPs due to climate warming. (Table S7). At the monthly scale (Fig. 14a–h), the SR / TR was the largest in June during the reference period, while it gradually shifted to May after the near term under SSP585 and SSP126. In short, the abovementioned results indicate that change in the SR under all of the SSPs is not the primary reason for the increase in the annual total runoff, but its contribution cannot be ignored.

5 Discussion

5.1 Impacts of snow changes on runoff

According to the climate projections of the two SSPs, the future climate of the USR basin will become warmer and wetter, leading to a continuous reduction in the snowfall and accelerated melting of the SC during 2021–2100. In addition, the reduction in the annual snowfall and SC will be greater under SSP585 than under SSP126. Moreover, it was found that there are inconsistencies in the seasonal fluctuations of the simulated result in the different periods under SSP126 and SSP585. For example, the summer and autumn SC increases in all of the periods under SSP126 compared with the reference period, but it significantly decreases under SSP585. Moreover, a snow-free summer would likely occur in the long term under SSP585. The above phenomena were also reflected in the snow storage, snowmelt, and snowmelt runoff. This is mainly due to the differences in the GCM models – that is, the warming rate of SSP58 is stronger than that of SSP126 during 1995–2100 (Fig. 8). Although SSP126 includes warming in the different periods, it has very stable changes in temperature and precipitation. SSP585, in contrast, is accompanied by strong warming and significant precipitation from the near term to the long term, which prompts a shift from more solid precipitation to liquid precipitation in each season, resulting in dramatic decreases in the snowfall, SC, and snow storage (Figs. 9, 10; Table 1).

In addition, significant differences in the response of the SC to climate warming are predicted in the different elevation ranges, and it was found that the SC in the low-elevation regions (<4500 m) will be more susceptible to climate warming than that at middle–high elevations. This is because the temperature and precipitation variations are mainly dependent on an elevation gradient. Under the same temperature increase, the warming rate would be faster in the low-elevation areas than in the middle- to high-elevation areas, particularly in areas below 4000 m a.s.l. on the TP. This is more likely to cause a large amount of snowfall to be converted into rainfall and increase the melting, which would result in a significant decrease in the amount of snowfall (Hock et al., 2019; Marty et al., 2017; Kapnick et al., 2014), especially in spring when more snowfall would occur. Moreover, the decrease in the SC in the low-elevation areas will also increase the absorption of solar radiation by the surface, increasing the surface temperature and further melting the snow (Scherrer et al., 2012). However, the response to the increase in temperature will be slower in the high-elevation areas than in the low-elevation areas, which would result in more precipitation in the high-elevation areas due to the effect of terrain uplift; thus, there would be a correspondingly greater amount of snowfall (Hock et al., 2019). This change pattern would remain at lower elevations until the end of the century; however, in the long term, under SSP585, the SC in high-elevation areas would also be greatly reduced under the effect of the large increase in temperature, and the ablation time of the SC would also be 1–2 months earlier than during the reference period. This is because the warming rate would slow down by the end of the century under SSP126, whereas the sharp temperature increase would be more intense in the long term under SSP585 (Fig. 8, Table S1).

For both SSPs, the increase in the temperature leads to a reduction in the snow storage over the projected period, resulting in a continued reduction in the total snowmelt and further causing a reduction in the snowmelt runoff, particularly in the cold season (Figs. 13, S4, S5). We found that, in the projected period, the spring and winter snowmelt and snowmelt runoff increase compared with the reference period, especially in May (Fig. S6), Moreover, the reduction in the spring snowfall leads to a reduction in the amount of snowpack that can be stored in the spring. These factors are the reasons why the meltwater and runoff in the summer are greatly reduced during the projected periods. However, the increased snowmelt in spring would cause an increase in the snowmelt runoff, which would make up for the reduction in the total runoff in spring to a certain extent and would play an important role in alleviating the drought before the monsoon period (Table S7, S8). Although the precipitation increases in autumn, the amount of snowfall is very limited due to the influence of the increase in temperature, which results in a decrease in the snowmelt and snowmelt runoff in autumn during the projected period (a significant decrease under SSP585 and a weak balance under SSP126) (Table S8, Fig. S6). Hock et al. (2019) and Nepal et al. (2021) pointed out a similar pattern of snowpack meltwater in alpine regions, and they attributed this pattern to increased rainfall due to the increase in temperature. During the projected period, the peak snowmelt in the USR shifts from July to June, but the pattern of the total runoff differs from that of the snowmelt. The total runoff still peaks in July (Figs. 14a–d and 15a–d), which is mainly influenced by the monsoon precipitation (Figs. 9 and 10c–d). Thus, the annual hydrological curve remains unchanged in the future. This finding is consistent with that of Su et al. (2016), who reported that the runoff patterns of rivers in monsoon regions would remain stable in the future. We also found that the proportion of snowmelt runoff to snowmelt would continue to decrease. In addition to climate warming, another possible explanation for this is that most of the snowmelt infiltrates into the soil and is stored in the snowpack, and it can also evaporate during runoff according to the calculations of the WEB-DHM-sf model. Under future climate warming and a continuous increase in precipitation, in view of the very small glacial area in the basin (<3 %), the influence of the snow on the hydrological processes in the basin would weaken, and precipitation would be the dominant factor affecting the runoff changes in the future (Fig. 10c, d).

5.2 Comparison with other studies

Although there are some differences in the driving data, models, and periods used in previous studies, the overall trends and patterns in the snow variables and runoff in the USR are still comparable and consistent. In this study, it was found that the total runoff exhibited an insignificant decreasing trend (−0.13 mm yr⁻¹) during the reference period, which is consistent with the recent findings of Yang et al. (2021). However, there are inconsistencies with several previous studies based on site observations that have reported a decrease in the total runoff, which is mainly due to the difference in the study periods. For example, our results also show that the total runoff increased during 1981–2018 (Fig. S7). The total runoff would also increase, which is consistent with the results of Lutz et al. (2014), Su et al. (2016), Zhao et al. (2019), and Khanal et al. (2021) for the projected period. A similar distribution pattern and seasonal changes in the SC were reported by Kraaijenbrink et al. (2021), who also found that the SC would decrease the most at low elevations in summer and winter. They argued that the shallow SC in elevated regions is more susceptible to warming, and the rapid reduction in the snowpack would lead to a shorter snow season, which would also cause a decrease in the albedo and would eventually form an elevation-dependent warming feedback cycle. The inter-annual and seasonal changes in the future snowmelt runoff are also very similar to the predictions of previous studies (Su et al., 2016; Zhao et al., 2019; Khanal et al., 2021; Kraaijenbrink et al., 2021). All of these studies found that the snowmelt would decrease and advance to spring, although there may be differences in a certain period. This can mainly be explained by the shift in the snow and rain as well as the increased snow melting rate due to warming as well as the decrease in the snowfall in the USR. The contribution of the snowmelt runoff to the total runoff during the reference period was about 17.6 %, which is between the values of 13.4 % (Yang et al., 2021) and 28.3 % (Lutz et al., 2014) reported in previous studies. This is largely due to the differences in the forcing data and hydrological model used to describe the physical snow-related processes. The contribution of the snowmelt runoff to the total runoff under SSP585 is consistent with previous studies that predicted significant decreases (Su et al., 2016; Zhao et al., 2019; Khanal et al., 2021). We also concluded that the snowmelt runoff would be important but would not be the main factor controlling the changes in the total runoff in the basin.

5.3 Uncertainties and limitations

The main uncertainties and limitations of this study were as follows. Very limited observed runoff data (for only 6 years) were used in the calibration and validation of the model. The lack of long-term runoff observations may cause uncertainties in the simulation results for the other periods. The SSP245 scenario, which is closer to the current state of development, was not considered here due to data limitations. Although ERA5 precipitation data were better than the other products in terms of the temporal and spatial distributions, most studies have pointed out that this product overestimates the precipitation during the monsoon period in high mountains (Yang et al., 2021; Khanal et al., 2021). The ERA5 product was evaluated at a single point, and the historical forcing data were not corrected. Although this effect could be reduced through model parameter calibration, it would inevitably lead to the overestimation of variables such as snow cover and runoff. For example, Khanal et al. (2021) pointed out that the cold bias in the ERA5 (even after bias correction) is the main cause of the overestimation of the winter snow cover, but an insufficient description of the snow sublimation and other physical mechanisms in the model may also lead to excessive snow cover. In addition, although the climate model data selected in the study can reflect the uncertainty range of all of the GCM datasets, there is still some uncertainty at the basin scale. Moreover, a simple bias correction method was used to conduct a secondary correction of the GCM based on the ERA5, which only ensured the consistency of the relative change trend; it did not improve the accuracy of the future predictions of the meteorological variables, such as the precipitation and temperature frequency distribution and seasonal variations, which may cause some uncertainty in the simulation results (Khanal et al., 2021; Zhao et al., 2019; Su et al., 2016). Therefore, to improve our understanding of the effects of the synergistic changes in the cryospheric components under climate change on the hydrological processes in high mountainous areas, it is necessary to enhance our understanding of the related hydrological processes through in situ experiments and to calibrate more reliable parameters to further improve the physical processes in the model.

In addition, the spatial resolution of the model has an important effect on the snow hydrological simulation in high mountainous regions with complex environments (Etchevers et al., 2001). Here, the accuracy of snow hydrological simulation between the model resolutions of 3 and 5 km was compared. The NSE, R², KGE, RB, and RMSE values between the simulated and observed discharge at the basin outlet during 1981–1987 based on the 3 km model resolution were 0.82, 0.8, 0.87, −7.93 %, and 346.12 m³ s⁻¹, respectively (Fig. S8a). It is shown that the accuracy of simulated discharge based on the 3 km resolution is slightly improved at the basin outlet compared with that of the 5 km model resolution. Meanwhile, the 8d and monthly FSCA of former was slightly improved at the basin scale compared with that of the latter according to the R², which is improved from 0.73 (8d) and 0.56 (monthly) to 0.74 and 0.59, respectively (Fig. S9a and c and Fig. 5a and c). Overall, the improvement was not obvious between the 3 km resolution and the 5 km resolution, which may be due to the following reasons: the model (WEB-DHM-sf) used in this study considers sub-grid parameterizations – that is, a model grid consists of a set of geometrically symmetrical hillslopes that are the basic hydrological units (BHUs) of the model. The topography parameters (slope and length) of the hillslope in each model grid can be calculated by a fine DEM with a 90 m resolution (Wang et al., 2009a). The sub-grid parameterization scheme can enable more detailed characterization of basin topography information such as slope, aspect, sky view factor, and elevation (Etchevers et al., 2001), thereby retaining the important information of topography-induced spatial variation in forcing inputs (e.g., the solar radiation, albedo, precipitation, and wind), which would reduce simulation bias caused by model resolution (Winstral et al., 2014; Sohrabi et al., 2019). Furthermore, the original resolution of the forcing data input to the model was above 10 km. Downscaling from low resolution to high resolution (3 and 5 km) through interpolation cannot improve the accuracy of the variables that affect the snow hydrological simulation, as coarse-scale data cannot transfer more effective information to the interpolated fine-scale data by downscaling the scale (Sohrabi et al., 2019). Nevertheless, this does not mean that it is not important to improve the model resolution, and the model resolution suitable for a specific watershed should be selected according to the characteristics of the watershed, the snow hydrological model used, and the accuracy of the forcing data.

6 Conclusions

In this study, we used multisource reanalysis data and four GCMs under two different SSPs to drive a validated cryosphere–hydrological model (WEB-DHM-sf) in order to analyze the impact of future snow changes on runoff in the USR basin. The main conclusions of this study are as follows:

1. From 1995 to 2100, the annual average temperature in the USR is projected to significantly increase under SSP126 (0.2 ^∘C per decade) and SSP585 (0.7 ^∘C per decade). The annual precipitation is also projected to significantly increase under SSP126 (5 mm per decade) and SSP585 (27.8 mm per decade). According to these simulations, the climate of the USR will become warmer and wetter in the future. On the seasonal scale, the temperature will significantly increase in all of the seasons, while the precipitation would increase in summer and autumn and decrease in the winter and spring. In the future, the USR is projected to become warmer and wetter during the monsoon season and warmer and drier during the non-monsoon season.

2. From the perspective of the entire study period, the annual and seasonal snowfall in the basin will significantly decrease, i.e., by 13 %–16 % in the near term and 16 %–49 % in the long term. Overall, the snowfall will also decrease in all of the seasons, with the greatest decrease occurring in summer and autumn, i.e., by 44 %–85 % and 21 %–59 % in the abovementioned respective seasons by the end of the century. The decrease in snowfall will directly affect the changes in the snow cover. The snow cover in the low-elevation areas will decrease significantly due to the warming, especially in winter and spring. Melting will also accelerate in the middle- to high-elevation areas due to the increase in temperature under SSP585, especially in the long term, and the ablation of snow will shift from May to March. A large amount of snowpack is melting in the USR, which will further lead to a continuous reduction in snow storage. Compared with the reference period, the snow storage in the USR in the future will exhibit negative growth in terms of both annual average (a >30 % decrease under SSP126 and an 85 % decrease under SSP585) and seasonal (a >30 % decrease in summer and autumn under SSP126 and a 90 % decrease under SSP585) changes. Thus, based on our simulations, climate warming will lead to unfavorable development of snow in the basin.

3. During the reference period, the contribution of the snowmelt runoff to the total runoff was about 17.6 %. Under the large temperature increase scenario, this contribution will become continuously smaller in the future and be about 2 % by the end of the century, leading to changes in the total runoff from a snow- and rain-dominated pattern to a rain-dominated pattern. Under the SSPs, the snow–rain runoff pattern would be maintained. The annual total runoff in the USR will increase significantly in the future, which will also increase the availability of water resources in the basin. At the seasonal scale, the total runoff will decrease in winter and spring and increase in summer and autumn, but the total annual hydrograph will remain unchanged. The increasing peak runoff may increase the risk of flooding in the future. The snowmelt runoff will significantly decrease, except in the long term under SSP126, and the meltwater peak will advance to June, with the largest proportion occurring in May. The increase in the spring snowmelt will make up for the reduction in the spring total runoff caused by the reduction in rainfall, thereby ensuring the availability of water resources in the basin during the growing season in the spring and alleviating the spring drought to a certain extent. The advance in the snowmelt will also lead to a decline in snowmelt in summer and change the pattern of snowmelt runoff in this season. These changes will have some impacts on the availability of water resources, ecosystem, and agriculture in high mountainous regions. As such, to reduce the risk of summer rain floods and spring droughts, it would be necessary to actively adjust the water resource allocation scheme to adapt to the impacts of climate change.