Runo ﬀ Changes from Urumqi Glacier No. 1 over the Past 60 Years, Eastern Tianshan, Central Asia

: Glaciers are vital to water resources in the arid land of central Asia. Long-term runo ﬀ records in the glacierized area are particularly valuable in terms of evaluating glacier recession and water resource change on both a regional and global scale. The runo ﬀ records of streams draining basins with 46% current glacier cover, located at the Urumqi Glacier No. 1 in the source area of the Urumqi River in eastern Tianshan, central Asia, were examined for the purpose of assessing climatic and glacial inﬂuences on temporal patterns of streamﬂow for the period 1959–2018. Results suggest that runo ﬀ from the catchment correlates well with temperature and associated precipitation data. During the period 1993–2018, it increased by 114.39 × 10 4 m 3 , which was 1.7 times the average runo ﬀ during the period 1959–1992. A simple water balance model is introduced to calculate the di ﬀ erent components of the runo ﬀ , including precipitation runo ﬀ from glacier surface and from nonglacial areas, glacier mass balance and glacial runo ﬀ . Thus, the long-term change of each component and its response to climate change are revealed. We found that the period 1997–2018 is likely to be the “peak water” (tipping point) of the glacial runo ﬀ resulting from shrinkage of glacier area. A simple model glacier surface ( R ) and from nonglacial ( R ng ), glacier mass MB and glacial runo ﬀ R g . The long-term changes of R and its components, as well as their response to climate change are revealed.


Introduction
Glaciers are vital to water resources, particularly in the arid land in Central Asia [1][2][3]. Most glaciers on Earth have retreated due to global atmospheric warming during the 20th century. Mountain glaciers around the globe have responded strongly to recent climate change and are expected to experience continued mass loss and retreat throughout the 21st century [4][5][6][7][8], which inevitably affects the water supply on their downstream communities. Glacier evolution and its impacts on hydrology is an urgent and critical question in practice as well as for scientific research, and presently, worldwide studies on glacial runoff are based on detailed observations from a rather small number of glacierized catchments or glaciers [9][10][11][12][13][14][15][16][17][18][19]. Consequently, detailed individual glacial runoff observation is imperative for the evaluation of glacier recession and water resource change on both a regional and global scale.

Water Balance Model
In order to present the observed long-term fluctuation of streamflow from the UG1 catchment and to analyze its interaction with temperature and precipitation, we introduce a simple water balance model to separate the different components of the streamflow. This allows us to characterize the temporal variability and governing factors of the different contributors. This model was previously partially employed and described for the UG1HMS catchment by Yang [20] and Li et al. [25]. By ignoring the spatial variation in the catchment and the temporal variability within a year for the variables in the model, this model is able to focus only on the total amount of the variables over the catchment in a yearly time scale. The annual water balance during a water balance year or a glacier mass balance year (from 1 September in previous year to 31 August) at the UG1HMS can be described by where Rg is glacial runoff originated from glacierized area, R is total runoff measured at the gauging station and Rng is the nonglacial runoff derived from precipitation over the nonglacierized area in the catchment. In a simplified condition, Rng can be calculated from = × − × ng (2) where P is precipitation in the catchment, Ac and Ag are the areas of the catchment and of the glacier, respectively. The variable αng is the runoff coefficient over the nonglacierized area (ratio of observed runoff and precipitation over the nonglacierized area).

Water Balance Model
In order to present the observed long-term fluctuation of streamflow from the UG1 catchment and to analyze its interaction with temperature and precipitation, we introduce a simple water balance model to separate the different components of the streamflow. This allows us to characterize the temporal variability and governing factors of the different contributors. This model was previously partially employed and described for the UG1HMS catchment by Yang [20] and Li et al. [25]. By ignoring the spatial variation in the catchment and the temporal variability within a year for the variables in the model, this model is able to focus only on the total amount of the variables over the catchment in a yearly time scale. The annual water balance during a water balance year or a glacier mass balance year (from 1 September in previous year to 31 August) at the UG1HMS can be described by where R g is glacial runoff originated from glacierized area, R is total runoff measured at the gauging station and R ng is the nonglacial runoff derived from precipitation over the nonglacierized area in the catchment. In a simplified condition, R ng can be calculated from where P is precipitation in the catchment, A c and A g are the areas of the catchment and of the glacier, respectively. The variable α ng is the runoff coefficient over the nonglacierized area (ratio of observed runoff and precipitation over the nonglacierized area). Because glacial runoff R g is defined as the runoff solely from the glacier, which includes precipitation runoff derived from the glacierized area (R pg ) and glacier mass budget (balance) (MB), then R g also can be described by where R pg can be calculated from where α g is the runoff coefficient for the glacier (ratio of observed runoff and precipitation on the glacial surface). According to Equations (1)-(3), MB can be calculated from In our case, MB is also obtained from in situ measurement on the UG1 by using the glaciological method (stake/snowpit method), which provides an independent parameter to validate this model.

