Recent Changes in Water Discharge in Snow and Glacier Melt-Dominated Rivers in the Tienshan Mountains, Central Asia

Global warming has generally led to changes in river runoffs fed by snow and glacier meltwater in mountain ranges. The runoff of the Aksu River, which originates in the Southern Tienshan Mountains, exhibited a positive trend during 1979–2002, but this trend reversed during 2002–2015. Through a comprehensive analysis, this study aims to estimate potential reasons for changes in the runoff of its two contrasting headwaters: the Toxkan and Kumalak Rivers, based on climatic data, the altitude of the 0 °C isotherm, glacier mass balance (GMB), snow cover area (SCA), snow depth (SD) and the sensitivity model. For the Toxkan River, the decrease in spring runoff mainly resulted from reductions in precipitation, whereas the decrease in summer runoff was mainly caused by early snowmelt in spring and a much-reduced snow meltwater supply in summer. In addition, the obvious glacier area reduction in the catchment (decreased to less than 4%) also contributed to the reduced summer runoff. For the Kumalak River, a sharp decrease rate of 10.21 × 108 m3/decade in runoff was detected due to summertime cooling of both surface and upper air temperatures. Reduced summer temperatures with a positive trend in precipitation not only inhibited glacier melting but also dropped the 0 °C layer altitude, resulting in a significant increase in summertime SCA and SD, a slowing of the glacier negative mass balance, and a lowering of the snow-line altitude.


Introduction
Mountain snow and glacier are an important component of the cryosphere and are considered as natural indicators of climatic change [1]. In recent years, a warming climate globally has generally caused significant changes in river runoff and regional hydrological cycles [2,3], which generally involves the water interaction between the atmosphere and earth, as well as water movement and storage on and under the surface of the earth [4]. Changes in the hydrological cycle have driven the alpine water resources, especially for rivers fed predominantly by snow and glacier meltwater [3, 5,6]. Snow and glacier meltwater play an important role in runoff formation, as their ablation and accumulation directly affect the hydrological processes and could serve as a regionally important buffer against drought stress [7][8][9][10]. However, under the warming climate scenario, the vulnerability and uncertainty of water resources from mountains are increasing, and their changes may cause a major consequence for water supplies downstream [5,11]. cycle in mountain areas has intensified and the rivers have been in a situation of high variability [60,69]. Rivers that are predominantly recharged by glacier meltwater have shown a noticeable increasing trend in runoff from June to September [70]. In addition, the magnitude and frequency of flooding have also increased [71], and the floods are often accompanied by an extreme glacier and snow-related events, such as a sharp increase in meltwater and glacier lake outburst floods [72][73][74].
In Central Asia, the Aksu River is a transboundary river originating in Kyrgyzstan and flowing into China. It provides about 70-80% of water to the Tarim River [24,71], which is the largest inland river in China, and thus plays a critical role in the socio-economic development and ecological security of the Tarim River basin [6,75,76]. The Aksu River has two headwaters, i.e., the Toxkan and Kumalak Rivers. The catchments of these two water bodies are covered by large areas of snow or glacier and have shown remarkably different hydrological responses to climate change [8]. Over the past 40 years, an obvious warming trend has been detected in the Aksu River catchment at a rate of 0.34~0.38 • C/decade which is notably higher than other regions of Northwestern China (0.32 • C/decade) [77] or Central Asia (0.31 • C/decade) [78]. Our results detect that from 1979-2002, runoff from the Aksu River has been increasing at a rate of 12.34 × 10 8 m 3 /decade. For the Toxkan River, the annual runoff has increased at a rate of 5.71 × 10 8 m 3 /decade, while for the Kumalak River, runoff has increased 6.63 × 10 8 m 3 /decade. It is worth noting that, compared to 1979-2002, the runoffs of both the Toxkan River and Kumalak River have been decreasing since 2002 at rates of −10.17 × 10 8 m 3 /decade and −10.21 × 10 8 m 3 /decade, respectively. This is a very interesting phenomenon and leads to the question of why this is happening under a warming climate.
To find comprehensive trend attributions for factors affecting runoff changes in the Toxkan and Kumalak Rivers during 1979-2002 and 2002-2015, we analyzed several variations in climate-related drivers. These include temperature, precipitation, and altitude of the 0 • C isotherm, snow cover area (SCA) and snow depth (SD), and glacier mass balance (GMB) according to statistical mechanics and the maximum entropy principle (SMMEP) model, In addition, we discuss the impacts of these drivers on hydrological processes, look for potential links, and examine responses to the drivers' variabilities. Ultimately, this investigation and comprehensive analysis will provide scientific evidence for the allocation and management of water resources in this area of the Tienshan Mountains, Central Asia.

Research Area
The Aksu River is a typical mid-latitude alpine transboundary river originating in Kyrgyzstan and flowing through the middle of the Southern Tienshan Mountains (Figure 1a). The river's basin covers an area measuring about 5.0 × l0 4 km 2 , with an average elevation of 2233 m a.s.l. The study region comprises two headwater catchments, namely the Toxkan River catchment and the Kumalak River catchment ( Figure 1a; Table 1). The Toxkan River originates in the Artbash Mountains (with the highest peak about 5897 m a.s.l.), is 457 km long, and covers a catchment area of 1.84 × 10 4 km 2 , with an average altitude of 3567 m a.s.l. In the 2010s, this area had a glacier cover area of approx. 3.66% and an average single glacier area of about 0.69 km 2 [79]. The Kumalak River, which originates in the Hantengri Mountains, is 293 km long and covers a catchment area of 1.28 × 10 4 km 2 . This catchment features the Tienshan Mountains' highest peaks-Tomor Peak (7439 m a.s.l.), which is also known as the Tomor Feng in China. A considerable portion of the Kumalak River catchment is glacierized (Figure 1a) and the average altitude is 3730 m a.s.l. Above this altitude, large glaciers (about 1.36 km 2 ) abound, covering 16.34% of the area in the 2010s (Table 1). Figure 1b,c shows the distribution area for different altitudes bands in the Toxkan and Kumalak catchments. These two branches-the Toxkan and the Kumalak Rivers-merge into the Aksu River at Xiaoxiake station, flowing from there to the Tarim River.
The Aksu River is a typical river in the Tienshan Mountains, in that the primary water sources of the river's runoff are mainly supplied by snow and glacier meltwater [26,80]. The Aksu river's runoff is mainly concentrated in summer, accounting for 65% of the river's total discharge. The contribution of snow-glacier meltwater to runoff is up to 45%, with precipitation accounting for about 33.1% [80]. For the Toxkan River, glacier meltwater accounts for only about 21.3% of the total annual runoff [81],

