Changes in Snow Phenology from 1979 to 2016 over the Tianshan Mountains , Central Asia

Snowmelt from the Tianshan Mountains (TS) is a major contributor to the water resources of the Central Asian region. Thus, changes in snow phenology over the TS have significant implications for regional water supplies and ecosystem services. However, the characteristics of changes in snow phenology and their influences on the climate are poorly understood throughout the entire TS due to the lack of in situ observations, limitations of optical remote sensing due to clouds, and decentralized political landscapes. Using passive microwave remote sensing snow data from 1979 to 2016 across the TS, this study investigates the spatiotemporal variations of snow phenology and their attributes and implications. The results show that the mean snow onset day (Do), snow end day (De), snow cover duration days (Dd), and maximum snow depth (SDmax) from 1979 to 2016 were the 78.2nd day of hydrological year (DOY), 222.4th DOY, 146.2 days, and 16.1 cm over the TS, respectively. Dd exhibited a spatial distribution of days with a temperature of <0 ◦C derived from meteorological station observations. Anomalies of snow phenology displayed the regional diversities over the TS, with shortened Dd in high-altitude regions and the Fergana Valley but increased Dd in the Ili Valley and upper reaches of the Chu and Aksu Rivers. Increased SDmax was exhibited in the central part of the TS, and decreased SDmax was observed in the western and eastern parts of the TS. Changes in Dd were dominated by earlier De, which was caused by increased melt-season temperatures (Tm). Earlier De with increased accumulation of seasonal precipitation (Pa) influenced the hydrological processes in the snowmelt recharge basin, increasing runoff and earlier peak runoff in the spring, which intensified the regional water crisis.


