Snow-Cover Area and Runo ﬀ Variation under Climate Change in the West Kunlun Mountains

: In recent years, the climate in the arid region of Northwest China has become warmer and wetter; however, glaciers in the north slope of the West Kunlun Mountains (NSWKM) show no obvious recession, and river ﬂow is decreasing or stable. This contrasts with the prevalent response of glaciers to climate change, which is recession and initial increase in glacier discharge followed by decline as retreat continues. We comparatively analyzed multi-timescale variation in temperature–precipitation–snow cover-runo ﬀ in the Yarkant River Basin (YRK), Karakax River Basin (KRK), Yurungkax River Basin (YUK), and Keriya River Basin (KRY) in the NSWKM. The Mann–Kendall trend and the mutation–detection method were applied to data obtained from an observation station over the last 60 years (1957–2017) and MODIS snow data (2001–2016). NSWKM temperature and precipitation have continued to increase for nearly 60 years at a mean rate of 0.26 ◦ C / decade and 5.50 mm / decade, respectively, with the most obvious trend (R 2 > 0.82) attributed to the KRK and YUK. Regarding changes in the average snow-cover fraction (SCF): YUK (SCF = 44.14%) > YRK (SCF = 38.73%) > KRY (SCF = 33.42%) > KRK (SCF = 33.40%). Between them, the YRK and YUK had decreasing SCA values (slope < − 15.39), while the KRK and KRY had increasing SCA values (slope > 1.87). In seasonal variation, the SCF of the three of the basins reaches the maximum value in spring, with the most signiﬁcant performance in YUK (SCF = 26.4%), except for YRK where SCF in spring was lower than that in winter ( − 2.6%). The runo ﬀ depth of all river basins presented an increasing trend, with the greatest value appearing in the YRK (5.78 mm / decade), and the least value in the YUK (1.58 mm / decade). With the runo ﬀ response to climate change, temperature was the main inﬂuencing factor of annual and monthly (summer) runo ﬀ variations in the YRK, which is consistent with the runo ﬀ -generation rule of rivers in arid areas, which mainly rely on ice and snow melt for water supply. However, this rule was not consistent for the YUK and KRK, as it was disturbed by other factors (e.g., slope and slope direction) during runo ﬀ generation, resulting in disruptions of their relationship with runo ﬀ . This research promotes the study of the response of cold and arid alpine regions to global change and thus better serve regional water resources management.


