The Impact of Climate Change on Hydrological Processes of the Glacierized Watershed and Projections

Abstract

Under the influence of climate change, the hydrological processes of glaciers have undergone significant changes, a fact which is seriously affecting agricultural production in the downstream region of the Tianshan Mountains, China. In order to explore the intrinsic relationship between climate change and hydrological elements, we proposed an “evaluation-driving-prediction” system to study it. First, we constructed a glacier-enhanced soil and water assessment tool model (GE-SWAT) and used a two-stage calibration method to optimize the model parameters. Next, a scenario analysis was used to evaluate the driving factors of historical runoff changes. Finally, we projected future runoff changes using bias-corrected regional climate model (RCM) outputs. The results of the case study on the Jinghe River Basin in the Tianshan Mountains show that from 1963 to 2016, total runoff increased by 13.3%, 17.7% of which was due to increasing precipitation and 1.8% of which was negated by rising temperatures. The glacier runoff increased by 14.5%, mainly due to the rising temperatures. A 3.4% reduction in snowmelt was caused by a lower snowfall/precipitation ratio, which significantly reduced the snowfall from June to August. The RCM projection indicated that the warming and humidification phenomenon in the study area will continue at least through to the mid-21st century. A consistent increase in glacier runoff and total runoff is projected, but the contribution rate of the glacier runoff will have little to no change under the RCP4.5 and RCP8.5 emission scenarios. Our research demonstrates the simulation performance of the GE-SWAT model in a basin with moderate glacier cover. This method is shown to be efficient in quantifying the impact of climate change on glacier hydrological processes and predicting future streamflow changes, providing a good research reference for similar regions.

1. Introduction

As the source catchment area of many rivers in Central Asia, the Tianshan Mountains incorporate 10,778 glaciers (with an area of approximately 13,567 km²) [1], supplying valuable freshwater resources to oasis agriculture and livestock downstream, and nurturing nearly 50 million people from Kazakhstan, Kyrgyzstan, and China [2]. However, the warming characteristics of the climate in the Tianshan Mountains have been very obvious for more than half a century. The warming rate in Central Tianshan and Eastern Tianshan is up to 0.45 °C/10a, which is 2.6 times the global average [3], and has been fluctuating at a high level since 1998 [4]. Since the end of the Little Ice Age, the increasing negative glacier mass balance in the Tianshan Mountains has provided increased meltwater (approximately 97.5% of the glaciers are in a state of retreat [1]), and the degradation of permafrost has increased infiltration, which ultimately affects the quantity and seasonal distributions of the streamflow [5]. In the past few decades, the peak value of the discharge hydrograph over the course of a year has shifted to an earlier time, the flood season runoff has increased, and the total runoff has increased [6]; however, it will eventually decrease, due to the reduction in glacier area [7], further exacerbating the instability of water resources and the contradiction between supply and demand. Due to the critical role of glaciers in downstream water supplies, systematic assessments of past, present, and future changes are required to develop imperative and efficient measures to adapt to these changes in climate.

Numerous studies have been carried out, focusing on predicting runoff trends by coupling hydrological models with global or regional climate models, but little attention has been paid to the influence mechanism and contribution rate of historical climate change on hydrological processes [8]. Meanwhile, there has been strong debate on the changing trend of future river runoff due to a lack of verification data for future periods [9]. Unger-Shayesteh et al. summarized approximately 100 studies on past changes in climate, glaciers, and runoff in Central Asian headwater catchments, and suggested that there is a need for sound attribution studies linking the observed hydrological changes to climatic and cryospheric changes in order to carry out more reasonable forecasts of future runoff [10]. The Jinghe River Basin (JRB), located in the central Tianshan Mountains, has a glacier coverage rate of 6.2%, with small glaciers (with a size less than 1 km²) accounting for 68.1% [11]. With climate warming, the decrease in the total glacier area in the JRB is accelerating and affecting the future. We found that the current literature is more focused on the response of glaciers and runoff in watersheds with high glacier coverage (>10%) to climate change, such as the Kumalak River [8,12], Manas River [13], Muzati River [14], and Urumqi River [15]. However, the effects of glacial meltwater have also been monitored in watersheds with moderate glacier coverage (<10%) [16].

A glacier module calculates the processes of glacier melt, sublimation/evaporation, and accumulation [17]. Glacier melt algorithms are the key elements in a glacier module, and have considerable impact on the calculation of the future glacier area and runoff, mainly due to the energy balance approach and the temperature-index approach [18]. The energy balance approach uses many parameters and a complex calculation process, which is suitable for limited areas with substantial observations [19]. Due to a lack of observational data, especially in high elevation zones, the temperature-index approach is still widely used. Despite the simplicity, this approach often performs remarkably well [20]. The melting process of a glacier is calculated using the temperature-index approach. The glacier module then adjusts the glacier area through the empirical glacier response models, and feeds back the response of the glacier to climate change (advance and retreat) into the model output. The empirical glacier response models are based on a specified relationship between glacier volume and area, including the volume-area (V-A) method [21], the accumulation area ratio (AAR) method [22], and the Δh-parameterization method [23]. To date, the glacier module has been basically completed.