Data Sets
In the water balance model described, the runoff (R) measurements at the UG1HMS began in 1959 and were interrupted during 1967-1979, but the data were subsequently reconstructed using meteorological data [20]. The observations at the gauging stations are carried out from May to September each year, where the observed water level records are converted to discharges based on rating curves, and then converted into total daily, seasonal and annual volumes in cubic meters or annual runoff depths (dividing annual volumes by catchment areas) in mm. Over 95% of the annual runoff at the stations occurs during the observation period, whereas for the rest of the year, the stream is mostly frozen. The runoff data have been compiled and internally published in annual reports of the Tianshan Glacier Station from 1980-2018 [25].
Because the amount of precipitation in the UG1HMS catchment (P) is crucial for the water balance model, an intensive study was carried out to convert the measurements at Daxigou meteorological station to the real amount of precipitation in the catchment [34,35]. By comparing the precipitation volume from different rainfall gauges mounted at Daxigou and within the catchment, a correction coefficient (1.26) of that observed precipitation at Daxigou was determined by the study, which indicates the annual precipitation in the catchment (P) equals 1.26 times that observed at Daxigou. This correction coefficient was affirmed and introduced by subsequent research [20,25]. It is also adopted in our study.
The area of UG1 (A g ), which was continuously shrinking throughout the entire period of observation, was measured from time-to-time by a variety of techniques such as theodolite, ground-based stereophotogrammetry, total station, GPS-RTK, remote sensing cartography and laser scanning. Based on these measurements, 14 topographic maps and the area data of the glacier were generated from 1962 throughout 2018. In our study, a regression equation (R 2 = 0.98, n = 14, p < 0.01) was created according to the area data shown in Figure 2, which was used to estimate the annual values of A g . Figure 2 indicates the shrinkage of the area is accelerating, which results from glacier mass loss and geometric adjustment to the mass loss through a glacier dynamic process [26]. The runoff coefficients for the area on the glacier (αg) and over the nonglacierized area (αng) are imperative input parameters in the water balance model. Based on long-term observation and a comprehensive analysis, Yang et al. [20] found that the values of these parameters varied in a relatively wide range for individual precipitation, but were relatively constant on interannual scale, ranging between 0.51-2.16, with an average value of 0.9 from 1959 to 1987. The best values of αg and αng for the UG1HMS catchment were thus estimated to be 0.85 and 0.7, respectively [20]. According to previous studies [36][37][38][39], the runoff coefficients are mainly determined by the magnitude of evapotranspiration. The interannual variations in terrestrial water storage and evapotranspiration are insignificant in a basin in a cold environment, like the UG1HMS catchment, which is covered by bare rock with little vegetation and has extensive snow and glacier cover. This might explain the relatively stable and high runoff coefficient of the UG1HMS catchment.
However, some more recent research suggested that the recent climate warming would have enhanced surface evaporation and thus decreased the values of the runoff coefficients [25,40]. To assess these runoff coefficients, test runs of the model were conducted in our study. We found that when the known pair of runoff coefficients is applied to the model, the modeled glacier mass balance agrees well with the glaciological mass balance ( Figure 3). This result suggests that this known pair of runoff coefficients is still suitable for the UG1HMS catchment.  The runoff coefficients for the area on the glacier (α g ) and over the nonglacierized area (α ng ) are imperative input parameters in the water balance model. Based on long-term observation and a comprehensive analysis, Yang et al. [20] found that the values of these parameters varied in a relatively wide range for individual precipitation, but were relatively constant on interannual scale, ranging between 0.51-2.16, with an average value of 0.9 from 1959 to 1987. The best values of α g and α ng for the UG1HMS catchment were thus estimated to be 0.85 and 0.7, respectively [20]. According to previous studies [36][37][38][39], the runoff coefficients are mainly determined by the magnitude of evapotranspiration. The interannual variations in terrestrial water storage and evapotranspiration are insignificant in a basin in a cold environment, like the UG1HMS catchment, which is covered by bare rock with little vegetation and has extensive snow and glacier cover. This might explain the relatively stable and high runoff coefficient of the UG1HMS catchment.
However, some more recent research suggested that the recent climate warming would have enhanced surface evaporation and thus decreased the values of the runoff coefficients [25,40]. To assess these runoff coefficients, test runs of the model were conducted in our study. We found that when the known pair of runoff coefficients is applied to the model, the modeled glacier mass balance agrees well with the glaciological mass balance ( Figure 3). This result suggests that this known pair of runoff coefficients is still suitable for the UG1HMS catchment. The runoff coefficients for the area on the glacier (αg) and over the nonglacierized area (αng) are imperative input parameters in the water balance model. Based on long-term observation and a comprehensive analysis, Yang et al. [20] found that the values of these parameters varied in a relatively wide range for individual precipitation, but were relatively constant on interannual scale, ranging between 0.51-2.16, with an average value of 0.9 from 1959 to 1987. The best values of αg and αng for the UG1HMS catchment were thus estimated to be 0.85 and 0.7, respectively [20]. According to previous studies [36][37][38][39], the runoff coefficients are mainly determined by the magnitude of evapotranspiration. The interannual variations in terrestrial water storage and evapotranspiration are insignificant in a basin in a cold environment, like the UG1HMS catchment, which is covered by bare rock with little vegetation and has extensive snow and glacier cover. This might explain the relatively stable and high runoff coefficient of the UG1HMS catchment.
However, some more recent research suggested that the recent climate warming would have enhanced surface evaporation and thus decreased the values of the runoff coefficients [25,40]. To assess these runoff coefficients, test runs of the model were conducted in our study. We found that when the known pair of runoff coefficients is applied to the model, the modeled glacier mass balance agrees well with the glaciological mass balance ( Figure 3). This result suggests that this known pair of runoff coefficients is still suitable for the UG1HMS catchment.  The mass balance of Urumqi Glacier No. 1 has been observed since 1959. The mass balance data are available from 1959 to 2018, except for 1967-1979. The missing data from 1966 to 1979 was reconstructed by meteorological data [41]. The single-point specific mass balance has been measured on a monthly basis from 1 May to 31 August by the stake/snow pit method. The stake network consists of not less than 42 stakes drilled into the glacier surface and evenly distributed with elevation. Using the data of specific mass balance and additional snow pits, glacier-wide mass balance, annual net accumulation and ablation during the period from September 1 in previous year to August 31 are calculated by contour maps of accumulation and ablation [21,22,42,43]. Glaciological data have been internally published in annual reports of the Tianshan Glaciological Station from 1980-2018, and in the Glacier Mass Balance Bulletin and Global Glacier Change Bulletin [44].