Introduction
The water cycle of a basin is affected by regional climate, topography, and human activities [1][2][3][4]. Global warming has accelerated the regional water cycle and has changed the rainfall-runoff process in mountainous areas, resulting in changes in runoff recharge from glacier meltwater and precipitation [5]. For extremely arid areas with sparse precipitation, meltwater is the main recharging source of rivers [6][7][8]. Water from the accelerated melting of snow and ice caused by climate warming can contribute to runoff in a certain period of time, which leads to an increase in runoff [6,9]. However, continued warming also accelerates snow covered area and glacier area loss, leading to changes in runoff processes and water resources [10,11].
In the context of global warming [12], temperature in the Tibetan Plateau (TP) is rising at a rate of 0.28 • C/decade, which leads to permafrost degradation, widespread glacial retreat, intensified snow melting, frequent occurrence of glacial-lake disasters, and the associated deterioration of the ecological environment, so this area can be regarded as sensitive to global climate change [13,14]. As the main glacier distribution area of the TP [15], the North Slope of the West Kunlun Mountains (NSWKM) comprises the Yarkant River Basin (YRK), Karakax River Basin (KRK), Yurungkax River Basin (YUK), and Keriya River Basin (KRY). Between them, the rivers Yarkant, Karakax and Yurungkax are the major tributaries of the longest continental river in China, the Tarim River (TR), and discharge mainly comes from the meltwater of glaciers and seasonal snow in the NSWKM [6,16]. Recent studies have shown that the glacier area of the TR Basin decreased by 3.3% from 1960 to 2000, but different sub-basins had very different decreases in glacier area, ranging from 0.7% to 7.9% [17]. Theoretically, more glacier mass loss is typically associated with more increased stream discharge, which is mainly contributed by glacier meltwater, but the phenomenon that regional runoff changes are not proportional to glacier melt occurs in parts of the Himalayan region [18]. Despite the vital importance of meltwater for water-resource management in these basins, little is known about the hydrological roles of snowmelt runoff under the changing climate [19]. Therefore, it is necessary to accurately understand the controlling factors for flow variability in the past several decades for both practical and scientific needs.
Currently, many studies on spatial-and temporal-variation characteristics, such as the average snow-cover area (SCA) [20][21][22], snow depth [23][24][25], and snow-cover days [26,27], have been carried out in various basins. With the development of remote-sensing technology, snow-cover data have also been developed from site-scale to remote-sensing products on the basis of Landsat, SPOT, and MODIS. From them, MODIS snow products are the most widely used [28]. However, these studies mostly focus on the investigation of snow cover in wet areas and rarely include systematic analysis and quantification of the relationship between alpine snow cover, climate, and runoff in arid areas. In a study of the variation of snow cover and runoff in the West Kunlun Mountains, Yasuda and Furuya [29] found that the annual average of SCA in the YUK increased by 1.4 km 2 from 1970 to 1989. However, it decreased by 3.5 km 2 from 1989 to 2001. Gardelle et al. [30] found that some basins of the Kunlun Mountains had significantly different SCA values, and most of them had increasing SCA values. With the intensification of global warming, the climate in most parts of arid areas of Northwest China, especially in Xinjiang, is projected to change toward a warmer-wetter state [31][32][33], and runoff is projected to increase significantly in most rivers [31]. However, some basins (such as the YUK and KRK) of the NSWKM experience no significant change in runoff, or may even experience a decreasing trend [32,34]. The difference in snow cover and runoff change in the NSWKM has received much attention by researchers [30,[35][36][37][38][39]. In addition, the YRK, YUK, KRK, and KRY sit on the western frontier of China and border Pakistan and India. Historically, however, the allocation of water resources of the Upper Indus River Basin (UIB) is one of the problems leading to intense conflict between India and Pakistan [40,41]. Therefore, the study of the impact of climate change on runoff and snow cover in this region can aide in the regional planning and management of water resources.
The objectives of this study are to (1) improve our understanding of the variation characteristics of the NSWKM climate, (2) quantitatively assess runoff dynamics and snow-cover variation at different altitudes and seasons, and (3) analyze the rule of runoff generation and investigate the controlling factors for runoff changes in the NSWKM under climate change. This research will promote the study of the response of cold and arid alpine regions to global change and thus better serve regional water resources management.

Study Area
Located in the southwestern margin of the TR Basin, the northern slope (74 • 20 -83 • E, 34 • 30 -38 • 30 N) of the West Kunlun Mountains includes rivers Yarkant, Tiznafu, Pishan, Karakax, Yurungkax, Qira, and Keriya from west to east [42][43][44]. Between them, the Yarkant, Karakax, and Yurungkax are the three largest rivers in the NSWKM [45]. The Keriya is a typical glacial meltwater-recharged river; thus, its water volume directly restricts the economic and social development of the Keriya oasis [16]. In this study, the YRK, KRK, YUK, and KRY were selected as the main basins to be studied. The basic information on these basins is shown in Table 1.  In this study, the missing data, including temperature and precipitation data, from some stations (in several years) were mainly interpolated by the data from the nearby stations that had complete data. In turn, the missing runoff data from some hydrological stations in the dry season were mainly interpolated by the recession-curve method [49,50]; and the missing summer runoff data were mainly interpolated by runoff data from upstream and downstream hydrological stations. In addition, the arithmetic average of temperature and precipitation from meteorological and hydrological stations in the same period was used to reflect overall climate conditions in the NSWKM.

Observation Data
The basic meteorological elements of high-altitude mountainous areas, such as temperature and precipitation, were observed by three Automatic Weather Stations (AWSs) in the study area in 2010-2014 ( Figure 1, Table 1). In this study, the daily temperature observed by the AWS was used to test the correlation with daily temperature from adjacent hydrological and meteorological stations in lowaltitude mountainous areas ( Figure 2). Since the results showed a very significant positive temperature correlation with correlation coefficients between 0.802 and 0.966 at the significant level of 0.001, it can be stated that long-term temperature data of the hydrological and meteorological stations in low-altitude mountainous areas were representative of temperature changes in highaltitude mountainous areas. In this study, the missing data, including temperature and precipitation data, from some stations (in several years) were mainly interpolated by the data from the nearby stations that had complete data. In turn, the missing runoff data from some hydrological stations in the dry season were mainly interpolated by the recession-curve method [49,50]; and the missing summer runoff data were mainly interpolated by runoff data from upstream and downstream hydrological stations. In addition, the arithmetic average of temperature and precipitation from meteorological and hydrological stations in the same period was used to reflect overall climate conditions in the NSWKM.

Observation Data
The basic meteorological elements of high-altitude mountainous areas, such as temperature and precipitation, were observed by three Automatic Weather Stations (AWSs) in the study area in 2010-2014 ( Figure 1, Table 1). In this study, the daily temperature observed by the AWS was used to test the correlation with daily temperature from adjacent hydrological and meteorological stations in low-altitude mountainous areas ( Figure 2). Since the results showed a very significant positive temperature correlation with correlation coefficients between 0.802 and 0.966 at the significant level of 0.001, it can be stated that long-term temperature data of the hydrological and meteorological stations in low-altitude mountainous areas were representative of temperature changes in high-altitude mountainous areas. There are few hydrological and meteorological stations in the mountainous study area. However, precipitation is significantly affected by landforms [51]. Therefore, we tested the correlation of monthly precipitation between hydrological and meteorological stations at high and low altitudes in summer ( Figure 3). The results showed that, although monthly precipitation had a lower coefficient of correlation between all stations than temperature did (Figure 2), they all passed the 99.99% significance test, showing significant positive correlation, i.e., long-term precipitation data from relatively low-altitude hydrological and meteorological stations could be used as proxy records to study precipitation in the mountainous area of the Tibetan Plateau.

Snow-Cover Data
Due to 8-day data synthesis, 8-day MODIS composite snow-cover products partly eliminated cloud obscuration, and their accuracy was improved [52]. Therefore, in this work, MOD10A2 (acquired from Terra satellite in 2001-2016) and MYD10A2 (acquired from Aqua satellite in 2002-2016) products were chosen to analyze snow dynamics. As high-runoff and high-precipitation zones in the basin, alpine areas, especially in glacier areas, were much more likely to be disturbed by clouds than those in nonglacial areas despite the remarkable cloud-elimination effect [36,53]. Therefore, we further optimized the cloud-elimination effect by combining the satellite images of the MOD10A2 and MYD10A2 snow products and introducing glacier boundaries to improve the usability of the satellite images [54]. The specific operation process was as follows: (1) Images of the Terra and Aqua snow covers were projected and formatted by MODIS Reprojection Tool (MRT); (2) Images of Terra and Aqua in the same period were synthesized by extracting snow-cover images of the basin and recoding the original classification codes [55]; There are few hydrological and meteorological stations in the mountainous study area. However, precipitation is significantly affected by landforms [51]. Therefore, we tested the correlation of monthly precipitation between hydrological and meteorological stations at high and low altitudes in summer ( Figure 3). The results showed that, although monthly precipitation had a lower coefficient of correlation between all stations than temperature did (Figure 2), they all passed the 99.99% significance test, showing significant positive correlation, i.e., long-term precipitation data from relatively low-altitude hydrological and meteorological stations could be used as proxy records to study precipitation in the mountainous area of the Tibetan Plateau.

Snow-Cover Data
Due to 8-day data synthesis, 8-day MODIS composite snow-cover products partly eliminated cloud obscuration, and their accuracy was improved [52]. Therefore, in this work, MOD10A2 (acquired from Terra satellite in 2001-2016) and MYD10A2 (acquired from Aqua satellite in 2002-2016) products were chosen to analyze snow dynamics. As high-runoff and high-precipitation zones in the basin, alpine areas, especially in glacier areas, were much more likely to be disturbed by clouds than those in nonglacial areas despite the remarkable cloud-elimination effect [36,53]. Therefore, we further optimized the cloud-elimination effect by combining the satellite images of the MOD10A2 and MYD10A2 snow products and introducing glacier boundaries to improve the usability of the satellite images [54]. The specific operation process was as follows: (1) Images of the Terra and Aqua snow covers were projected and formatted by MODIS Reprojection Tool (MRT); (2) Images of Terra and Aqua in the same period were synthesized by extracting snow-cover images of the basin and recoding the original classification codes [55]; (3) On the basis of the watershed glacier boundary, the cloud pixels of the glacier areas were identified as snow pixels to obtain the final sequence of the snow image data ( Figure 4).
Water 2019, 11, x FOR PEER REVIEW 6 of 23 (3) On the basis of the watershed glacier boundary, the cloud pixels of the glacier areas were identified as snow pixels to obtain the final sequence of the snow image data ( Figure 4).   In this paper, the method proposed by Paudel and Andersen [20] was used to evaluate the effectiveness of the cloud-elimination method ( Table 2). After two-step cloud removal, 77.23% of cloud pixels were eliminated on average. Among them, the two satellite-image synthesis methods eliminated 48.61% of cloud pixels on average, underestimated 3.66% snow-cover pixels, and In this paper, the method proposed by Paudel and Andersen [20] was used to evaluate the effectiveness of the cloud-elimination method ( Table 2). After two-step cloud removal, 77.23% of cloud pixels were eliminated on average. Among them, the two satellite-image synthesis methods eliminated 48.61% of cloud pixels on average, underestimated 3.66% snow-cover pixels, and overestimated 0.89% snow-cover pixels. After the introduction of the glacier-snow-cover area, there was only slight overestimation, but no underestimation. Therefore, it can be seen that the method of removing cloud pixels by introducing a glacier-snow-cover area was effective. Statistical analysis of the changes in the content of cloud pixels was performed before and after the removal of 648 images. From the original images, 80 images contained more than 16% of cloud pixels, while 56 images contained 10% to 16% of cloud pixels. After synthesis of the two sets of satellite images, there were only 38 images with a cloud-pixel content of more than 16%, and 42 images with a cloud-pixel content of between 10% and 16%. When it was introduced into the glacier area, only 3 of them exceeded 16%, while 14 of them ranged from 10% to 16%. After 2 steps of cloud removal, only 2.6% of the total number of images contained more than 10% of the cloud pixels and most of images contained less than 2% of the cloud pixels. Thus, it was concluded that the method of eliminating cloud pixels in the SCA was very effective in this area [54].

Statistical Analysis
The Mann-Kendall method, a popular nonparametric test method [56,57], has been widely used to detect changes in meteorological and hydrological time series [47,58,59]. Given a statistical time series S k from i = 1, 2, . . . , k, for each variable (x i , x j ), if j < i, x i > x j , an upward trend is stated; r i is the total number of dual values with a growth trend.
The mean and variance of S k are computed as follows: Water 2019, 11, 2246 8 of 21 The sequential values of the UF k statistic are then calculated as When UF > 0, this indicates that the sequence has an upward trend. In contrast, UF < 0 indicates a decreasing trend. The change point is determined as the intersection point of UF and UB within the confidence interval (value of Zc = ±1.96). Thus, a significant increase or decrease can be identified at a 95% confidence level. UB is the backward trend sequence generated with the reverse data series of UF and can also be calculated by the function of Equation (4). The intersection of UF and UB indicates an abrupt change at the corresponding time [60].

Variation Characteristics of Climate on North Slope of West Kunlun Mountains
Analysis of the annual average temperature of the NSWKM from 1957 to 2017 ( Figure 5) showed that all stations had an increasing annual average temperature. Among them, the Hotan and Kaqun stations had the highest (0.41 • C/decade) and lowest warming rate (0.12 • C/decade), respectively. The NSWKM had a warming rate of 0.26 • C/decade (average for data from meteorological and hydrological stations). Except for Tashkurgan station, which had a lower average annual temperature than the NSWKM, the other stations had a similar average annual temperature to that of the NSWKM. Among them, Wuluwati station had the most significance (R 2 = 0.86), and therefore further reflected the characteristics of temperature changes in the entire region.
The NSWKM had an average annual precipitation of 77 mm, with a positive trend in time of 5.50 mm/decade (increased by 7.14% per decade). All stations showed increasing annual average precipitation ( Figure 6). Kulukelangan station had the most significant increasing average annual precipitation with an increasing rate of 10.5 mm/decade (about 9.15% per decade), while Tongguziluoke station had the most correlative annual average precipitation with the NSWKM (R 2 = 0.82) and was therefore considered as representative. NSWKM had a warming rate of 0.26 °C/decade (average for data from meteorological and hydrological stations). Except for Tashkurgan station, which had a lower average annual temperature than the NSWKM, the other stations had a similar average annual temperature to that of the NSWKM. Among them, Wuluwati station had the most significance (R 2 = 0.86), and therefore further reflected the characteristics of temperature changes in the entire region. The NSWKM had an average annual precipitation of 77 mm, with a positive trend in time of 5.50 mm/decade (increased by 7.14% per decade). All stations showed increasing annual average precipitation ( Figure 6). Kulukelangan station had the most significant increasing average annual precipitation with an increasing rate of 10.5 mm/decade (about 9.15% per decade), while To study SCA change regulation in the four basins at different altitudes, we divided their altitudes according to elevation difference, using 500 m as a reference, and drew the annual average snow-cover-fraction (SCF) change curves for different altitudes (Figure 8). Results showed that the YRK, KRK, YUK, and KRY had average SCAs of 17,977.8 km 2 (SCF = 38.60%), 6652.8 km 2 (SCF = 33.39%), 6514.0 km 2 (SCF = 43.95%), and 2799.3 km 2 (SCF = 33.29%), respectively. During seasonal change, the four basins had maximum SCF values in spring. For the whole basin (Figure 8, black curve), annual change curves were M-shaped, i.e., the largest SCF was in spring and autumn, and the smallest SCF was in summer and midwinter (December-January). The maximum SCF of the different basins occurred in different periods of time and in different forms. From these, the YRK had a maximum SCF mainly in mid-March and mid-late October, the KRK from mid-March to late May and mid-October, the YUK from mid-March to late May and from late September to late October, and the KRY in early May and late September. The four basins showed the same change in annual SCF at each altitude, i.e., a single peak below 4000 m, mainly in winter; a significant double peak between 4000 and 6000 m; and a longer peak duration at higher altitudes; additionally, higher glacier coverage above 6000 m resulted in little change in SCF. This study area showed bimodal change in annual SCF and a smaller SCA in winter than in spring and autumn, possibly because of lower precipitation in winter, smaller snow supply, lower temperature, and higher wind speed in mountainous areas in winter leading to smaller critical wind speeds of wind-blown snow, easier occurrence of wind-blown snow, and snow redistribution and sublimation [61]. Tongguziluoke station had the most correlative annual average precipitation with the NSWKM (R 2 = 0.82) and was therefore considered as representative.  Figure 1).

