The Impact of Glacial Shrinkage on Future Streamflow in the Urumqi River Source Region of Eastern Tien Shan, Central Asia

Abstract

Understanding changes in runoff due to climate variations in glacier-dominated headwaters is key to managing water resources and dryland watersheds effectively and rationally. The continuous glacier shrinkage caused by climate warming has significantly impacted the water supply and ecological systems in the vast arid regions of Central Asia, attracting extensive public concern. The study results indicate an increase in total runoff at the Urumqi River source region during both the baseline (1997–2016) and mid-century (2040–2059) periods, encompassing rain, glacier meltwater, and snowmelt components. Compared to the baseline period, the temperature increases by the mid-century under the three climate scenarios (SSP1−26, SSP2−45, and SSP5−85) range from 0.98 to 1.48 °C. In this region, during the period from 1997 to 2016, glacier meltwater was the dominant component of runoff, comprising 42.10–43.79% of the total, followed by snowmelt at 29.64–30.40% and rainfall contributions of 26.56–27.49%. Additionally, glacier storage in this typical catchment responds quickly to temperature fluctuations, significantly impacting runoff. The Urumqi River source region’s runoff exhibits heightened sensitivity to these temperature shifts compared to precipitation effects. We hypothesized three glacier coverage scenarios: unchanged at 100% glaciation, reduced by half to 50%, and fully retreated to 0% glaciation. Analysis of these scenarios demonstrated that glaciers are pivotal in runoff formation. Under the SSP1−26, SSP2−45, and SSP5−85 climate scenarios, glaciers contributed additional runoff increases of 51.61%, 57.64%, and 62.07%, respectively. Generally, glaciers play a critical role in supplying water in dry areas. Thus, accurately forecasting future water resource shifts in high-altitude glacier regions is crucial for downstream water resource management and utilization.

1. Introduction

Tien Shan, often referred to as the “water tower of Asia”, is the main water supply for Central Asia’s rivers [1,2]. Around 100 million inhabitants depend on its glaciers for vital freshwater, essential to the region’s mountain hydrology [3,4,5,6]. With the ongoing effects of climate change, glaciers worldwide, particularly smaller ones, are increasingly retreating due to their sensitivity to weather changes [7,8,9,10]. In the Tianshan Mountains, known for their small- to medium-sized glaciers, this retreat greatly affects the water availability downstream. Understanding the interplay between climate change and glacier dynamics is crucial for grasping how hydrological processes in Tianshan’s arid zones respond, affecting runoff from glaciated areas [11,12,13].

In glacial basins, climate change and glacier evolution under different hydroclimatic conditions can significantly affect water flow in different basins, and the resulting runoff changes are more complex. Therefore, it is crucial to integrate hydrologic modeling with climate forecasting in order to determine the extent of runoff changes in differently sized watersheds. Many studies evaluating climate change’s impact on runoff use general circulation models (GCMs) to drive hydrological models. Nevertheless, the resolution of these models is frequently inadequate for their directly application to small-scale watershed hydrological studies [14,15]. While numerous downscaling methods have been applied to hydrological studies, their capability to capture complex spatial precipitation is comparatively limited, presenting challenges to regional hydrological simulation and prediction. Consequently, climate models with relatively fine resolution and suitable downscaling methods are imperative for studying small- to medium-scale watersheds [16,17,18]. Regional climate models (RCMs) provide a coarse resolution that, when integrated with digital elevation models and land use data, facilitates the examination of geographic characteristics in small areas. This integration is essential for simulating and predicting outcomes with hydrological models [19,20]. Although climate models alone, even with relatively high accuracy, may not suffice for hydrological simulation predictions in some small watersheds, applying statistical downscaling to the output of regional climate models helps bridge this gap [21,22]. The percentile method in statistical downscaling is frequently employed to analyze extreme meteorological conditions post-downscaling, and it is also utilized for the revision of precipitation data [23]. The method effectively characterizes diverse precipitation variations across different levels and over time. It can be utilized in ensemble studies of various emission scenarios for distinct GCMs, effectively replicating stochastic climates for driving impact assessment models and evaluating the uncertainties in climate change projections.

Runoff from the Urumqi River source region is critical for the middle and lower reaches and for Urumqi City, serving as a vital water resource that supports socio-economic growth and ecological stability within the basin. Located in a typical glaciated alpine catchment of the northwest, this area was selected for our study. Additionally, analyzing changes in this region’s runoff is essential for aiding mass balance studies and shaping regional policies within the Urumqi River basin [24,25]. Earlier research in this region was largely devoted to simulating historical runoff, with few studies exploring the impact of future climate shifts on basin runoff [26,27,28]. Additionally, statistical analysis has been applied to examine the relationship between runoff and meteorological changes using historically measured data [27,29]. This method relies on the availability of high-quality, long-term data series, limiting its ability to establish causality and often overlooking the influences of other factors, such as glaciers. To quantify different hydrological components within the glacier basin, methods such as chemical tracers, artificial neural networks, and various balance and glacio-hydrological models have been employed [15,30,31]. Notably, glacio-hydrological models are extensively utilized for reconstructing mass balance and runoff sequences and are integrated with glacier evolution and various meteorological models to project changes in future runoff [26,32,33].

This study investigates how climate and glacier interactions affect runoff variations in a typical arid region. Characterized by a glacial environment, the area is modeled using the Hydrologiska Byråns Vattenavdelning (HBV) hydrological model [34]. Initially, it simulates glacial runoff for 1997–2016 using observed weather data. Subsequently, the model projects future runoff through bias correction and statistical downscaling of the Coupled Model Intercomparison Project Phase 6 (CMIP6) meteorological forecasts under multiple models and scenarios. The objectives include: (1) Assessing runoff sensitivity to climate shifts; (2) determining runoff component ratios; and (3) forecasting future runoff under varied climate and glacier coverage scenarios.

