Modeling Glacier Mass Balance and Runoff in the Koxkar River Basin on the South Slope of the Tianshan Mountains, China, from 1959 to 2009

Abstract

Water resources provided by alpine glaciers are an important pillar for people in the arid regions of west China. In this study, the Hydrologiska Byrans Vattenavdelning (HBV) light model was applied to simulate glacier mass balance (GMB) and runoff in the Koxkar River Basin (KRB) on the south slope of Mount Tumur, in the western Tianshan Mountains. Daily temperature and precipitation were calculated by multiple linear regressions and gradient-inverse distance weighting, respectively, based on in-situ observed data by automatic weather stations (AWSs) in the Koxkar River Basin (KRB; 2007–2009) and four meteorological stations neighboring the basin (1959–2009). Observed daily air temperature and precipitation were input into HBV model. The runoff data in 2007/2008 and 2008/2009 were used to calibrate and validate the model in 2009/2010 and 2010/2011. Generally, the model simulated runoff very well. The annual glacier mass balance and runoff were calculated by the HBV model and were driven by interpolated meteorological data between 1959 and 2009. The calculated glacier mass balances were reasonable, and were compared with nearby glaciers. The results indicate the decreasing trend of mass balance in the Koxkar Glacier, with an average value of ablation of −370.4 mm·a⁻¹ between 1959 and 2009. The annual runoff showed an increasing trend (5.51 mm·a⁻¹). Further analysis showed that the runoff is more sensitive to temperature than precipitation in KRB.

1. Introduction

Glaciers concentrated in mountain regions are important water resources in arid and semi-arid areas in west China [1,2,3]. Rising temperatures have resulted in the loss of more than half of the glaciers in Asia since 1900 [4]. In the northwest arid region of China, the importance of meltwater is particularly prominent. Meltwater accounts for 22% of the runoff from mountains. Melted water from glaciers from the Tianshan Mountains supplied 38.5% of the total runoff in the Tarim River basin [5]. Since the 1980s, runoff has significantly increased, flowing out from the Tianshan Mountains. The glacial retreat has major impact on water resources in arid regions [6].

Hydrology in the arid region has become an important topic [7]. The hydrological model can simulate hydrological processes; it has practical significance and scientific value for exploring new ideas about management of water resources in the river basin [8]. However, due to the extremely high altitude and difficult environment of the Tianshan Mountains, hydro-meteorological data are rare in this areas, especially in the glaciated watersheds. Additionally, there have been few applications of hydrological models to the south slope of Mount Tumur, Tianshan Mountains. There are more than 2.52 million people and 1.33 million acres of irrigated agricultural land in the lower reaches of the river basin which depend on the runoff coming from high mountains. Therefore, it is urgent to investigate water resources from glaciers in this area. It will also help understand the changes in water resources under the climate change.

The HBV (Hydrologiska Byrans Vattenavdelning) model was developed by the Swedish Meteorological and Hydrological Institute (SMHI) for cold regions in 1973. The HBV model, with its strong practical aspects and improved simulation results, has been applied successfully to forecast runoff for more than 200 river basins in Sweden and has been successfully applied in more than 40 countries of the world to study hydrological processes under climate change, changes of permafrost, and changes of land use [9,10,11,12,13,14]. The HBV model is a semi-distributed conceptual model; it produces better results of simulation when less data input is required [15]. It has also been successfully used in different watersheds in the cold regions of west China. Kang et al. [16] established a conceptual hydrological model for upstream of Hei river in Qilian Mountains of west China based on the principle of the HBV model. The model is accordance with the characteristics of formation of mountain runoff in study area, where the basin was divided into two basic landscape zones, alpine tundra-snow and mountain vegetation [16]. In addition, the model has been applied to glaciated watersheds in the Tibetan Plateau, where it also performed well; the HBV light model includes a glacier module which is more applicable in mountain regions where glaciers are distributed [17].

