Impact of Snowpack on the Land Surface Phenology in the Tianshan Mountains, Central Asia

: The accumulation and ablation processes of seasonal snow significantly affect the land surface phenology in a mountainous ecosystem. However, the ability of snow to regulate the alpine land surface phenology in the arid regions is not well described in the context of climate change. The impact of snowpack changes on land surface phenology and its driving factors were investigated in the Tianshan Mountains using the land surface phenology metrics derived from satellited products and a snow dataset from downscaled regional climate model simulations covering the period from 1983 to 2015. The results demonstrated that the annual mean start of growing season (SOS) and length of growing season (LOS) experienced a significant ( p < 0.05) decrease and increase with a rate of − 2.45 days/decade and 2.98 days/decade, respectively. The significantly advanced SOS and increased LOS were mainly seen in the Western Tianshan Mountains and Ili Valley regions with elevations from 2500 to 3500 m a.s.l and below 3000 m a.s.l, respectively. During the early spring, the significant decline in snow cover fraction (SCF) could advance the SOS. In contrast, snowmelt amount and annual maximum snow water equivalent (SWE) have an almost equally substantial positive correlation with annual maximum vegetation greenness. In particular, the SOS of grassland was the most sensitive to variations of snow cover fraction during early spring than that of other vegetation types, and their strong relationship was mainly located at elevations from 1500 to 2500 m a.s.l. Its greenness was significantly controlled by the annual maximum snow water equivalent in all elevation bands. Both decreased SCF and increased temperature in the early spring caused a significant advance of the SOS, consequently prolonging the LOS. Meanwhile, more SWE and snowmelt amount could significantly promote vegetation greenness by regulating the soil moisture. The results can improve the understanding of the snow ecosystem services in the alpine regions under climate change.


Introduction
Vegetation has a critical effect on Earth's ecosystem, water-heat exchange, and carbon cycles [1][2][3]. The variability of land surface phenology (the start of growing season (SOS), the end of growing season (EOS), and the length of growing season (LOS)) are controlled by climatic conditions and partly human activities [4,5]. The alpine environment is vulnerable to climate change, and is facing the warming faster than low-elevation regions [6,7]. A warming climate may increase the amount of energy received by seasonally snowcovered areas [8], benefiting the vegetation greening in the high latitude and altitude surface phenology in the entire TS. The results of this study will improve our understanding of snow ecosystem service.

Study Area
The TS are the largest mountain system of Central Asia (Figure 1a), spanning from Xinjiang (China), southeastern Kazakhstan and Uzbekistan to Kyrgyzstan. They have a length of approximately 2500 km and measure a width of 250-350 km, and the average altitude is close to 4000 m above sea level (a.s.l). The TS are geographically divided into four parts: (1) the Western Tianshan Mountains (WTS), (2) the Northern Tianshan Mountains (NTS), (3) the Eastern Tianshan Mountains (ETS), and (4) the Southern Tianshan Mountains (STS) [44]. Affected by the westerlies and the terrain obstruction, the TS exhibit a distinct continentality gradient with an increasing temperature and precipitation from southeast to northwest. The annual total precipitation and mean temperature amount to 329.3 mm and 4.6 • C, respectively [35]. The annual total precipitation exceeds 500 mm along the northern slope of the TS and the WTS, but amounts are below 100 mm around the Tarim Basin [45]. The precipitation in the ETS and NTS mainly occurs in spring and early summer, which is earlier than that in the STS (summer) but later than the WTS (later winter to early spring) [44]. The abundant precipitation in the mountainous regions forms a "wet island" in the arid land and benefits to shape the richly glaciated and snow-covered areas [44,46], which provides the freshwater for the mountain-oasis-desert ecosystem [33]. In addition, the annual maximum snow mass storage and mean snowfall in the TS reaches about 97.85 gigatons and 84.53 mm, respectively [37,47]. The TS exhibit a rich vegetation diversity ecosystem (Figure 1b), and grassland is the dominant land cover type (Figure 1c). related snow metrics (snow cover fraction (SCF), annual maximum snow water equivalent (SWEmax), and snowmelt amount); and (2) to explore the relationships between snow changes and vegetation dynamic. This is the first attempt to reveal how the snow changes modulate land surface phenology in the entire TS. The results of this study will improve our understanding of snow ecosystem service.