2. Data and Methodology

2.1. Study Area

Located on the northern flank of the middle Tien Shan in northwest China’s arid territory, the Urumqi River source region spans 28.9 km² and contains seven glaciers, together covering 5.2 km². Roughly 80% of this terrain consists of bare soil, rocks, and alpine meadows. Urumqi Glacier No. 1, a summer-accumulation type within this continental climate, experiences simultaneous accumulation and ablation during summer months [35]. Owing to persistent shrinkage, it split into two smaller entities in 1993. Since establishing the Tianshan Glacier Station (TGS) in 1959, climatic, glacier, and hydrological data have been recorded. At 3405 m, the Zongkong hydrological station, located at the outlet in the Urumqi River source region, monitors precipitation and runoff. Additionally, the Daxigou Meteorological Station, operational since 1958 at 3539 a.s.l., falls under the administration of the China Meteorological Administration (CMA) as a national meteorological station (Figure 1). The Urumqi River source region has a typical continental alpine climate, and according to the data from the Daxigou meteorological station, its precipitation is mainly concentrated from May to September, accounting for 88 per cent of the total annual precipitation. The region’s average annual temperature stands at −5.1 °C. January, the coldest month, averages −15 °C, while July, the warmest, reaches 5 °C. Influences from Atlantic and Arctic Ocean currents primarily determine the water vapor sources. The colder months from October to March see scant rainfall, in contrast to the summer period from June to August, when most precipitation, predominantly snow, occurs.

2.2. Model Description

Developed by the Swedish Meteorological and Hydrological Institute (SMHI) [34], the HBV model stands as a semi-distributed hydrological framework with a robust physical basis. It partitions the study area into various elevation and vegetation zones and distinct sub-watersheds for simulations and calculations. Key components of the model include snow accumulation and melt, soil moisture accounting, groundwater management, and simple routing routines [36]. Thanks to its straightforward design and the ease of gathering input data, the HBV-Light version is widely utilized across diverse terrains such as High Mountain Asia [37,38], Europe [39,40], and North America [41].

2.2.1. Snow and Glacier Routine

Within this module, rainfall is distinguished between rain and snow based on the temperature exceeding or falling below a threshold temperature, TT (°C). Additionally, the snow correction factor, SFCF (−) (where “(−)” indicates the parameter is dimensionless and unitless), is employed to rectify systematic errors arising from snowfall measurement and unaccounted snow evaporation in the model. Snowmelt, $M_{snow}$ (mm d⁻¹), is computed by a degree day method using the degree day factor, ${CFMAX}_{snow}$ (mm d⁻¹ °C⁻¹) (Equation (1)).

$$M_{snow}={CFMAX}_{snow}·\left(T\left(t\right)-TT\right)$$

(1)

Glacier melt, $M_{ice}$ (mm d⁻¹), is calculated by the same method. “t” represents the time step of the model simulation, with the unit of t being days. In addition, in order to eliminate the errors caused by snow and glacier ice albedo and topography, a glacier correction factor $C_g$ (−) and a topography correction factor $C_a$ (−) are introduced as shown in Equation (2) [42,43].

$$M_{ice}={CFMAX}_{ice}·\left(T\left(t\right)-TT\right)·C_g·C_a$$

(2)

When the temperature $T\left(t\right)$ decreases to less than $TT$, the snow melting process is interrupted, the liquid water is re-frozen, and the freezing amount $R$ (mm d⁻¹) is calculated by the re-freezing coefficient $CFR$ (−) (Equation (3)).

$$R=CFR·M_{ice}$$

(3)

For the model run, glacier coverage data spanning the base period (1997–2016) are essential. These data, derived from continuous time series of glacier coverage information obtained through remote sensing imagery, serve as the foundation. The model incorporates three distinct glacier coverage scenarios to simulate future runoff. Firstly, two extreme scenarios are considered: one where glacier area remains intact (100%) and another where glaciers vanish entirely (0%). Although these extremes are unlikely in reality, they are included to evaluate how glaciers’ presence or absence influences future runoff modeling. Additionally, a third scenario is adopted, wherein the glacier area is halved compared to the end of the base period (50%). This scenario is based on scholarly simulations of glaciers in the Urumqi River source region [44]. Δh−parameterization is an empirical approach used to simulate glacier retreat [45]. By calculating changes in glacier thickness and area over time, this method updates the glacier’s geometry, thereby simulating the impact of glacier retreat on hydrological processes. Additionally, the study area is divided into 25 elevation bands, each with an interval of 41.68 m. Based on field surveys, land use cover is categorized into two types (grassland and bare land), and their area proportions in each elevation band are calculated and used as input data. To enhance the accuracy of glacier change simulations, the study area is further subdivided into 100 elevation bands using Δh−parameterization. The glacier proportion is calculated in each elevation band of the study area, as well as the water equivalent. The water equivalent of the glacier can be determined using the thickness data of the glacier.

2.2.2. Soil Routine

In this procedure, groundwater recharge and actual evaporation are modeled as functions of actual water storage. The ratio of total amount to soil ($Recharge$, mm d⁻¹) and groundwater flux ($P\left(t\right)$, mm d⁻¹) is represented by a power function relationship involving the current water content of the soil box ($SM\left(t\right)$, mm) and its maximum value ($FC$, mm) (Equation (4)). BETA determines the relative contribution of rain or snowmelt to runoff.

$$\frac{Recharge}{P\left(t\right)}={\left\{\frac{SM\left(t\right)}{FC}\right\}}^{BETA}$$

(4)

If $SM\left(t\right)$ is less than $FC·LP$ (LP, −), the actual evaporation ($E_{act}$) equals the potential evaporation ($E_{pot}$) multiplied by the ratio of SM(t) to (FC∙LP) (Equation (5)). LP is the soil moisture threshold for actual evapotranspiration (AET) to reach potential evapotranspiration (PET). Conversely, actual evaporation ($E_{act}$) is equivalent to the potential evaporation ($E_{pot}$) if $SM\left(t\right)$ is above $FC·LP$.