Meteorological and Streamflow Data
In the study, the daily mean surface precipitation and temperature from five meteorological stations (Aheqi, Aksu, Tuergate, Baicheng, and Bayinbuluk, 1104~3504 m a.s.l.) in China during 1979-2015 were provided by the National Meteorological Information Center (Table 2). Considering the scarcity of meteorological stations in this region, very few of the stations have a long time series. Hence, they can't adequately represent the precipitation changes in the alpine regions spatially and comprehensively. Consequently, the current study used monthly precipitation from the GPCC monthly precipitation product V.2018 (V8, available in 0.25 • resolution) for the period 1979-2015. This product is based on data from more than 79,000 stations spread across the global land surface. The gridded precipitation dataset is obtained from the Global Precipitation Climatology Centre (GPCC), which is operated by Deutscher Wetterdienst under the auspices of the World Meteorological Organization (WMO). It has higher accuracy with its spatial resolution of 0.25 • and should, therefore, perform better with regard to real amounts of precipitation in mountain regions with complex topography [85]. Other data include monthly gridded surface temperature (two-meter air temperature) from the ERA-Interim dataset, which is downloaded from the European Center for Medium-Range Weather Forecasts. It provides a 0.25 • horizontal resolution for the period 1979-2015, which was calculated using daily mean temperatures. Compared to other gridded datasets, ERA-Interim exhibits higher spatial resolution and better reflects the spatiotemporal distribution of air temperature in the Tienshan Mountains [86]. Therefore, the gridded GPCC precipitation and ERA-Interim temperature data have been adopted in this study to analyze climate variations in the Aksu River basin.
To further estimate the relationship between the variations of 0 • C isotherm and runoff in the Kumalak River during the period 1979-2015, the 12-h data (pressures, temperatures, and heights) from the Kuqa radiosonde station (code 51644) over the summer months (June to August) for the period 1979-2015 were obtained from the University of Wyoming. The monthly cumulative streamflow data from 1979-2015 for the Shaliguilanke hydrological station (Toxkan River) and Xiehela hydrological station (Kumalak River) used in this study were obtained from the Aksu River Management Bureau.