In this study, the HBV light model and the principle of water balance were employed to simulate the change in depth of runoff, evaporation, and glacier mass balance (GMB) in the Koxkar River Basin (KRB) between 1959 and 2009. The work aims to reconstruct climate change in the past 50 years in the region in order to understand changes in glacier mass balance and runoff in response to climate change. This work may provide knowledge for water resource management in the arid regions of Central Asia.

2. Study Area and Data Collection

2.1. Study Area

The KRB is located in the south of Mount Tumur of Western Tianshan Mountains; it is one of the source regions of the Akesu River (Figure 1). The area of KRB is 116.5 km², the area of glacier is 89.6 km² (76.9%), and its elevation ranges from 2960 to 6342 m.a.s.l. The length of Koxkar glacier is 25.1 km, the ice volume is 15.79 km³. Annual precipitation is about 533 mm in the KRB, which is supplied mainly by moisture from the Atlantic and Arctic. The precipitation which occurs from June to August accounts for 50%, 70% of precipitation occurs from May to September, and about 30% occurs during the cold season. Annual precipitation near the snow line is between 750 and 850 mm [18,19,20]. There is one automatic weather station (AWS) and one hydrological station located in the terminal of glacier (Base Camp) (2974 m.a.s.l.) (Figure 1), and four other AWSs located in the glacier with elevations of 3200, 3400, 3700, and 4200 m.a.s.l., respectively (Figure 1).

2.2. Geographical Data

The HBV light model requires geographical data including a digital elevation model (DEM), aspect, and land cover. The 30-m-resolution Global DEM (http://datamirror.csdb.cn) and a digital vector map of modern glaciers were obtained from the Chinese Glacier Inventory [20]. The basin was divided by altitude zones, aspect zones, glacier zones and nonglacier zones (Figure 2). The land cover data were taken from the Science Data Center for Cold and Arid Regions (http://westdc.westgis.ac.cn/). The types of land cover include bare land, shrub, grassland, snow cover, glacier, and water (

Table 1. Types of vegetation in the Koxkar River Basin.

Cover Type	Area (km2)	Proportion
Bare Land	9.3	8.0%
Shrub	1.8	1.5%
Grassland	6.15	5.3%
Snow Cover	6.9	5.9%
Glacier	89.6	76.9%
Water	2.7	2.32%

).

2.3. Observations and Meteorological Data

The HBV light model requires daily air temperature, daily precipitation, gradients of temperature and precipitation data for input. Five AWSs have been established in KRB since 2007, they have operated until present (Figure 1). Precipitation was monitored by both manual rain gauge and automatic snow/rain gauge (T-200B weighting rain gauge). This study used air temperature and precipitation from the AWS at the Base Camp (Camp-AWS) as input data, they are of higher quality. Based on the altitudes of the five AWSs, we calculated the gradients of temperature and precipitation: they are 0.66 °C (100 m)⁻¹ and 4.6% (100 m)⁻¹, respectively. Since meteorological observations started in 2007 in the KRB, the long-term air temperature and precipitation (1959–2009) data were reconstructed using four national meteorological stations near to KRB (NMS´s) (Figure 1,

Table 2. Details of the Camp-AWS and the four national meteorology stations (NMS) around the Koxkar River Basin. AWS: automatic weather station.

Station	Longitude (°E)	Latitude (°N)	Altitude (m)	Years
Akesu	80.233	41.167	1103	1959–2009
Baichen	81.900	41.783	1229	1959–2009
Aheqi	78.450	40.930	1985	1959–2009
Zhaosu	81.133	43.150	1851	1959–2009
Koxkar	80.167	41.700	2974	2007–2009

) (for the method see Section 4).

2.4. Hydrological and Glacial Data

The hydrological station was established at the base camp (BC) in 2007 (Figure 1). Water stage was measured by an automatic HOBO Water Level Logger (Onset, Bourne, MA, USA) with an interval of 30 min; velocity of flow was measured manually by a hydrodynamometer. Through several measurements at different water stages, the relationship between water stage and discharge was established. Then, the discharge was calculated at the BC from 2007 to 2011.