$$E_{act}=E_{pot}·min\left\{\frac{SM\left(t\right)}{FC·LP},1\right\}$$

(5)

2.2.3. Response and Routing Routine

In the HBV-Light model, the soil is partitioned into upper and lower reservoirs, with PERC (mm d⁻¹) specifying the percolation rate from the upper soil reservoir (SUZ, mm) to the lower soil reservoir (SLZ, mm). The runoff from the soil reservoir ($Q_{GW}$) (mm d⁻¹) is calculated by a linear equation (Equation (6)), where $K_0$ (d⁻¹), $K_1$ (d⁻¹), and $K_2$(d⁻¹) are all regression coefficients.

$$Q_{GW}\left(n\right)=K_2·SLZ+K_1·SUZ+K_0·\max \left(SUZ-SLZ,0\right)$$

(6)

Finally, this runoff is converted to the simulated runoff ($Q_{sim}$) (mm d⁻¹) using a triangular weighting function defined by the parameter $MAXBAS$ (−) (Equations (7) and (8)). During the model simulation, $u$ represents the time step.

$$Q_{sim}\left(n\right)=\sum _{i=1}^{MAXBAS}c\left(i\right)·Q_{GW}\left(n-i+1\right)$$

(7)

$$c\left(i\right)=\underset{i-1}{\overset{i}{\int }}\frac{2}{MAXBAS}-\left|u-\frac{MAXBAS}{2}\right|·\frac{4}{{MAXBAS}^2}d(u)$$

(8)

The Nash–Sutcliffe coefficient of efficiency ($R_{eff}$) and the coefficient of determination ($R^2$) are used to evaluate the simulation results (Equations (9) and (10)). The closer $R_{eff}$ and $R^2$ are to 1, the better the simulation results of the model are.

$$R_{eff}=1-\frac{\sum {\left(Q_o-Q_s\right)}^2}{\sum {\left(Q_o-\overline{Q_s}\right)}^2}$$

(9)

$$R^2=\frac{{\left(\left(\sum \left(Q_o-\overline{Q_o}\right)·\left(Q_s-\overline{Q_s}\right)\right)\right)}^2}{{\sum \left(Q_o-\overline{Q_o}\right)}^2·{\sum \left(Q_s-\overline{Q_s}\right)}^2}$$

(10)

In the formula, $Q_s$ is the simulated runoff; $Q_o$ is the observed runoff; and $\overline{Q_s}$ and $\overline{Q_o}$ are the mean of the observed and simulated values, respectively.

2.3. Data Description

2.3.1. CMIP6 Climate Scenario

Presently, the Coupled Models Working Group (WGCM) of the World Climate Research Program (WCRP) has initiated CMIP6. This is the most extensive climate model comparison program to date, featuring the highest number of models, the most diverse numerical experiments, and the largest volume of simulation data submitted over the past 20 years [46]. For this study, five CMIP6 datasets were selected: CESM2, CESM2-WACCM, CanESM5, BCC-CSM2-MR, and NorESM2-LM [47]. These models were chosen for their exceptional performance in the Tien Shan region of Central Asia, where their Taylor Skill Scores (TSS) were 0.49, 0.56, 0.40, 0.63, and 0.59, respectively [48]. Furthermore, the scenarios SSP1−26 (sustainability), SSP2−45 (middle of the road), and SSP5−85 (fossil fuel) were utilized to run the HBV-Light model [49].

2.3.2. Percentile Statistical Downscaling

Based on the measured meteorological data, the daily precipitation from CMIP6 is adjusted using the percentile scaling method. This approach achieves precipitation correction by establishing a mathematical function relationship between the CMIP6−simulated data and the observed daily precipitation during the reference period, as depicted in Equation (11) [50].

$$P_{adj,d}=P_{CMIP6,d}·\left(\frac{P_{obs,Q}}{P_{CMIP6,Q}}\right)$$

(11)

$P_{adj,d}$ is the corrected daily precipitation, and $P_{CMIPs,d}$ is the daily CMIP6−simulated precipitation. $P_{CMIP6,Q}$ and $P_{obs,Q}$ are the monthly percentiles of the CMIP6−simulated precipitation and the observed precipitation, respectively. For enhanced precision, percentiles such as the 10th, 20th, 30th, 40th, 50th, 60th, 70th, 80th, 85th, 90th, 93rd, 95th, 96th, 97th, 98th, 99th, and 100th are computed for each month in this study. Drawing from the distribution of daily precipitation across an extended time series [51,52], we employ a refined scheme with sparse percentile values corresponding to small precipitation and dense percentile values corresponding to large precipitation.

2.3.3. Measured Data

Meteorological data, including daily temperature and precipitation, were sourced from the Daxigou meteorological station (

Table 1. Information notes on Daxigou meteorological station.

Meteorological Station	Elevation (m a.s.l.)	Latitude (°N)	Longitude (°E)	Air Pressure	Air Temperature	Precipitation	Radiation	Relative Humidity	Snow Depth	Wind Speed and Direction
Daxigou	3539	43.113	86.843	Vaisala PTB 110	Campbell Scientific CS215	Geonor T200B	Kipp&Zonen CNR4	Campbell Scientific CS215	Campbell Scientific SR50A	Young 05103

). The station is equipped with a Geonor T-200B auto-weighing precipitation gauge (Campbell, Mississauga, ON, Canada), a Young 05103 wind monitor (R.M. Young, Traverse, MI, USA), a CNR4 pyranometer (Kipp & Zonen, Utrecht, The Netherlands), an HC2-S3 temperature and humidity sensor (Campbell), and an SR50A snow depth sensor (Campbell Scientific, USA). It was designed in cooperation with the Norwegian Meteorological Institute and the Norwegian Geotechnical Institute. Daily runoff data from the Zongkong hydrological gauging station (

Table 2. Information notes on Zongkong hydrological gauging station.

Hydrometeorological Stations	Elevation (m a.s.l.)	Latitude (°N)	Longitude (°E)	Catchment Area (km2)	Glaciation (%)	Measuring Principle
Zongkong	3405	43.117	86.870	28.9	12	Planar microstrip array antenna CW + PCR

) were used to calibrate and validate the HBV-Light model simulations. The entire catchment area is 18% covered by glaciers, with alpine meadows and barren land as the underlying surfaces. Runoff data were measured using a TTC-RD-600S radar flow meter.