Results and Discussion
Based on Equations (1), (2), (4) and (5), glacial runoff (R g ), nonglacial runoff (R ng ), precipitation runoff derived from the glacier surface (R pg ) and glacier mass balance (MB) at the UG1HMS are computed. The results show that the contribution of glacial runoff (R g ) to the total runoff (R) over the past 60 years is dominating, accounting for 70%, and the remaining 30% comes from nonglacial runoff (R ng ) derived from precipitation over the nonglacierized area in the catchment. The glacial runoff (R g ) is composed of the glacier mass balance (MB) and the precipitation runoff formed from the glacier surface (R pg ), which, respectively, account for 26% and 44% of the R g on average for the past 60 years. Figure 4 shows the variations of R, R g and MB during the period of 1959-2018. To explore the drivers behind the fluctuation, the average temperature during the maximum glacier-melting period from July to August and the annual precipitation are also plotted in the same figure. The temporal variations of R pg and R ng are shown in Figure 5. Because both precipitation runoff derived from the glacierized area (R pg ) and from nonglacierized areas (R ng ) are related to precipitation, the annual precipitation series is also included in this figure.
Water 2020, 12, 1286 6 of 15 The mass balance of Urumqi Glacier No. 1 has been observed since 1959. The mass balance data are available from 1959 to 2018, except for 1967-1979. The missing data from 1966 to 1979 was reconstructed by meteorological data [41]. The single-point specific mass balance has been measured on a monthly basis from 1 May to 31 August by the stake/snow pit method. The stake network consists of not less than 42 stakes drilled into the glacier surface and evenly distributed with elevation. Using the data of specific mass balance and additional snow pits, glacier-wide mass balance, annual net accumulation and ablation during the period from September 1 in previous year to August 31 are calculated by contour maps of accumulation and ablation [21,22,42,43]. Glaciological data have been internally published in annual reports of the Tianshan Glaciological Station from 1980-2018, and in the Glacier Mass Balance Bulletin and Global Glacier Change Bulletin [44].