3. Hydrologiska Byrans Vattenavdelning (HBV) Model

The HBV model is a semi-distributed model, which means that a catchment can be separated into different elevation and vegetation zones as well as into different subcatchments [21,22]. The HBV light is a conceptual mode which is based on the HBV model [23]. It was coupled with a more detailed snow and glacier melt subroutine, employing the degree-day approach [24]. The subroutine has a physical basis. For this study, the HBV light model, which takes into account various aspect classes in each elevation belt and couples Monte Carlo and genetic algorithms to calibrate the parameters automatically [25,26], was applied to model runoff depth, glacier runoff, and GMB in the KRB. The model is composed mainly of the input data, which were processed as semi-distributed by four modules: the glacier ablation module, soil module, response module, and flow concentration module (Figure 3) [22].

The main formulas of each module are described as follows. The model simulates daily discharge using daily precipitation and temperature as input data. Precipitation is simulated to be either snow or rain depending on whether the temperature is above or below a threshold temperature, TT (°C). All precipitation simulated to be snow, i.e., falling when the temperature is below TT, is multiplied by a snowfall correction factor, SFCF. In the snow and glacier routine, the amount of melting water is calculated by the degree-day method [27], which uses positive degree-day (while daily mean temperature is above 0 °C) multiplied by a factor. This model also considers the influence of different aspects and different melt intensity of glacier and snow [25] as shown in Equation (2). Meltwater and precipitation are retained within the snowpack until the amount exceeds a certain fraction, CWH, of the water equivalent of the snow. Refreezing of meltwater is also considered in this model [28,29,30]. Liquid water within the snowpack refreezes according to Equation (3). Rainfall and snowmelt (P) are divided into water filling the soil box and groundwater recharge depending on the relation between water content of the soil box (SM, in mm) and its largest value (FC, in mm) as shown in Equation (4). Actual evaporation from the soil box equals the potential evaporation if SM/FC is above LP (the threshold of Soil water content and Field water holding capacity), whereas a linear reduction is used when SM/FC is below LP, as shown in Equation (5) [23]. Groundwater recharge is added to the upper groundwater box (SUZ, in mm). PERC (mm·day⁻¹) defines the maximum percolation rate from the upper to the lower groundwater box (SLZ, in mm). Runoff from the groundwater boxes is computed as the sum of two or three linear outflow equations depending on whether or not SUZ is above a threshold value, UZL (mm), as shown in Equation (6). This runoff is finally transformed by a triangular weighting function defined by the parameter MAXBAS, as shown in Equation (7) to give the simulated runoff (mm·day⁻¹). If different elevation zones are used, the changes in precipitation and temperature with elevation are calculated using the two parameters P_CALT (%/100 m) and T_CALT (°C/100 m), as shown in Equations (8) and (9). The long-term mean of the potential evaporation, E_pot, M, for a certain day of the year can be corrected to its value at day t, E_pot(t), by using the deviations of the temperature, T(t), from its long-term mean, TM, and a correction factor, CET (°C⁻¹), as shown in Equation (10) [31,32].

$$\begin{array}{ll}{P}_L(h)=P_t(h)& T_t\left(h\right)>TT;\\ P_S(h)=P_t(h)SFCF& T_t\left(h\right)<TT\end{array}$$

(1)

$$Melt=DDF(T(t)-TT)$$

(2)

$$Refreezing=CFRC{F}_{MAX}(TT-T(t))$$

(3)

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

(4)

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

(5)

$$Q_{GW}(t)=K_{2}SLZ+K_{1}SUZ+K_{0}max(SUZ-UZL,0)$$

(6)