2.3.4. Spatial Data

To fulfill the input requirements for glacier coverage during the baseline period (1997–2016), remotely sensed imagery and pre-existing datasets were selected as primary data sources. Glacier and snow information within the study area was extracted through supervised classification based on remotely sensed imagery. Land use data were obtained from the China Annual Land Cover Dataset, with a resolution of 30 m [53]. This dataset classifies land use into nine categories: cropland, forest, shrub, grassland, barren, impervious, water, snow/ice, and wetland. Remote sensing images primarily originated from Landsat 8 OLI_TIRS and Landsat 4–5 TM satellites (

Table 3. Detailed information on the satellites Landsat 8 OLI_TIRS and Landsat 4-5 TM.

Satellite	Sensor	Launch Date	Resolution	Spectral Range	Bit Depth	Repeat Cycle	Data Format	Sampling Method
Landsat 8	OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor)	11 February 2013	30 m	OLI: 0.433–2.355 μm (11 bands), TIRS: 10.60–12.51 μm (2 bands)	12 bits	16 days	GeoTIFF, JPEG	triple convolution
Landsat 4–5 TM	TM (Thematic Mapper)	Landsat 4: 16 July 1982; Landsat 5: 1 March 1984	30 m	TM: 0.45–12.50 μm (7 visible and infrared bands), TIRS: 10.40–12.50 μm (1 band)	8 bits (multispectral) or 6 bits (thermal)	16 days	GeoTIFF, JPEG	triple convolution

), with long-term continuous imagery from both satellites downloaded from the USGS (https://www.usgs.gov/) (last access: 6 October 2023). The image acquisition period was from 1997 to 2016, selecting images from June to September each year (when glaciers are in the ablation season and the boundaries are clearer), with cloud cover less than 5%.

3. Model Preparation and Application

3.1. Meteorological Data Downscaling

The study selected the years 2007–2016 to demonstrate the reliability of percentile correction, using 2007–2011 for calibration and 2012–2016 for validation. Percentiles ranging from the 10th to the 100th were utilized for monthly precipitation in this analysis. The observed actual precipitation and the percentile-corrected precipitation during both periods exhibited good agreement, with R² values of 0.99 and 0.97, respectively. Specific results are depicted in Figure 2. The line in the graph represents the total precipitation over the years 2007–2011 and 2012–2016.

Studies on the spatio-temporal distribution of air temperature and precipitation in the CMIP6 model have shown that the predicted values are consistent with meteorological observations in the Tianshan Mountains. However, the subgrid-scale features and dynamics need further attention. The climate model’s meteorological outputs were not directly used in our hydrological model due to their lack of representation of local subgrid-scale features and dynamics [54]. Therefore, a relatively simple statistical downscaling method, nonlinear regression, was employed to establish the relationship between observations from the Daxigou weather station and predicted climate variables from CMIP6 at four neighboring grid points for the period of 1961–2000. Daily temperature and precipitation were calculated using the approach proposed in ref. [55].

$$y={\beta }_1\times x_1+{\beta }_2\times x_2+{\beta }_3\times x_3+{\beta }_4\times x_4+b$$

(12)

The temperature and precipitation after statistical downscaling, denoted as y, were derived from the projected temperature or precipitation time series using CMIP6, represented by $x_i$ (i = 1, 2, 3, 4). The values of meteorological variables (daily temperature (°C) precipitation (mm)) represented by $x_i$ were obtained from the four grid points adjacent to the Daxigou weather station within the CMIP6 climate model. The corresponding parameters ${\beta }_i$ (i = 1, 2, 3, 4) are detailed in

Table 4. Parameters defining relationships between observed and projected data. T: monthly mean air temperature (°C a⁻¹); P: monthly precipitation (mm a⁻¹); R²: coefficient of determination.

Variable	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{\beta }}_3$	${\mathit{\beta }}_4$	b	R2
T (°C)	0.14	−0.16	0.21	0.81	−0.48	0.963
P (mm)	1.70	−0.68	−0.64	0.24	0.88	0.996

. Given the presence of spatial correlation, future predictions of temperature and precipitation were significantly influenced by the surrounding topographic conditions. The temperature and precipitation time series for 2012–2016, obtained through CMIP6 downscaling, closely matched observations from the Daxigou meteorological station, with correlation coefficients of 0.96 and 0.99, respectively. This suggests the effective performance of the statistical relationships. A comparison between the downscaled temperature and precipitation data and the observed meteorological data is illustrated in Figure 3. In summary, statistical downscaling of meteorological model outputs can make the temperature and precipitation data more consistent with regional characteristics and dynamic changes, providing an important data basis for subsequent watershed runoff simulation.

3.2. Changes in Projected Climate

Figure 4 and Figure 5, respectively, depict the multi-year monthly averages for temperature and precipitation during 2040–2059, as well as the long-term monthly temperature and precipitation variation curves for 2040–2059. Compared to the baseline period (1997–2016), all three climate scenarios showed upward trends in both metrics. The average temperature increase ranged from 0.98 to 1.48 °C, with precipitation increases observed in SSP1−26, SSP2−45, and SSP5−85, respectively. The SSP5−85 scenario recorded the highest summer temperatures, followed by SSP2−45 and SSP1−26. Both the spring and summer precipitation were higher in the mid-century compared to the baseline period.