In hydrological simulation, the adaptation and rationality of parameter values are extremely important. However, in cases where the observed discharge is the only available data for parameter optimization [24], this may lead to “equifinality for different parameters”; i.e., if there are multiple parameter sets to simulate sufficient runoff in spite of an error in the glacier melt simulation [25]. This problem is more prominent in watersheds with moderate glacier coverage, where the glacier module may be less sensitive to glacier melt. Some scholars have stated that, when simulating the hydrological processes of glacierized areas, the relevant parameters of the model should be optimized based on glacier information data to improve the “authenticity” of the simulation results [26]. Nevertheless, most alpine mountains lack glacier-monitoring projects due to the harsh climate and complex terrain, and only a few studies have used long-term monitoring data to verify the simulation experiments of hydrological processes in the Alps [27], Himalayas [28,29], Tianshan Mountains [30], and Karakoram [31]. In view of the scarcity of glacier observation data, the first and second Chinese Glacier Inventories (CGI), released by the Chinese Academy of Sciences, are particularly important. At present, the use of the CGI is concentrated in three aspects: (1) data input—integrate the CGI and a land-use/land-cover map to drive the model with an accurate distribution map of the glaciers [32]; (2) formula improvement—calibrate the glacier volume-area scaling ratio to reduce the error estimation of the glacier volume changes [33]; and (3) result verification—calibrate and validate the parameters, including the degree-day factor [9].

In this study, we constructed a glacier-enhanced soil and water assessment tool (GE-SWAT) model, tested its performance in the JRB located in the Tianshan Mountains, analyzed the response of glacier hydrological processes to climate change, and attempted to provide research references for similar watersheds with moderate glacier coverage. First, two performance evaluation criteria, the observed runoff and Chinese Glacier Inventories (CGI), were used for parameter adjustment to ensure an optimal model performance. Then, based on the warming and humidification over nearly 60 years, the possible causes of the changes in the glacier hydrological processes were explored through scenario analysis. Finally, we forced the GE-SWAT model with the bias-corrected regional climate model (RCM) outputs to analyze the characteristics of future climate change and to project the evolution trend of the runoff components.

2. Study Area

The JRB, located in the northern slope of the Eastern Tianshan, covers a source area of 1419 km², and the elevation ranges between 674 and 4342 m a.s.l., with a generally increasing trend from the northwest to the southeast [34]. The Jinghe River is fed by two major tributaries, the Dongdujinghe River to the south and the Wutujinghe River to the north. The JRB has a typical temperate continental climate, with the water vapor transport mainly coming from westerly winds from Siberian anticyclones. The long-term average annual temperature, precipitation, and runoff observed from 1961 to 2016 were −1.19 °C (based on an altitude of 2150 m a.s.l.), 522 mm, and 4.66 × 10⁸ m³, respectively. According to the first Chinese Glacier Inventory for the JRB, the fraction of glacier coverage is 6.2%, and the glacier area is 87.38 km², with the largest one extending over 10 km². There are 129 glaciers in the JRB, including 27 glaciers in the altitude range of 3000–3600 m a.s.l., 79 glaciers in the 3000–3600 m a.s.l. range, and 23 glaciers above the 3800 m a.s.l. altitude. The abundant alpine glacier meltwater is a major source of streamflow, as is the Ebinur Lake basin.

3. Data and Methodology

3.1. Data Description

3.1.1. Underlying GIS-Referenced Data