Variation Characteristics of Snow Cover on North Slope of West Kunlun Mountains
Comparative analysis of the variation characteristics of the average, maximum, and minimum values of SCA in four typical basins on the NSWKM from 2001 to 2016 (Figure 7) indicated that, except for the YUK, the three other basins showed insignificant increase in minimum SCA. The YRK and YUK showed insignificant decrease in annual average SCA, while the KRK and KRY showed  To study SCA change regulation in the four basins at different altitudes, we divided their altitudes according to elevation difference, using 500 m as a reference, and drew the annual average snow-cover-fraction (SCF) change curves for different altitudes (Figure 8). Results showed that the YRK, KRK, YUK, and KRY had average SCAs of 17,977.8 km 2 (SCF = 38.60%), 6652.8 km 2 (SCF = 33.39%), 6514.0 km 2 (SCF = 43.95%), and 2799.3 km 2 (SCF = 33.29%), respectively. During seasonal change, the four basins had maximum SCF values in spring. For the whole basin ( Figure 8, black curve), annual change curves were M-shaped, i.e., the largest SCF was in spring and autumn, and the smallest SCF was in summer and midwinter (December-January). The maximum SCF of the different basins occurred in different periods of time and in different forms. From these, the YRK had a maximum SCF mainly in mid-March and mid-late October, the KRK from mid-March to late May and mid-October, the YUK from mid-March to late May and from late September to late October, and the KRY in early May and late September. The four basins showed the same change in annual SCF at each altitude, i.e., a single peak below 4000 m, mainly in winter; a significant double peak between 4000 and 6000 m; and a longer peak duration at higher altitudes; additionally, higher glacier coverage above 6000 m resulted in little change in SCF. This study area showed bimodal change in annual SCF and a smaller SCA in winter than in spring and autumn, possibly because of lower precipitation in winter, smaller snow supply, lower temperature, and higher wind speed in mountainous areas in winter leading to smaller critical wind speeds of wind-blown snow, easier occurrence of wind-blown snow, and snow redistribution and sublimation [61].