Figure 6 illustrates the temperature and precipitation changes projected by CMIP6 under three climate scenarios (SSP1−26, SSP2−45, and SSP5−85) relative to the baseline (1997–2016). Generally, both temperature and precipitation exhibited increasing trends. A comparison of the monthly changes showed that the largest increases in temperature and precipitation occurred in the spring and the smallest in the autumn. These two trends remained relatively consistent. In terms of temperature, winter showed the largest increase, with January’s increase exceeding 3 °C in all three climate scenarios compared to the lowest increases in summer and autumn. In comparison to the baseline period, the temperature in September exhibited declines of −1.29 °C, −0.93 °C, and −0.46 °C under the three climate scenarios, respectively. As for precipitation, there was an overall increase in precipitation across almost all months, with the most significant increase occurring in summer, with an average increase of 90.6%. All in all, compared with the baseline period (1997–2016), the temperature and precipitation in the three climate scenarios under CMIP6 showed increasing trends. Although the temperature also showed an increasing trend, the difference in the melting period was small, but the difference in winter was larger. The general rise in temperature and precipitation aligned with findings from previous studies [37,56].

3.3. Calibration and Validation

The study utilized daily runoff data from 1997 to 2006 for calibrating the HBV-Light model, with the subsequent selection of 2007 to 2016 as the model validation period. Figure 7 illustrates the simulated and observed runoff for both the calibration period (1997–2006) and the validation period (2007–2016). As depicted in the figure, the daily runoff trends aligned closely with precipitation variability during both the calibration and validation periods. Additionally, given precipitation’s crucial role as a runoff component, the annual peak of runoff corresponded to the peak of precipitation. The peak runoff was based on daily values. In addition, daily runoff is concentrated in May to September each year, accounting for more than 95% of the annual runoff, and the remaining months of the runoff account for a very small proportion.

Figure 7. Daily runoff from observations and simulations in the calibration period (1997–2006) (a) and validation period (2007–2016) (b) are shown by red and black lines, and precipitation is shown by gray lines. The performance of the HBV-Light model during calibration and validation was assessed using the Nash–Sutcliffe efficiency and the coefficient of determination. Both metrics ranged from 0.68 to 0.79 for the respective periods. Specifically, the Nash–Sutcliffe efficiency values were 0.69 and 0.68 during calibration and 0.79 and 0.78 during validation. Detailed results are shown in

Table 5. Modeling performance for calibration and validation periods: observed and simulated runoff, Nash–Sutcliffe efficiency (Reff), and coefficient of determination (R²).

Calibration					Validation
Period	Qobs	Qsim	R2	Reff	Period	Qobs	Qsim	R2	Reff
1997–2006	1734.53	1651.16	0.69	0.68	2007–2016	1966.93	1762.31	0.79	0.78

. These results indicate that the model can effectively simulate daily runoff variations in the Urumqi River source region during both the calibration and validation periods, demonstrating good applicability of the HBV-Light model in glacier areas of small watersheds.

Figure 7. Daily runoff from observations and simulations in the calibration period (1997–2006) (a) and validation period (2007–2016) (b) are shown by red and black lines, and precipitation is shown by gray lines. The performance of the HBV-Light model during calibration and validation was assessed using the Nash–Sutcliffe efficiency and the coefficient of determination. Both metrics ranged from 0.68 to 0.79 for the respective periods. Specifically, the Nash–Sutcliffe efficiency values were 0.69 and 0.68 during calibration and 0.79 and 0.78 during validation. Detailed results are shown in Table 5. These results indicate that the model can effectively simulate daily runoff variations in the Urumqi River source region during both the calibration and validation periods, demonstrating good applicability of the HBV-Light model in glacier areas of small watersheds.

Determining parameter values is a crucial aspect of hydrological modeling and serves as a pivotal factor in constructing a hydrological model that accurately fits the data [57]. The sensitivity of the parameters utilized in the HBV-Light model was analyzed, and suitable values for the main parameters essential for model operation were obtained. This determination took into account previous studies conducted in the study area and the expertise gained from parameter calibration, as depicted in

Table 6. Description of the optimization parameter set of HBV-Light and its meaning.

Parameter	Description	Value	Units
Temperature and precipitation gradients
PCALT	Precipitation gradient	3	%/100 m
TCALT	Temperature gradient	0.7	°C/100 m
Snow routine
TT	threshold temperature	2	°C
CFMAX	degree-day factor of snow	5.19	mm·(°C·d) −1
SP	seasonal variability	0.1	−
SFCF	snowfall correction factor (−)	0.738	−
CFR	refreezing coefficient (−)	0.05	−
CWH	water holding capacity (−)	0.0002	−
Soil routine
FC	maximum soil moisture storage	500	mm
LP	soil moisture threshold for AET to reach PET	0.22	−
BETA	parameter determining runoff contribution from rain or snowmelt	2.40	−
Glacier routine
CFGlacier	glacier correction factor	0.48	−
CFSlope	slope correction factor	2.45	−
KSI	snow to ice conversion factor	0.001	d−1
KGmin	minimum outflow coefficient	0.18	d−1
dKG	outflow coefficient range	0.0001	d−1
AG	calibration parameter	0.75	mm−1
Response routine
PERC	threshold parameter	0.75	mm/d
UZL	threshold parameter	50	mm
K0	storage coefficient 0	0.19	d−1
K1	storage coefficient 1	0.19	d−1
K2	storage coefficient 2	0.0001	d−1
Routing routine
MAXBAS	length of triangular weighting function	1	−

.