Results and Discussion
Based on Equation (1), (2), (4) and (5), glacial runoff (Rg), nonglacial runoff (Rng), precipitation runoff derived from the glacier surface (Rpg) and glacier mass balance (MB) at the UG1HMS are computed. The results show that the contribution of glacial runoff (Rg) to the total runoff (R) over the past 60 years is dominating, accounting for 70%, and the remaining 30% comes from nonglacial runoff (Rng) derived from precipitation over the nonglacierized area in the catchment. The glacial runoff (Rg) is composed of the glacier mass balance (MB) and the precipitation runoff formed from the glacier surface (Rpg), which, respectively, account for 26% and 44% of the Rg on average for the past 60 years. Figure 4 shows the variations of R, Rg and MB during the period of 1959-2018. To explore the drivers behind the fluctuation, the average temperature during the maximum glacier-melting period from July to August and the annual precipitation are also plotted in the same figure. The temporal variations of Rpg and Rng are shown in Figure 5. Because both precipitation runoff derived from the glacierized area (Rpg) and from nonglacierized areas (Rng) are related to precipitation, the annual precipitation series is also included in this figure.

Total Runoff (R)
R demonstrates a significant amplification from 1959 to 2018 with large interannual fluctuation. It can be roughly divided into two stages ( Figure 4). During the first stage from 1959 to 1992, the value of R was relatively low, with an average of 157.39 × 10 4 m 3 . In the second stage from 1993 to 2018, R rises dramatically, particularly in its baseline. The value reaches as high as 271.78 × 10 4 m 3 on average, an increase of 114.39 × 10 4 m 3 , which is 1.7 times compared with the first stage. Today, R appears at a high level with a large interannual variability.
The long-term change of temperature is overall consistent and well correlated with total runoff (R 2 = 0.68, n = 60, p < 0.01) as shown in Figure 4. It also exhibits two stages: rising strongly after 1992 as R does. The average value of the temperature during the period of 1993-2018 was 1.0 °C higher than that in 1959-1992. Because air temperature is well correlated with meltwater in mountains areas [45,46], this indicates that the escalating runoff largely originates from the increased meltwater production as a result of the mass loss of the glacier.
The relationship between precipitation and runoff is more complex. On the one hand, precipitation increases significantly during the period of 1996-2018, with an average annual value of 517 mm, 80 mm higher than that of 1959-1995. This is highly consistent with the overall increase of runoff. On the other hand, the interannual correlation between runoff and precipitation from 1959-2018 is relatively weak (R 2 = 0.35, n = 60, p < 0.01). The reason is perhaps that while enhanced precipitation increases the overall quantity of water entering the watershed, it also diminishes glacier melting on an annual timescale.
Our study indicates that R depends on both the timing and magnitude of precipitation and temperature. It rapidly rises to a high level after 1992 coinciding with the double-stepped increases in temperature and precipitation. The increase in precipitation resulted in the increase of the total amount of water entering the catchment, and the increase in temperature aggravated the melting of snow and ice in the glacierized area.

Precipitation Runoff Derived from the Nonglacierized Area (Rng) and Glacier Surface (Rpg)
The values of precipitation runoff derived from the nonglacierized area (Rng) calculated by Equation (2) ranges between 42.34 × 10 4 m 3 and 91.23 × 10 4 m 3 with an average of 60.51 × 10 4 m 3 over the past 60 years, which accounts for 30% of the total runoff (R) on average. According to Equation (2), Rng is determined by the amount of precipitation, the runoff coefficient of the nonglacierized area

Total Runoff (R)
R demonstrates a significant amplification from 1959 to 2018 with large interannual fluctuation. It can be roughly divided into two stages (Figure 4). During the first stage from 1959 to 1992, the value of R was relatively low, with an average of 157.39 × 10 4 m 3 . In the second stage from 1993 to 2018, R rises dramatically, particularly in its baseline. The value reaches as high as 271.78 × 10 4 m 3 on average, an increase of 114.39 × 10 4 m 3 , which is 1.7 times compared with the first stage. Today, R appears at a high level with a large interannual variability.
The long-term change of temperature is overall consistent and well correlated with total runoff (R 2 = 0.68, n = 60, p < 0.01) as shown in Figure 4. It also exhibits two stages: rising strongly after 1992 as R does. The average value of the temperature during the period of 1993-2018 was 1.0 • C higher than that in 1959-1992. Because air temperature is well correlated with meltwater in mountains areas [45,46], this indicates that the escalating runoff largely originates from the increased meltwater production as a result of the mass loss of the glacier.
The relationship between precipitation and runoff is more complex. On the one hand, precipitation increases significantly during the period of 1996-2018, with an average annual value of 517 mm, 80 mm higher than that of 1959-1995. This is highly consistent with the overall increase of runoff. On the other hand, the interannual correlation between runoff and precipitation from 1959-2018 is relatively weak (R 2 = 0.35, n = 60, p < 0.01). The reason is perhaps that while enhanced precipitation increases the overall quantity of water entering the watershed, it also diminishes glacier melting on an annual timescale.
Our study indicates that R depends on both the timing and magnitude of precipitation and temperature. It rapidly rises to a high level after 1992 coinciding with the double-stepped increases in temperature and precipitation. The increase in precipitation resulted in the increase of the total amount of water entering the catchment, and the increase in temperature aggravated the melting of snow and ice in the glacierized area.