Study Area
The TS are the largest mountain system of Central Asia (Figure 1a), spanning from Xinjiang (China), southeastern Kazakhstan and Uzbekistan to Kyrgyzstan. They have a length of approximately 2500 km and measure a width of 250-350 km, and the average altitude is close to 4000 m above sea level (a.s.l). The TS are geographically divided into four parts: (1) the Western Tianshan Mountains (WTS), (2) the Northern Tianshan Mountains (NTS), (3) the Eastern Tianshan Mountains (ETS), and (4) the Southern Tianshan Mountains (STS) [44]. Affected by the westerlies and the terrain obstruction, the TS exhibit a distinct continentality gradient with an increasing temperature and precipitation from southeast to northwest. The annual total precipitation and mean temperature amount to 329.3 mm and 4.6 °C, respectively [35]. The annual total precipitation exceeds 500 mm along the northern slope of the TS and the WTS, but amounts are below 100 mm around the Tarim Basin [45]. The precipitation in the ETS and NTS mainly occurs in spring and early summer, which is earlier than that in the STS (summer) but later than the WTS (later winter to early spring) [44]. The abundant precipitation in the mountainous regions forms a "wet island" in the arid land and benefits to shape the richly glaciated and snow-covered areas [44,46], which provides the freshwater for the mountain-oasis-desert ecosystem [33]. In addition, the annual maximum snow mass storage and mean snowfall in the TS reaches about 97.85 gigatons and 84.53 mm, respectively [37,47]. The TS exhibit a rich vegetation diversity ecosystem (Figure 1b), and grassland is the dominant land cover type (Figure 1c).    control process before its release [48]. GIMMS NDVI 3g has been widely applied for the long-time trend analysis of vegetation biomass and its phenology dynamics in the Tibet Plateau surrounding areas due to its good temporal consistency and fine quality [48,49]. The land cover type was collected from the European Space Agency (ESA) Climate Change Initiative (CCI) project (https://www.esa-landcover-cci.org/, accessed on 20 February 2020), which consists of a series of global land cover maps from 1992 to the present and the overall accuracy exceeds 75.4% [50], and it also exhibited highest overall accuracy in China compared with other six coarse-resolution land cover datasets [51]. In addition, the vegetation parameter has a significant effect on the land surface water-heat process simulation. Thus, to keep the consistency of land cover in the previous land surface simulation and later land surface phenology analysis, the basis land cover was derived from the WRF basis land over (Figure 1b, CCI-LC 2000 map) for further relationship analysis about the phenology from different land covers and the output from WRF model. [37].

Climate Data
The daily gridded SCF, SWE, surface air temperature (T2), surface albedo, snowmelt amount, and soil moisture content (SMOIS) with a 9 km spatial resolution from 1982 to 2015 were obtained from the outputs of WRF-AWR 4.01 [52] coupled with Noah-MP model [53] one-way double-nested dynamic downscaling simulation, which was forced by 6 h interval ERA5 reanalysis from October 1982 to September 2018. Compared with in situ and remote sensing observations, the snow characteristics, precipitation, and temperature from WRF downscaled simulation showed high accuracy, especially the correlation coefficient of daily snow depth and SWE overall exceeded 0.85 [37]. The detailed physical parameterizations configuration and set up of WRF simulation was shown in [37]. Since the annual maximum snow mass storage in the sub-regions of TS occur in February, March, and April (FMA) [37], the snow-related climate metrics during February, March, and April (SCF FMA , Snowmelt FMA , T2 FMA , SMOIS FMA , and Albedo FMA ) were chosen to analyze the snow effect on early spring land surface phenology.

Land Surface Phenology Method
The thresholds methods and change detection methods were widely applied to estimate land surface phenology metrics [54]. Previous studies demonstrated that the dynamic threshold method is simple and could effectively derive land surface phenology metrics in Central Asia and Tibet Plateau [41]. Considering the validation data are still scarce in the TS, the NDVI ratio of 20% has a good accuracy in SOS estimation in surrounding regions of the TS [41] and the EOS was determined by NDVI ratio of 60% in Tibet Plateau [55,56], which has a similar land cover and alpine environment. Therefore, the dynamic thresholds method (NDVI ratio of 20% and 60%) was used to determine the SOS and EOS, respectively [57,58]. The 15-day NDVI data were interpolated to a daily time scale, and the Savitzky Golay filter was used to smooth out the noise of the NDVI time series [59]. The NDVI ratio could be calculated as: where NDVI t represents the NDVI value at the time step t (Day of year, DOY), and NDVI min and NDVI max demonstrate the minimum and maximum NDVI values in a year, respectively. The length of LOS is equal to EOS minus SOS. In addition, the annual mean NDVI value (1983-2015) < 0.1 was used to exclude the background values (bare land) from sparse vegetation ( Figure 1b) [60]. The GIMMS NDVI3g product was regridded to a 9 km spatial resolution in order to have the same spatial resolution with WRF outputs before land surface phenology calculation.

Statistical Analyses
The trend magnitudes of land surface phenology and climate factors, and their significant levels were estimated by the Sen's method [61]  test [62,63], respectively. In addition, the relationships between the land surface phenology and climate factors were investigated by the Pearson Correlation coefficient (r).

Spatiotemporal Distribution of Land Surface Phenology
The annual mean SOS, LOS, and NDVI max were 106.30 (Day of year, DOY), 162.24 days, and 0.41 in the TS from 1983 to 2015, respectively. The SOS showed an increase from the low-elevation regions to high-elevation regions (Figure 2a), in which the values mainly varied from 70 to 135. In contrast, the LOS exhibited an opposite pattern (Figure 2c), in which a longer growing season was found in the low-elevation region (elevation < 3000 m a.s.l) compared with the high-elevation regions (elevation ≥ 3000 m a.s.l). The high NDVI max values prevailed in the WTS and the low-elevation regions of the Ili Valley, and the corresponding values exceeded 0.6 ( Figure 2e). sparse vegetation (Figure 1b) [60]. The GIMMS NDVI3g product was regridded to a 9 km spatial resolution in order to have the same spatial resolution with WRF outputs before land surface phenology calculation.

Statistical Analyses
The trend magnitudes of land surface phenology and climate factors, and their significant levels were estimated by the Sen's method [61] and Mann-Kendall (M-K) trend test [62,63], respectively. In addition, the relationships between the land surface phenology and climate factors were investigated by the Pearson Correlation coefficient (r).

Spatiotemporal Distribution of Land Surface Phenology
The annual mean SOS, LOS, and NDVImax were 106. 30 (Figure 3). The LOS experienced a significant increase (p < 0.05) with a rate of 2.98 days/decade. Similarly, most significantly increasing LOS values were noticed in the lowelevation regions of the Ili valley and WTS (Figure 2d), a few significantly decreasing values were found in the low-elevation regions of the northern slope the TS. The NDVImax showed a slight increase with a rate of 0.004/decade. The significantly increasing NDVImax values were mainly noticed in the WTS and the intersection of NTS, STS, and ETS ( Figure  2f).

Spatiotemporal Variability of SWE, SCF, and Snowmelt
The variations of SWE, SCF, and snowmelt in the TS are shown in Figure 4. The SWEmax showed a slight increase) with a rate of 2.31 mm/decade in the TS. It exhibited a consistent augmentation in most regions of TS except for the STS (Figure 4a). The significantly increasing values of SWEmax were mainly observed in the ETS. In contrast, the SCFFMA experienced a slight decline with a rate of −1.23%/decade. Almost all regions of the TS showed a consistent decrease (Figure 4b), and a significant decrease in SCFFMA was seen in the STS. Conversely, the SnowmeltFMA exhibited a significant augmentation (p < 0.05) with a rate of 9.41 mm/decade. The significantly increased values were noticed in the ETS and high-elevation regions of the Ili valley and WTS (Figure 4c), but the decreased values prevailed in the STS.