The parameters PCALT, TCALT, TT, CFMAX, and CFGlacier were determined by referencing existing studies [58,59,60], while other parameters were established through a combination of manually defining parameter ranges and employing the Monte Carlo calibrate determination method. The calibration process involved initially setting parameter ranges based on comprehensive watershed observation data and prior research, followed by optimizing parameter combinations using the Monte Carlo calibrate determination method [61]. During the correction process, TT, CFMAX, CFGlacier, SFCF, CWH, and CFR exhibited the highest sensitivity, followed by FC, LP, and CFSlope, with BETA and MAXBAS being the least sensitive. Furthermore, PCALT, TCALT, SFCF, FC, and CFMAX primarily influenced flow production, while K₀, K_1, and K₂ predominantly affected the flow process, with higher values leading to increased peak flow and faster recession. MAXBAS values reflected the speed of both fast and slow runoff responses.

The model still has shortcomings in the simulation of peak runoff in individual years. If the measured value is inconsistent with the simulated value, the analysis indicates that the main reason is that the parameter values used in the model are fixed values and their dynamic changes are not considered. With global warming, changes in the underlying surface conditions in the study area will also affect the parameters used in the model. In addition, the relevant parameters listed in Table 6 will change with climate changes, such as CFMAX and TT. The conceptual hydrological model has some defects in the process of regional hydrological simulation. It does not fully consider the spatial heterogeneity of precipitation or the underlying surface elements (land, vegetation, soil and topography, etc.), and cannot adequately reflect the regional hydrological process and runoff changes under the background of climate change. Moreover, under the influence of climate change, the landscape type of the study area is likely to change.

4. Result and Discussion

4.1. Components of Runoff

In the study area, runoff is classified into glacier runoff, snowmelt runoff, and rain runoff based on their formation processes, with baseflow representing a form of mutual transformation and recharge among these three types of runoffs. Generally, glacial meltwater constitutes nearly half of the total runoff, with rain runoff and snowmelt runoff following. There is minimal variation in each runoff component between the calibration and validation periods. As illustrated in Figure 8, runoff in the Urumqi River source region is primarily concentrated between May and September. From 2007 to 2016, rain runoff from May to June accounts for the largest proportion at 64.34%, followed by snowmelt runoff at 19.57% and glacier meltwater runoff at 16.07%. From July to August, glacier meltwater increases, accounting for 52.50% of the total runoff, with snowmelt runoff and rain runoff contributing 25.63% and 21.87%, respectively. Notably, glacier meltwater runoff peaks in August, while July sees the highest proportions of snowmelt runoff and rain runoff. The overall runoff in the basin steadily increases from May to August, followed by a sharp decline in September. Before May and after September, runoff is mainly contributed by snowmelt runoff and rain runoff due to the influence of westerly airflow. Compared to summer, the runoff contribution in winter is much smaller because the study area experiences maximum snowfall in summer, with only a small amount of snowfall in winter.

Table 7. The share of the three runoff components (rain, glacier meltwater, snowmelt) in the total runoff during the hydrological simulation.

Period	Contribution (%)
	Rain	Glacier Meltwater	Snowmelt
1997–2006 (calibration period)	27.49%	42.10%	30.40%
2007–2016 (validation period)	26.56%	43.79%	29.64%

presents the runoff composition in the Urumqi River source region from 1997 to 2016, encompassing both the calibration and validation periods. The table further analyzes the recharge ratio of glacial meltwater, snowmelt, and rain runoff. Notably, the proportion of glacier meltwater runoff was the largest, constituting 42.10% and 43.79% for the two respective periods. Snowmelt runoff followed, accounting for 30.40% and 29.64% of the total runoff during these periods. Meanwhile, the proportion of rain runoff was comparatively smaller, at 27.49% and 26.56%, respectively.

Figure 9 illustrates the changes in total runoff, rain runoff, snowmelt runoff, and glacier meltwater in the Urumqi River source region during the model simulation period (1997–2016). The glaciers in the study area are typical representatives of mountain glaciers, with high summer precipitation and glacier ablation occurring in July and August. The majority of the runoff (95%) occurs during the spring (March–May) and summer (June–August) months, with minimal runoff observed in the autumn (September–November) and winter (December–February). The presence of autumn runoff is limited to the years 2008–2010 and 2013. Spring precipitation is dominated by solid precipitation, with snowmelt runoff representing the primary contributor to runoff. This runoff peaks in the watershed as time progresses into the summer months, when precipitation gradually increases and glacial ablation increases. Glacial runoff plays a crucial role in the overall runoff volume. In years characterized by abundant precipitation but insufficient thermal energy, glacial melt is minimal, leading to increased glacier accumulation and a subsequent reduction in glacial meltwater. Conversely, during years of limited precipitation, the available thermal energy for glacier melt intensifies, resulting in enhanced glacial melting and the release of substantial amounts of glacial meltwater. This regulatory function of glaciers on runoff has earned them the designation of “solid reservoirs”. During arid and rainless summers, the presence of glacial meltwater ensures the continuity of river flow, preventing river channels from drying up. Additionally, in the event of heavy rainfall, a portion of the precipitation is stored on the glacier in solid form, thereby reducing the risk of flooding.

4.2. Future Runoff Cycle

Figure 10 depicts the monthly average simulation results for the baseline period (1997–2016) and future climate scenarios (2040–2059) under three glacier cover scenarios relative to the glacier area in 2016. The glacier cover scenarios consist of: no change in glacier area (100% glacierization), a 50% retreat of glacier area (50% glacierization), and complete loss of glacier (0% glacierization).