Precipitation Runoff Derived from the Nonglacierized Area (R ng ) and Glacier Surface (R pg )
The values of precipitation runoff derived from the nonglacierized area (R ng ) calculated by Equation (2) ranges between 42.34 × 10 4 m 3 and 91.23 × 10 4 m 3 with an average of 60.51 × 10 4 m 3 over the past 60 years, which accounts for 30% of the total runoff (R) on average. According to Equation (2), R ng is determined by the amount of precipitation, the runoff coefficient of the nonglacierized area (α ng ) and the value of nonglacierized area (A c -A g ) in the catchment. Consistent with the variation of precipitation, the long-term change of R ng can be roughly divided into two periods: from 1959 to 1995 and from 1996 to 2018. During the period of 1996-2018, R ng averaged 72.7 × 10 4 m 3 , compared with 52.9 × 10 4 m 3 during the period of 1959-1995, thus showing an increase of 19.8 × 10 4 m 3 or 27% ( Figure 5).
With respect to α ng in Equation (2), it is constant over time and not temperature dependent in the research catchment because the nonglacierized area mainly includes bare hillside, with sparse vegetation coverage, the loss of precipitation results mainly from evaporation and snow sublimation on stable rock and soil surface with an elevation between 3659 m and 4480 m. Rising temperature may elevate the rainfall ratio in liquid precipitation and enhances the melting of snowpack in the catchment; however, it does not change the annual value of R ng , as all snowpack in the nonglacierized area eventually melts away at the end of summer and supplies to R ng .
The precipitation runoff derived from the glacier surface R pg is calculated by using Equation (4), ranging between 67.13 × 10 4 m 3 and 125.38 × 10 4 m 3 with an average value of 88.58 × 10 4 m 3 over the past 60 years, and accounts for 44% of the total runoff (R). According to Equation (2), R pg is determined by the amount of precipitation, the runoff coefficient of the glacier (α g ) and the glacier area (A g ) in the catchment. The temporal variation of R pg ( Figure 5) is generally consistent with precipitation, but unlike the R ng , the average value in the period of 1996-2018 is 91.69 × 10 4 m 3 , compared with that of 86.64 × 10 4 m 3 in the period of 1959-1995, thus showing just a little increase.
Regarding α g in Equation (4), it is found to be influenced by multiple factors, because the transformation of precipitation to runoff on a glacier surface is more complicated than that on the nonglacierized area. A glacier usually can be divided into accumulation and ablation zones. The dividing line is called the equilibrium line. The surface of the accumulation zone is composed of firn (snow that lasts for more than a year), and the surface of the ablation zone is usually bare ice in summer. All precipitation occurring in the ablation zone is able to convert into runoff at an annual scale, thus making α g a constant value in the ablation zone. However, in the accumulation zone, a portion of snowfall is turned into glacier through the processes of granulation, densification and dynamic metamorphism, which may take several years [47]. The other portion is melted and penetrates into the firn layer on the glacier, which either refreezes and becomes part of the glacier, or forms runoff at the interface between firn and glacial ice and transforms into R pg [48]. The α g in the ablation zone is obviously greater than in the accumulation zone where it is positive, depending on temperature. Furthermore, higher temperature can enlarge the ablation zone and reduce the accumulation zone by raising the equilibrium line (average value is 4073 m a.s.l.), which also increases α g . The α g value (0.85) used in this study is a multiyear observed average value. On an annual scale, α g varies greatly and largely depends on temperature.
According to Equations (2) and (4), the shrinkage of glacier area A g resulted in a reduction of R g and an amplification of R ng . This explains the larger increase in R ng than that in R g during the period of 1996-2018. Regarding the total runoff R, it is reduced by the shrinkage of A g . This is because α g is greater than α ng , and the increase in R ng was not able to compensate for the decrease in R g .