Runoff Variation Characteristics on North Slope of West Kunlun Mountains
From 1957 to 2017, the four basins of the NSWKM had increasing runoff depth (Figure 9). The YRK had the most significant increasing runoff depth (5.78 mm/decade), while the YUK had the most insignificant (1.58 mm/decade). The four basins had a minimum monthly runoff depth from December to March and, apart from the YRK in August, a maximum monthly runoff depth in July ( Figure 10).
Statistical analysis of the probability of minimum and maximum monthly mean runoff occurrence in the four basins ( Table 3) also showed that the YRK was very different in terms of its maximum monthly runoff depth compared to the three other basins, possibly because glaciers blocked the lakes upstream of the YRK, and flash floods mainly appeared in August [54].

Runoff Variation Characteristics on North Slope of West Kunlun Mountains
From 1957 to 2017, the four basins of the NSWKM had increasing runoff depth (Figure 9). The YRK had the most significant increasing runoff depth (5.78 mm/decade), while the YUK had the most insignificant (1.58 mm/decade). The four basins had a minimum monthly runoff depth from December to March and, apart from the YRK in August, a maximum monthly runoff depth in July ( Figure 10).

Dynamics of Abrupt Change in Runoff Depth
The Mann-Kendall (M-K) trend test of annual runoff depth showed that the Zc values of the annual average runoff depth of the YRK and KRY passed the 0.05 significance level tests and significantly increased ( Figure 11) with increasing rates of 5.78 and 4.27 mm/decade, respectively. The annual runoff depth of the KRK and YUK increased insignificantly, with increasing rates of 1.71 and 1.58 mm/decade, respectively. This result was qualitatively consistent with the M-K trend test results obtained by Xu et al. [8], who analyzed the annual runoff changes of the YRK, KRK, and YUK from 1957 to 2008, i.e., showing that the annual runoff of the YRK significantly increased, while that of the KRK and YUK did not. The M-K trend test of monthly runoff depth showed that the monthly runoff depth of the YRK increased; this increase was significant from January to June, and from September to December, but it was insignificant from July to August. The result is due to runoff timing not changing in the main glacier melt months of July and August. Moreover, the runoff depth of the KRY in June and July did not pass the 95% significance tests and, in the other months, Statistical analysis of the probability of minimum and maximum monthly mean runoff occurrence in the four basins ( Table 3) also showed that the YRK was very different in terms of its maximum monthly runoff depth compared to the three other basins, possibly because glaciers blocked the lakes upstream of the YRK, and flash floods mainly appeared in August [54].