Compared to the baseline glacier area in 2016, watershed runoff shows an increasing trend from 2040 to 2059 when assuming a constant glacier area. The largest increase is observed under the SSP2−45 climate scenario, with a rise of 34.83%, which aligns with expectations in a warming climate with rising temperatures and precipitation. However, this scenario is hypothetical in this study and is unlikely to occur in the future. Similarly, with 50% glacier coverage, the overall basin runoff trend decreases. Only the SSP2−45 climate scenario shows an increasing trend of 4.97%, while the other two climate scenarios, SSP1−26 and SSP5−85, exhibit decreasing trends of −15.70% and −19.44%, respectively. Comparing the 100% and 50% glacier cover scenarios, it becomes evident that changes in the glacier area due to climate change have a more substantial impact on runoff when keeping the climate scenario constant, especially when the percentage of glaciers in the study area is relatively large. However, 100% glacier glaciation is only a hypothetical result of our study and may not actually occur in the future. This is because, in the wake of such a drastic climate change, the glacier area could shrink dramatically. If the glacier area completely retreats from 2040 to 2059, glacier runoff decreases in all climate scenarios, with a 50.03% reduction in runoff in the SSP5−85 climate scenario and the least reduction in the SSP2−45 climate scenario, as well as a 23.26% decrease in watershed runoff. The specific results are presented in

Table 8. Runoff variations relative to the baseline period (1997–2016) for the future period (2040–2059) under three SSP and Glacier scenarios.

Climate Scenarios	Glacier Scenarios
	100% Glaciation	50% Glaciation	0% Glaciation
SSP1−26	10.58%	−15.70%	−41.03%
SSP2−45	34.38%	4.97%	−23.26%
SSP5−85	12.04%	−19.44%	−50.03%

.

The Urumqi River source region experiences the majority of its runoff from May to September. The runoff peaks almost always occur in July each year, coinciding with intense glacier melt and the highest precipitation levels. Considering the baseline runoff (1997–2016), we find that its peak runoff is comparable to the peak runoff under SSP1−26 (100% glaciation), SSP2−45 (50% glaciation), and SSP5−85 (100% glaciation). By comparison, it can be observed that the runoff under SSP1−26 and SSP2−45 is primarily concentrated in June and July, while the runoff under SSP5−85 is mainly concentrated in July and August. Furthermore, the differences in runoff volumes under the three glaciation scenarios are quite significant, indicating that the presence of glaciers greatly influences the downstream runoff volumes in glacier-fed watersheds. Figure 11 illustrates the variations in total runoff and its three components in the study area from 2040 to 2059. The analysis reveals that, due to glacier retreat, rain runoff constitutes the majority of the total runoff during this period, while snowmelt runoff and glacier meltwater runoff contribute smaller proportions. This observation is consistent with the trend of increasing warmth and humidity in the arid regions of Northwest China as a result of global warming.

Compared with the two assumptions of 100% glaciation and 50% glaciation, the glacier runoff under the three climate scenarios of 0% glaciation decreases sharply. In the context of global climate warming, the trend of glacier melting is inevitable. As the source of the Urumqi River Basin, the Urumqi River source region’s total runoff is closely related to climate change and glacier area changes. This relationship is also important in the Urumqi River Basin [62,63]. In order to verify the reliability and completeness of the simulation results, the results of the future runoff simulation of this study were compared with the results of Rounce 2022 [64] (Figure 12). The latter simulated global glacier changes using a coupled PyGEM (python glacier evolution model)-OGGM (The Open Global Glacier Model) model. The reduction in runoff due to glacier retreat was offset by recharge due to climate change compared to glacier retreat. A comparative analysis revealed that the total runoff in the Urumqi River source region exhibited an increasing trend from 2040 to 2059 under the three climate scenarios (SSP1−26, Sustainability; SSP2−45, Middle of the road; SSP5−85, Fossil-fuelled). However, the latter study indicated a contrasting decreasing trend in total runoff. Additionally, studies conducted by scholars in the research area indicate that runoff changes peaked around 2013 [65]. Combining these findings with the simulation results, it can be concluded that, despite the runoff in the Urumqi River source region having already reached its peak, future runoff is expected to continue increasing due to the impacts of climate change. The continuous retreat of glaciers will not offset the trend of increasing runoff caused by increased precipitation. By the mid-century (2040–2059) period, total runoff in the study area will be primarily composed of rain runoff, with snowmelt runoff and glacier meltwater runoff accounting for a smaller proportion. This result is consistent with findings from other studies [65,66]. At the headwaters of glaciated river basins, the combination of increased glacier melts and lengthening of the melt season provides an important basis for production development of downstream water uses. Urumqi has more than 4 million residents, and more than 70% of its total water consumption relies on the water source of the Urumqi River. Analysis of runoff components and future runoff changes in the source areas of rivers in the Tianshan region is crucial to water resources and disaster management.

4.3. Sensitivity Analysis

To evaluate climate change’s impact on hydrology, we performed a sensitivity analysis on hydrological simulations, focusing on variations in temperature and precipitation. Leveraging insights from the IPCC’s Sixth Assessment Report and prior studies [67,68], we created simulations with various temperature and precipitation parameters: T + 1.5 °C, T − 1.5 °C, P + 20%, P − 20%, T + 1.5 °C and P + 20%, T − 1.5 °C and P + 20%, T + 1.5 °C and P − 20%, and T − 1.5 °C and P − 20%, totaling eight scenarios.

As shown in Figure 13a,b, both temperature and precipitation positively correlate with runoff. Additionally, the combined scenario analysis indicates that temperature predominantly influences runoff variation (Figure 13c). Specifically, watersheds with glaciers as the primary recharge source are most responsive to temperature changes, with precipitation being the secondary factor.

The detailed variations are presented in

Table 9. Runoff change percentages in various sensitivity tests.

Parameter	Temperature		Precipitation		Temperature and Precipitation
Experiment	T + 1.5 °C	T − 1.5 °C	P + 20%	P − 20%	T + 1.5 °C P + 20%	T − 1.5 °C P + 20%	T + 1.5 °C P − 20%	T − 1.5 °C P − 20%
Runoff change (%)	33.10%	−41.14%	1.81%	−15.50%	45.67%	−35.81%	22.06%	−46.50%