$$\begin{array}{l}{Q}_{sim}(t)={\displaystyle \sum _{i=1}^{MAXBAS}c(i)QGW(t-i+1)}\\ \mathrm{where}c(i)={\displaystyle \underset{i-1}{\overset{i}{\int }}\frac{2}{MAXBAS}}-|u-\frac{MAXBAS}{2}|\frac{4}{MAXBA{S}^2}du\end{array}$$

(7)

$$P(h)=P_0(1+\frac{P{C}_{ALT}(h-h_0)}{10000})$$

(8)

$$T(h)=T_0-\frac{T{C}_{ALT}(h-h_0)}{100}$$

(9)

$$\begin{array}{l}{E}_{pot}(t)=(1+C_{ET}(T(t)-T_M))E_{pot,M}\\ \mathrm{but}0\le E_{pot}(t)\le 2E_{pot,M}\end{array}$$

(10)

The R_eff and the coefficient of determination (r²) are used for assessment of simulations [33]. The formulae are as follows (Equations (11) and (12)).

$$R_{eff}=1-\frac{{\displaystyle \sum (Q_{sim}(t)-Q_{obs}(t)}{)}^2}{{\displaystyle \sum {(Q_{obs}(t)-\overline{Q_{obs}})}^2}}$$

(11)

$$r^2=\frac{{({\displaystyle \sum (Q_{obs}-\overline{Q_{obs}})(Q_{sim}-\overline{Q_{sim}})})}^2}{{\displaystyle \sum {(Q_{obs}-\overline{Q_{obs}})}^2{\displaystyle \sum {(Q_{sim}-\overline{Q_{sim}})}^2}}}$$

(12)

4. Reconstruction of Temperature and Precipitation Data

Understanding the climate, discharge, and GMB in past decades is important for predicting future trends. In this work, we used the HBV hydrological model to simulate runoff and GMB from 1959 to 2009. Long-term observed data are unavailable for the basin, so it was necessary to reconstruct the historical runoff and GMB. The daily air temperature and precipitation in the KRB were reconstructed from 1959 to 2009 based on in-situ meteorological observations in the KRB between 2007 and 2009 and four NWSs around the KRB for the period 1959–2009.

4.1. Reconstruction of Temperature

There are many methods to interpolate temperature, e.g., inverse distance weight, trend surface, kriging interpolation, and cokring interpolation; these methods are usually used to interpolate the spatial distribution of temperature. However, they require a high number of data points for calculation. In this study, we had only four NMSs for interpolation (Table 2). Thus, the previous temperature data were interpolated by the statistical method of multiple regression analysis, which depended only on the temperature recorded at the four stations. Two equations were formulated because the four stations had records of various lengths.

$$T_0=0.3098\times T_1-0.2314\times T_2-0.0304\times T_3+0.7211\times T_4, R^2=0.942$$

(13)

$$T_0=1.2292\times T_1-0.3211\times T_2-0.1725\times T_3, R^2=0.918$$

(14)

T₀ is the daily air temperature at Camp-AWS, T₁, T₂, T₃, and T₄ are the daily air temperatures of Akesu, Baichen, Aheqi, and Zhaosu, respectively. The result calculated by Equation (13) (R² = 0.942) was better than that calculated by Equation (14) (R² = 0.918), so we used Equation (13) to reconstruct the daily temperature. The reconstructed daily temperatures fit well to the observation ones from 2007 to 2009 (R² = 0.91; p < 0.01; Figure 4). The temperature increased 0.32 °C (10a⁻¹; p < 0.01) in the past semi-century in the KRB (Figure 5), which is a little higher than the average change in the southern slope of the Tianshan Mountains. (0.30 °C, 10a⁻¹) [34].

4.2. Reconstruction of Precipitation

The precipitation, which has large temporal variation, was influenced by the terrain in the Tianshan Mountains. The reconstruction of precipitation was conducted as follows. We used stepwise multiple regression analysis to eliminate stations that recorded precipitation farther from the Camp-AWS precipitation in a step-by-step manner. Stations were eliminated in the order of Akesu, Baichen, Aheqi, and Zhaosu, respectively. The correlation between Camp-AWS and Akesu was not significant due to Akesu being located in the lowland city. Therefore, we used data from other three stations located in the mountain regions to interpolate the Camp-AWS precipitation.