Spatiotemporal Variability of SWE, SCF, and Snowmelt
The variations of SWE, SCF, and snowmelt in the TS are shown in Figure 4. The SWE max showed a slight increase) with a rate of 2.31 mm/decade in the TS. It exhibited a consistent augmentation in most regions of TS except for the STS (Figure 4a). The significantly increasing values of SWE max were mainly observed in the ETS. In contrast, the SCF FMA experienced a slight decline with a rate of −1.23%/decade. Almost all regions of the TS showed a consistent decrease (Figure 4b), and a significant decrease in SCF FMA was seen in the STS. Conversely, the Snowmelt FMA exhibited a significant augmentation   Figure 4 indicate the trend at a 0.05 significant level.

Response of Land Surface Phenology to Snow
The relationships between land surface phenology metrics and snow during 1983-2015 are shown in Figure 5. The SCFFMA had a significant positive effect (r = 0.64, p < 0.01) on SOS in the entire TS, especially in the WTS and low-elevation regions of the Ili Valley (Figure 5a). In contrast, the significant positive areas influenced by the SWEmax (Figure 5b)   Figure 4 indicate the trend at a 0.05 significant level.

Response of Land Surface Phenology to Snow
The relationships between land surface phenology metrics and snow during 1983-2015 are shown in Figure 5. The SCFFMA had a significant positive effect (r = 0.64, p < 0.01) on SOS in the entire TS, especially in the WTS and low-elevation regions of the Ili Valley (Figure 5a). In contrast, the significant positive areas influenced by the SWEmax (Figure 5b) and SnowmeltFMA (Figure 5c) were smaller and concentrated in the WTS and part lowelevation regions of the intersection of WTS and NTS, respectively. Conversely, the Snow-meltFMA in the high-elevation regions had a significant negative effect on SOS, particularly in high-elevation regions of the Ili Valley. Compared with a mixed pattern of the relationships between LOS and SCFFMA or SnowmeltFMA (Figure 5d,f) Figure 5 means no data, and the black dots in Figure 5 indicate the correlation coefficient at a 0.05 significant level.