Glacier Mass Balance (MB)
The glacier mainly gains mass (solid precipitation) in the accumulation zone and loses mass through melting in the ablation zone [49]. MB is defined as the net change in mass of a glacier over a balance year (from 1 September in previous year to 31 August). Thus, MB presents the runoff derived from glacier melting itself. A positive value of MB indicates the surplus of glacier mass gained from snowfall. On the contrary, a negative value indicates the loss of glacier mass released to glacial runoff. Over the past 60 years, MB calculated by the water balance model (Equation (4)) shows a strong negative growth tendency (Figure 4), ranging between 66.82 × 10 4 and −213.71 × 10 4 m 3 , with an average of −57.87 × 10 4 m 3 (negative value means mass loses of the glacier). These values were verified by the in situ observation: the mass balance measured by using stake/snowpit method at the same period of time is consistent with calculated glacier mass balance series (Figure 3), ranging from 72.01 × 10 4 to −216.21 × 10 4 m 3 , with an average of −56.12 × 10 4 m 3 . According to the value of MB, glacier melting releases 54.88 × 10 4 m 3 on average to the catchment annually, which accounts for 26% of R in the past 60 years.
The long-term variation of MB is inversely correlated with air temperature (R 2 = − 0.67, n = 60, p < 0.01). According to glacier physics, glacier mass balance is determined by the energy budget (balance) of the glacier surface. The energy comes mainly from solar radiation. As one of the factors in the energy budget, air temperature plays a vital role in glacier melting [26]. If the ice is covered by snow, the energy is first used to melt snow. If there is no snow cover, the ice starts to melt, generating glacial runoff [50]. Actually, to explain mass balance and its interaction with meteorological elements (here referring to temperature and precipitation), a mass balance model coupled with a surface energy balance model on the glacier is required. From a previous study based on a mass balance model, three mechanisms are found to be responsible for an enhanced mass loss of UG1 [26]. The first is air temperature rise during the melting season as described by positive degree days (PDD, the sum of daily mean air temperatures above the melting point in a year) [24,26], which directly enhances the melting rate of the glacier. Due to the rapid ascent of air temperature, the simultaneous increase in precipitation did not prevent UG1 from the acceleration of mass loss in this case. This is because the warming temperature elevates the ratio of rainfall to snowfall precipitation and enhances glacier melting [26]. The second is the ice temperature rise of UG1, which reduces the heat in raising the ice surface to the melting point and in the refreezing of percolation of meltwater as well. As a result, it increases the sensitivity of the glacier to climate warming, making it an important contributor to the accelerated mass loss. The third is the albedo reduction of the glacier surface due to the continued increase of the ablation area. The positive feedback of glacier recession on the ablation area has most likely resulted in an accelerated mass loss. The low albedo in the ablation was found to be caused by mineral dust, which is mainly derived from the surrounding Asian deserts, as well as organic matter, such as living cyanobacteria on the dust surface, which is produced quickly as temperature rises.
The relationship between MB and precipitation is complicated. A higher snowfall can moderate glacier melting thus making MB more positive, and oppositely, a higher rainfall can accelerate glacier melting thus making MB more negative. Additionally, the shrinkage of the glacier area (A g ) results in a reduction in the total amount of glacier melt and leads to an extra positive MB value.