LSDPC Snow Depth and MODIS Snow Cover Data
To include a long time series of observation of SCA and SD change in our study, the Long-Term Snow Depth Product for China (LSDPC), provided by the Cold and Arid Regions Sciences Data Center at Lanzhou during 1979-2015, was used here, with a spatial resolution of 0.25 • . This original dataset was retrieved from the scanning multichannel microwave radiometer (SMMR, 1979(SMMR, -1987, special sensor microwave/imager (SSM/I, 1987(SSM/I, -2007, and special sensor microwave imager/sounder (SSMI/S, 2008-2015), which were provided by the National Snow and Ice Data Center (NSIDC). However, because these sensors are mounted on different platforms, there may be inconsistencies in the data, cross-calibration of the light temperature may improve the consistency. Compared to the SSM/I and SSMI/SSD datasets, this dataset has higher accuracy and lower biases [87], which are adjusted and validated using the meteorological stations, regarding frost-cover, the liquid water content in snow layer, vegetation, and surface water bodies. The accuracy of the LSDPC has been greatly improved by cross-calibrating the brightness temperatures of the different sensors; it has also been improved through applying a modified Chang algorithm that uses the decision-tree method [87,88]. By comparing the MODIS SCA in winter, the mean accuracy of this passive microwave remote-sensing data is 86% [89], with relative biases decreasing from 66.18% to −1.5% [87]. This approach has been widely used in the Tibetan Plateau, the Tienshan Mountains, and the plains and mountain regions of China [59,75,85,90].
To further detect the SCA and snowline changes in a higher resolution, the MODIS/Terra Snow Cover eight-day (MOD10A2) product was utilized for 2002-2015, as it is generated by reading eight days of 500 m resolution MOD10A1 tiles. If snow is observed in a cell on any given day in the eight-day period, the cell is mapped as snow. In contrast, if no snow is found, the cell is filled with the clear-view observation that occurred most often (e.g., snow-free land, lake, etc.), which can greatly reduce the impact of clouds. It has an absolute accuracy of~93% [91] and has been widely used to detect SCA in topographically complex Tienshan Mountains and Karakoram ranges [27,60,78]. The study area is composed of two MOD10A2 tiles (h23v04, h24v04). We mosaicked these two tiles and reprojected them into an Albers equal-area conic projection using the MODIS Reprojection Tool (MRT). Finally, the mosaicked images were saved as a Geotiff file format, covering the 14-year period from 2002 to 2015. To reduce the influence from the larger cloud cover, we eliminated the records of these mosaicked images when the cloud area of the Aksu River catchment exceeded 20%. The eliminated data were replaced by adjacent images in space (note there were only a few relevant images, but these may have been highly influential). The values of the eliminated and nonexistent images were calculated from the average values of the adjacent images.

ASTER DEM
The Advanced Spaceborne Thermal Emission and Reflection Radiometer-Digital Elevation Model (ASTER DEMV2) was selected to carry out our investigation, which provides a resolution of 30 m. This dataset used in our study was obtained from China's Geospatial Data Cloud Site.
All of the datasets used in this study are summarized in Table 2.

Trend Analysis
In this study, a Mann-Kendall trend test was adopted to identify the trends of temperature, precipitation, SCA, SD, and streamflow [92,93]. This is a non-parametric test and considered to be more objective for the trends of time series, a method that has been widely used in the trend analysis of hydrological and climatic variables [85,94]. Additionally, linear regression was applied to analyze changing rates in temperature, precipitation, SCA, SD, and streamflow during each period in 1979-2015.

Summer Altitude of 0 • C Isotherm Calculations
More than 60% of the runoff in the Aksu River is determined by summer discharge. The percentage is higher in the Kumalak River (about 70%), which mainly depends on a large proportion of glacier meltwater. Temperature is identified as the dominant driver for runoff changes in the Kumalak River [26], so it is necessary to detect temperature changes that include both the land surface and the upper air.
To explore the runoff changes in the Kumalak River that respond to upper air temperatures, it is important to detect the summer 0 • C isotherm changes in the upper air temperatures of the catchment. Summer runoff originates from the mountain ranges and is strongly associated with the upper air temperature. Chen et al. [95] reported that similar change trends were found between the summer runoff and 0 • C level height in Northwest China. Li et al. [12] also found that the summer runoff of the Hotan River from the Kunlun Mountains showed a strong association with summertime upper air temperature. To accurately explore the effects of upper-air temperature on river runoff, a typical long-term radiosonde station (Kuqa) located in the Southern Tienshan Mountains was chosen in our study. In the troposphere, when the height rises with pressure drops, temperature generally presents a positive trend. Thus, in our study, the pressures, heights, and temperatures of the Kuqa radiosonde station data from 1979 to 2015 were used to calculate the altitude's 0 • C isotherm. We discovered that the summer 0 • C layer in the station is basically between the 500 hPa and the 850 hPa standard pressure level. Note that there are several 0 • C isotherms, but we chose the lowest one for our purposes.
In carrying out our investigations, we first determined, the nearest upper and lower standard pressure level at the position of the 0 • C layer at 0 and 12 o'clock noon each day. Next, the twelve-monthly altitudes of the 0 • C isotherm at 0 o'clock (midnight) and 12 o'clock (noon) during 1979-2015 were calculated by using the linear interpolation method, assuming that the temperature changes uniformly in the vertical direction of the two nearest standard pressure levels. The average monthly values from June-August were then calculated. Finally, the altitude of the 0 • C isotherm in the summertime was formulated as follows: where T is the temperature ( • C), H is the altitude (m), j is the 0 • C layer, and q and p are the lower and upper standard pressure levels of the 0 • C isotherm, respectively.

Glacier Mass Balance Calculations
Glacier mass balance (GMB) is considered a key indicator of climate change and is directly related to changes in hydrological processes. We want to estimate the runoff changes based on GMB, especially for rivers dominated by glacier meltwater. However, the shortage of monitoring GMB data in mountain regions, particularly for long time series, made this impractical without resorting to model building. Previous work found that runoff and precipitation in the mountains have close negative-exponential relations to their areas, with maximum precipitation and runoff occurring in glacierized areas [86]. According to the observed runoff and precipitation datasets, GMB in this study during 1979-2015 was reconstructed by using statistical mechanics and a maximum entropy principle (SMMEP) model, which has already been applied in the Tienshan Mountains [13]. Specifically, this model's reliability and accuracy have been validated in research on the Urumqi River [86].

Snow Cover Mapping Calculations
Snow is considered a regionally important buffer in the hydrological system against drought and could mitigate the speeds of a warming climate. Hence, the response of SCA to climate change has a greater influence on the regional hydrological process. The length of snow cover area duration (SCAD) is an indicator of snow phenology and has either a greater or lesser effect on the regulation of the hydrological system [85]. In the present study, normal SCAD is recorded as the number of days Remote Sens. 2020, 12, 2704 9 of 28 annually (1st January to 31st December) when SD exceeded 1.0 cm [94] from 1979-2015. This can be expressed as shown in Equations (2) and (3): where SCAD j is the snow cover area duration in year j (d); m is the number of days beginning with January 1st and ending with December 31st for any given year within the study period; S i is the pixel value, such that if S i < 1, we give it a new value of 0 (indicating no snow cover), whereas if S i > 1, we give it a new value of 1 (indicating snow cover); SCA i is the snow cover area (%) in day i of the study area; TNP i is the total number of these pixels (S i > 1) in day i; N is the total pixel number of the study area.

Sensitivity and Correlation Analysis between Snow Cover Area, Snow Depth, and Runoff
In this study, to accurately present the relationship between temperature, precipitation, SCA, SD, and runoff in the study area, the Pearson correlation test was used to detect the degree of correlation between them. Both SCA and SD have a direct effect on these rivers fed by snow or glacier meltwater in mountain regions. As the fastest-warming occurs in wintertime and springtime in the Tienshan Mountains, both SCA and SD have decreased at the same time as snowmelt runoff in spring increased. A positive precipitation increase in summer was observed to accompany SCA and SD increases. However, we need to find a way to quantify the response of runoff to SCA and SD for our research. To quantify the influence of SCA and SD changes on river runoff in the two catchments under study, a sensitivity model was employed. This sensitivity analysis theory, as recommended by Zheng et al. [96], was used to evaluate the sensitivity of monthly runoff to SCA/SD in the headwaters of the Yellow River Basins. Kour et al. [97] also used the model to detect the monthly river discharge to SCA in the Chenab Basin, Western Himalayas. The equation is given as: where C i is the SCA/SD, ε is the sensitivity coefficient, Q i is the monthly runoff, and C and Q are average SCA/SD and runoff elements during 2002-2015, respectively. The ε indicates that SCA/SD elements that have changed by 1% can result in an ε% change in river runoff.

Changes in River Runoff
From 1979 to 2015, the runoffs of the Kumalak and Toxkan Rivers both have increased substantially, with increasing rates being 0.56 × 10 8 m 3 /decade and 1.83 × 10 8 m 3 /decade (Figure 2a), respectively. For the first period , the runoffs both showed a significant upward trend (p < 0.01), particularly for the Toxkan River's runoff in spring, summer, and autumn, which accounted for 16.52%, 53.15%, and 21.32% of the total annual runoff increase (Table 3), respectively. Comparatively, the runoff of the Kumalak River in autumn and summer showed an obvious increase, accounting for 18.80% and 79.13% of the total annual runoff increase. However, since 2002, the runoffs of both rivers showed a significant negative trend (p < 0.05). In the Toxkan River, spring and summer runoffs decreased by 58.57% and 39.21%, respectively, compared to the annual runoff decrease during 2002-2015, while in the Kumalak River, the summer runoff showed an even more obvious decline, accounting for 81.96% of the annual runoff decrease (Figure 2c and Table 3).
Mean monthly runoff in the Toxkan and Kumalak rives during the periods 1979-2002 and 2002-2015 are shown in Figure 2f. As can be seen, the runoffs from these rivers began to increase in April and peaked in July/August. Runoff variations in the Toxkan River have a maximum contribution (19.41%) of the total average annual runoff in July, followed by August (18.90%) and June (16.04%), while for the Kumalak River, the maximum contribution occurs in August (27.74%), followed by July (26.68%) and June (14.12%). By contrasting the variations in monthly runoffs in the two periods, we observe a notable average runoff increase in April (76.94%) for the Toxkan River during 2002-2015, which accounts for about 49.14% of the total average annual runoff increase, while no obvious change is observed in the Kumalak River. Furthermore, we observe a runoff decrease in both rivers for July and August, especially for the Kumalak River, whose runoffs decreased by 8.64% and 8.19%, respectively. Runoff decreases of 9.24% and 5.22% are observed in the Toxkan River.  Figure 2f. As can be seen, the runoffs from these rivers began to increase in April and peaked in July/August. Runoff variations in the Toxkan River have a maximum contribution (19.41%) of the total average annual runoff in July, followed by August (18.90%) and June (16.04%), while for the Kumalak River, the maximum contribution occurs in August (27.74%), followed by July  Note: FSC/AC means a fraction of seasonal variables (runoff, temperature, precipitation, SCA, and SD) change relative to the total annual variables. The annual negative or positive symbols show the negative or positive trends in annual variables, while the seasonal negative or positive symbols indicate that they will inhibit or promote the annual trends. The physically meaningful percentages indicate the relative contribution of the seasonal to the annual trend in variables.
Since 1979, seasonal runoff contributed 20.51% (spring), 57.11% (summer), 16.98% (autumn), and 5.40% (winter) to the annual runoff in the Toxkan River, whereas in the Kumalak River, the contributions were 8.91%, 68.87%, 17.26%, and 4.96%, respectively. In considering the seasonal runoff changes during these two periods, we can see that changes in spring and summer runoff have the greatest influence on annual runoff in the Toxkan River, while variations in summer runoff have the greatest influence on the Kumalak. Therefore, in our study, we propose the hypothesis that the potential causes for annual runoff changes in the Toxkan and Kumalak Rivers are based on the characteristics of these two contrasting rivers. This study hypothesizes that the spring runoff decrease in the Toxkan River since 2002 is the result of precipitation decrease or slower snowmelt, while the negative trends in summer runoff of the Toxkan and Kumalek Rivers are caused by positive precipitation or slower glacier melt.

Relationship between Temperature, Precipitation, and Runoff
Obvious warming in the Aksu River catchment was observed at a rate of 0.021 • C/year during 1979-2015 (p < 0.05). We found a greater temperature increase (0.026 • C/year, p < 0.01) in the Toxkan River catchment, whereas a slower temperature increase was observed in the Kumalak River catchment (0.013 • C/year) and a negative temperature trend in the northeastern parts ( Figure S1a). These rises were especially notable in spring, at rates of 0.055 • C/year and 0.042 • C/year, respectively (p < 0.01 and p < 0.05). Summer temperatures showed only a slight increase at a rate of 0.013 • C/year for the Toxkan River catchment, while no summer temperature increase was detected in the Kumalak River catchment (Figure 3b). Annual precipitation rates for both the Toxkan and Kumalak Rivers showed a remarkable increase of 2.02 mm/year and 2.04 mm/year, respectively, while an astounding rate of increase of 21.72 mm/decade was found in the Tuergate station at 3504 m a.s.l. An analysis of annual precipitation from 1979 to 2015 revealed that precipitation in the catchments mainly occurred in summer and contributed about 50% of total annual precipitation. During the past 40 years, summertime precipitation increased by about 1.12 mm/year and 0.73 mm/year, which accounted for 55% and 36% of the total annual precipitation increase, respectively.
By contrasting temperature and precipitation changes in the Toxkan and Kumalak River catchments Based on our hypothesis, in this part, we hypothesized that spring and summer runoff in the Toxkan River and summer runoff in the Kumalak River are associated with precipitation decreases. A close relation was found between spring runoff and precipitation in the Toxkan River (R = 0.51, p < 0.01), but no obvious correlation was noted in the Kumalak River (R = 0.12). Therefore, the obvious decrease in the runoff in the Toxkan in spring since 2002 was related to the significant decrease in springtime precipitation (a sharp decrease rate of 16.55 mm/decade, which was about 3.23 times stronger than the total annual precipitation increase). For summer runoff, a positive correlation (R = 0.28) was found between runoff and precipitation in the Toxkan River, while in contrast, a significant negative correlation (R = −0.34, p < 0.05) was found between the same elements in the Kumalak River. The decrease in runoff despite positive trends in precipitation during summer confirms that precipitation changes are not the cause for summer runoff decreases in either the Toxkan or Kumalak Rivers.
A close relation (R = 0.77, p < 0.01) between the standard summer runoff of the Kumalak River and the altitude of the 0 • C layer height isotherm was found for the 1979-2015 study period ( Figure 5 (Figure 3f), which also showed a strong correlation (R = 0.64, p < 0.01) with summer runoff. Therefore, both the summer surface temperature and upper-air temperature and runoff have consistent trends, which indicates that the changes in summer runoff are inhibited by decreases in summertime temperatures and accelerated by increases.
River catchment, whereas a slower temperature increase was observed in the Kumalak River catchment (0.013 °C/year) and a negative temperature trend in the northeastern parts ( Figure S1a). These rises were especially notable in spring, at rates of 0.055 °C/year and 0.042 °C/year, respectively (p < 0.01 and p < 0.05). Summer temperatures showed only a slight increase at a rate of 0.013 °C/year for the Toxkan River catchment, while no summer temperature increase was detected in the Kumalak River catchment (Figure 3b). Annual precipitation rates for both the Toxkan and Kumalak Rivers showed a remarkable increase of 2.02 mm/year and 2.04 mm/year, respectively, while an astounding rate of increase of 21.72 mm/decade was found in the Tuergate station at 3504 m a.s.l. An analysis of annual precipitation from 1979 to 2015 revealed that precipitation in the catchments mainly occurred in summer and contributed about 50% of total annual precipitation. During the past 40 years, summertime precipitation increased by about 1.12 mm/year and 0.73 mm/year, which accounted for 55% and 36% of the total annual precipitation increase, respectively.  m/year and −0.13/year, which has a strong correlation (R = 0.73, p < 0.01). Meanwhile, an 0.1 °C average summer surface temperature decrease was detected in the Kumalak River catchment during 2002-2015 compared to 1979-2002 (Figure 3f), which also showed a strong correlation (R = 0.64, p < 0.01) with summer runoff. Therefore, both the summer surface temperature and upper-air temperature and runoff have consistent trends, which indicates that the changes in summer runoff are inhibited by decreases in summertime temperatures and accelerated by increases.

Relationship between Glacier Mass Balance and Runoff
Using the SMMEP model to calculate the GMB (Figure 6), we found that the mean annual GMB and cumulative GMB in the Toxkan catchment were about −0.32 and −11.96 m from 1979 to 2015, respectively. Larger negative GMB was detected in the Kumalak catchment, with mean annual and cumulative GMB reaching about −0.50 and −18.36 m, respectively, which was consistent with Duethmann et al.'s [26] findings. Compared to different periods, larger negative GMB was found in the Toxkan and Kumalak catchments during 1979-2002, with mean annual GMB of −0.36 and −0.57

Relationship between Glacier Mass Balance and Runoff
Using the SMMEP model to calculate the GMB (

Relationship between Snow Cover Area, Snow Depth, and Runoff
Approximately 64.66% of the Aksu River's total catchment area is covered by permanent or seasonal snow. The average annual SCA is about 57.13% in the Toxkan River catchment and 75.13% of the Kumalak River catchment. Spatially, the SCA is mainly distributed in the Central and Eastern parts of the Toxkan River catchment and the Eastern and Northwestern part of the Kumalak River catchment. The winter SCA in the Toxkan River catchment was about 90.17%, but decreased to 15.71% in summer, while the winter SCA of the Kumalak River catchment reached above 95% even in the warmest summer months, with the permanent snow in the catchment maintaining an average SCA above 35%.
Compared to different periods, SCA and SD of the Toxkan River catchment during 1979-2015 both showed a negative trend, especially from 2002 onward (Figure 7a,c), while an obvious positive trend in SCA and SD in the Kumalak River catchment was exhibited from 2002 to 2015, at rates of 5.23%/decade and 2.48 mm/decade, respectively (Figure 7b,d).
Because the Toxkan River's annual runoff is mainly influenced by snow meltwater in spring, we propose a hypothesis that the spring runoff changes were mainly caused by the snow meltwater. We found that the winter/spring SCA and SD strongly affect the spring runoff (as witnessed by significance test levels of 0.01 and 0.05, respectively) in the Toxkan River (Table 4). For example, Pearson's correlation coefficient values showed no obvious negative correlation between annual GMB and summer runoff in the Toxkan River. For the Kumalak River, the annual GMB showed a significantly negative correlation with annual runoff (R = −0.52, p < 0.01) and a strong negative correlation with summer runoff (R = −0.56, p < 0.01). The close relationship between GMB and summer runoff and significant (p < 0.05) negative GMB increase in the Kumalak River since 2002 have confirmed our hypothesis, which suggests that changes in GMB have directly affected the runoff dynamics in the summer melt season.

Relationship between Snow Cover Area, Snow Depth, and Runoff
Approximately 64.66% of the Aksu River's total catchment area is covered by permanent or seasonal snow. The average annual SCA is about 57.13% in the Toxkan River catchment and 75.13% of the Kumalak River catchment. Spatially, the SCA is mainly distributed in the Central and Eastern parts of the Toxkan River catchment and the Eastern and Northwestern part of the Kumalak River catchment. The winter SCA in the Toxkan River catchment was about 90.17%, but decreased to 15.71% in summer, while the winter SCA of the Kumalak River catchment reached above 95% even in the warmest summer months, with the permanent snow in the catchment maintaining an average SCA above 35%.
Compared to different periods, SCA and SD of the Toxkan River catchment during 1979-2015 both showed a negative trend, especially from 2002 onward (Figure 7a,c), while an obvious positive trend in SCA and SD in the Kumalak River catchment was exhibited from 2002 to 2015, at rates of 5.23%/decade and 2.48 mm/decade, respectively (Figure 7b,d).  To clarify the influence of SCA and SD changes on the Toxkan River and Kumalak River runoff, we further calculated the sensitivity coefficients of monthly runoff to SCA and SD (Figure 8). We found substantial differences between the two rivers, with the Toxkan River presenting a much higher sensitivity coefficient. Specifically, the sensitivity coefficients of the monthly runoff to SCA and SD in the Toxkan are −0.13~2.49 and −0.12~0.95, respectively. The sensitivity coefficient in April indicates that if SCA or SD increase by 1%, then the runoff will increase by 2.49% or 0.95%.
Meanwhile, negative sensitivity coefficients were observed in the Kumalak River with the εSCA(SCA, Q) and εSD(SD, Q) being about -[−1.84~−0.16] and [−0.3~−0.14%], respectively. The sensitivity coefficients were for March to September, which means that there is an increase in SCA or SD by 1%, then the runoff of the Kumalak River will decrease by 0.16~1.84% or 0.14~0.3%. The Because the Toxkan River's annual runoff is mainly influenced by snow meltwater in spring, we propose a hypothesis that the spring runoff changes were mainly caused by the snow meltwater. We found that the winter/spring SCA and SD strongly affect the spring runoff (as witnessed by significance test levels of 0.01 and 0.05, respectively) in the Toxkan River (Table 4). For example, during 1979-2002, the increasing rates of winter and spring SCA in the Toxkan River catchment were 6.51 and 2.91 times higher than the annual increase, respectively, and the increasing rate of winter SD was 2.13 times higher than the annual rate (Table 3). During this period, SCA/SD changes agreed well with the substantial increase in spring runoff attributed to the increasing meltwater from winter/spring SCA and SD. However, since 2002, both SCA and SD in winter/spring showed negative trends that are 1.75 and 1.58 times higher than the rates of annual decrease. As a result, a large runoff decrease was observed in spring, which accounts for 58.57% of the total annual runoff decrease (Table 3). In contrast, a significant inverse relevance (with significance levels of 0.05 and 0.01, respectively) is indicated between the runoff and SCA/SD in the summertime for the Kumalak River (Table 4). From 2002 onwards, the summer SCA and SD of the Kumalak River catchment showed a positive trend (Table 3), while the Kumalak River's runoff exhibited a notable decrease (Figure 2c,f). To clarify the influence of SCA and SD changes on the Toxkan River and Kumalak River runoff, we further calculated the sensitivity coefficients of monthly runoff to SCA and SD (Figure 8). We found substantial differences between the two rivers, with the Toxkan River presenting a much higher sensitivity coefficient. Specifically, the sensitivity coefficients of the monthly runoff to SCA and SD in the Toxkan are −0.13~2.49 and −0.12~0.95, respectively. The sensitivity coefficient in April indicates that if SCA or SD increase by 1%, then the runoff will increase by 2.49% or 0.95%.

Mechanism Analysis
The cause mechanisms differ with regard to runoff changes in the Toxkan and Kumalak Rivers during 1979-2002 and 2002-2015. Thus, a systematic analysis of runoff changes in these water bodies was conducted based on the multiple driving factors shown in Figure 9. During 1979-2002, as mentioned in the conclusions from the hypotheses above, the apparent runoff increase for the Toxkan River was directly caused by positive trends in precipitation and snow meltwater, i.e., the largest precipitation increase was in summer, autumn and spring, accompanied by the greatest runoff increase, which was attributed to 91% of the annual runoff increase. For increases in snow meltwater, we found a positive trend in winter/spring SCA and SD with a positive trend in temperature. This increase resulted in a greater volume of snow meltwater in spring, while WSCT displayed an advancing trend of −0.29 d/year ( Figure S2). Furthermore, since 2002, we detected a negative trend in annual runoff that was mainly caused by the precipitation decrease in spring (p < 0.05). At the same time, obvious negative trends in SCA and SD in winter/spring contributed to the negative trend in snow meltwater, despite the advancing trend in WSCT. The sensitivity coefficients were for March to September, which means that there is an increase in SCA or SD by 1%, then the runoff of the Kumalak River will decrease by 0.16~1.84% or 0.14~0.3%. The sensitivity is also used to reflect the different sources of meltwater feeding the river runoff. For example, the Toxkan River's runoff, which is made up of snow meltwater, shows higher sensitivity than the Kumalak River, which showed a negative response to the increase in SCA and SD.
According to our hypothesis, we conclude that the notable decrease in winter/spring SCA and SD since 2002 have greatly influenced the decrease in spring runoff of the Toxkan River. Furthermore, the notable increase in summer SCA and SD provide evidence that the expansion of the cryosphere with negative GMB has resulted in a decrease in summer runoff in the Kumalak River.

Mechanism Analysis
The cause mechanisms differ with regard to runoff changes in the Toxkan and Kumalak Rivers during 1979-2002 and 2002-2015. Thus, a systematic analysis of runoff changes in these water bodies was conducted based on the multiple driving factors shown in Figure 9. During 1979-2002, as mentioned in the conclusions from the hypotheses above, the apparent runoff increase for the Toxkan River was directly caused by positive trends in precipitation and snow meltwater, i.e., the largest precipitation increase was in summer, autumn and spring, accompanied by the greatest runoff increase, which was attributed to 91% of the annual runoff increase. For increases in snow meltwater, we found a positive trend in winter/spring SCA and SD with a positive trend in temperature. This increase resulted in a greater volume of snow meltwater in spring, while WSCT displayed an advancing trend of −0.29 d/year ( Figure S2). Furthermore, since 2002, we detected a negative trend in annual runoff that was mainly caused by the precipitation decrease in spring (p < 0.05). At the same time, obvious negative trends in SCA and SD in winter/spring contributed to the negative trend in snow meltwater, despite the advancing trend in WSCT. For the Kumalak River, a negative trend in summer precipitation was observed during 1979-2002. The summer runoff showed a significant increase (p < 0.05) under rising temperatures, contributing to 79.13% of the annual runoff increase. Since 2002, however, a negative trend was observed in summer runoff, which contributed to 81.96% of the annual runoff decrease. As a result, the positive trend in summer precipitation did not inhibit this trend of annual runoff decrease. For the Kumalak River, a negative trend in summer precipitation was observed during 1979-2002. The summer runoff showed a significant increase (p < 0.05) under rising temperatures, contributing to 79.13% of the annual runoff increase. Since 2002, however, a negative trend was observed in summer runoff, which contributed to 81.96% of the annual runoff decrease. As a result, the positive trend in summer precipitation did not inhibit this trend of annual runoff decrease. Meanwhile, from 1979-2002, obvious positive trends in summer SCA and SD were found, accounting for 56.87% and 90.05% of the annual SCA and SD increase. An average temperature decrease was also found compared to the period 1979-2002. Therefore, we considered that the significant decrease in the Kumalak River's summer runoff since 2002 directly resulted from the decrease in meltwater from glacier and snow in summer, owing to the decline in summer temperatures and substantial precipitation increase. Both of these factors have been confirmed by the negative trend in GMB, with positive trends in SCA and SD.
Furthermore, we analyzed the SCA and snow-line variations by using MOD10A2 images from 2002 onwards. These findings were consistent with our results, in that SCA showed a negative trend in the Toxkan River catchment, with a decreasing rate of 78.11 km 2 /year since 2002. In addition, the snow-line in the Toxkan River catchment showed a rising trend, at a rate of 4.16 m/year. This trend was especially obvious in summer, at an increasing rate of 2.74 m/year. In contrast, the snowline of the Kumalak River catchment showed a continuous decline at a rate of −9.08 m/year, though slightly higher in summer (about −10.67 m/year). The continuous decrease in the Kumalak River catchment's snowline is partly owed to the decreasing summertime temperatures and increasing precipitation.
To further explore the reasons for the reductions in runoff in the Kumalak River since 2002, we analyzed SCA changes for each elevation by dividing the study region into twelve zones at every 500 m (Table 5). We found that the SCA above 4000 m in the Toxkan River catchment exhibited a positive trend in summer (month with the highest temperature), which was mirrored by the Kumalak River catchment's SCA above 3500 m. Moreover, as the glacier termini of the Kumalak River catchment are mainly concentrated between 3900 and 4300 m, the expanding summer SCA in these altitudes is likely caused by the increase in precipitation and a decrease in temperature. By analyzing the variations in summer precipitation and temperature at different altitudes during 2002-2015 ( Figure 10), we can see that precipitation at these altitudes (3500-4500 m) in the Kumalak River catchment showed an increasing rate with decreasing temperatures, which means that the cryosphere in these areas (above 3500 m) where glaciers developed showed an expansion trend. These changes have caused obvious increases in SCA and SD as well as a lowering of the snowline. Additionally, there has been a slowdown in glacier retreat, resulting in reduced runoff of glacier meltwater to the Kumalak River.
Note: ↑ and ↓ denotes the increase and decrease in the snow cover area.
In contrast, we found a positive trend in temperature but a negative trend in precipitation at 3500-4000 m (and 3900-4300 m) in the Toxkan River catchment during summertime, which also caused a negative trend in SCA (Table 5). For the high mountain regions, the increased precipitation Remote Sens. 2020, 12, 2704 20 of 28 with decreased temperatures expanded the cryosphere. This trend is especially notable in summer, which further clarifies why the Kumalak Rives' runoff has presented a downward trend in the summertime since 2002.

Discussion
Under global warming, positive trends in alpine river runoff have generally been observed [6,14,98]. However, not all rivers have experienced increasing trends in the runoff, for example, a declining trend in annual runoff was observed in the Hotan River (northern slope of the Kunlun Mountains) and the Hunza River (the Karakoram Range) [12,27]. Significant precipitation increases in alpine regions in summer is often accompanied by a cooling process [99], particularly in the higher mountains with greater precipitation. In these regions, when there is a significant rise in summer precipitation, the change from solid to liquid or water vapor will absorb more heat and thus decrease the temperature, leading to an increase in SCA and SD if the temperature remains below the melting point [100,101].
In China, glaciers are primarily cold glaciers that accumulate from summer precipitation. In our study of the Southern Tienshan Mountains, accumulation and ablation both mainly occur in summer [102], as the photograph in Figure 9c shows. An increase in summertime SCA and SD will then lead to a negative trend in glacier mass loss [103,104] and reduce runoff from glacial meltwater. Recently, several authors [7,105] have confirmed that SCA and SD are increasing in the mountain regions and that this increase may have positive influences on glacier mass balance. The decreasing runoff trend from increased SCA was also observed in the Chenab Basin of the Himalayas [97]. With the positive trend in precipitation, summertime temperatures in the Kumalak River basin have shown a downward trend at a rate of −0.58 K/decade during 1998-2009 [100]. Wang et al. [106] also found an apparent temperature decrease in the Kumalak River catchment since 2000. This may be the reason for the significant slowing of the Southern Inylchek Glacier's mass loss during 1997-2009 [107] and for the Gregoriev ice cap [108] in the study area, which ultimately resulted in a decrease in Kumalak River runoff since 2002.
It appears that basins with a higher fraction of glacierized area generally experience increasing runoff trends in a warming climate (e.g., Middle Tienshan Mountains), while basins with less or no glacierization showed large variations in runoff changes (e.g., East Tienshan Mountains) [5,60,109]. For example, runoff in the Karatal River basin in the North Tienshan Mountains, which has a glacier area of more than 5%, showed a generally positive trend, while the runoffs in basins with glacier areas of less than 2% exhibited a negative trend [66,109]. Runoffs in the Toudaogou and Erdaogou River basins in the East Tien Shan showed a considerable negative trend, decreasing by about 12.5% and 6.6%, respectively, since 1998 as a consequence of the sharp decrease in glacier area (less than 1%) [110]. Ragettli et al. [5] found runoff in the Juncal catchment in Chile is projected to sharply decrease because of the total ice melt already being beyond its tipping point. Elsewhere, the Tarfala catchment, which has a 30% glacierized area, saw a significant increasing trend in floods, despite having similar climate conditions as sub-arctic northern Sweden. There was also a significant decrease in summer runoff in the Abisko catchment, whose glacierized area covers no more than 1% [96]. In our study, the increasing rate of the Kumalak River runoff (more than 16% glacierized area) was notably higher than that of the Toxkan River (with a glacier area of less than 4% of its catchment area) during 1979-2002 (Table 1). The Kumalak River runoff is projected to increase until the 2050s, in large part due to the sizeable fraction of its catchment area being covered by glaciers [84].
A systematic analysis of runoff changes was adopted in this study, despite some uncertainties in the datasets. For example, regarding the consistency and limitations of datasets used in our research, none of the meteorological stations have location changes or missing data during the study period. As a GPCC product, the gauge-based gridded precipitation dataset is influenced by a number of stations over time. Climatological infilling with good data coverage occurs after 1951, whereas prior to 1901, the data are relatively poor. Compared to other GPCC near real-time products, GPCC V.2018 (V8) has a higher accuracy, which is why it has been recommended for hydrometeorological model verification and water cycle studies. In our research, GPCC precipitation is consistent with precipitation fluctuations at both annual and seasonal scales at the local stations, with R 2 approaching 0.72-1 and absolute percentage bias 1.45-9.09%. Temperatures in the EARA-Interim dataset were used here, giving the highest consistent trend in annual and monthly temperature; specifically, R 2 approached about 0.50-0.85 (higher for the high-altitude stations) and 0.98-0.99, respectively. No location change information about these two hydrological stations were found during the study period, although local changes in equipment may have affected the consistency of the run-off data sets.
Furthermore, hydrological facilities in the upstream of the hydrological station were also considered, which may have influenced the hydrological process. In our study, no reservoirs were situated upstream of the hydrological stations, but two hydropower stations were found upstream of the Xiehela hydrological-station on the Kumalak River. One of the hydropower stations ("Dashixia") was built starting in October 2019, which puts it beyond our research period. The other hydropower station ("Xiaoshixia") began construction in December 2009 and started operations in 2012. Although this infrastructure might change the natural runoff in the downstream, its influence between 2012 and 2015 is relatively short, and therefore it would have a greater influence on the daily runoff, but less influence on monthly or annual runoff. In other words, the significant consistent trend continued to be upper air temperatures and runoff in the summertime, which persists in having a significant correlation throughout these two (p < 0.01). For the Toxkan River, whose stations are mainly distributed downstream, only one hydropower station ("Biedilike") was found in the upper area, and it has only been operational since March 2012. Given the influence of complex terrain, elevation, underlying surface, alpine climate, and a lower spatial resolution, misclassification in extracting SCA might be caused by the relatively coarse spatial resolution. Even though the MOD10A2 product has a higher resolution of 500 m and therefore was adopted for analyzing the SCA and snow-line variations at different altitudes, this dataset was only available after March 2000, and only for the period from 2002 to 2015.
Compared to the observed GMB data, there are some uncertainties in the physical glacier mass model used here, e.g., altitudes of the meteorological and hydrological stations. While previous studies have further confirmed our results, the area loss in the Qingbingran No. 72 Glacier in the Tomor region was 0.02 km 2 /year during 2008-2009, which is much lower than during 1964-2008 (0.03 km 2 /year) [111].
Zhao et al. [84] also found that since 2000, the shrinkage rate of the glacier area in the Aksu River basin shows a lower rate of 0.67%/year compared to 1990-2000, when the reduction rate was 0.89%/year. Meanwhile, Wang et al. [112] discovered that Qingbingtan Glacier No. 74 had the highest shrinkage rate, showing a 21.5% decrease in glacier area from 1964 to 2009, while a slower glacier area loss of −3.7 ± 2.7% from 1990 to 2010 was found in the Kumalak catchment [100]. In addition, both glacier areas (Glacier No. 72  Interestingly, some glacier areas have even increased since 1999. Examples include the Samoilowitsch Glacier, whose area increased by more than 10% [85], and the Northern Inylchek Glacier, which increased by~6.5 km 2 [113]. Compared to the features of the loss area, GMB is considered a direct indicator of climate change and is related to changes in hydrological processes. Pieczonka et al. [114] reported that the Aksu River's runoff increased by about 20% during~1975-1999, mainly due to the glacier's negative mass loss of about −0.35 ± 0.34 m/year. A remarkable glacier mass loss (with a rate of −0.51 ± 0.36 m/year) was found in the outer ranges of the Ak-Shrink massif (northwestern part of the Kumalak River catchment). In addition, Pieczonka et al. [76] showed that the glacier mass loss in the Aksu River catchment slowed during 1999-2009, at a loss rate of −0. 23 Figure 11). Figure 11. Glacier mass budgets for (a) Aksu catchment [76], (b) Qingbingtan Glacier No. 74 [112], (c) Southern Inylchek Glacier [107], and (d) South of Tomur Peak [114].
Meanwhile, as shown in Figure 11c, the Southern Inylchek Glacier (largest glacier in the Tomor region) displayed a smaller retreat from 1999-2007 (0.28 ± 0.46 m/year) compared to ~1975-1999 (0.43 ± 0.10 m/year) [107]. As well, Farinotti et al. [114] found that the glacier area in the Tienshan Mountains had decreased by 18% and glacier mass by about 27% over the past 50 years, especially from 1970 to 1990. During 1999-2009, the Tomor region experienced only a moderate loss in mass. Therefore, we believe that the decrease in the Kumalak River runoff since 2002 is closely related to the slower decrease in glacier mass loss. The slowing of glacier mass loss in the Kumalak River basin in recent years appears to be a key factor in the considerable decrease in the runoff.
With a continuously changing climate, the hydrological process in mountain ranges will become more vulnerable to climate change. Thus, combining information from glaciers, snow, and meteorological and hydrological data by using the glacial-hydrological model to predict future water resources is our next work.
Meanwhile, as shown in Figure 11c, the Southern Inylchek Glacier (largest glacier in the Tomor region) displayed a smaller retreat from 1999-2007 (0.28 ± 0.46 m/year) compared to~1975-1999 (0.43 ± 0.10 m/year) [107]. As well, Farinotti et al. [114] found that the glacier area in the Tienshan Mountains had decreased by 18% and glacier mass by about 27% over the past 50 years, especially from 1970 to 1990. During 1999-2009, the Tomor region experienced only a moderate loss in mass. Therefore, we believe that the decrease in the Kumalak River runoff since 2002 is closely related to the slower decrease in glacier mass loss. The slowing of glacier mass loss in the Kumalak River basin in recent years appears to be a key factor in the considerable decrease in the runoff.
With a continuously changing climate, the hydrological process in mountain ranges will become more vulnerable to climate change. Thus, combining information from glaciers, snow, and meteorological and hydrological data by using the glacial-hydrological model to predict future water resources is our next work.

Conclusions
This study investigated the impact of climate, GMB, SCA, and SD changes on water resources of two contrasting rives-the Toxkan and the Kumalak-originating in the Tienshan Mountains, based on multiple data sets from 1979 to 2015. The main conclusions are given below.
Over the past few decades, these two rivers that are recharged by snow and glacier meltwater have experienced an increasing trend in runoff under a warming climate. However, the runoffs of the Toxkan and Kumalak Rivers revealed different trends. Although exhibiting a strongly positive trend at rates of 5.71 × 10 8 m 3 /decade and 6.63 × 10 8 m 3 /decade from 1979 to 2002, the upward trend was reversed to downward at rates of 10.17 × 10 8 m 3 /decade and 10.21 × 10 8 m 3 /decade, respectively, post-2002. In total, the changes in annual runoff of the Toxkan River are mostly influenced by the spring and summer runoff, whereas the changes in annual runoff of the Kumalak River are mainly controlled by the summer runoff.
For the Toxkan River, the increased spring and summer runoff during 1979-2002 was mainly caused by the positive trends in spring snow meltwater and summer precipitation, whereas the decreased spring runoff since 2002 was directly related to the significant decrease in spring precipitation, at a rate of −16.55 mm/decade (3.23 times higher than the rate of annual decrease). Meanwhile, the tremendous decrease in wintertime/springtime SCA and SD (with rates 1.76 and 1.58 times higher than the annual decrease, respectively) also contributed to the decrease in spring runoff. The reduction in summer runoff resulted from the decrease in snow meltwater, owning to the earlier snowmelt period during spring as well as to the sharp shrinkage in the glacier area in the catchment.
For the Kumalak River, increasing temperature with a related increase in glacier meltwater was the dominant driver for the summer runoff increase during 1979-2002. However, the substantial decrease in summer runoff since 2002 was related to the drop in summertime temperatures (both surface and upper-air temperatures). Summer temperature reductions have slowed the melting of glaciers, with mean annual GMB decreasing to −0.37 m; they have also increased summer SCA and SD, accounting for 56.87% and 90.05% of the annual increase, respectively. In addition, falling average temperatures while the positive trend in precipitation and SCA at higher elevations (3500-4500 m) where glaciers develop. We also found that the average snowline of summertime has decreased by about 300 m.