Dynamics of Abrupt Change in Runoff Depth
The Mann-Kendall (M-K) trend test of annual runoff depth showed that the Z c values of the annual average runoff depth of the YRK and KRY passed the 0.05 significance level tests and significantly increased ( Figure 11) with increasing rates of 5.78 and 4.27 mm/decade, respectively. The annual runoff depth of the KRK and YUK increased insignificantly, with increasing rates of 1.71 and 1.58 mm/decade, respectively. This result was qualitatively consistent with the M-K trend test results obtained by Xu et al. [8], who analyzed the annual runoff changes of the YRK, KRK, and YUK from 1957 to 2008, i.e., showing that the annual runoff of the YRK significantly increased, while that of the KRK and YUK did not. The M-K trend test of monthly runoff depth showed that the monthly runoff depth of the YRK increased; this increase was significant from January to June, and from September to December, but it was insignificant from July to August. The result is due to runoff timing not changing in the main glacier melt months of July and August. Moreover, the runoff depth of the KRY in June and July did not pass the 95% significance tests and, in the other months, significantly increased. Additionally, the runoff depth of the YUK increased significantly from January to May and from September to December, but insignificantly in June and August; however, the average monthly runoff depth of the YUK decreased in July. The KRK had similar interannual change in monthly runoff depth to that in the YUK but a much higher decrease rate of runoff depth in August than that in the YUK. This result might be because the YUK had much larger runoff depth in summer than in the other months, resulting in an insignificant increase in annual runoff depth [62]. However, the runoff depth of the KRK decreased in summer and significantly increased in winter, possibly because of river regulation due to the newly built Wuluwati Reservoir. As a result, the annual runoff depth at the Hydrological Station (HS) of Wuluwati was subject to change [63]. In addition, temperature rising in winter led to the wide degradation of permafrost in the study area, resulting in surface-water infiltration into groundwater in the early stage, an increase in groundwater storage capacity in the basin [64,65], and an increase in basin runoff in winter. resulting in surface-water infiltration into groundwater in the early stage, an increase in groundwater storage capacity in the basin [64,65], and an increase in basin runoff in winter. An M-K abrupt runoff change test in the four basins was carried out (Figure 12). Results showed that the annual runoff depth of the YRK and KRY changed abruptly in 1999, which was nearly 10 years later than the M-K abrupt-change results of the annual mean temperature in the two basins obtained by Ling et al. [66]. However, the YUK and KRK had no abrupt change in annual runoff depth. The M-K abrupt change test of monthly runoff change in the four river basins showed that the YRK underwent an abrupt increase of runoff depth in a total of nine months. The abrupt change in runoff depth occurred the earliest in May 1990 and the latest in March and April 2005, and mainly occurred in the other months in the late 1990s. The KRK suffered an abrupt change in months other than July, an abrupt decrease of runoff depth in August, especially in 1995, and an abrupt increase in other months, the earliest of which was in February 1978 and the latest of which was in April 2007. The YUK underwent an abrupt increase in runoff depth from January to May and from September to December, especially in January, February, November, and December in the late 1990s, and in the An M-K abrupt runoff change test in the four basins was carried out (Figure 12). Results showed that the annual runoff depth of the YRK and KRY changed abruptly in 1999, which was nearly 10 years later than the M-K abrupt-change results of the annual mean temperature in the two basins obtained by Ling et al. [66]. However, the YUK and KRK had no abrupt change in annual runoff depth. The M-K abrupt change test of monthly runoff change in the four river basins showed that the YRK underwent an abrupt increase of runoff depth in a total of nine months. The abrupt change in runoff depth occurred the earliest in May 1990 and the latest in March and April 2005, and mainly occurred in the other months in the late 1990s. The KRK suffered an abrupt change in months other than July, an abrupt decrease of runoff depth in August, especially in 1995, and an abrupt increase in other months, the earliest of which was in February 1978 and the latest of which was in April 2007. The YUK underwent an abrupt increase in runoff depth from January to May and from September to December, especially in January, February, November, and December in the late 1990s, and in the other months after 2000. The KRY suffered an abrupt change eight months after 2000, i.e., in November and December in 2001 and later from February to April after 2006.