The gradient of precipitation and inverse distance weighting (IDW) methods were combined to reconstruct Camp-AWS precipitation from 1959 to 2009. The gradients of precipitation were 2.05 mm/100 m, 3.54 mm/100 m, and 1.6 mm/100 m between Camp-AWS and Baichen, Aheqi, and Zhaosu, respectively, from 2007 to 2009. Based on the three precipitation gradients, we calculated the precipitations at the same elevation as Camp-AWS from 1959 to 2009. Then, we used IDW to calculate the precipitation at Camp-AWS from 1959 to 2009, as shown in Equation (15).

$$P_0=\frac{S_1}{S_1+S_2+S_3}\times \frac{P_{0\mathrm{annual}}}{P_{1\mathrm{annual}}}\times P_1+\frac{S_2}{S_1+S_2+S_3}\times \frac{P_{0\mathrm{annual}}}{P_{2\mathrm{annual}}}\times P_2+\frac{S_3}{S_1+S_2+S_3}\times \frac{P_{0\mathrm{annual}}}{P_{3\mathrm{annual}}}\times P_3$$

(15)

In Equation (15), P₀ is the daily precipitation of Camp-AWS and P₁, P₂, and P₃ are the daily precipitations of the NMSs at Baichen, Aheqi, and Zhaosu, respectively. P_{0 annual}, P_{1 annual}, P_{2 annual}, and P_{3 annual} are the average annual precipitation of Camp-AWS, Baichen, Aheqi, and Zhaosu, respectively, between 2007 and 2009. S₁, S₂, and S₃ are the distances between Camp-AWS and Baichen, Aheqi, and Zhaosu, respectively.

Comparison of reconstructed with observed monthly precipitation is shown in Figure 6 (R² reached 0.91; p < 0.01). Because precipitation is very important for the simulation of runoff and GMB, we needed to obtain a more vigorous validation for the reconstructed result. Thus, the coefficient of efficiency (R_eff); Equation (11) was used to evaluate the reconstructed daily result. R_eff was 0.88 from 2007 to 2009, which is fairly satisfactory considering the great varying of precipitation in the mountainous regions (Figure 6). By applying the above method, the annual precipitation at Camp-AWS from 1959 to 2009 was reconstructed (Figure 7). The precipitation increased by 4.8 mm·a⁻¹ (p < 0.05), consistent with other studies in the regions [34].

5. Results and Discussion

5.1. Model Calibration and Validation

Except for meteorological data, there are 13 parameters in the HBV model that need to be calibrated (

Table 3. Definition of and set of optimized parameters of HBV for the Koxkar River Basin.

Parameter	Description	Unit	Range	Value	Method of Estimation
Meteorological data
PCalt	Gradient of precipitation	%/100 m	4.6	4.6	Observation
TCalt	Gradient of temperature	°C/100 m	0.66	0.66	Observation
Snow and glacier routine
TT	Threshold temperature	°C	−1–2.7	2.7	[20]
DDF	Degree-day factor of snow	mm·°C−1·day−1	9.5	9.5	Observation
CFR	Refreezing coefficient	-	0.2	0.2	[19] and calibration
SFCF	Snowfall correction factor		0–1	0.05
Cg	Factor for increased melt of ice	-	1.3	1.3	[30] and calibration
Ca	Factor for increased melt of south slope to north slope	-	1.1–2	1.1	Calibration
Soil routine
FC	Maximum soil box (SM)	mm	100–450	300	[36] and calibration
LP	SM threshold for reduction of evaporation	-	0.3–1	0.9	[36] and calibration
BETA	Shape coefficient	-	1–5	1.5	[36] and calibration
Response routine
K1	Recession coefficient	day−1	0.01–0.2	0.01	[36] and calibration
K2	Recession coefficient	day−1	0.001–0.1	0.05	[36] and calibration
PERC	Maximal flow from upper to lower GW-box	mm·day−1	0–6	2	[36] and calibration
Routing routine
MAXBAS	Routing, length of weighting function	day	1–5	1.9	[36] and calibration

). To obtain a proper set of parameters, we gave every parameter a reasonable range according to the literature. A set of parameters representing the characteristics of the KRB was then automatically obtained by observation and by employing the Monte Carlo method [35,36,37]. We selected hydrological data from 2007/2008 and 2008/2009 for calibration. In addition, data for 2009/2010 and 2010/2011 were used for validation.