Introduction
Snow plays a critical role in regional and global water cycles, as well as in climate systems [1][2][3].Snow cover influences land surface energy budgets and atmospheric circulation patterns due to its high surface albedo and good thermal insulation [4][5][6][7][8].Snow phenology, i.e., the snow onset day, snow end day, and snow cover duration, is essential in the representation of snow variability and has direct implications on the growth of vegetation, snowmelt timing, freshwater supply, and irrigation in snowmelt-dominated basins [5,[9][10][11][12][13][14].Therefore, it is necessary and meaningful to investigate the spatiotemporal variability of snow phenology.
Compared with the sparse observations of meteorological stations, remote sensing is an effective way to monitor snow dynamics of past decades on regional and global scales [4,8,[15][16][17].Optical remote sensing uses the normalized difference snow index (NDSI) to identify the snow cover extent with a high resolution [18,19].Snow phenology conditions are mainly influenced by the variability of the snow cover extent [8,20,21].Increasing temperatures have shortened the snow cover duration (D d ) [6,22], delayed the snow onset day (D o ) [23], and resulted in an earlier snow end day (D e ) [24].In contrast, enhanced mid-latitude westerlies have caused an enlarged snow cover extent in High Asia and Tibet with prolonged D d , earlier D o , and later D e [25][26][27].However, the cloud mask and lack of snow depth detection methods limit the application of snow products from optical remote sensing [28][29][30][31].Passive microwave remote sensing based on the microwave spectral gradient method provides a potential way to monitor the spatiotemporal variations of snow depth and snow cover under cloudy conditions [32][33][34].The snow depth products retrieved from passive microwave sensors have been successfully applied in vegetation phenology [35][36][37], snow climatology, and snow hydrology [12,31,[38][39][40].
Situated in inland Eurasia and far away from the oceans, the Tianshan Mountains (TS) are called the water tower of Central Asia.The main rivers that are recharged by glacier/snow melt water (e.g., the Ili River, Syr Darya River, Amu Darya River, Tarim River, and Chu River) originate from the TS, forming one of the largest irrigated zones in the world [41][42][43].Previous studies indicated that the TS have experienced a significant warming trend over the past few decades [44][45][46].The snowfall/precipitation ratio showed a decreasing trend accompanied by increasing temperatures in the TS [47], which intensified glacier and snow melt [2], leading to increasing runoff and earlier peak runoff from glacier/snow melt in these recharging river basins [45,[48][49][50].Several studies have investigated the impact of glacier change on water resources in the TS [42,[51][52][53]], but only a few studies have addressed the snow cover changes with optical remote sensing [54,55].In addition, different moisture sources cause diverse climate zones in the TS [56,57], but most studies have focused on only a part of the TS, especially on the areas within China [43,58].Research on the snow phenology concerning the entire TS and its response to the different climate types in the sub-regions is lacking.Therefore, it is necessary to survey the spatiotemporal variability of snow phenology and its potential influencing factors throughout the entire TS.
This study investigated the variability of snow phenology from 1979 to 2016 based on passive microwave snow depth data from the TS.The objectives of this study were as follows: (1) to derive snow phenology, including D o , D e , D d , and SD max from daily snow depth data; (2) to analyze the spatiotemporal distributions and trends of snow phenology, as well as their potential influencing factors; and (3) to discuss the potential impact of the changes of snow phenology on water resources.This is the first study to describe snow phenology using passive microwave snow depth data for the entire TS.The results will help to improve the understanding of snow phenology and lead to better management of the regional water resources.

Study Area
As the largest mountain system in Central Asia, the TS are approximately 250-350 km wide and over 2500 km long, spanning from Uzbekistan to Kyrgyzstan, southeastern Kazakhstan, and Xinjiang (China) with 800,000 km 2 (Figure 1a) [59].The TS have abundant precipitation due to westerly circulation and unique topography, exhibiting heavily glaciated and snow-covered regions [42,60].The people living in the surrounding areas are heavily dependent on glacier/snow melt water for their fresh water supply [41,61,62].The TS are geographically divided into four parts: I. the Western Tianshan Mountains (WTS), II.The Northern Tianshan Mountains (NTS), III.The Central Tianshan Mountains (CTS), and IV.The Eastern Tianshan Mountains (ETS) [41,43].The annual mean precipitation across the TS is 329.3 mm, which is concentrated between April and June, and the annual mean temperature is 4.6 • C (Figure 1c).The WTS and NTS are characterized by a relatively humid climate, while the CTS and ETS have a typical continental climate [63].The maximum precipitation in the NTS and ETS occurs in the spring and early summer, which is later than that in the WTS (from late winter to early spring) but earlier than that in the CTS (summer) [59].Mountains (CTS), and IV.The Eastern Tianshan Mountains (ETS) [41,43].The annual mean precipitation across the TS is 329.3 mm, which is concentrated between April and June, and the annual mean temperature is 4.6 °C (Figure 1c).The WTS and NTS are characterized by a relatively humid climate, while the CTS and ETS have a typical continental climate [63].The maximum precipitation in the NTS and ETS occurs in the spring and early summer, which is later than that in the WTS (from late winter to early spring) but earlier than that in the CTS (summer) [59].

Precipitation and Temperature Datasets
The daily surface air temperature and precipitation data from 50 meteorological stations in the TS for 1979-2016 were collected from the China Meteorological Administration (CMA) (http://data.cma.cn) and National Climate Data Center (NCDC) (https://www.ncdc.noaa.gov).The

Remote Sensing Snow Depth Products
The passive microwave snow depth dataset, within the range from 60-140 • E and 15-55 • N at 25 km spatial resolution (http://westdc.westgis.ac.cn), which is derived from the brightness temperature data from the National Snow and Ice Data Center (NSIDC), including the scanning multichannel radiometer (SMMR) (1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987), special sensor microwave/imager (SSM/I) , and special sensor microwave imager/sounder (SSMI/S) (2008-2016), was employed in this study.The time consistency of this dataset was improved by cross calibration of the light temperatures of different sensors [64].This daily database was developed through a modified Chang algorithm, which was based on the in situ snow data of China, taking the impact of forest cover, liquid water content in snow layer, and surface water body into consideration [11,30,64,65].Compared with the snow depth from SSM/I and SSMI/S, this dataset has high accuracy and low biases in China [65] and allows the successful investigation of snow depth in China, the Tibetan Plateau, northeastern China, and Xinjiang [11,27,30,31].

Precipitation and Temperature Datasets
The daily surface air temperature and precipitation data from 50 meteorological stations in the TS for 1979-2016 were collected from the China Meteorological Administration (CMA) (http://data.cma.cn) and National Climate Data Center (NCDC) (https://www.ncdc.noaa.gov).The data from meteorological stations were applied to analyze the number of days during which the temperature was below 0 • C in the study area.Monthly gridded surface air temperature  was derived from the Global Historical Climatology Network 2 and Climate Anomaly Monitoring System (GHCN_CAMS, http://www.esrl.noaa.gov/psd/).The gridded 2 m temperature datasets were merged from two large individual datasets of station observations with 0.5 • × 0.5 • spatial resolution [66].The gridded surface air temperature data were resampled to 0.25 • × 0.25 • spatial resolution for resolution consistency.The Global Precipitation Climatology Centre products (GPCC, https://www.dwd.de/EN/ourservices/gpcc/gpcc.html) were more suitable for precipitation research in Central Asia when compared with the climatic research unit (CRU) and Willmott and Matsuura (WM) precipitation data [67].Therefore, monthly precipitation gridded data at 0.25 • × 0.25 • spatial resolution from 1979 to 2016 were derived from the GPCC (V2018), which were based on data from more than 79,000 stations [68].The gridded data concerning air temperature and precipitation were used to investigate the relationship between snow phenology changes and climate variations in the TS.As a typical snowmelt-dominated river in the TS, the Kaidu River was chosen to discuss the variability of snow phenology impacts on water resources.The daily runoff data in the Kaidu River from 1979 to 2011 were collected from the Xinjiang Hydrological Bureau.

Snow Phenology Calculation
To investigate the variability of snow phenology over the TS, this study defined the period from September to the following August as a snow hydrological year.The snow cover onset day (D o ) for the hydrological year was defined as the first day of the first five consecutive snow cover days on each grid to reduce the influence of instantaneous snow on the snow phenology computations [27].The snow cover end day (D e ) for the hydrological year was defined as the last day of the last five consecutive snow cover days on each grid.The snow cover duration (D d ) for the hydrological year was calculated from D o to D e .The period from September to the following February was defined as the accumulation season, and the period from March to August was defined as the melt season, with the maximum snow depth (SD max ) peaking around March (Figure 1d).

Sensitivity Analysis of Snow Phenology
The mean monthly air temperature (T a ) and total monthly precipitation (P a ) during the accumulation season, maximum snow depth (SD max ), mean monthly air temperature (T m ), and total monthly precipitation (P m ) during the melt season were applied to quantify their contributions to snow phenology over the TS.According to Peng et al. [1,69] and Chen et al. [1,69], a multiple linear regression model (Equation ( 1)) can be used to analyze the variability of D o , which was mainly driven by T a and P a .
where the regression coefficients β 1 and β 2 were defined as the sensitivity of D o to T a and the sensitivity of D o to P a , respectively.D e was mainly determined by the T m and SD max , and the SD max in the snow season was dependent on T a and P a [1,69].Thus, Equations ( 2) and (3) were employed in the sensitivity analysis of SD max and D e , respectively.
The regression coefficients β 4 and β 5 in Equation ( 2) were defined as the sensitivity of SD max to T a and SD max to P a .Similarly, the regression coefficients β 7 and β 8 in Equation (3) were defined as the sensitivity of D e to T m and D e to SD max .β 3 , β 6 , and β 9 are the residual error.Therefore, the sensitivity analysis was performed based on Equations ( 1)-( 3) to detect the attributes of the snow phenology changes over the TS.

Trend Analysis and Statistical Analysis
The trends in snow phenology, air temperature, and precipitation were calculated by the Mann-Kendall (M-K) test [70,71], which is a nonparametric method of monotonic trends that has successfully been applied to detect trends in a time series [72,73].The Sen method was used to estimate the slope of the trend [74].The significance of a trend was assumed if the Z (M-K) value was unequal to zero with a significance level of less than 0.05.The variations in snow phenology and its relationship with air temperature and precipitation were investigated by correlation analysis and linear regression.3c).The SD max in the upper reaches of the Kaidu, Ili, and Syr Darya Rivers, and the eastern ETS showed a significant downward trend (Figure 3d).The areas with a significant increase in SD max were located in the Ili Valley and the upper reaches of the Aksu and Chu Rivers.

Changes in Snow Phenology over the TS
De and Dd in the TS experienced a significant downward trend at the 0.05 significance level (Figure 4a).Spatial distributions of temporal trends in mean Do, De, Dd, and SDmax over the TS during 1979-2016 are displayed in Figure 3. Do was earlier in the Ili Valley and the upper reaches of the Chu and Syr Darya Rivers but delayed on the southern slope of the TS; in particular, it was significantly later in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the ETS, and the Fergana Valley (Figure 3a).Compared with the changes in Do, De was earlier across almost all of the TS; in particular, it was significantly earlier in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the ETS, and the Fergana Valley but delayed in sporadic regions (Figure 3b).Combined with the variations of Do and De, Dd decreased in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the northern slope of the ETS, and the Fergana Valley but increased in the upper reaches of the Chu River and the Ili Valley (Figure 3c).The SDmax in the upper reaches of the Kaidu, Ili, and Syr Darya Rivers, and the eastern ETS showed a significant downward trend (Figure 3d).The areas with a significant increase in SDmax were located in the Ili Valley and the upper reaches of the Aksu and Chu Rivers.

Spatiotemporal Variations of Temperature and Precipitation over the TS
This study calculated the number of days with mean air temperatures below 0 °C from meteorological station observations over the TS from 1979 to 2016 to compare the results of Dd derived from passive microwave instruments.Small values for days with a mean air temperature <0 °C occurred (<120 days) at the southern edge of the TS, in the Ili Valley, and in the WTS (Figure 5a).This pattern was generally consistent with the spatial distribution of Dd (Figure 2c).Days with a mean air temperature <0 °C at all the stations showed a downward trend, and the significantly decreasing (significance level at 0.05) stations were located in the Aksu and Kaidu Rivers, Ili Valley, and ETS (Figure 5b).
As displayed in Figure 5c,d, the Ta and Tm experienced a significant increase in almost all of the regions over the TS (significance level at 0.05).Moreover, both the rising rates of Ta and Tm in the upper reaches of the Syr Darya, Aksu River, NTS, and ETS reached 0.07 °C per year.The increasing rate in Tm was larger than that in Ta over the TS.Pa and Pm in all regions over the TS from 1979-2016 showed an upward trend, except for the upper reaches of the Syr Darya River (Figure 5e,f).Pa with significant increase was observed in the upper reaches of the Aksu River, Ili Valley, and the northern WTS, where the increase rate was above 1 mm per year.In contrast to Pa, Pm showed significant increase in the western area of the ETS and the upper reaches of the Aksu, Ili, and Kaidu Rivers.Although significantly increased winter temperatures over the TS were reported [43,44,47], strengthened westerlies brought more air moisture flux [25,76,77] and increased snowfall, which was beneficial to snow cover accumulation [47,78].The increasing rate of precipitation was more profound in the west TS during winter and spring and in the east TS during summer [43].

Spatiotemporal Variations of Temperature and Precipitation over the TS
This study calculated the number of days with mean air temperatures below 0 • C from meteorological station observations over the TS from 1979 to 2016 to compare the results of D d derived from passive microwave instruments.Small values for days with a mean air temperature <0 • C occurred (<120 days) at the southern edge of the TS, in the Ili Valley, and in the WTS (Figure 5a).This pattern was generally consistent with the spatial distribution of D d (Figure 2c).Days with a mean air temperature <0 • C at all the stations showed a downward trend, and the significantly decreasing (significance level at 0.05) stations were located in the Aksu and Kaidu Rivers, Ili Valley, and ETS (Figure 5b).  5 indicate that the trends were significant (significance level at 0.05).

Sensitivity Analysis of Snow Phenology over the TS
Sensitivity analysis was performed based on Equations ( 1)-( 3) to detect the attributes of the snow phenology changes over the TS, and the results are shown in Figure 6.The results showed a positive sensitivity of Do to Ta at 0.91 days °C−1 across the TS from 1979 to 2016, which indicated that the Do would be delayed 0.91 days, with Ta increased by 1 °C.In contrast, the significant negative sensitivity of Do to Pa was −0.11 days mm −1 (significance level at 0.01), which indicated that Do would be advanced 0.11 days, with Pa increased by 1 mm.The results indicated that the increased Ta would hinder the increase of SDmax (−0.67 cm °C−1 ), but the increased Pa would promote the significant increase of SDmax (0.04 cm mm −1 , significance level at 0.05).Compared with the sensitivity of De to SDmax (0.29 days cm −1 ), the magnitude of sensitivity of De to Tm (−5.36 days °C−1 , significance level at 0.01) was significantly larger (Figure 6e,f).As displayed in Figure 5c,d, the T a and T m experienced a significant increase in almost all of the regions over the TS (significance level at 0.05).Moreover, both the rising rates of T a and T m in the upper reaches of the Syr Darya, Aksu River, NTS, and ETS reached 0.07 • C per year.The increasing rate in T m was larger than that in T a over the TS.P a and P m in all regions over the TS from 1979-2016 showed an upward trend, except for the upper reaches of the Syr Darya River (Figure 5e,f).P a with significant increase was observed in the upper reaches of the Aksu River, Ili Valley, and the northern WTS, where the increase rate was above 1 mm per year.In contrast to P a , P m showed significant increase in the western area of the ETS and the upper reaches of the Aksu, Ili, and Kaidu Rivers.Although significantly increased winter temperatures over the TS were reported [43,44,47], strengthened westerlies brought more air moisture flux [25,76,77] and increased snowfall, which was beneficial to snow cover accumulation [47,78].The increasing rate of precipitation was more profound in the west TS during winter and spring and in the east TS during summer [43].

Sensitivity Analysis of Snow Phenology over the TS
Sensitivity analysis was performed based on Equations ( 1)-( 3) to detect the attributes of the snow phenology changes over the TS, and the results are shown in Figure 6.The results showed a positive sensitivity of D o to T a at 0.91 days • C −1 across the TS from 1979 to 2016, which indicated that the D o would be delayed 0.91 days, with T a increased by 1 • C. In contrast, the significant negative sensitivity of D o to P a was −0.11 days mm −1 (significance level at 0.01), which indicated that D o would be advanced 0.11 days, with P a increased by 1 mm.The results indicated that the increased T a would hinder the increase of SD max (−0.67 cm • C −1 ), but the increased P a would promote the significant increase of SD max (0.04 cm mm −1 , significance level at 0.05).Compared with the sensitivity of D e to SD max (0.29 days cm −1 ), the magnitude of sensitivity of D e to T m (−5.36 days • C −1 , significance level at 0.01) was significantly larger (Figure 6e,f).

Factors Driving Snow Phenology Changes over the TS
Changes in air temperature and precipitation played a key role in the variations in the snow phenology over the TS.Do was more sensitive to Ta than Pa, which could be explained by the delayed Do over the TS (Figure 4a).Do was earlier in the Ili Valley and the upper reaches of the Chu and Syr Darya Rivers due to the significantly increased Pa and negative sensitivity of Do to Pa (Figures 3a, 5e,  and 6b).The significantly increased Ta and highly positive sensitivity of Do to Ta caused the significantly delayed Do in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the ETS, and the Fergana Valley (Figures 3a, 5e, and 6a).The significantly increased Tm and highly negative sensitivity of De to Tm resulted in the significantly earlier De over the TS, especially in high altitude regions of the Aksu, Kaidu, and Ili Rivers, the northern ETS, and the Fergana Valley.The significantly increased SDmax and highly positive sensitivity of De to SDmax promoted the delayed De.Previous studies proved that high net radiation with thin snow cover in spring across the Northern Hemisphere created a high sensitivity of De to temperature, which resulted in earlier De over Eurasia [6,69,79].
Although the magnitude of sensitivity of SDmax to Ta was larger than Pa (Figure 6c,d), the increased rate of Pa was larger than Ta (Figure 5c,e).For example, areas with significantly increased SDmax, such as the central part of the ETS and the east of the WTS, as well as the west of the NTS and CTS, were caused by the significantly increased Pa and highly positive sensitivity of SDmax to Pa (Figures 3d, 5e, and 6d).The significantly increased Ta and highly negative sensitivity of SDmax to Ta were attributed to the significantly decreased SDmax in the eastern ETS and the upper reaches of the Kaidu, Ili, and Syr Darya Rivers (Figures 3d, 5c, and 6c).
The increased SDmax and earlier Do caused longer Dd in the upper reaches of the Chu River and the Ili Valley (Figure 3c).The significantly increased Tm advanced the De and subsequently shortened Dd in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the northern ETS, and the Fergana Valley.De dominated the snow season variations due to the higher Tm across the TS, which was coincidental with the studies in the Northern Hemisphere [6,69,79].However, the Do in the Tibetan Plateau controlled the variability of snow phenology [27].In contrast, decreased Ta, which was driven

Factors Driving Snow Phenology Changes over the TS
Changes in air temperature and precipitation played a key role in the variations in the snow phenology over the TS.D o was more sensitive to T a than P a , which could be explained by the delayed D o over the TS (Figure 4a).D o was earlier in the Ili Valley and the upper reaches of the Chu and Syr Darya Rivers due to the significantly increased P a and negative sensitivity of D o to P a (Figures 3a,  5e and 6b).The significantly increased T a and highly positive sensitivity of D o to T a caused the significantly delayed D o in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the ETS, and the Fergana Valley (Figures 3a, 5e and 6a).The significantly increased T m and highly negative sensitivity of D e to T m resulted in the significantly earlier D e over the TS, especially in high altitude regions of the Aksu, Kaidu, and Ili Rivers, the northern ETS, and the Fergana Valley.The significantly increased SD max and highly positive sensitivity of D e to SD max promoted the delayed D e .Previous studies proved that high net radiation with thin snow cover in spring across the Northern Hemisphere created a high sensitivity of D e to temperature, which resulted in earlier D e over Eurasia [6,69,79].
Although the magnitude of sensitivity of SD max to T a was larger than P a (Figure 6c,d), the increased rate of P a was larger than T a (Figure 5c,e).For example, areas with significantly increased SD max , such as the central part of the ETS and the east of the WTS, as well as the west of the NTS and CTS, were caused by the significantly increased P a and highly positive sensitivity of SD max to P a (Figures 3d, 5e and 6d).The significantly increased T a and highly negative sensitivity of SD max to T a were attributed to the significantly decreased SD max in the eastern ETS and the upper reaches of the Kaidu, Ili, and Syr Darya Rivers (Figures 3d, 5c and 6c).
The increased SD max and earlier D o caused longer D d in the upper reaches of the Chu River and the Ili Valley (Figure 3c).The significantly increased T m advanced the D e and subsequently shortened D d in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the northern ETS, and the Fergana Valley.D e dominated the snow season variations due to the higher T m across the TS, which was coincidental with the studies in the Northern Hemisphere [6,69,79].However, the D o in the Tibetan Plateau controlled the variability of snow phenology [27].In contrast, decreased T a , which was driven by the Pacific Ocean's surface cooling in winter and amplified Arctic warming effects [80][81][82], resulted in earlier D o in the Tibetan Plateau (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) [27,79].

Potential Impact of Snow Phenology Changes on Water Resources
Snow melt is an important recharge source for the TS rivers [83], which means that the variability of snow phenology significantly impacted the local runoff regimes.The Kaidu River is a typical snowmelt runoff river in the TS [48,84], where spring runoff from the Bayanbulak and Dashankou stations increased by 36.5% and 11.21%, respectively, and the runoff peak advanced between the 1980s and 2000s (Figure 7a,b).Although D d and SD max showed downward trends in the Kaidu River, the increased P m and advanced D e caused runoff changes.The advanced D e , with increased precipitation in the upper reaches of the Aksu, Syr Darya, Kaidu, and the Ili Rivers contributed to increased snowmelt earlier in spring.The maximum runoff peak in the snowmelt-dominated rivers over the TS showed an advanced trend, but the runoff in summer changed insignificantly [45].The runoff pattern of the Syr Darya River shifted from spring and early summer to late winter and early spring due to climate warming [45].In the Toshkan River, an important tributary of the Aksu River, snowmelt runoff in the cold season increased by 65% and 56% due to the rising temperatures during the periods from 1960 to 1997 and from 1998 to 2015 [45,85], respectively.In addition, some studies indicated that the snowfall/precipitation ratio was expected to decrease in the 21st century over the TS, which could reduce the snow cover area and aggravate snow melt [1,86].The runoff shift in snowmelt-dominated rivers has a big impact on water resource management, which may increase the risk of spring flooding and cause a time mismatch between irrigation watering and crop growth, considering population growth and climate warming.

Potential Impact of Snow Phenology Changes on Water Resources
Snow melt is an important recharge source for the TS rivers [83], which means that the variability of snow phenology significantly impacted the local runoff regimes.The Kaidu River is a typical snowmelt runoff river in the TS [48,84], where spring runoff from the Bayanbulak and Dashankou stations increased by 36.5% and 11.21%, respectively, and the runoff peak advanced between the 1980s and 2000s (Figure 7a,b).Although Dd and SDmax showed downward trends in the Kaidu River, the increased Pm and advanced De caused runoff changes.The advanced De, with increased precipitation in the upper reaches of the Aksu, Syr Darya, Kaidu, and the Ili Rivers contributed to increased snowmelt earlier in spring.The maximum runoff peak in the snowmelt-dominated rivers over the TS showed an advanced trend, but the runoff in summer changed insignificantly [45].The runoff pattern of the Syr Darya River shifted from spring and early summer to late winter and early spring due to climate warming [45].In the Toshkan River, an important tributary of the Aksu River, snowmelt runoff in the cold season increased by 65% and 56% due to the rising temperatures during the periods from 1960 to 1997 and from 1998 to 2015 [45,85], respectively.In addition, some studies indicated that the snowfall/precipitation ratio was expected to decrease in the 21st century over the TS, which could reduce the snow cover area and aggravate snow melt [1,86].The runoff shift in snowmelt-dominated rivers has a big impact on water resource management, which may increase the risk of spring flooding and cause a time mismatch between irrigation watering and crop growth, considering population growth and climate warming.

Limitation and Outlook
The passive microwave data could provide effective and long series of daily snow depth data to monitor the snow phenology over the TS from 1979 to 2016.However, there is uncertainty in snow depth retrieval at high-altitude regions in the study area where in situ observations are lacking.This uncertainty also leads to limitations in the separation of the snow accumulation season and snow melt season by the maximum snow depth in high-altitude regions.Several studies pointed out that passive microwave remote sensing has a limited ability to detect wet snow during the snow melt season, which may underestimate the Dd [87][88][89][90].Meanwhile, misclassification and error in deriving snow cover were attributed to relatively coarse spatial resolution, as well as the complexity of snow characteristics and topography [91][92][93][94].Combined with optical remote sensing, passive microwave remote sensing and a land surface model can effectively improve the monitoring accuracy of snow phenology and snow depth [95].In addition, five consecutive snow cover days in each grid for

Limitation and Outlook
The passive microwave data could provide effective and long series of daily snow depth data to monitor the snow phenology over the TS from 1979 to 2016.However, there is uncertainty in snow depth retrieval at high-altitude regions in the study area where in situ observations are lacking.This uncertainty also leads to limitations in the separation of the snow accumulation season and snow melt season by the maximum snow depth in high-altitude regions.Several studies pointed out that passive microwave remote sensing has a limited ability to detect wet snow during the snow melt season, which may underestimate the D d [87][88][89][90].Meanwhile, misclassification and error in deriving snow cover were attributed to relatively coarse spatial resolution, as well as the complexity of snow characteristics and topography [91][92][93][94].Combined with optical remote sensing, passive microwave remote sensing and a land surface model can effectively improve the monitoring accuracy of snow phenology and snow depth [95].In addition, five consecutive snow cover days in each grid for calculating D o and D e can reduce the influence of instantaneous snow on snow phenology computations.However, this method may also cause later D o and earlier D e , consequently underestimating the D d .It neglects the spatial heterogeneity of snow by using a fixed day to define the snow melt season and snow accumulation season in each grid, bringing uncertainty into the sensitivity analysis of the snow phenology.Sensitivity analysis based on the monthly precipitation and air temperature cannot discriminate between the effect of snow and rain on snow phenology due to the insufficient daily in situ observations.A comprehensive equation, which considers the change rate for sensitivity analysis based on the daily data, will be more quantitatively capable of describing the change in the snow phenology parameters over the years.The daily data from the climate model simulations and satellite products can be used for further study.Other climate variables (such as relative humidity and wind) and complex topography also have an important impact on snow, which require intensive observation networks and large amounts of field survey data.Alpine vegetation plays a crucial role in a mountain ecosystem, so the variability of snow phenology and snow depth have been reported to significantly influence vegetation growth [24,96,97].For example, a shorter D d would lead to an earlier start and longer growing season on the Tibetan Plateau [10,37].The snow water equivalent and D d play an equal role in the growth of grassland and sparse vegetation [98].Moreover, the mean snow depth significantly regulates desert vegetation growth by persistently impacting soil moisture [36].The TS has plenty of vegetation and well-developed animal husbandry [99].Therefore, the relationships between the variations in snow phenology and vegetation growth over the TS need further research.

Conclusions
This study investigated the spatiotemporal variability of snow phenology and snow depth over the TS using a passive microwave daily snow depth dataset from 1979 to 2016, as well as exploring the impacts and attributes of snow changes.The main findings of the study are as follows: 1.
The snow end day and snow cover duration across the TS experienced a significant decrease, and the snow end day dominated the variability of snow season.The snow end day was earlier across the TS with increasing air temperatures during the melt season, especially in the high-altitude regions and the Fergana Valley.

2.
Increasing precipitation during the accumulation season may result in increased maximum snow depth in the Ili Valley and the upper reaches of the Aksu and Chu Rivers, subsequently causing a longer snow cover duration.The snow cover duration was shortened in the upper reaches of the Kaidu River basin, high-altitude regions of the Aksu and Ili Rivers, northern ETS, and Fergana Valley due to the fact that the increasing temperatures delayed the snow onset day and caused the occurrence of an earlier snow end day.

3.
The earlier snow end day with increased precipitation during the accumulation season in the snowmelt-dominated river basins contributed to increased snowmelt and advanced the runoff peak to earlier in spring.The increased snowmelt will increase the risk of floods.The shifted runoff peak will introduce temporal mismatch between the water supply and the crop-growing season, increasing the difficulties of the regional reservoir and agriculture management and intensifying contradictions for the planning and utilization of water in neighboring countries.

Figure 1 .
Figure 1.(a) Location of the Tianshan Mountains (TS) and the distribution of the meteorological stations.(b) Annual precipitation over the TS from 1979 to 2016 (Global Precipitation Climatology Centre products (GPCC)).(c) Monthly distribution of precipitation and temperature in the TS during 1979-2016 (GPCC and Global Historical Climatology Network 2 and Climate Anomaly Monitoring System (GHCN_CAMS)).(d) Averaged snow depth time series over the TS during 1979-2016.

Figure 1 .
Figure 1.(a) Location of the Tianshan Mountains (TS) and the distribution of the meteorological stations.(b) Annual precipitation over the TS from 1979 to 2016 (Global Precipitation Climatology Centre products (GPCC)).(c) Monthly distribution of precipitation and temperature in the TS during 1979-2016 (GPCC and Global Historical Climatology Network 2 and Climate Anomaly Monitoring System (GHCN_CAMS)).(d) Averaged snow depth time series over the TS during 1979-2016.

2. 3 . 4 .
Contributions of D o and D e to D d The standardized z-score was used to quantify the contributions of D o and D e to D d anomalies.The standard score of a raw score x [75] can be defined as follows: z = x − µ σ (4) where µ is the mean of D o , D e , and D d over the TS from 1979 to 2016 and σ is the standard deviation of D o , D e , and D d in the corresponding period.By regressing the z-scores of D d against the z-scores of D o and D e , the contributions of D o and D e to D d anomalies can be derived from the regression coefficients.
Phenology over the TS Spatial distributions of D o , D e , D d , and SD max over the TS during 1979-2016 are shown in Figure 2. The mean D o , D e , D d , and SD max across the TS during 1979-2016 were the 78.2nd day of the hydrological year (DOY; i.e., 17 November), 222.4thDOY (i.e., 10 April), 146.2 days, and 16.1 cm, respectively.Snow cover appeared earlier in the upper reaches of the Syr Darya, Aksu, Kaidu, and Ili Rivers and later in the Ili Valley, Fergana Valley, and edge of the TS (Figure 2a).On the other hand, the spatial distribution of D e was opposite to the pattern of D o (Figure 2b): later D e corresponded to earlier D o , and vice versa.Combining the patterns of D o and D e , the spatial distribution of D d displayed a longer period in the upper reaches of the Aksu and the Kaidu Rivers, but a shorter period in the Ili Valley, Fergana Valley, and edge of the TS (Figure2c).The large values in SD max occurred in the upper reaches of the Syr Darya, Aksu, and Kaidu Rivers, as well as the Ili Valley, while the small values occurred in the Fergana Valley, southern edge of the TS, eastern ETS, and lower part of the Ili Valley (Figure2d).earlierDo, and vice versa.Combining the patterns of Do and De, the spatial distribution of Dd displayed a longer period in the upper reaches of the Aksu and the Kaidu Rivers, but a shorter period in the Ili Valley, Fergana Valley, and edge of the TS (Figure2c).The large values in SDmax occurred in the upper reaches of the Syr Darya, Aksu, and Kaidu Rivers, as well as the Ili Valley, while the small values occurred in the Fergana Valley, southern edge of the TS, eastern ETS, and lower part of the Ili Valley (Figure2d).

Figure 2 .
Figure 2. Averaged snow phenology over the period of 1979-2016: (a) snow onset day D o , (b) snow end day D e , (c) snow cover duration D d , and (d) maximum snow depth SD max .

3. 2 .
Changes in Snow Phenology over the TS D e and D d in the TS experienced a significant downward trend at the 0.05 significance level (Figure 4a).Spatial distributions of temporal trends in mean D o , D e , D d , and SD max over the TS during 1979-2016 are displayed in Figure 3. D o was earlier in the Ili Valley and the upper reaches of the Chu and Syr Darya Rivers but delayed on the southern slope of the TS; in particular, it was significantly later in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the ETS, and the Fergana Valley (Figure 3a).Compared with the changes in D o , D e was earlier across almost all of the TS; in particular, it was significantly earlier in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the ETS, and the Fergana Valley but delayed in sporadic regions (Figure 3b).Combined with the variations of D o and D e , D d decreased in the high-altitude regions of the Aksu, Kaidu, and Ili Rivers, the northern slope of the ETS, and the Fergana Valley but increased in the upper reaches of the Chu River and the Ili Valley (Figure Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 16

Figure 2 .
Figure 2. Averaged snow phenology over the period of 1979-2016: (a) snow onset day Do, (b) snow end day De, (c) snow cover duration Dd, and (d) maximum snow depth SDmax.

Figure 3 .
Figure 3. Temporal trends in snow phenology: (a) snow onset day Do, (b) snow end day De, (c) snow cover duration Dd, and (d) maximum snow depth SDmax.The black dots in Figure 3 indicate that the trends were significant (significance level at 0.05).
The contributions of Do and De to Dd were estimated by their normalized values over the TS from 1979 to 2016 (Figure 4).Changes in De accounted for 52.8% of changes in Dd over the TS during 1979-2016 (Figure 4b), while changes in Do accounted for 47.2% of the changes in Dd.In addition, the correlation analysis between Do, De and Dd indicated that both Do and De had a significant (at 0.01 significance level) relationship with Dd (Figure 4c,d) during the period of 1979-2016.However, this relationship between De and Dd (R 2 = 0.68) was more significant than that between Do and Dd (R 2 = 0.39).The above results revealed that changes in De caused the most variations in the snow season over the TS during 1979-2016.

Figure 3 . 16 Figure 4 .
Figure 3. Temporal trends in snow phenology: (a) snow onset day D o , (b) snow end day D e , (c) snow cover duration D d , and (d) maximum snow depth SD max .The black dots in Figure 3 indicate that the trends were significant (significance level at 0.05).3.3.Contributions of D o and D e to D d The contributions of D o and D e to D d were estimated by their normalized values over the TS from 1979 to 2016 (Figure 4).Changes in D e accounted for 52.8% of changes in D d over the TS during 1979-2016 (Figure 4b), while changes in D o accounted for 47.2% of the changes in D d .

Figure 4 .
Figure 4. (a) Annual variations of snow onset day D o , snow end day D e , and snow cover durations D d from 1979 to 2016 over the TS; (b) Contributions from D o and D e to D d across the TS during 1979-2016 (normalized values); (c) Relationship between D o and D d ; (d) Relationship between D e and D d .R 2 is the coefficient of determination.

Figure 5 .
Figure 5. Spatial distribution of changes in temperature and precipitation over the TS from 1979 to 2016.The spatial pattern of (a) days with a mean temperature of <0 °C and (b) changes from meteorological observations.The variations of (c) annual mean Ta, (d) Tm, (e) Pa, and (f) Pm.The black crosses and dots in Figure 5 indicate that the trends were significant (significance level at 0.05).

Figure 5 .
Figure 5. Spatial distribution of changes in temperature and precipitation over the TS from 1979 to 2016.The spatial pattern of (a) days with a mean temperature of <0 • C and (b) changes from meteorological observations.The variations of (c) annual mean T a , (d) T m , (e) P a , and (f) P m .The black crosses and dots in Figure 5 indicate that the trends were significant (significance level at 0.05).

Figure 6 .
Figure 6.(a) Sensitivity of Do to Ta (days °C−1 ), (b) sensitivity of Do to Pa (days mm −1 ), (c) sensitivity of SDmax to Ta (cm °C−1 ), (d) sensitivity of SDmax to Pa (cm mm −1 ), (e) sensitivity of De to Tm (days °C−1 ), and (f) sensitivity of De to SDmax (days cm −1 ).The black dots in Figure 6 indicate that the correlation was significant (significance level at 0.05).

Figure 6 .
Figure 6.(a) Sensitivity of D o to T a (days • C −1 ), (b) sensitivity of D o to P a (days mm −1 ), (c) sensitivity of SD max to T a (cm • C −1 ), (d) sensitivity of SD max to P a (cm mm −1 ), (e) sensitivity of D e to T m (days • C −1 ), and (f) sensitivity of D e to SD max (days cm −1 ).The black dots in Figure 6 indicate that the correlation was significant (significance level at 0.05).