Climate Change Impact on Runoff Depth
The Pearson correlation of annual runoff depth with the annual average temperature and precipitation in the YRK, KRK, YUK, and KRY was analyzed. Rivers mainly supplied by ice and snow meltwater usually have good correlation between temperature and precipitation. Therefore, we also analyzed the partial correlation between temperature and runoff depth and between precipitation and runoff depth. The analysis results are shown in Table 4. The YRK had a significant positive correlation between the runoff depth and temperature (single correlation coefficient was 0.360, and partial correlation coefficient was 0.356, both of which passed the test of significance at the 0.05 level); however, there was no correlation with the runoff depth and annual precipitation (both the single and partial correlation coefficients were close to 0), which was consistent with the results of Ling et al. [67], who showed that annual average temperature played a dominant role in the influence of annual runoff depth in the YRK. Although the KRK had a significant correlation between the annual runoff depth and annual precipitation (single correlation coefficient was −0.367 and partial correlation coefficient was −0.400, both of which passed the 0.05 test of significance), both the KRK and the YRK had insignificant positive correlation between annual runoff depth and annual average temperature and between annual runoff depth and annual precipitation, which disagreed with the runoff generation law of rivers that are mainly supplied by ice and snow meltwater [6,31]. This result may be due to other factors, such as earthquakes, regional geological structures, and glacier-blocked lakes, generating Glacier Lake Outburst Floods (GLOFs) that affected the impact of climatic elements on runoff [54,68,69]. The annual runoff depth was positively correlated with the annual average temperature and annual precipitation in the KRY. The annual precipitation influenced the annual runoff depth more greatly (single correlation coefficient was 0.411, and partial correlation coefficient was 0.426, both of which passed the 0.01 significance level test) than it influenced the annual average temperature (single correlation coefficient was 0.250, and the partial correlation coefficient was 0.277, both of which passed the 0.1 significance level test), which was in agreement with earlier results of Xu et al. [32], who showed that annual runoff depth was greatly influenced by both annual precipitation and average temperature. However, it was also found in this study that precipitation influenced runoff depth more greatly than temperature did in the KRY. The contribution of glaciers and snow melting to river flow is an important parameter for analysis of the impact of temperature and precipitation changes on river flow [70]. Singh et al. [18], studying the total flow of the Himalayan