For the period of calibration, R_eff and R² are 0.77 and 0.76 for 2007/2008, and 0.87 and 0.84 for 2008/2009, respectively. For the period of validation, R_eff and R² are 0.89 and 0.89 for 2009/2010 and 0.80 and 0.81 for 2010/2011, respectively (Figure 8). The results are acceptable on the daily scale, especially considering the uncertainty in both observation and simulation data in such a mountainous area. The modeling result reflects the actual runoff well (Figure 9). All R_eff values of the daily runoff are greater than 0.76. R_eff is even better (0.95) for the monthly result (Figure 10), and there is an only 6% difference between the observed and simulated annual runoff. Considering the dominant role of glacier runoff in the KRB, the reason for the above differences may be that increasing the temporal resolution led to decreases in the accuracy of the degree-day model [38].

5.2. Simulation of Runoff and GMB for 1959–2009

We assumed that the set of parameters did not change from 1959 to 2009, and the depth of runoff and evaporation were simulated (Figure 6a,b) using the reconstructed data. The runoff depth increased 5.5 mm·a⁻¹ (p < 0.01), and it increased more significantly after the 1990s than before the 1990s. The simulated result is consistent with the results for other watersheds in which runoff is observed in the Tianshan Mountains. [39]. Evaporation increased by 3.5 mm·a⁻¹ in the KRB. On the basis of the water balance principle shown in Equation (16), the change in water volume in the basin (ΔS) equals the total precipitation (P) minus the runoff (R), soil evaporation and sublimation (E_s), and evaporation at the glacier surface (E_g).

$$\Delta S=P-R-E_s-E_g$$

(16)

$$\Delta S=\Delta S_s+B$$

(17)

$$B=P-R-E_s-E_g$$

(18)

ΔS is the sum of the soil moisture change (ΔS_S) and the GMB, as shown in Equation (17). The long-term GMB can be calculated according to Equation (18). Because there is no permafrost in the KRB, it can be assumed that ΔS_S did not change. The contribution of the change in ice content in permafrost to the water balance is less than 1% on a yearly scale [40]. Precipitation was calculated by considering the precipitation gradient. E_S was simulated, and E_g was regarded as constant at 137.2 mm (Figure 11c) [41]. Because there is no long-term observation of GMB in the KRB, we chose the two nearest glaciers (No. 1, in the headquarter of the Urumuqi river, and Tuyuksu) in the Tianshan Mountains, which have been observed in situ, for comparison from 1959 to 2009. As can be seen from

Table 4. Comparison of simulated and observed glacier mass balance accumulation between Koxkar Glacier and adjacent glaciers.

Glacier	Koxkar (Simulation)	No. 1 Glacier (Observation)	Tuyuksu Glacier (Observation)
1959–2009 glacier mass balance accumulation	−18,530 mm	−15,200 mm	−20,384 mm

, the value of GMB in Koxkar Glacier is between the values for the No. 1 and Tuyuksu glaciers, thus the simulated result is reasonable and credible. Glacial ablation in the KRB increased at a rate of 4.2 mm·a⁻¹ (p < 0.01) from 1959 to 2009, and the average elevation of glacier decreased by 18.53 m (16.68 m w. Equation), or at a rate of 370 mm·a⁻¹ (Figure 11c).