Glacial Runoff (R g )
Originating in the glacier, the glacial runoff R g of UG1 is composed of mass balance (MB) and precipitation runoff (R pg ), and is calculated by subtracting nontotal runoff R ng from total runoff R (Equation (1)). As result, the values of glacial runoff over the past 60 years fluctuate significantly between 25.61 × 10 4 and 304.35 × 10 4 m 3 with an average of 146.45 × 10 4 m 3 , which accounts for 70% of R. This indicates that R is largely determined by R g .
The long-term change of R g and the five-year moving average curve demonstrate an overall rising tendency (Figures 4 and 6). To look into the details, we performed a trend test for glacial runoff (Figure 7). Results from the Mann-Kendall (M-K) mutation test showed that the glacial runoff had an abrupt change in 1992. The only point of intersection between UF and UB (two confidence lines of α = 0.05)occurred in 1992 and the intersection point was within two confidence lines of α = 0.05, indicating that an abrupt change of the glacial runoff occurred in 1992. It was also confirmed by the mathematical analysis of the curve change point. The UF was higher than 1.96 after 1997, and the tendency of UF to continue to increase was obvious (Figure 7a). The minimum flow of its cumulative anomaly curve appeared around 1992 and experienced a fluctuation from low to high (Figure 7b). Both results indicate there was a dividing time line around 1992, and this result passed the significance test of α = 0.01. As a result, similar to R, the glacial runoff can be characterized by two stages: one from 1959 to 1992 and the other from 1993 to 2018. During the first stage, R g was relatively low with an average value of 105.24 × 10 4 m 3 . In the second stage, however, R g experienced sustained growth.
The average value, reaching 200.33 × 10 4 m 3 , increased by 94.64 × 10 4 m 3 and corresponded to a value that was 1.904 times the first stage.  The maximum value of 266 × 10 4 m 3 of the first stage appeared in 1986. A previous study examined the hydrological process of this specific year, in which the mechanism was clearly revealed [51]. The extremely high Rg in 1986 was attributed to record high temperatures that occurred during July-August. The average value of temperature in July-August was 1.1 °C higher than that observed in the past years. As a result, the average monthly runoff was over 2.65 times the monthly average for prior years. Due to the remarkably high temperature, the ablation zone of UG1 in 1986 was largely expanded into the former accumulation zone, which dramatically enhanced the glacier melting. Subsequently, a huge mass of UG1 was lost and transformed into Rg.
During The value of Rg equal the algebraic sum of mass balance (MB) and precipitation runoff derived from glacier surface (Rpg) (Equation (3)), therefore, it should be closely related to air temperature and  The maximum value of 266 × 10 4 m 3 of the first stage appeared in 1986. A previous study examined the hydrological process of this specific year, in which the mechanism was clearly revealed [51]. The extremely high Rg in 1986 was attributed to record high temperatures that occurred during July-August. The average value of temperature in July-August was 1.1 °C higher than that observed in the past years. As a result, the average monthly runoff was over 2.65 times the monthly average for prior years. Due to the remarkably high temperature, the ablation zone of UG1 in 1986 was largely expanded into the former accumulation zone, which dramatically enhanced the glacier melting. Subsequently, a huge mass of UG1 was lost and transformed into Rg.
During The value of Rg equal the algebraic sum of mass balance (MB) and precipitation runoff derived from glacier surface (Rpg) (Equation (3)), therefore, it should be closely related to air temperature and The maximum value of 266 × 10 4 m 3 of the first stage appeared in 1986. A previous study examined the hydrological process of this specific year, in which the mechanism was clearly revealed [51]. The extremely high R g in 1986 was attributed to record high temperatures that occurred during July-August. The average value of temperature in July-August was 1.1 • C higher than that observed in the past years. As a result, the average monthly runoff was over 2.65 times the monthly average for prior years. Due to the remarkably high temperature, the ablation zone of UG1 in 1986 was largely expanded into the former accumulation zone, which dramatically enhanced the glacier melting. Subsequently, a huge mass of UG1 was lost and transformed into R g .
During the second stage from 1993 to 2018, the R g continuously rose and gradually reached its maximum period of 1997-2007 with an average value of 213.02 × 10 4 m 3 , and then demonstrated a slight decreasing tendency with great annual fluctuations from 2008 to 2018, as shown in Figure 4. The average value in this period of time declined to 189.07 × 10 4 m 3 .
The value of R g equal the algebraic sum of mass balance (MB) and precipitation runoff derived from glacier surface (R pg ) (Equation (3)), therefore, it should be closely related to air temperature and the precipitation. Over the past 60 years, R g was found to be well correlated with temperature (R 2 = 0.65, n = 60, p < 0.01) but poorly correlated with precipitation (R 2 = 0.2, n = 60, p > 0.05).
As from previous results, a higher temperature was beneficial to producing both MB and R pg , but as for the precipitation, it was more complex. The effect of snowfall and rainfall in precipitation were different. Firstly, in a year with abundant snowfall, both R pg and R ng had higher values, thus contributing more water to R. However, at the same time, because a higher snowfall can protect the glacier from melting, the value of MB turns to be more positive (the glacier gains more mass), so that the glacier generates less water to R g , so as to R. This is opposite in a year with less snowfall. Therefore, the opposite effect of snowfall on R pg and MB over a glacier can actually stabilize R g , and thereby R. Because of the adjustment capacity to streamflow, the glacier is regarded as "solid reservoir," which is extremely advantageous to the water supply in arid land of central Asia.
Secondly, the rainfall in the precipitation over the glacier can bring heat to the surface and accelerate the glacier melting, so that the effect of rainfall on MB is opposite to that of snowfall. The ratio between snowfall and rainfall is thus significant to R g , which is sensitive to air temperature. The precipitation in the study area mainly occurs in summer. The average value of precipitation from May through September accounts for 87% of the whole year. Therefore, a higher temperature can reduce the value of the ratio of snowfall to rainfall, and therefore promote glacier melting. This phenomenon is quite different from the Alps in Europe, where more precipitation is brought by westerly current and occurs during winter as snowfall, which can efficiently protect the glacier from melting [52,53].
As the result, the amplification of R g during the period from 1992 to 2018 is mainly attributed to the negative increases in MB caused by an elevated temperature, and also the rise in R pg triggered by precipitation increase. In addition, the reduction of the ratio between snowfall and rainfall caused by a rising temperature is believed to play a significant positive role.
Here, special attention must be paid to the period of 2008-2018 in Figure 6. During this time period, contrary to a continuous increase of temperature (as shown in Figure 4), R g shows an unexpected slightly downward trend. Meanwhile, precipitation is fluctuating with an overall descending trend. Two reasons are found to be possibly responsible for the decline of R g . The first is the general decrease in precipitation, which diminishes R pg . The second must be the reduction of the glacier area A g . The shrinkage of A g not only leads to a decrease in R pg according to Equation (4), but also reduces the total amount of melt on the glacier (making MB more positive in value), both of which can be found from the temporal changes of R pg and MB. Due to the reduction of A g , a much larger growth in R ng than that in R pg occurs, as shown in Figure 5. The MB runoff is found to demonstrate a slightly decreasing trend during the period 2008-2018, as shown in Figure 4, which is evidently because of the shrinkage of A g . This process was studied and described by previous research [16,54]. According to the studies, as glaciers recede, water is released from long-term glacial storage (MB). Thus, annual R g volume typically increases until a maximum is reached, often referred to as "peak water," beyond which R g decreases because the reduced glacier area can no longer support rising meltwater volumes. As a matter of fact, the area of UG1 showed a decrease of 21%, from 1.95 km 2 in 1962 to 1.54 km 2 in 2017. If the negative effect of R g was enough to offset the positive effect of rising temperature on R g , the peak water (tipping point) of the glacial runoff already appeared during the period of 1997-2007. To confirm this finding, however, longer observation and more in-depth quantitative analysis are further required. Apparently, the decreasing trend of glacial runoff after 2007 resulted from both diminishing precipitation and decreasing glacier area.