. When the temperature increases by 1.5 °C (T + 1.5 °C), runoff increases (+22.06%), even with a 20% decrease in precipitation (P − 20%). Conversely, when the temperature decreases by 1.5 °C (T−1.5 °C), runoff decreases (−35.81%), even with a 20% increase in precipitation (P + 20%). Notably, the two scenarios resulting in the most significant increase and decrease in runoff, T + 1.5 °C and P + 20% and T − 1.5 °C and P − 20%, exhibit runoff changes of 45.67% and 46.50%, respectively. Consequently, the quality of the temperature dataset has the most substantial impact on the uncertainty of watershed runoff simulation when using this model. An increase in temperature accelerates glacier melt, resulting in more glacier meltwater runoff, and also increases liquid precipitation, which further contributes to runoff. Conversely, a decrease in temperature slows down glacier melt, lowers the equilibrium line altitude, puts the glacier in a state of overall mass accumulation, reduces glacier meltwater runoff, and increases solid precipitation, ultimately leading to a reduction in basin runoff. When predicting future runoff changes, the parameters used in models are determined based on historical data. With climate warming, the underlying surface conditions within the basin will change, causing the parameters used in the models (e.g., TT, CFMAX, PCALT, and TCALT) to change as well, thereby affecting the simulation results. Furthermore, the study assumes three glaciation scenarios (100% glaciation, 50% glaciation, and 0% glaciation) to analyze the impact of glacier presence on basin runoff. However, glacier melt and retreat are processes involving complex dynamics. Given the small size of the study area, the Urumqi River Source Region, and the limited variety of underlying surface types, this study categorizes the underlying surface types into two categories: glacier areas and bare rock areas. This categorization overlooks the evolutionary process of the underlying surface in the study area over time, which requires further research and discussion.

Temperature and precipitation are important factors affecting the variation in glaciers [69,70]. The Urumqi River source region primarily receives solid precipitation, which significantly supplements the glaciers. Due to their relatively small areas, these glaciers are highly sensitive to changes in temperature and precipitation. Approximately 40% of the total runoff in this region is derived from glacial meltwater. To analyze the sensitivity of glacier runoff to variations in temperature and precipitation, different change scenarios were constructed. The assumed change scenarios included a temperature increase of 1.5 °C, a temperature decrease of 1.5 °C, a precipitation increase of 20%, and a precipitation decrease of 20%. Among them, the temperature rise was more consistent with the future climate change in the Tianshan region. Comparing the glacier runoff under different scenarios, it was found that the glaciers in this region are particularly sensitive to temperature changes compared with precipitation changes. The main reason for this is that the volume of the glacier is relatively small, and although the glacier has a certain lag in response to climate change, a slight increase in temperature will cause the height of the mass balance line of the glacier to rise, and the line of the points on the glacier surface where the sum of the accumulation and melting amount is zero in a certain period of time (generally one year). This can divide the glacier into an accumulation area and a melting area. In this scenario, the melting area of the glacier expands and the accumulation area shrinks, which eventually leads to retreat of the end of the glacier and thinning of the thickness of the glacier in the melting area. Similarly, when the regional temperature decreases, the height of the mass balance line of the glacier decreases, the overall mass balance of the glacier is positive, and the glacier is in an accumulation state.

In evaluating the sensitivity of runoff in the study area, we established control experiments by assigning incremental and decremental changes to temperature and precipitation, thereby assessing the differential sensitivity of runoff to these two factors. Indeed, this method allows for evaluation of the response relationship between runoff and climatic factors. However, changes in climatic factors are actually a gradual process. Simply increasing or decreasing the values of climatic factors to analyze their relationship can overlook many details in the interaction processes between temperature, precipitation, and the underlying surface in the study area, such as glacier melt and precipitation recharge to glaciers. The coefficient of variation (CV) is commonly used to characterize runoff variability. Based on the simulation results, the CV values of runoff during the melt season (June, July, and August) in the study area from 1997 to 2016 were calculated to be 0.25, 0.08, and 0.2, respectively. The CV is smallest in July, which is the period of the most intense glacier melt and the highest precipitation in the study area. During this period, the variability of glacier runoff is higher.

5. Conclusions

In this study, percentile correction and statistical downscaling were applied to process the output from the CMIP6 regional climate model, aiming to provide more accurate meteorological data for the HBV-Light hydrological model. Glacial meltwater accounted for the largest share of the runoff component, ranging from 42.10% to 43.79%, followed by snowmelt runoff at 29.64% to 30.40%, and rain runoff constituted 26.56% to 27.49%. In the period from 2040 to 2059, the runoff in the Urumqi River source region is projected to increase. During this simulation phase, the total runoff in the study area will mainly be composed of rain runoff. Despite the continuous retreat of glaciers due to climatic influences, the increase in precipitation is expected to gradually compensate for the reduction in runoff caused by glacier retreat. By designing control experiments with different variations in temperature and precipitation, the sensitivity of the HBV model’s simulated runoff to these two meteorological factors was evaluated. The results indicated that temperature changes have a significant dominant effect on runoff in the study area. Under the T + 1.5 °C P − 20% scenario, runoff still increased by 22.06%. Additionally, under the T − 1.5 °C P + 20% scenario, runoff still decreased by −35.81%. The primary reason for this is that the study area is located in a mid-to-high-altitude region. Temperature changes not only affect glacier melt, but also influence the form in which precipitation occurs (solid, liquid, or a mixture of both), thereby impacting total runoff. The rational evaluation and analysis of the composition of river runoff, which is mainly fed by glaciers, is very important for areas highly dependent on their output water resources, and the simulation and evaluation of large-scale watershed runoff will also be the main goal of later work.