Influence of Snow on Land Surface Phenology
The spring land surface phenology is strongly associated with snow cover dynamics in the alpine regions [17,64]. Previous studies suggested that the spring land surface phenology is principally controlled by air temperature [12,65]. Changes in SOS were also primarily influenced by the air temperature during the early spring in the TS (Figure 6), which also confirmed the above views. Meanwhile, the dramatic warming in spring could cause significantly earlier snowmelt dates and the reduction of the SCF FMA [66], which was recognized as the most ecologically relevant evidence for spring phenology [28,67]. Variation of SCF FMA has a significant positive effect on the SOS. A significant advanced snow end date and shortened snow cover duration have been reported in the TS due to a warming climate [35,36], which could result in a decrease in Albedo FMA and feedback to the atmosphere, aggravating warming in spring [68]. Therefore, a decline in SCF FMA could trigger a widely earlier SOS in the TS due to increases in soil temperature and air temperature during the early spring, especially in the low-elevation regions of the Ili Valley and WTS (Figure 5a), which might contribute to the vegetation germination [69].
ation of SCFFMA has a significant positive effect on the SOS. A significant advanced snow end date and shortened snow cover duration have been reported in the TS due to a warming climate [35,36], which could result in a decrease in AlbedoFMA and feedback to the atmosphere, aggravating warming in spring [68]. Therefore, a decline in SCFFMA could trigger a widely earlier SOS in the TS due to increases in soil temperature and air temperature during the early spring, especially in the low-elevation regions of the Ili Valley and WTS (Figure 5a), which might contribute to the vegetation germination [69]. Compared with SCFFMA, SWEmax, and SnowmeltFMA could significantly regulate the annual maximum vegetation greenness in the TS (Figure 6), which could be partly attributed to the interannual variability of temperature and soil moisture during the early growing season [28]. The vegetation growth is limited by water in the low-elevation regions, whereas in alpine and high-latitude regions it is temperature limited [65]. The degree of dependency on snow for vegetation growth is determined by the legacy effect of snow on soil moisture and the demand for moisture conditions during the growing season [28]. Vegetation greenness is greatly sensitive and vulnerable to the water conditions in the TS [38]. A high SWEmax value is always associated with a large snow depth and longer snow cover duration, raising winter soil temperatures and reducing the risk of soil exposure to freezing temperatures [17]. Meanwhile, more SnowmeltFMA could increase soil moisture contents and promote vegetation greening in the subsequent growing season until the summer [27]. Hence, a widespread significant increase in SnowmeltFMA is a factor Compared with SCF FMA , SWE max , and Snowmelt FMA could significantly regulate the annual maximum vegetation greenness in the TS (Figure 6), which could be partly attributed to the interannual variability of temperature and soil moisture during the early growing season [28]. The vegetation growth is limited by water in the low-elevation regions, whereas in alpine and high-latitude regions it is temperature limited [65]. The degree of dependency on snow for vegetation growth is determined by the legacy effect of snow on soil moisture and the demand for moisture conditions during the growing season [28]. Vegetation greenness is greatly sensitive and vulnerable to the water conditions in the TS [38]. A high SWE max value is always associated with a large snow depth and longer snow cover duration, raising winter soil temperatures and reducing the risk of soil exposure to freezing temperatures [17]. Meanwhile, more Snowmelt FMA could increase soil moisture contents and promote vegetation greening in the subsequent growing season until the summer [27]. Hence, a widespread significant increase in Snowmelt FMA is a factor that contributed to a rise in NDVI max . The wetting and warming tendency might also benefit vegetation greenness in the TS.