Conclusions and Outlook
The Urumqi River headwater region has the longest runoff records of any glacier in China. This study analyses streamflow at Urumqi Glacier No. 1 gauging station (R). A simple water balance model is used to separate the different components of R, including precipitation runoff from the glacier surface (R pg ) and from nonglacial areas (R ng ), glacier mass balance MB and glacial runoff R g . The long-term changes of R and its components, as well as their response to climate change are revealed.

•
Average total runoff at the UG1HMS was about 206.96 × 10 4 m 3 over the past 60 years, of which the average R g accounted for 70%, among them R pg and MB accounted for 44% and 26%, respectively. The rest was precipitation runoff in nonglacierized areas. R demonstrated a significant increase in the period of 1993-2018, with an increase of 114.39 × 10 4 m 3 , corresponding to 1.7 times R during the period of 1959-1992. This increase was correlated with temperature and associated with precipitation, indicating that both the mass reduction of the glacier and the elevated precipitation were contributors. At present, R and its components are characterized by high values with great annual fluctuation.

•
As important components, precipitation runoff from the glacier surface (R pg ) and nonglacial areas (R ng ) were determined by the amount of precipitation in the catchment, showing a step increase after 1995. A higher temperature can elevate the runoff coefficient on the glacier (α g ) resulting in an enhanced R pg , and vice versa. The recession of the glacier area, which reduces R g and increases R ng , had an overall negative impact on the R value.

•
The modeled glacier mass balance MB was inversely correlated with air temperature and showed a strong negative growth trend over the past 60 years. For precipitation, A higher solid precipitation can moderate glacier melting thus making MB more positive, and oppositely, a higher liquid precipitation can accelerate glacier melting thus making MB more negative.

•
Over the past 60 years, the long-term change of glacial runoff showed an increasing trend, particularly after 1992. After reaching its maximum during the period 1997-2007, it decreased slightly from 1998 to 2018. We found that the general reduction of precipitation and particularly the shrinkage of the glacier area were responsible for the downward trend. The period of 1997-2007 is likely to be the "peak water" (tipping point) of the glacial runoff, as described by Huss et al. (2018) [16]. However, to verify it, longer observation and a more in-depth quantitative analysis are required.
In this study, to reveal a long-term change of runoff over the glacierized area and minimize the inherent temporal noise of the high-resolution time series, we analyzed the annual periods. However, the relationship and the interactive mechanism between the runoff or its components and meteorological elements are better described using a higher temporal resolution scale. In addition, the division of time period in this study is inevitably subjective and arbitrary. For this reason, a higher temporal resolution of the hydrological model needs to be introduced in future studies.