Two periods of land-use data (1990 and 2018) were obtained from the Resource and Environment Data Center of the Chinese Academy of Sciences (http://www.resdc.cn, accessed on 20 August 2021), with a spatial resolution of 100 m. The digital elevation model (grid size of 30 × 30 m) and the soil texture data were taken from ASTER-GDEM data and the Harmonized World Soil Database, respectively. The digital vector map of glaciers was extracted from the first and second Chinese Glacier Inventories (FCGI and SCGI) (http://www.ncdc.ac.cn, accessed on 20 August 2021).

3.1.2. Streamflow Monitoring and Meteorological Data

Monthly streamflow records for the Jinghe hydrological station were collected from the Hydrology and Water Resource Survey Bureau in Bole, China for the 1961–2010 period, and from the Water Resources Bulletin in Bozhou, China for the 2010–2016 period. The National Meteorological Centre provided the 0.5-degree gridded datasets (http://data.cma.cn/, accessed on 20 August 2021), including daily precipitation and temperature data from 1961 to 2016. These datasets were obtained by interpolating the observation data of 2474 national meteorological stations. These grid points can be regarded as the virtual meteorological stations, which well-reflect the spatial characteristics of precipitation in complex terrain in the Tianshan Mountains [17,35].

3.1.3. RCM Data and Bias Correction

In this study, the RCMs generated in phase 1 of the Coordinated Regional Climate Downscaling Experiment for East Asia (CORDEX-EA, https://cordex.org/, accessed on 20 August 2021) were used as future climate forcing data to predict future runoff changes under the two emission scenarios (RCP4.5 and RCP8.5). CORDEX-EA used the 1979–2005 historical meteorological conditions to simulate climate change for 2006–2049. In order to reduce the uncertainty of a single model, it was necessary to use the ensemble average of multiple RCM results [36,37]. We selected four RCMs with better performances based on the research results of Yin et al. [38], as shown in

Table 1. The four regional climate models used in this study and their characteristics.

Model Name	Spatial Resolution	Time Series	Institute	Published Time
YSU-RCM [39]	50 km	1980–2050	Yonsei University	2013
SNU-MM5 [40]	50 km	1979–2049	Seoul National University	2009
SNU-WRF [41]	50 km	1979–2049	Seoul National University	2008
RegCM4 [42]	50 km	1979–2050	Kongju National University	2012

.

While the RCMs were considered to have high accuracy, the important regional factors (e.g., terrain, vegetation, and cloud cover) that affect the local climate could not be fully considered, especially for the high-altitude area [29,43]. Hence, we used the multivariate bias correction (MBC) approach to correct the bias in the RCM projections. MBC, which was developed by Mehrotra et al. [44], is a user-friendly technology that can simultaneously adjust the deviations of multiple variables on multiple time scales. For different time scales, the MBC adopts a three-step correcting procedure to adjust the univariate and multivariate deviations one by one. First, it calculates the statistical distribution of each variable, such as the mean value, standard deviation, and the lag-0 and lag-1 auto cross-correlation coefficients. Next, the mean and variance of each variable are corrected step by step. Finally, a multivariate autoregressive model is used to correct the time and cross-variable dependence bias. The bias-corrected time series are aggregated or averaged to the next time scale, and the same process is repeated.

The standard multiple autoregressive first-order model used for observed and RCM data is as follows:

$${\hat{Z}}_i^h=C{\hat{Z}}_{i-1}^h+D{\epsilon }_i$$

(1)

$${\hat{Z}}_i^g=E{\hat{Z}}_{i-1}^g+F{\epsilon }_i$$

(2)

where ${\hat{Z}}_i^h$ is monthly observations, ${\hat{Z}}_i^g$ represents monthly RCM outputs, i indicates time steps, C and D are the lag-1 and lag-0 cross-correlations of the observations, E and F represent the corresponding coefficients matrices for the RCM outputs, and ${\epsilon }_i$ is the normalized vector obtained from the sequence.

To obtain the adjusted data for the month i of year t, the calibration model is revised as:

$$Z_{t,i}^{’g}=C_i{Z}_{t,i-1}^{’g}+D_i{F}_i^{-1}{\hat{Z}}_{t,i}^g-D_i{F}_i^{-1}{E}_i{\hat{Z}}_{t,i-1}^g$$

(3)

where $Z_{t,i}^{’g}$ is the adjusted time series of the previous month in year t. After correction, the time series $Z^{’g}$ is readjusted through the average and standard deviation of the observations to obtain the final corrected time series ${\overline{Z}}^g$. A detailed description of the MBC can be found in Mehrotra et al. [45,46].

3.2. GE-SWAT Model Description

3.2.1. Brief Description of the SWAT Model

The soil and water assessment tool (SWAT) generates a hydrological response unit (HRU) by considering the comprehensive effects of soil properties, terrain, vegetation, and human land management to flexibly handle various complex conditions [47]. HRUs are the basic computational units of the SWAT model. During the modeling process, SWAT first assumes that the area of the HRUs remains unchanged and simulates the hydrological cycle inside a single HRU separately, then connects the HRUs through the river network to obtain the outlet sectional flow [48]. In this study, we delineated the JRB into 71 subbasins (Figure 1) and 1983 HRUs, of which each subbasin was further divided into 8 elevation zones based on elevation difference, and each glacier was considered as an HRU.

The SWAT model currently uses the conventional temperature-index method to estimate snowmelt, and uses a linear two-reservoir method to simulate baseflow, but it is deficient in the simulation of glaciers [8]. Therefore, the SWAT model coupled with a glacier module was constructed to simulate the glacier hydrological processes. It should be noted that, in this study, glacier runoff refers to the runoff generated by glacier melting, and excludes the runoff generated by snow melting.

3.2.2. Glacier Module

A glacier module contains three main algorithms: glacier melt algorithm, glacier volume-area scaling law, and glacier accumulation rate. First, the initial glacier area is used to calculate the glacier meltwater; meanwhile, considering that the meltwater refreezes and the glacier accumulates, the glacier area is updated using the V-A method. A detailed description of the glacier accumulation rate can be found in Luo et al. [13] and Yin et al. [17]. The SWAT model contains the source code for calculating HRU evaporation, so we did not calculate the sublimation/evaporation of the glacier regions separately here.

An enhanced temperature-index approach [49,50] correlates the daily glacier melt with temperature and terrain, including a consideration of the influence of potential solar direct radiation, thus accounting for the enormous spatial difference in the degree-day factors (DDFs) on steep slopes. Here, daily temperature data in each elevation zone were computed by a constant temperature lapse rate (TLAPS), and potential solar radiation was calculated based on the research results of Garnier et al. [51]:

$$M=\{\begin{array}{c}\left(F_M+R_{ice}{I}_{pot}\right)\times \left(T-T_{mlt,ice}\right), T>T_{mlt,ice}\\ 0, T\le T_{mlt,ice}\end{array}$$

(4)

where $M$ is the daily glacier meltwater, $T_{mlt,ice}$ represents the threshold temperature of the glacier melt, $T$ illustrates the daily mean temperature, below $T_{mlt,ice}$ means that no melting has occurred, $F_M$ is the ice temperature melt factor, $R_{ice}$ denotes the ice radiation melt factor, the default value of which is 1.44 × 10⁻⁵ [52], and $I_{pot}$ represents the potential solar direct radiation.

There is a lack of glacier-monitoring projects in the JRB. Thus, we updated the glacier area based on the V-A formula introduced by Chen and Ohmura [53] to simulate the glacier dynamic variation:

$$V_g=c\times {\left(A_g\right)}^{\gamma }$$

(5)

where $V_g$ and $A_g$ are the volume and surface area of the glacier, respectively. The constant $c$= 0.04 and the dimensionless scaling coefficient $\gamma$= 1.35 were obtained from ice-penetrating radar thickness measurements calculated by Liu et al. [54] in the Qilian and Tianshan Mountains in China.

The relationship between the volume, area, and equivalent water depth of glaciers can be expressed as:

$$V_g=\frac{W_g\times A_g}{{\rho }_i}$$

(6)

where ${\rho }_i$ is the bulk density of the ice, the value is 900 kg × m⁻³, and $W_g$ is the equivalent water depth of the ice.

3.3. Model Performance Evaluation

In this study, the CGI and runoff data were used for model calibration using a two-stage calibration method. First, we described the FCGI and SCGI as one special land-cover class, and merged them into the land-use maps of two periods, respectively, to use accurate spatial distribution data of glaciers to drive the model. Second, we used the initial glacier input data from the FCGI to simulate the glacier mass balance in each subbasin, and transformed them to obtain the glacier area change in the subbasins, further correcting the DDFs with reference to the updated glacier distribution datasets of the SCGI. Finally, the runoff was used to calibrate and validate other parameters of the model.

The Nash-Sutcliffe efficiency (NSE) coefficient, the percent bias (PBIAS), and the coefficient of determination (R²) were selected to evaluate the goodness of fit of the hydrological process. The NSE represents the degree of curve matching between the simulated and natural runoff values, which is between −∞ and 1.0. The closer the value is to 1.0, the better the curve-matching degree. The PBIAS measures the average tendency of the simulated runoff to be larger or smaller than the natural runoff [55]. Positive PBIAS values indicate that an overestimation of the model results relative to the measurement conditions, negative values indicate a model underestimation bias, and the optimal value is 0.0. The R² value describes the proportion of the variance between the simulated and natural runoff values. The closer R² is to 1.0, the better the curve fit between the simulated and natural runoff values. For monthly runoff, it is generally considered that when 0.75 ≤ NSE ≤ 1.0 and PBIAS ≤ ±10%, the simulation yielded “very good” results [56].

3.4. Attribution Analysis for Runoff Change

In this study, the Mann-Kendall method was employed to determine the variation trend of the annual runoff series, and the Pettitt test and sliding t-test (step size was 7) were employed to detect the change point, so as to estimate the hydrometeorological change degree before and after the change point. The change point occurred in 1997 (p = 0.102; t = −2.323). The entire period of streamflow was then subdivided into two historical periods—Period 1 (1961–1997) and Period 2 (1998–2016).

We simulated the water cycle process in Period 1 using the GE-SWAT model. While keeping the model parameters unchanged, we regenerated the climate datasets for both scenarios and used them as the driving data of the model. Finally, we obtained the runoff series under the influence of changes in temperature and precipitation. In Scenario 1, precipitation remains at the Period 1 level and the temperature undergoes natural variations, which investigates the impact of temperature changes on runoff ($\Delta$Q_t). In Scenario 2, the temperature remains at the Period 1 level and precipitation undergoes natural variations, which investigates the impact of precipitation changes on runoff ($\Delta$Q_p) [12]. The natural scenario means that both precipitation and temperature are maintained at the level of Period 1:

$$\Delta Q_{tot}=\Delta Q_{lucc}+\Delta Q_{cl}=\Delta Q_{lucc}+\Delta Q_t+\Delta Q_p$$

(7)

where $\Delta$Q_lucc and $\Delta$Q_cl represent the land-use-change-induced and climate-change-induced runoff changes, respectively, and $\Delta$Q_tot is the total runoff change.

4. Results

Next, we evaluate the model performance based on the historical data of the glacier area and runoff, and quantitatively calculate the impact of precipitation and temperature changes on hydrology. We also forecast the trend and quantity of future runoff using a model forced by the bias-corrected RCM data under RCP4.5 and RCP8.5.

4.1. GE-SWAT Model Performance

4.1.1. Model Parameterization

A two-stage calibration method was used to calibrate the parameters of the glacier regions and the total catchment step by step.

Table 2. The optimal parameter values for the Jinghe River Basin during Period 1.

Parameters	Units	Chosen Value	Definitions
ALPHA_BF	days	0.01	Baseflow alpha factor
CN2	—	40.00	Initial SCS runoff curve number for moisture condition II
LAT_TIME	days	22.00	Lateral flow travel time
ESCO	—	0.10	Soil evaporation compensation factor
PLAPS	mm km−1	44.50	Precipitation lapse rate
TLAPS	°C km−1	−4.10	Temperature lapse rate
SMFMX	mm °C−1 d−1	2.00	Melt factor for snow, 21 June
SMFMN	mm °C−1 d−1	1.00	Melt factor for snow, 21 December
SFTMP	°C	2.00	Snowfall temperature
SMTMP	°C	1.50	Snowmelt base temperature
Bgmlt,12	mm °C−1 d−1	12.00	Glacier melt degree-day factor, 21 December
Bgmlt,6	mm °C−1 d−1	5.50	Glacier melt degree-day factor, 21 June
gmlt_tmp	°C	0.50	Glacier melt base temperature

shows the optimal parameter values selected for the JRB.

(1) Glacier regions

The parameters that have a higher correlation with glacier mass balance are DDFs and TLAPS [25]. Here, the starting point FCGI and ending point SCGI were used to verify the simulated value of glacier area change in each subbasin, further correcting the DDFs. The initial values of TLAPS and precipitation gradient (PLAPS) in the JRB were calculated, respectively, based on the statistical relationship between the elevation and multiyear average temperature and precipitation of 18 virtual meteorological stations on the northern slope of the Tianshan Mountains, China. PLAPS was further optimized according to glacier area change and the water balance of the JRB.

(2) Total catchment

The parameter calibration principle was as follows: first ensure the water balance, then correct the hydrological process. The initial parameter estimates of the model were derived from the physically acceptable values provided by the SWAT model database and previous research results for the vicinity [57]. We compared the simulated and actual hydrological processes based on the evaluation indicators NSE, PBIAS, and R², and used a trial-and-error method to adjust the parameters in the management, groundwater, soil component, subbasin, HRU, and other modules.

4.1.2. Glacier Area Validation

Figure 2 shows the observed and simulated changes in glacier area of each subbasin. The observed glacier area based on the FCGI and SGCI decreased by 27.3% (23.82 km²) from 1959 to 2007, and the simulated glacier area lost 31.4% during the same period. We also calculated the NSE and PBIAS of the glacier area changes in each subbasin, which were 0.87 and 10.79%, respectively, indicating that the simulated changes in glacier area were similar to the observed values. In addition, glacier meltwater accounted for 20.2% of the total runoff in the JRB according to the simulation results, which was relatively consistent with the glacier investigation results of Wang et al. [58]. In their survey, the proportion of glacier meltwater in the JRB was 20.6% between 1964 and 2004.

4.1.3. Simulation Runoff Validation

Figure 3 compared the observed and simulated monthly discharge curves, and the two curves were closely aligned. The monthly flow is characterized by high flows of 15–65 m³/s from June to August, and low flows of 1–10 m³/s from November to April, which were also reconstructed in the hydrological process simulation. The NSE, PBIAS, and R² evaluation indices for Period 1 are listed in

Table 3. Evaluation indices for assessing performance of the simulation.

	Calibration Period			Validation Period			Period 1
	NSE	PBIAS	R2	NSE	PBIAS	R2	NSE	PBIAS	R2
The GE-SWAT model	0.84	7.6%	0.89	0.84	−3.7%	0.86	0.85	−0.7%	0.87
Traditional SWAT model	0.83	−11.1%	0.84	0.72	−23.7%	0.81	0.77	−19.7%	0.81

. We chose 1961–1962 as the warm-up period, 1963–1980 as the calibration period, and 1981–1997 as the verification period. GE-SWAT outperformed the traditional SWAT model in the simulation results with a good ability to capture the change process of the runoff. The NSE, PBIAS, and R² for the GE-SWAT were 0.85, −0.7%, and 0.87, respectively. However, the NSE, PBIAS, and R² for the traditional SWAT model were 0.77, −19.7%, and 0.81, respectively. The simulation yielded “very good” results, which demonstrated that the structure of the GE-SWAT was credible, that the selection of parameters was reasonable, and that further analysis could be performed based on the output results of the model.

4.2. Attribution of Streamflow Changes

4.2.1. Characteristic of Streamflow Changes

Table 4. The changing characteristics of hydrometeorological elements in a seasonal scale.

Season	Runoff (mm/10a)	Glacier Runoff (mm/10a)	Precipitation (mm/10a)	Temperature (°C /10a)
Spring	0.30	0.00	3.21 *	0.27
Summer	4.33 *	1.48	1.76	0.24 *
Autumn	4.56	2.50	5.06	0.30 *
Winter	2.37 *	0.31 *	5.33 *	0.26 *
Year	11.56 *	4.29 *	15.36 *	0.27 *

* Indicates a significant trend (significance level α = 0.05).

summarizes the variation characteristics of total runoff, glacier runoff, and meteorological elements during the 1963–2016 period.

(1) In the past 60 years, the hydrometeorological elements in the JRB have shown an increasing trend. From 1963 to 2016, the temperature increased at a rate of 0.27 °C/10a, and precipitation at a rate of 15.36 mm/10a. Affected by climate warming and humidification, the total runoff in the JRB climbed at a rate of 11.56 mm/10a, with an increased change point occurring in 1997, i.e., an annual average increase of 13.3% from1998 to 2016. During the same period, glacier runoff also increased year by year, increasing by 14.5%. The contribution of glacier meltwater to total runoff has changed significantly (Figure 4), from 11.2% to 32.8%, with an average of approximately 20.2%, suggesting that it is an important source for the JRB.

(2) The increasing range of hydrological elements varies greatly in each season, corresponding to the seasonal variation characteristics in precipitation and temperature. From 1963 to 2016, the autumn increase rate of total runoff was up to 4.56 mm/10a, dominating the increase degree of annual total runoff, and was related to the substantial increase in autumn precipitation (5.06 mm/10a). Additionally, the increase in glacier meltwater is also most obvious in autumn (2.50 mm/10a), which corresponds to the most significant warming trend in autumn (0.30 °C/10a). In the next section, we will quantitatively identify the effects of temperature and precipitation changes on hydrological processes over the past 60 years.

4.2.2. Attribution of the Effect of Temperature Changes on Runoff

Figure 5 plots the seasonal variation in hydrological elements under three scenarios, and Figure 6 plots their annual variation. Glacier and snowmelt are the most sensitive indicators of climate change. The snowmelt was mainly concentrated in April to July, with the peak time in June, and the glacier melt reached a maximum in July and August. After August, when the air temperature dropped to zero and a new snowfall season began, the snowmelt increased again and the contribution of glacier meltwater began to decline. On average, glacier melt continued to accelerate under the warming mode of scenario 1, with a total increase of 14.6%. In addition, the annual snowmelt decreased by 13.5%, due mainly to the decreasing snowfall/precipitation ratio from June to August caused by the higher temperature (Figure 5b); in particular, the snowfall decreased significantly in June. Overall, in the mutual offset of water changes including glacier runoff, snowmelt, evaporation, and soil water content, the rising temperature caused a slight decrease in the total runoff, with a contribution rate of −1.8%. It is noteworthy that agricultural water consumption accounts for more than 90% of the total water consumption in the JRB, and water demand is concentrated from May to August [59]. Glacier meltwater, in particular, supplies the river from May to August, providing a good supplement (18–32%) to the water demand in the JRB. However, glacier ablation is projected to occur earlier as the temperature rises, leading to a reduction in water supply during the water demand period, thereby affecting the agricultural activities in the middle and downstream reaches.

4.2.3. Attribution of the Effect of Precipitation Changes on Runoff

For Scenario 2, the amount of snowmelt increased by 10.1%, and the increase was particularly significant from April to August. On one hand, this was attributed to the increased snowfall in the cold season, which provided good storage conditions for snowmelt from April to August. On the other hand, the precipitation between May and August largely occurred as rain, so the heat energy it possessed accelerated the snow melting process. Under Scenario 2, the annual hydrograph of glacier melt remained practically unchanged relative to the hydrograph of the natural scenario (Figure 6a), indicating that the change in temperature largely governs the melting process of glaciers. The increase in precipitation only led to a 0.1% decrease in glacier melt, which may be explained by the increase in snow cover area between October and March affecting the exposure time of the glacier. Moreover, the abundant precipitation supplemented the soil water content and increased the evaporation in the JRB. In general, the main driver of the increase in runoff was the increasing precipitation, with a contribution rate of 17.7%.

The causes of runoff changes in historical periods were evaluated and quantified, as shown in

Table 5. Contribution of various factors to runoff changes in the Jinghe River Basin.

Factor Contribution	$\Delta$Qtot	$\Delta$Qlucc	$\Delta$Qcl	$\Delta$Qt	$\Delta$Qp
in mm	41.75	−9.86	51.61	−5.72	57.33
in %	13.3	−2.6	15.9	−1.8	17.7

. In addition, we input the land-use data for 2018 into the GE-SWAT model under the natural scenario climatic conditions to obtain the runoff sequence under the land-use change scenario. The simulation results indicated that the changes in land-cover class caused the runoff to decrease by 8.29 mm, which was consistent with the calculation results of Equation (7) (−9.86 mm). Therefore, we concluded that the method used to assess the contributions of the relevant environmental factors to the changes in runoff was effective [60]. The total runoff in the JRB increased by 13.3%, of which the temperature rise and precipitation increase contributed −1.8% and 17.7% of the increase in total runoff, respectively, and the land-use change contributed −2.6%. Thus, the total runoff change in the JRB is mainly derived from the precipitation changes.

4.3. Prediction of Runoff Changes

4.3.1. RCM Bias Correction

The correction of the RCM datasets was necessary because the original RCM datasets had lower precipitation and higher temperatures than the observed data (Figure 7a,c,e). After the bias correction using the MBC method, the NSE coefficients of the daily maximum temperature and daily minimum temperature were both higher than 0.74 from 1961 to 2016. Daily precipitation is highly random with many uncertain factors, so we only calculated the NSE of the monthly precipitation series [61], and the NSE was greater than 0.71. Figure 7b,d,f show a good match between the bias-corrected RCM results and the observed data for the same historical period. In this study, we used the RCM ensemble means as future climate conditions to predict the evolution trend of the meteorological elements.

4.3.2. Projected Climate Change

The bias-corrected RCM results predict that the warming and humidification of the JRB will continue at least until the mid-21st century (Figure 8).

There is a significant increase in temperature for each emission scenario, and the RCP8.5 scenario is generally warmer than the RCP4.5 scenario. The temperature increase in the JRB is ongoing and is expected to reach a peak in the 2030s under the two emission scenarios, after which the rate of temperature increase will be gentle. The RCP4.5 scenario predicts that the temperature rise rate will be 0.42 °C/10a (1961–2049), and that the annual average temperature in the 2040s will rise by 2.36 °C compared with the baseline period (1961–2016). Under the RCP8.5 scenario, the changes in temperature for the JRB show significant warming rates of 0.53 °C/10a, and the temperature will be 3.07 °C higher than the baseline period by the 2040s (

Table 6. Climate change for future periods under two different RCPs relative to 1961–2016.

	Temperature (°C)			Precipitation (%)
	2020s	2030s	2040s	2020s	2030s	2040s
RCP4.5	1.80	2.49	2.36	6.45	8.88	10.43
RCP8.5	2.59	3.11	3.07	10.49	6.23	14.48

).

As time passes, precipitation in the JRB is projected to rise amidst fluctuations. Under the RCP4.5 scenario, the annual precipitation is projected to increase by 10.43% compared with the baseline period, with an increase rate of 11.0 mm/10a. Under the RCP8.5 scenario, the annual average precipitation is projected to increase by 14.48% by the 2040s, which is equivalent to an increase in precipitation of 12.2 mm/10a. It should be noted that precipitation in the 2030s has obvious signs of decline under the RCP8.5 scenario, which will bring definite uncertainty to the hydrological process.

4.3.3. Future Hydrology

Before simulating the future, we optimized the model parameters for the runoff sequence in Period 2. The NSE, PBIAS, and R² values of the runoff series were 0.82, −1.3%, and 0.85, respectively. The parameters for Period 2 were used to predict the future. The interdecadal evolution characteristics of the total runoff, glacier runoff, and glacier runoff contribution are shown in Figure 9. The RCP4.5 and RCP8.5 scenarios suggest that there will be a continuous wet period for the JRB during the forecast period, due to the comprehensive effects of rising temperatures and increased precipitation. The analysis of simulation results based on the bias-corrected RCM data shows that the total runoff has increased at an average rate of 13.96 mm/10a under the RCP4.5 scenario, and 23.82 mm/10a under the RCP8.5 scenario. We also observed that the increase rate of runoff is higher than the change rate of 11.56 mm/10a in the historical period under the two emission scenarios. Compared with the historical period (1961–2016), the total runoff in the future period (2017–2049) will increase by 27.5% and 48.5% under the RCP4.5 and RCP8.5 scenarios, respectively.

Under the two emission scenarios, the glacier runoff will slowly climb at a rate of 2.55 mm/10a (RCP4.5) and 6.22 mm/10a (RCP8.5), and “glacier inflexion” will not yet be reached in the JRB during the forecast period. It is observed that the glacier runoff contribution will not change greatly, as the rainfall-induced runoff still largely governs the total amount of runoff. Under the RCP4.5 scenario, the contribution of glacier meltwater to total runoff decreases slightly (contributes 19.0% of runoff), and rainfall runoff contributes a remarkable portion of the total runoff, which weakens the regulation effect of glaciers on the streamflow. The annual glacier melt is greater under RCP8.5 compared with RCP4.5. Under the RCP8.5 scenario, the role of glacier runoff in the regulation of rivers slowly increases, with an annual contribution of 20.2% in the historical period and 21.3% in the future period. However, the sharp reduction in glacier area is likely to aggravate the future decrease in glacier meltwater, leading to an increasing variability of water resources at yearly time scales in the long run. This projection result is not conducive to the development and utilization of water by water-resource managers.

5. Discussion

5.1. Uncertainty of Input Data

(1) The 0.5-degree gridded datasets

The GE-SWAT model has always shown good performance in simulating glacier hydrological processes, except that runoff in the 1970s was generally overestimated, which may be attributed to the uncertainty of meteorological forcing data in alpine watersheds. In view of the scarcity of observed climate data in high elevation zones, atmospheric reanalysis datasets (such as TRMM [62], APHRODITE [8], and PGMFD [63]) are commonly used as a substitute for observed data. The 0.5-degree gridded datasets were used as the model input for historical periods in this study, and the meteorological input data of different elevation zones were then adjusted to make the physical meaning of the model clearer by setting the elevation zone, PLAPS, and TLAPS. The initial PLAPS was calculated from the datasets of 18 virtual meteorological stations inside and adjacent to the basin. However, 12 of the selected stations are concentrated in the west, which may increase the error of precipitation in the east as it is affected by topographic differences and atmospheric circulation [64]. Subsequently, the PLAPS was optimized by combining the NSE and PBIAS coefficients according to Luo et al. [13] and Wang et al. [8]. As the PLAPS increases, the NSE coefficient first increases and then decreases, while the PBIAS keeps increasing (Figure 10). The PLAPS was finally selected as 44.5 mm km⁻¹, and at this time, the index NSE was 0.85 and the PBIAS was −0.7%, meaning that the simulation results are only slightly underestimated.

(2) RCM outputs

Generally, the uncertainty of RCM outputs and the MBC method is much greater than that of a hydrological model. In this study, the ensemble average of the four RCMs was used to reduce systematic errors and biases. When evaluating the accuracy of daily precipitation after bias correction, Teng et al. [65] and Chiew et al. [66] pointed out that the effect of simulating historical hydrological processes with the corrected precipitation data should be comprehensively considered. Gan et al. [67] compared the effects of runoff simulation using bias-corrected and simultaneous observed precipitation data to investigate whether the bias-corrected data affects the simulation effect of the model. Their results suggested that very little bias was found between the two runoff sequences. Prior to simulating the hydrological responses under future climate change in this study, both the daily bias-corrected and observed precipitation data from the period between 1979 and 2005 were used to drive the GE-SWAT model, thus judging the accuracy of the prediction. The NSE, PBIAS, and R² values were 0.82, −2.21%, and 0.83, respectively. We believe that the bias-corrected RCM outputs have the ability to reproduce historical weather conditions, enabling us to further predict future hydrological changes.

The Tianshan Mountains are highly ridged, with a wide 1000–5000 m a.s.l. elevation range [13]. There are currently 39 ground weather stations in the Tianshan Mountains, China, 15 of which are located from 1000 to 2000 m a.s.l., and only 4 are located above 2000 m a.s.l. elevation. The scarcity of meteorological observations makes the study of changes in river runoff and its impact more challenging. In the alpine watersheds, adding ground observation sites, enhancing remote sensing observation technology, and improving the comprehensive observation network of water cycles from upstream glaciers to downstream valleys are urgently needed [27].

5.2. Uncertainty of DDF Factors

Glacier ablation is the transformation and release process of glacier water resources from solid to liquid, and the mutual compensation between it and precipitation may trigger uncertainties in basin water balance [68]. Immerzeel et al. [69] indicated that river flows in the upper Yangtze and Brahmaputra Rivers would perhaps decrease by mid-century because of the reduced glacier meltwater. Zhao et al. [9] reported the opposite results; i.e., a sustained increase in total runoff from the Yangtze and Brahmaputra Rivers was expected, mainly due to increased rainfall runoff that compensated for the loss of snow and ice melting. However, Su et al. [70] had different views on the causes of the runoff change. They believe that the accelerated melting of glaciers is the main reason for the total runoff increase for the Brahmaputra River. When predicting the possible response of glaciers to future climate change, models must be able to simulate glacier mass balance with reasonable accuracy, not just flow.

As with all models, the application of the GE-SWAT depends on parameter settings [71]. Figure 11 shows the sensitivity of the GE-SWAT model to DDFs. When the DDFs change by 0.5 unit (mm °C⁻¹ d⁻¹), the glacier runoff will change by approximately 32%, and the contribution of glacier runoff will change by approximately 4.5%. The DDFs largely govern the simulation accuracy of glacier hydrological processes, which can lead to significant changes in the contribution of glacier meltwater to total runoff and its future trend. The conventional temperature-index method is based on the empirical statistical relationship between air temperature and glacier melt, and the same DDF setting means that regional differences in climate control the amount of glacier melt [71]. The enhanced temperature-index approach can well-represent the inhomogeneity of DDFs in space and time. We then adopted the two-stage method to carry out multiparameter calibration, which reduced the uncertainty caused by the lack of future test data, and were very confident in the authenticity of the hydrological process reflected by the model.

5.3. Potential Risks and Responses of Runoff Change

Glacier water resources play an important role in socio-economic development and ecological environment restoration. Agricultural production is the main body of the local economy, with its water consumption accounting for more than 90% of the total water, and the peak period of water demand is concentrated from May to August. Glacier meltwater has a good supplementary effect during the peak water consumption period, contributing 18-32%, even though the fraction of glacier coverage in the JRB is only 6.2%. Ebinur Lake, the terminal lake of the JRB, is not only a wetland nature reserve, but also the second largest lake in Xinjiang, China. In recent years, the contradiction between supply and demand has become increasingly prominent in the JRB, as the increasing irrigation area and population in the downstream lowland have intensified the water consumption [34]. Although the upstream water continues to increase, the water entering the lake from the estuary of the Jinghe River is still decreasing, resulting in a nearly 57% reduction in the surface area of Ebinur Lake since the 1950s [72,73]. Despite a series of local ecological environmental protection measures and water-saving engineering projects being implemented to improve the ecological environment of Ebinur Lake, real adaptation policies need to take into account the impact of future climate change on water resources, not just current climate variability, if they are to be measures that deliver sustainable benefits to society [74].

Current long-term planning for the JRB, used, for example, in the comprehensive planning of water resources, rarely considers the impact of glacier meltwater. As “solid reservoirs”, glaciers stabilize and regulate streamflow by adding surface runoff during arid periods [75]. Both future total runoff and glacier runoff in the JRB will increase, according to the prediction results of RCMs. It is noteworthy that under the RCP4.5 scenario, the regulatory effect of glaciers on the hydrological regime is weakened as a result of increased rainfall runoff. Therefore, changes in runoff fluctuations should be investigated to ensure a more balanced water resource allocation. Despite certain errors in the future projections, it is a fact that glaciers are shrinking under the trend of global warming. The annual retreat rate of the glacier area in the JRB (0.57%) is much higher than that of other basins in the Tianshan Mountains, such as the Manas River Basin [13], the Sary-Djaz-Kumaric River Basin [8], and the Kuitun River Basin [6], which are 0.29%, 0.14%, and 0.38%, respectively. In the short term, the accelerated glacier melting may not only change the hydrology and its cycle characteristics, but also bring a series of disasters [76], such as glacier floods, glacier lake outbursts, glacier debris flows, etc. In the long run, the glaciers will potentially disappear as the effect of climate warming continues. If this occurs, how will the glacierized basins ensure the security of the water resources? Generally, rising temperatures increase the demand for agricultural water. In the JRB, only one Xiatianji reservoir has been constructed, the water system connectivity project is lacking, and the ability to resist drought and flood is weak. Once the glacier meltwater replenishment disappears, the food security of the downstream regions and the ecosystem security of Ebinur Lake will be threatened. It is necessary for water managers to comprehend the impact of glacial retreat on water supply, and the fact that they may disappear with climate warming, even in watersheds with moderate glacier coverage.

6. Conclusions

In this work, we initially established an “evaluation-driving-prediction” system. This system allowed us to conduct a systematic analysis of glacier hydrological processes under past, present, and future climatic conditions, with the past concentrating on quantitatively evaluating the impact of climate change on runoff and the future focusing on predicting an increase or decrease trend and runoff quantities. The results of the case study conducted for the JRB in the Tianshan Mountains are as follows:

(1)

A two-stage calibration method was applied to improve the performance of the glacier-enhanced SWAT model, which allowed it to have the ability to capture the change process of the runoff in the basin with a moderate glacier cover.

(2)

From 1963 to 2016, the mountain river showed an increasing trend; in particular, a significant increase in runoff of 13.3% was seen in the last 20 years, −1.8% of which was caused by temperature rises and 17.7% by an increase in precipitation. Specifically, with temperatures rising, glacier runoff increased by 14.6%, and snowmelt decreased by 13.5%. Increased precipitation had a negligible effect on glacier melt but provided a larger snow water equivalent, resulting in a 10.1% increase in snowmelt from April to August.

(3)

The RCM projections indicated that the warming and humidification phenomenon in the JRB will continue at least through to the mid-21st century. Under the RCP4.5 (RCP8.5) scenario, the annual mean temperature is projected to increase by 0.42 (0.53) °C/10a, precipitation by 11.0 (12.2) mm/10a, total runoff by 13.96 (23.82) mm/10a, and glacier runoff by 2.55 (6.22) mm/10a.

(4)

The glacier runoff contributed 20.2% of streamflow to only 6.2% of the catchment area. The glacier area retreat rate (0.57%/year) in the JRB was higher than that of other catchments in the Tianshan Mountains. There is a need for glacier water resources change adaptation measures to be discussed and implemented, even in catchments with a moderate glacier cover.