Climate Change Impact on Runoff Depth
The Pearson correlation of annual runoff depth with the annual average temperature and precipitation in the YRK, KRK, YUK, and KRY was analyzed. Rivers mainly supplied by ice and snow meltwater usually have good correlation between temperature and precipitation. Therefore, we also analyzed the partial correlation between temperature and runoff depth and between precipitation and runoff depth. The analysis results are shown in Table 4. The YRK had a significant positive correlation between the runoff depth and temperature (single correlation coefficient was 0.360, and partial correlation coefficient was 0.356, both of which passed the test of significance at the 0.05 level); however, there was no correlation with the runoff depth and annual precipitation (both the single and partial correlation coefficients were close to 0), which was consistent with the results of Ling et al. [67], who showed that annual average temperature played a dominant role in the influence of annual runoff depth in the YRK. Although the KRK had a significant correlation between the annual runoff depth and annual precipitation (single correlation coefficient was −0.367 and partial correlation coefficient was −0.400, both of which passed the 0.05 test of significance), both the KRK and the YRK had insignificant positive correlation between annual runoff depth and annual average temperature and between annual runoff depth and annual precipitation, which disagreed with the runoff generation law of rivers that are mainly supplied by ice and snow meltwater [6,31]. This result may be due to other factors, such as earthquakes, regional geological structures, and glacier-blocked lakes, generating Glacier Lake Outburst Floods (GLOFs) that affected the impact of climatic elements on runoff [54,68,69]. The annual runoff depth was positively correlated with the annual average temperature and annual precipitation in the KRY. The annual precipitation influenced the annual runoff depth more greatly (single correlation coefficient was 0.411, and partial correlation coefficient was 0.426, both of which passed the 0.01 significance level test) than it influenced the annual average temperature (single correlation coefficient was 0.250, and the partial correlation coefficient was 0.277, both of which passed the 0.1 significance level test), which was in agreement with earlier results of Xu et al. [32], who showed that annual runoff depth was greatly influenced by both annual precipitation and average temperature. However, it was also found in this study that precipitation influenced runoff depth more greatly than temperature did in the KRY. The contribution of glaciers and snow melting to river flow is an important parameter for analysis of the impact of temperature and precipitation changes on river flow [70]. Singh et al. [18], studying the total flow of the Himalayan rivers, found that the contribution of snow cover and glacier melting to the total flow of rivers increased with the increase of altitude. The difference of vertical gradient is the main factor affecting regional snow cover and glacier melting [43,46] because of the relatively low elevation and more precipitation of the Keriya River Basin, and lead to precipitation contributed more to runoff than temperature in the region. Table 4. Pearson correlation coefficients and partial correlation coefficients between annual runoff (R) and mean air temperature (T) and precipitation (P) in four basins. Note: ***, **, and *, significance levels of 0.001, 0.01, and 0.05, respectively. |, partial correlation analysis between R and removal of P (|P) or T (|T).

Basin/Correlation
The study area had high temperatures and abundant rainfall in summer (i.e., June, July, and August), and the runoff of the four basins was mainly recharged by ice and snowmelt water. Therefore, we mainly analyzed the relationship between runoff and temperature and between runoff and precipitation in summer ( Figure 13). According to the change in SCA (Figures 7 and 8), the YRK, KRK, and YUK had a positive regression coefficient between temperature and runoff depth from June to July, but the KRK and YUK had much smaller regression coefficients between temperature and runoff depth in August than in June and July. On the basis of the main months when peak temperature (Table 3) occurred in the two basins, the average temperature in August was close to that in July and much higher than that in June. Therefore, the two basins should have had significant correlation between runoff depth and temperature in August in principle, but this was not the case. Results also showed that the relationship between runoff depth and temperature in the KRK and YUK in August did not conform to runoff yield and the concentration law of rivers dominated by glacier recharge. In addition, the YRK, KRK, and YUK had a smaller regression coefficient between runoff depth and precipitation than between runoff depth and temperature in summer, showing that temperature greatly influenced runoff depth in the three basins. However, when examining runoff depth, temperature, and precipitation ( Figure 13(d1) or Figure 13(d2)), we found that the KRY had a larger regression coefficient between monthly precipitation and runoff depth than with temperature in summer, indicating that monthly runoff depth was mainly influenced by precipitation in the KRY in summer. Zhao, et al. [71] suggested that the difference between the KRK and YUK might be caused by a decrease in upper air temperature that led to glacier ablation. However, Zhou, et al. [72] argued that the difference in site conditions led to the formation of a microclimate between regions, resulting in the mitigation of ice and snow melt. In this study, not all the different weather stations reflect the importance of microclimates. However, it can provide a reference for better research in this region. Therefore, the effects of topography, sudden disasters, and human activities should also be comprehensively considered for the difference in runoff between the two basins in summer.