5.3. Sensitivity Analysis of Runoff

In this study, we used the HBV model to analyze the sensitivity of runoff to changes of temperature and precipitation in the KRB. Climate change in the area is manifested mainly as changes in temperature and precipitation, and it has had a significant impact on the hydrological processes of the watershed, which is partly covered by glacier. Since the required input of the HBV model is daily data and the time interval of predicted meteorological data is monthly, in order to quantitatively study the sensitivity of runoff to changes in temperatures and precipitation, we assumed that the set of parameters and the glacier did not change in the future. Thus, the temperature was increased by 1, 2, and 4 °C, the precipitation was increased by 10%, 20%, and 30%, and the runoff was simulated by the HBV model.

When the temperature increased by 1 °C to 4 °C (T + 1 °C, T + 2 °C, T + 4°C), the annual runoff increased by 27.4%, 56.2%, and 127%, respectively (Figure 12,

Table 5. Change of runoff under various changes of temperature.

	T + 1 °C	T + 2 °C	T + 4 °C
Changes of runoff (%)	27.4	56.2	127.2

). As can be seen from Figure 12, in the case of rising temperature, increases in the amount of runoff occur mainly from May to October. The month with the largest runoff is August, which shows that the runoff was supplied mainly by meltwater from the glacier. The distribution of internal annual runoff changes when temperature is higher (Figure 12). Significantly increased runoff normally appears in May, and significantly decreased runoff normally appears in October. However, when temperature is higher, the significantly increased and decreased runoff periods shift to earlier and later, respectively.

When precipitation increases by 10%, 20%, and 30% (P + 10%, P + 20%, P + 30%), monthly runoff increases by 6.1%, 11.3%, and 19.8%, respectively (Figure 13,

Table 6. Change of runoff under various changes of precipitation (P).

	P + 10%	P + 20%	P + 30%
Change of runoff (%)	6.1	11.3	19.8

). Increased amounts of runoff occur mainly in the ablation season (from May to October). However, when precipitation increases, the runoff does not increase very significantly. This is mainly because the runoff in the Tianshan Mountains is limited by temperature conditions. When only precipitation increases, runoff does not have a significant increase. Under the precipitation increases, increased runoff is concentrated mainly in June, July, and August. Within this period of time, radiation is sufficient and increased precipitation on the glacier leads to melt, which leads to increased runoff in the watershed. Increased precipitation causes a reduction in runoff from October to April (Figure 13). The increase of air temperature and precipitation has a great influence on the seasonal distribution of runoff. As the temperature increases, the increased degree of ablation decreases. The coefficient variations (Cv) of runoff increase at first and then decrease (Figure 14). However, the differences in Cv are not obvious when temperature increases 1 °C and 2 °C, indicating that the fluctuation of runoff is not significant in this situation. When the temperature increases 4 °C, Cv decreases, and the fluctuation of the runoff decreases. With the increase of precipitation, Cv becomes greater, and the fluctuation of runoff increases, but fluctuation is lower than that caused by the changes of temperature.

5.4. Discussion

Uncertainty of reconstructed meteorological data. Precipitation determines the accumulation of glacier, and temperature determines the melting. Additionally, the amount and distribution of interannual variability of precipitation affects the supply and activity of glacier [42]. In this study, we used polynomial regression and IDW based on gradients to reconstruct temperature and precipitation in the KRB from 1959 to 2009. The accuracy of the reconstructed temperature is higher than precipitation; this is mainly due to the temperature being influenced by single factors in the mountainous region [43]. The accuracy of the reconstructed monthly precipitation is better than that of daily precipitation. Precipitation in mountains is influenced by many factors, e.g., elevation, slope, aspect, surface conditions, etc. [44,45]. We considered only the influence of elevation in this paper and did not consider other factors. However, although the interpolated accuracy of precipitation was influenced, the result can reflect the long-term trend of precipitation.