Responses of Different Vegetation Types to Snow
The responses of different vegetation types' phenology to snow dynamics exhibited a distinct difference in the TS, which was consistent with results obtained from other regions such as the Tibetan Plateau and high-elevation regions of the Northern Hemisphere [65,69]. The SOS of grassland was more sensitive to SCF FMA variations than other vegetation types (Figure 7a). A similar result was also demonstrated in the low-elevation regions of the Ili Valley based on the MODIS snow and vegetation products [42]. Like the Alps [16], the SOS of grassland was influenced primarily by a decline in snow cover, secondarily by snowmelt amount, while the NDVI max was equally affected by SCF FMA and SWE max . Its SOS had a significant positive relationship with SCF FMA , but a negative relationship with Snowmelt FMA , indicating that snow plays a vital role in grassland growth onset. This strong positive relationship was mainly located at elevations from 1500 to 3000 m a.s.l (Figure 8a). Notably, although the SCF FMA was not significant with LOS of grassland, it could significantly affect the LOS below 3000 m a.s.l (Figure 8b). Moreover, Snowmelt FMA significantly affected the LOS of grassland (Figure 7b), and its greenness was significantly correlated with the SWE max amount in all elevation bands (Figures 7c and 8c). Overall, the Remote Sens. 2022, 14, 3462 9 of 14 relations between sparse vegetation and snow variations showed a similar pattern with grassland. The effect of SCF FMA on LOS of sparse vegetation was more significant than grassland except for elevation below 1000 m a.s.l. Some herb vegetation germination and growth, such as ephemeral plants, heavily depend on the snow conditions around the TS [29]. Large SWE max and Snowmelt FMA suggested that snow provided excellent thermal protection and essential moisture supply for shorter plant types [42,56]. Generally, earlier snowmelt with more snowmelt amount could significantly advance the SOS and prolong LOS for grassland, promoting the increase in NDVI max [16,27]. In contrast, although more Snowmelt FMA might also replenish the soil moisture, benefit the greenness, and prolong the LOS, earlier snowmelt might restrain the early growth and significantly reduce the NDVI max for shrubland and forest due to the exposure to the risk of lower soil temperature and wind damage [17,70,71]. It is noted that snow amount could not significantly promote the increase in NDVI max for forest compared with SCF FMA , because the air temperatures are more closely correlated with forest than with other vegetation types (Figure 7c) [72]. The main effect of snow on forest growth is to prevent frost damage [73]. The largest precipitation zone in the TS is near the forest line [45], and precipitation is the most important source of soil moisture for later forest growth. In addition, the terrain features (elevation, aspect, and slope) strongly affect alpine species richness, productivity, and nutrient throughout altering the snow metrics dynamics and thermal regime [17,71].

Limitations and Uncertainties
The performance of the snow dataset generated by a regional model depends on the forcing data accuracy, physical schemes and model parameterizations, and spatial resolution [43,74]. A cold bias, precipitation overestimation, and lack of data assimilation caused a large uncertainty of snow simulation in the high-elevation regions, particularly in the snow ablation period [37,75]. In addition, the spatial heterogeneity of snow could not be revealed in the complex terrains using a relatively coarse resolution (9 km) [76,77]. Similarly, the spatial resolution of GIMMS NDVI 3g is appropriate for landscape-scale research, but maybe not be sufficient for specific species or plant communities at a finer scale [17,78]. The synergistic analysis from the higher remote sensing data (such as MODIS, Landsat, Sentinel-2, Gaofen, and digital cameras) may provide a feasible way to investigate relationships between snow cover and vegetation and their interactions under complex terrain (i.e., ≤500 m) and micro-climate conditions [79][80][81]. The cloud cover, temporal resolution, and other production processes may also induce the NDVI product uncertainties [49].