Uncertainty Analysis
This study focused on investigating the annual and interannual snow-cover and runoff-change characteristics of four typical basins of the NSWKM and on analyzing their relationship with temperature and precipitation. However, snow melt was also affected by slope, slope direction, and type of surface cover. Additionally, snow was influenced by the distribution of wind speed and surface temperature. Although the data from ground-observation stations used in this study were representative of different basins, they were still unable to fully reflect the characteristics of climate change at different altitudes (e.g., precipitation was also affected by elevation, slope, and slope direction). In addition, we have not given sufficient consideration to the debris-covered area in our classification of glacier area. Although this area is small, it may affect our final results to some extent. In view of the above problems, we will adopt the new methods for further analysis in the following studies. Zhao, et al. [71] suggested that the difference between the KRK and YUK might be caused by a decrease in upper air temperature that led to glacier ablation. However, Zhou, et al. [72] argued that the difference in site conditions led to the formation of a microclimate between regions, resulting in the mitigation of ice and snow melt. In this study, not all the different weather stations reflect the importance of microclimates. However, it can provide a reference for better research in this region. Therefore, the effects of topography, sudden disasters, and human activities should also be comprehensively considered for the difference in runoff between the two basins in summer.

Uncertainty Analysis
This study focused on investigating the annual and interannual snow-cover and runoff-change characteristics of four typical basins of the NSWKM and on analyzing their relationship with temperature and precipitation. However, snow melt was also affected by slope, slope direction, and type of surface cover. Additionally, snow was influenced by the distribution of wind speed and surface temperature. Although the data from ground-observation stations used in this study were representative of different basins, they were still unable to fully reflect the characteristics of climate change at different altitudes (e.g., precipitation was also affected by elevation, slope, and slope direction). In addition, we have not given sufficient consideration to the debris-covered area in our classification of glacier area. Although this area is small, it may affect our final results to some extent. In view of the above problems, we will adopt the new methods for further analysis in the following studies.

Conclusions
As one of the key factors that restrict the socioeconomic development of oases in arid areas of Northwest China, water resources mainly come from ice and snow meltwater in high mountains. Understanding the runoff-change law and the runoff-generation mechanism in the arid regions of Northwest China is very significant for understanding the water cycle in cold and arid regions. On the basis of meteorological and hydrological data from the observation stations from the last 60 years (1957-2017) and snow-cover data collected by MODIS (2001MODIS ( -2016 in the study area, we used various mathematical and statistical methods to study the relationship between temperature, precipitation, snow cover, and runoff in typical basins (i.e., the YRK, KRK, YUK, and KRY) on the NSWKM. The NSWKM has had increasing temperature (0.26 • C/decade) and precipitation (5.50 mm/decade) in the past 60 years, especially in the KRK and YUK (R 2 > 0.82).
The average SCF of four typical basins in the NSWKM is: YUK (SCF = 44.14%) > YRK (SCF = 38.73%) > KRY (SCF = 33.42%) > KRK (SCF = 33.40%). From them, the YRK and YUK had decreasing SCA values (slope < −15.39), while the KRK and KRY had increasing SCA values (slope > 1.87). In seasonal variation, except for the SCF of YRK in spring being lower than that in winter (−2.6%), the SCF of the three other basins reached the maximum value in spring, with the most significant performance being in the YUK (SCF = 26.4%).
Variation in runoff depth in the four studied basins of the NSWKM showed an increasing trend. Between them, the YRK (5.78 mm/decade) had the most significant runoff depth, and the YUK (1.58 mm/decade) had the most insignificant increasing runoff depth. In the variation of monthly runoff depth, except for the maximum monthly runoff depth of the YRK that appeared in August, the maximum monthly runoff depth of the three other basins appeared in July. The cause was analyzed in combination with climate factors (e.g., temperature, precipitation) and snow-cover change. We found that temperature was the main factor affecting annual and monthly runoff variations in summer in the YRK, which was in agreement with the runoff-generation law for arid areas. However, temperature and precipitation were not the main factors affecting annual and monthly runoff variations in summer in the YUK and KRK. Therefore, runoff in the KRK and YUK is disturbed by other factors (e.g., slope and slope direction) during runoff generation, resulting in disruptions of the relationship between runoff and climate factors; hence, runoff changes in the KRK and YUK disagree with the runoff-generation law for rivers dominated by ice-and snow-meltwater recharge.