Temporal and spatial differences of gradients of temperature and precipitation. Gradients of precipitation and temperature are important input data for the HBV model, and these are fixed values in the model. However, gradients of precipitation and temperature have some differences under different altitudes and surface conditions [46]. Results show that the gradient of temperature has significant differences in areas of human activity, forest vegetation zones, and glaciated areas in each month, which can be up to 0.62 °C 100 m⁻¹ at the same altitude in different months [47]. In the KRB, there are only six types of land cover and there is no human activity. Additionally, the area of glacier accounts for 76.9% in the watershed, so we assume that differences in surface conditions have little influence on the gradient of temperature. In the Tianshan Mountains, the gradient of precipitation is more complex than the gradient of temperature. Differences in lapse rate of precipitation are not only in different aspects, in different aspects of the same month they are also significantly different. The gradient of precipitation in different months show differences at different elevations [45,48]. Because of limitations from observation conditions, the gradients of precipitation and temperature in the KRB are averages, which makes the simulated process of runoff uncertain at the daily scale.

Temporal and spatial differences of degree-day factor. The degree-day factor of glacier and snow was an average value in the HBV model when we simulated runoff in the KRB. The degree-day factor is a key parameter in the degree-day model, and, with improved time resolution, the accuracy of the degree-day factor decreases [38]. The degree-day factor is a sensitive parameter in the HBV model when runoff in a watershed includes glacier cover. It is influenced by temperature, altitude, slope, and aspect [49,50]. Previous studies have indicated that the degree-day factor shows temporal and spatial differences in the Tianshan Mountains and Tibetan Plateau [19,51,52]. The area of the KRB is 116.7 km², and the value of degree-day factor is 9.5 mm·°C⁻¹·day⁻¹ in the model. Therefore, the monthly runoff is more accurate than the daily runoff, which was simulated by the HBV model. The trend of the GMB results is consistent with increased glacier ablation in the Tianshan Mountains and is reasonable and credible at the annual scale.

Influence of debris on glacial runoff. Rock glaciers are expected to become increasingly important as current warming continues; they are an important part of the cryospheric water store. Research suggests that they could be important for local water management [53,54]. Glacier melt can be significantly altered by debris cover, but knowledge about the fraction and thickness of debris cover on the Tian Shan glaciers is still sparse [55]. The structure of glacier is complex, the storage and release of water within the large complex glacier will have a great impact on the runoff process. The debris can temporarily store large amounts of meltwater, it regulates runoff at the end of the glacier in monthly and seasonal scales. The accumulation and release of this water is a very complex process in glaciers, it is not linearly applied to a regulation factor that can be accurately described. Therefore, the process of simulated results cannot be accurately fitted to the observed runoff and this is an important factor in errors of modeling.

6. Conclusions

In this study, variations of runoff and glacial characteristics over the past 50 years were simulated and the sensitivities of runoff under different temperature and precipitation variations in the KRB were analyzed. Temperature and precipitation were found to increase by 0.32 °C (10a)⁻¹ and 4.8 mm·a⁻¹, respectively, from 1959 to 2009 at Camp-AWS in the KRB. Runoff depth increased at 5.6 mm·a⁻¹. The simulated GMB of the KG was approximately 370 mm·a⁻¹. The glaciers thinned 18.53 m from 1959 to 2009. When the temperature was increased by 1, 2, and 4 °C, the annual runoff increased by 27.4%, 56.2%, and 127%, respectively. When the precipitation increased by 10%, 20%, and 30%, the monthly runoff increased by 6.1%, 11.3%, and 19.8%, respectively. The runoff is more sensitive to variation of temperature than variation of precipitation. Increases in precipitation and temperature will lead to changes in the internal annual distribution of runoff. Temperatures, precipitation, and other climatic and topographical factors are not isolated, changes of these factors are closely related and affect the process of runoff. The drainage system and internal complex storage of large-scale glacier surfaces lead to the differences between simulated runoff and observed runoff.