Limitations and Uncertainties
The performance of the snow dataset generated by a regional model depends on the forcing data accuracy, physical schemes and model parameterizations, and spatial resolution [43,74]. A cold bias, precipitation overestimation, and lack of data assimilation caused a large uncertainty of snow simulation in the high-elevation regions, particularly in the snow ablation period [37,75]. In addition, the spatial heterogeneity of snow could not be revealed in the complex terrains using a relatively coarse resolution (9 km) [76,77]. Similarly, the spatial resolution of GIMMS NDVI 3g is appropriate for landscape-scale research, but maybe not be sufficient for specific species or plant communities at a finer scale [17,78]. The synergistic analysis from the higher remote sensing data (such as MODIS, Landsat, Sentinel-2, Gaofen, and digital cameras) may provide a feasible way to investigate relationships between snow cover and vegetation and their interactions under The dynamic threshold is a simple and efficient method to derive the land surface phenology in existing studies [57,58], but the lack of in situ observations to validate its performance increased the uncertainty of land surface phenology calculation in this study. Land cover experienced a dramatic change with a rapid increase in crop area around the TS during the past decades [82]; vegetation types based on the fixed CCI-LC 2000 map might not accurately reflect the land cover changes and their phenology. For example, large proportions of crop area were converted from grassland or sparse vegetation around the TS [83], which may increase the greenness and advance the SOS. Meanwhile, previous studies reported that the expansion of irrigated areas could affect the mountainous climate by increasing precipitation and reducing the temperature, and then modifying the regional snow status [84,85]. Similarly, the glacier areas are significantly shrinking due to climate warming, which increased glacier/snow melt water and benefited vegetation greening [71].

Conclusions
This study investigated the variability of land surface phenology metrics and their responses to the snow dynamics in the Tianshan Mountains from 1983 to 2015 using the GIMMS NDVI product and a snow dataset from a regional climate model dynamic downscaling simulation. The main findings of this study include the following: 1 The annual mean start of growing season was 106.30 (day of the year) and exhibited a significant decrease trend (p < 0.05, −2.45 days/decade), which concentrated in the Ili Valley and Western Tianshan Mountains with an elevation from 2500 to 3500 m a.s.l. In contrast, the length of growing season was 162.24 days and exhibited a significant increase trend (p < 0.05, 2.98 days/decade), which was mainly seen in low-elevation regions (elevations below 3000 m a.s.l) of the Ili Valley and Western Tianshan Mountains. 2 Snow cover fraction during February-April has a significant positive effect on the start of growing season. In contrast, snowmelt amount during February-April and annual maximum snow water equivalent have an almost equally significant positive correlation with NDVI max . In particular, the start of growing season of grassland was the most sensitive to variations of snow cover fraction during February-April than that of other vegetation types, and their strong relationship was mainly located at elevations from 1500 to 3000 m a.s.l. In addition, its greenness was significantly affected by the annual maximum snow water equivalent in all elevation bands.

3
Both decreased snow cover fraction and increased temperature in the early spring period caused the significant advance of the start date of vegetation growing season, consequently prolonging the length of vegetation growing season. Large annual maximum snow water equivalent and more snowmelt amount could significantly promote the increase in NDVI max by consistently regulating the soil moisture.
The variability of large-scale atmospheric circulations, such as El Nino and Atlantic multi-decadal oscillation (AMO), could significantly impact regional temperature and precipitation [86,87], and their relations with snow and vegetation need to be investigated further. Notably, although the significantly advanced SOS was found in the Tianshan Mountains and snow variability plays a critical role, it is necessary to quantify the contribution of snow, temperature, and precipitation to the change in the early spring vegetation growth. Meanwhile, the relationships between the spring land surface phenology and elevation-dependent warming should be considered using a finer resolution product in the next work.

Conflicts of Interest:
The authors declare no conflict of interest.