Relationship between Winter Snow Cover Dynamics , Climate and Spring Grassland Vegetation Phenology in Inner Mongolia , China

The onset date of spring phenology (SOS) is regarded as a key parameter for understanding and modeling vegetation–climate interactions. Inner Mongolia has a typical temperate grassland vegetation ecosystem, and has a rich snow cover during winter. Due to climate change, the winter snow cover has undergone significant changes that will inevitably affect the vegetation growth. Therefore, improving our ability to accurately describe the responses of spring grassland vegetation phenology to winter snow cover dynamics would enhance our understanding of changes in terrestrial ecosystems due to their responses to climate changes. In this study, we quantified the spatial-temporal change of SOS by using the Advanced Very High Resolution Radiometer (AVHRR) derived Normalized Difference Vegetation Index (NDVI) from 1982 to 2015, and explored the relationships between winter snow cover, climate, and SOS across different grassland vegetation types. The results showed that the SOS advanced significantly at a rate of 0.3 days/year. Winter snow cover dynamics presented a significant positive correlation with the SOS, except for the start date of snow cover. Moreover, the relationship with the increasing temperature and precipitation showed a significant negative correlation, except that increasing Tmax (maximum air temperature) and Tavg (average air temperature) would lead a delay in SOS for desert steppe ecosystems. Sunshine hours and relative humidity showed a weaker correlation.


Introduction
Spring phenology is an important indicator of vegetation dynamics, and an important contributor to vegetation activity and ecosystem functions.The onset date of spring phenology (SOS) can affect the carbon cycle of terrestrial ecosystems [1][2][3], hydrologic processes [4], nutrient cycling [5], and land-atmosphere interactions [6].The spatial-temporal dynamics of SOS are sensitive to changes in the climate and natural environment [7,8].Moreover, the El Niño-Southern Oscillation (ENSO) system is recognized as one of the dominant sources of climate variability at interannual timescales, with important effects on temperature, precipitation, and drought [9].Many studies have shown that there is a significant correlation between ENSO events and vegetation phenology change [10,11].Therefore, accurately monitoring SOS and exploring its responses to climate change are becoming important aspects of global change research [12,13].
Inner Mongolia, which is located along the northern border of China, is ecologically typified by mid-latitudinal and semi-arid temperate grassland.Furthermore, it is regarded as the typical terrestrial transect by the International Geosphere-Biosphere Programme for global change research [14].Due to climate change and overgrazing, the vegetation ecosystem of Inner Mongolia is very fragile, while land desertification and degradation have become very severe [15].In recent decades, the Chinese government has implemented a 'return-farmland-to-grassland' strategy to reverse the land degradation in pastures [16,17].However, the ecosystems of this region are inherently fragile and unstable, making them especially vulnerable to climate change [18,19].
In general, warmer temperatures or increasing precipitation would lead to advanced SOS.Recent studies showed that increasing the maximum air temperature (Tmax) in winter had a more significant influence on spring phenology than minimum air temperature (Tmin) [20,21].In contrast, Shen et al. [22] found that changes in the SOS were more strongly associated with a warmer winter Tmin rather than Tmax in the Tibetan Plateau.However, the radiation characteristics and thermodynamic properties of winter snow cover could affect the land (heat/cold and content of moisture), atmosphere, biological characteristics, and ecosystem to various degrees [23][24][25][26].The insulative properties of snowpack also could protect the over-winter survival of vegetation from wind and low air temperature [27].Furthermore, winter snow cover could provide fresh water for vegetation activities.With global warming, the timing of snowmelt advance would prompt vegetation growth.Meanwhile, global warming would also result in the increasing number of low-temperature events.The higher Tmin could alleviate frost damage, while the higher Tmax might exacerbate drought effects over the cold and dry areas.In temperate regions, as the vegetation usually experiences a long cold winter season, it does not begin to grow again without a certain cumulative temperature and precipitation.Although previous studies have investigated the effect of snow cover on the alpine vegetation in Tibetan Plateau, the responses of temperate grassland vegetation to changes of winter snow cover still need to be further investigated, as the relationships between SOS and climate vary across different vegetation types [28,29].In addition, the effects of spring phenology and climate change (including Tmin, Tmax, and frozen soil thawing) during the snow cover melt on SOS need to be further investigated.
Due to its long data record and continuity since 1981, the Normalized Difference Vegetation Index (NDVI) derived from the Advanced Very High Resolution Radiometer (AVHRR) has been widely used for monitoring long-term vegetation changes at global, continental, and regional scales [30].In this study, we used the logistic fitting method to retrieve the vegetation spring phenology from 1982 to 2015, and explored the relationship between winter snow cover dynamics, climate, and vegetation spring phenology in Inner Mongolia for the first time.The main aims of this present study are described as follows: (1) Quantify the spatial heterogeneity and temporal trends in vegetation spring phenology for different temperate grassland vegetation types in Inner Mongolia.(2) Explore the underlying mechanisms related to the effects of winter snow cover dynamics (including snow depth (SD), snow cover duration (SCD), snow cover onset date (SCOD), and snow cover end date (SCED)) on the spring phenology.(3) Investigate the responses of the spring phenology to climate factors (including Tmin, Tmax, Tavg, precipitation (PRE), onset date of frozen soil thawing (GST), relative humidity (RHU), and sunshine duration (SSD)) during SCED and SOS.
The findings of this study will improve our understanding of grassland vegetation dynamics and its connections with the climate, enhance our ability to predict the magnitude and direction of spring phenology, and determine changes in the structure and functioning of the grassland ecosystems.

Study Area
Inner Mongolia is located in the northwestern part of China (37 1a).It covers a grassland area of 792,000 km 2 , accounting for 67% of the total area.The elevation of Inner Mongolia is shown in Figure 1c.The elevation increases from the northeast to the southwest, ranging from 88 to 3365 m.The study area is a transition zone between the arid and semi-arid northwest inland region, which is affected by Asian monsoons.The climate in Inner Mongolia can be characterized by a cold, dry winter and warm, rainy summer [31].The annual rainfall decreases gradually from the east (~400 mm) to the west (~200 mm) [32], with higher rain use efficiency than most other areas in the semi-arid and arid regions of China.The annual mean temperature ranges from −1.6  [31].The annual rainfall decreases gradually from the east (~400 mm) to the west (~200 mm) [32], with higher rain use efficiency than most other areas in the semi-arid and arid regions of China.The annual mean temperature ranges from −1.6 °C to 9.2 °C [33].The annual average SOS ranges from April to July [34,35].The distribution of grassland vegetation types is shown in Figure 1b, which was obtained from a 1:1,000,000 map of vegetation in China (Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC), http://www.resdc.cn).Grassland vegetation types shown in this map include three different vegetation types: meadow steppe, typical steppe, and desert steppe.

Data
The meteorological data, including daily precipitation and temperature, were acquired from the China Meteorological Data Sharing Service System (http://www.escience.gov.cn).A total of 33 meteorological stations were selected to represent the major steppes in the Inner Mongolia (Figure 1c).
The primary data sets were the surface reflectance datasets from the AVHRR.These data sets were developed by the Long Term Data Record (LTDR) project and distributed through NASA's Goddard Space Flight Center (available at http://ltdr.nascom.nasa.gov).The AVHRR archival reflectance products (AVH09C1-version4), with a 0.05° spatial resolution from 1982 to 2015, are generated based on atmospheric correction using Rayleigh scattering, ozone, water vapor, and aerosol corrections in addition to masked out clouds and heavy aerosols [36].

Preprocessing of NDVI Time Series Data
Before NDVI calculation and data generation, pre-processing (including radiation correction, atmospheric correction, and coordinate transformation) was needed [37,38].After this, NDVI was calculated from the albedo of the first channel (visible light) and the second channel (near infrared) of the AVHRR sensor.The equation is calculated as follows: where VIS is the albedo of the visible light (0.58-0.68 μm) and NIR is the albedo of the near infrared channel (0.725-1.10 μm).The distribution of grassland vegetation types is shown in Figure 1b, which was obtained from a 1:1,000,000 map of vegetation in China (Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC), http://www.resdc.cn).Grassland vegetation types shown in this map include three different vegetation types: meadow steppe, typical steppe, and desert steppe.

Data
The meteorological data, including daily precipitation and temperature, were acquired from the China Meteorological Data Sharing Service System (http://www.escience.gov.cn).A total of 33 meteorological stations were selected to represent the major steppes in the Inner Mongolia (Figure 1c).
The primary data sets were the surface reflectance datasets from the AVHRR.These data sets were developed by the Long Term Data Record (LTDR) project and distributed through NASA's Goddard Space Flight Center (available at http://ltdr.nascom.nasa.gov).The AVHRR archival reflectance products (AVH09C1-version4), with a 0.05 • spatial resolution from 1982 to 2015, are generated based on atmospheric correction using Rayleigh scattering, ozone, water vapor, and aerosol corrections in addition to masked out clouds and heavy aerosols [36].

Preprocessing of NDVI Time Series Data
Before NDVI calculation and data generation, pre-processing (including radiation correction, atmospheric correction, and coordinate transformation) was needed [37,38].After this, NDVI was calculated from the albedo of the first channel (visible light) and the second channel (near infrared) of the AVHRR sensor.The equation is calculated as follows: where VIS is the albedo of the visible light (0.58-0.68 µm) and NIR is the albedo of the near infrared channel (0.725-1.10 µm).
Satellite-derived NDVI data are usually interrupted by contamination from the atmosphere (e.g., aerosol and dust), clouds, and snow [39].Although the maximum-value composite technique has been used to reduce such contamination [40], there is still residual noise in these time series data, which would greatly hinder further application of the data.The Savitzky−Golay filter [41] was employed to fill the gaps and smooth the time series in a local moving window, which provided a simplified least-squares-fit convolution for smoothing and computing derivatives of a set of consecutive values.This filter performs best in most situations when smoothing different vegetation types using various satellite data [42].The Savitzky-Golay filter is computed as follows: where Y is the original NDVI value, Y * is the resultant NDVI value, and C i is the coefficient for the ith NDVI of the smooth window.N is the number of convoluting integers and is equal to the smoothing window size (2m + 1).The running index of the ordinate data in the original data table is represented by j, and m represents the half width of the smoothing window.

Retrieval of Spring Phenology
In this study, the logistic fitting method was used to identify the SOS.The logistic model was employed to fit the temporal variation of the filtered NDVI data to an annual growth phase.Figure 2 shows the schematic diagram for the retrieval of spring phenology using the logistic model.The logistic model function has the following form [43]: where y(t) is NDVI at time t, parameters a and b control the shape of the curve, c + d is the amplitude values, and d represents the initial background NDVI value.After this, the SOS was retrieved from pre-processed NDVI data using the inflection point-based method (RCmax), which has been described and tested in diverse biomes in different climatic zones, and validated for robustness by extracting phenological metrics [44].The RCmax method is based on the rate of change (RC) of the fitted NDVI curve, which is the derivative of the curvature calculated with Equation ( 4): where z = e a+bt and a, b and c have the same definitions as in Equation ( 3).After this, we defined the SOS as the time when the rate of curvature change achieved the first local maximum.
ISPRS Int.J. Geo-Inf.2019, 8, x FOR PEER REVIEW 4 of 13 Satellite-derived NDVI data are usually interrupted by contamination from the atmosphere (e.g., aerosol and dust), clouds, and snow [39].Although the maximum-value composite technique has been used to reduce such contamination [40], there is still residual noise in these time series data, which would greatly hinder further application of the data.The Savitzky−Golay filter [41] was employed to fill the gaps and smooth the time series in a local moving window, which provided a simplified least-squares-fit convolution for smoothing and computing derivatives of a set of consecutive values.This filter performs best in most situations when smoothing different vegetation types using various satellite data [42].The Savitzky-Golay filter is computed as follows: where Y is the original NDVI value,  * is the resultant NDVI value, and  is the coefficient for the ith NDVI of the smooth window.N is the number of convoluting integers and is equal to the smoothing window size (2m + 1).The running index of the ordinate data in the original data table is represented by j, and m represents the half width of the smoothing window.

Retrieval of Spring Phenology
In this study, the logistic fitting method was used to identify the SOS.The logistic model was employed to fit the temporal variation of the filtered NDVI data to an annual growth phase.Figure 2 shows the schematic diagram for the retrieval of spring phenology using the logistic model.The logistic model function has the following form [43]: where y(t) is NDVI at time t, parameters a and b control the shape of the curve, c + d is the amplitude values, and d represents the initial background NDVI value.After this, the SOS was retrieved from pre-processed NDVI data using the inflection point-based method (RCmax), which has been described and tested in diverse biomes in different climatic zones, and validated for robustness by extracting phenological metrics [44].The RCmax method is based on the rate of change (RC) of the fitted NDVI curve, which is the derivative of the curvature calculated with Equation (4): where z  and a, b and c have the same definitions as in Equation ( 3).After this, we defined the SOS as the time when the rate of curvature change achieved the first local maximum.

Statistical Analyses
The least-squares method [46] was used to simulate the temporal trends in SOS during the study period.The trends in SOS for each pixel can be calculated using Equation ( 5): where i is the serial number from 1 to n, n is the total number of years, and X i is the parameter (SOS) of ith year.Furthermore, slope < 0 indicates a decreasing trend, while slope > 0 indicates an increasing trend.In addition, the statistical significance of the slope was assessed at the 95% confidence level.
The study used the latitude and longitude of each meteorological station to extract the corresponding SOS values from 1982 to 2015.Pearson correlation coefficients were calculated between SOS and climatic factors (e.g., temperature and precipitation) at each station during the study period.The correlations represented the sensitivity of SOS to climatic factors [47].

The Performances of Satellite-Based SOS
The SOS derived from satellite data required ground observations to determine their reliability and accuracy.Thus, we compared the satellite-derived SOS with ground-observed data from 17 phenology field observation stations in the Inner Mongolia grassland during 2004-2013.The mean absolute error between the satellite-derived and observed spring phenology was 8.9 days, and the root mean square error was 10.2 days.The results indicate that the retrieved SOS data from the AVHRR NDVI are reliable.

Spatial Variation in Vegetation Spring Phenology
By using long-term satellite NDVI data and the least-squares method, we simulated the temporal trends in SOS during the study period of 1982-2015.SOS showed a wide dynamic range, ranging from day of year (DOY) 100 to 150.Early SOS regions were mainly located in the eastern part of Inner Mongolia.A small region in the southwest also had an early SOS, i.e., between DOY 80 and 100.Regions with late SOS, e.g., after DOY 130, were found at the north borders of the study area and the west regions.For other regions, SOS generally fell in the range between DOY 106 and 135, which is generally consistent with previous results [20].
Furthermore, the inter-annual trends of SOS were calculated over the entire grassland vegetation types (Figure 3b).Most regions showed a significant advanced trend in SOS, accounting for the total area of 52.7%.There was a delaying trend in 34.30% of the study area.However, regions with a delayed SOS were mainly located in the east and west.Overall, the mean SOS over the entire study area significantly advanced at a rate of 0.3 ± 0.74 days/year (p < 0.05) from 1982 to 2015.Over 34 years, all three grassland vegetation types exhibited an earlier spring season, as the SOS significantly advanced by 0.27 ± 0.47 days/year (p < 0.05) for temperate meadow steppe, 0.32 ± 0.65 days/year (p < 0.05) for temperate typical steppe, and 0.47 ± 1.48 days/year (p < 0.05) for temperate desert steppe.
The statistical results of the SOS and the trends of three grassland vegetation types are shown in Figure 4. On average, the SOS occurs on DOY 120 for meadow steppe and typical steppe, and on DOY 130 in the desert steppe.In general, the SOS of meadow steppe varies from DOY 110 to 130, and between DOY 100 and 140 for typical steppe.The desert steppe shows a wider range, with the SOS ranging from DOY 90 to 150.In addition, meadow steppe shows an advanced trend at rates between −0.5 and 0.5 days/years.The desert steppe has similar advanced trends, with that of typical steppe within 1 day/year.A majority of the trends in desert steppe vary from −1 to −0.5, which indicates that the onset date of spring season is experiencing a faster advance in temperate China.

Temporal Variations in Vegetation Spring Phenology
The changes in the SOS regime are not only spatially complex, but also exhibit distinct temporal variations.The trends over a long time period of several decades conceal short-term fluctuations in the SOS.To assess the impact of the analysis timeframe on our regression results, we analyzed trends with different window sizes (i.e., the length of the study period) ranging from one year to 34 years, across all start-year and window-size combinations, which were averaged over the different grassland vegetation types (Figure 5).

Temporal Variations in Vegetation Spring Phenology
The changes in the SOS regime are not only spatially complex, but also exhibit distinct temporal variations.The trends over a long time period of several decades conceal short-term fluctuations in the SOS.To assess the impact of the analysis timeframe on our regression results, we analyzed trends with different window sizes (i.e., the length of the study period) ranging from one year to 34 years, across all start-year and window-size combinations, which were averaged over the different grassland vegetation types (Figure 5).Within the study period, the regional average SOS shows that there were significant negative trends universally (i.e., earlier spring phenology), regardless of the timeframe over which the regression was performed.For meadow steppe, the average SOS fluctuated but generally advanced from 2006 with some short positive trends (1986-1988, 1998-2000, and 2004-2005).The temporal variation of the SOS change for typical steppe vegetation indicated that there were two obvious receding trends (1982)(1983)(1984)(1996)(1997)(1998)(1999).After this time period, there was a rapid advance, especially after 2010.The desert steppe displayed different patterns.The average SOS showed significant receding trends during these periods (from 1982 to 1987 and 1998 to 2004).This implies that while the three decade trend in SOS has been negative, the trend has become more variable in the past decade.These anomalies may be attributed to the strong effects of the EI Niño and La Niña events [48].Therefore, the vegetation phenology can act as a comprehensive indicator of climate changes.

The Relationship between Winter Snow Cover, Climate, and Vegetation Spring Phenology
To further explore the relationship between winter snow cover dynamics and vegetation spring phenology in the temperate grasslands of Inner Mongolia, we defined three parameters to investigate the snow cover dynamics in each hydrological year (1 August to 31 July of the next year) [46].
(1) Snow cover duration (SCD): the number of snow covered days in a hydrological year; (2) Onset date of snow cover accumulation (SCOD): the start dates (Julian day) of the stable snow cover duration in a hydrological year [49,50].It is defined as the first day when snow cover is initially observed to exist for at least five consecutive days; (3) End date of snow cover (SCED): the ending date (Julian day) of the snow cover.It is defined as the last day when snow cover is last observed to exist for at least five consecutive days.
Although some researchers have investigated the relationship between the SOS, temperature, and precipitation in different months, they did not take the impact of climate change during the SCED and SOS into account [51,52].In this study, we not only explored the relationship between the winter snow cover and vegetation spring phenology, but also considered the various climate factors (e.g., Within the study period, the regional average SOS shows that there were significant negative trends universally (i.e., earlier spring phenology), regardless of the timeframe over which the regression was performed.For meadow steppe, the average SOS fluctuated but generally advanced from 2006 with some short positive trends (1986-1988, 1998-2000, and 2004-2005).The temporal variation of the SOS change for typical steppe vegetation indicated that there were two obvious receding trends (1982)(1983)(1984)(1996)(1997)(1998)(1999).After this time period, there was a rapid advance, especially after 2010.The desert steppe displayed different patterns.The average SOS showed significant receding trends during these periods (from 1982 to 1987 and 1998 to 2004).This implies that while the three decade trend in SOS has been negative, the trend has become more variable in the past decade.These anomalies may be attributed to the strong effects of the EI Niño and La Niña events [48].Therefore, the vegetation phenology can act as a comprehensive indicator of climate changes.

The Relationship between Winter Snow Cover, Climate, and Vegetation Spring Phenology
To further explore the relationship between winter snow cover dynamics and vegetation spring phenology in the temperate grasslands of Inner Mongolia, we defined three parameters to investigate the snow cover dynamics in each hydrological year (1 August to 31 July of the next year) [46].
(1) Snow cover duration (SCD): the number of snow covered days in a hydrological year; (2) Onset date of snow cover accumulation (SCOD): the start dates (Julian day) of the stable snow cover duration in a hydrological year [49,50].It is defined as the first day when snow cover is initially observed to exist for at least five consecutive days; (3) End date of snow cover (SCED): the ending date (Julian day) of the snow cover.It is defined as the last day when snow cover is last observed to exist for at least five consecutive days.
Although some researchers have investigated the relationship between the SOS, temperature, and precipitation in different months, they did not take the impact of climate change during the SCED and SOS into account [51,52].In this study, we not only explored the relationship between the winter snow cover and vegetation spring phenology, but also considered the various climate factors (e.g., GST, PRE, Tavg, Tmax, Tmin, RHU, and SSD) during SCED and SOS.We found that the mean winter snow depth presented a declining trend (0.01 cm/year) in the desert steppe, and an increasing trend (>0.01 cm/year) in both the meadow steppe and typical steppe during the study period (Figure 6).Additionally, the mean snow cover duration presented a significant decrease (p < 0.1) at rates of −0.342 days/year, −0.069 days/year and −0.038 days/year for meadow steppe, typical steppe, and desert steppe, respectively.The SCOD presented a significant delay in both meadow steppe and desert steppe.The advanced rate of SCED in typical steppe (−0.225 days/year) was larger than that in meadow steppe (−0.08 days/year) and desert steppe (−0.007 days/year).
ISPRS Int.J. Geo-Inf.2019, 8, x FOR PEER REVIEW 8 of 13 GST, PRE, Tavg, Tmax, Tmin, RHU, and SSD) during SCED and SOS.We found that the mean winter snow depth presented a declining trend (0.01 cm/year) in the desert steppe, and an increasing trend (>0.01 cm/year) in both the meadow steppe and typical steppe during the study period (Figure 6).Additionally, the mean snow cover duration presented a significant decrease (p < 0.1) at rates of −0.342 days/year, −0.069 days/year and −0.038 days/year for meadow steppe, typical steppe, and desert steppe, respectively.The SCOD presented a significant delay in both meadow steppe and desert steppe.The advanced rate of SCED in typical steppe (−0.225 days/year) was larger than that in meadow steppe (−0.08 days/year) and desert steppe (−0.007 days/year).Focusing on the climate factors during SCED and SOS, we also found that the temperatures (Tavg, Tmax, and Tmin) followed an increasing trend from 1982 to 2015.The precipitation in desert steppe presented a non-significant trend.The SSD presented a negative trend, which might be due to the increase in atmospheric anthropogenic aerosols [53].The RHU presented a positive trend at rates of 0.123, 0.011 and 0.188 for meadow steppe, typical steppe and desert steppe, respectively.The start date of frozen soil thawing (the first day when ground soil temperature above 0 °C) presented a slight but non-significant increase (from −0.257 to −0.07 days/year) for the whole study area.
For the three vegetation types, changes in the SOS were closely correlated with temperature, precipitation and winter snow cover (Figure 7).Many climate factors showed a positive relationship between the SOS, winter snow cover parameters (e.g., SD, SCD and SCED), and GST for three grassland vegetation types.This indicates that deeper snow depth (or later snowmelt, later frozen soil thawing, or longer duration of snow cover) might lead to later SOS.In contrast, later snowfall might lead to an earlier SOS.The temperate meadow steppe and the typical steppe were more weakly associated with the start of snow cover compared with other snow cover parameters.By contrast, the SOS and SCOD had a significant correlation for the desert steppe.Moreover, the SOS was more sensitive to variations of temperature and precipitation during the SCED and SOS.In contrast, the SOS showed a positive correlation with Tmax and Tavg for the desert steppe.In addition, it was also reported that the increase in relative humidity had a weak effect on advancing the SOS for the three vegetation types in Inner Mongolia.Focusing on the climate factors during SCED and SOS, we also found that the temperatures (Tavg, Tmax, and Tmin) followed an increasing trend from 1982 to 2015.The precipitation in desert steppe presented a non-significant trend.The SSD presented a negative trend, which might be due to the increase in atmospheric anthropogenic aerosols [53].The RHU presented a positive trend at rates of 0.123, 0.011 and 0.188 for meadow steppe, typical steppe and desert steppe, respectively.The start date of frozen soil thawing (the first day when ground soil temperature above 0 • C) presented a slight but non-significant increase (from −0.257 to −0.07 days/year) for the whole study area.
For the three vegetation types, changes in the SOS were closely correlated with temperature, precipitation and winter snow cover (Figure 7).Many climate factors showed a positive relationship between the SOS, winter snow cover parameters (e.g., SD, SCD and SCED), and GST for three grassland vegetation types.This indicates that deeper snow depth (or later snowmelt, later frozen soil thawing, or longer duration of snow cover) might lead to later SOS.In contrast, later snowfall might lead to an earlier SOS.The temperate meadow steppe and the typical steppe were more weakly associated with the start of snow cover compared with other snow cover parameters.By contrast, the SOS and SCOD had a significant correlation for the desert steppe.Moreover, the SOS was more sensitive to variations of temperature and precipitation during the SCED and SOS.In contrast, the SOS showed a positive correlation with Tmax and Tavg for the desert steppe.In addition, it was also reported that the increase in relative humidity had a weak effect on advancing the SOS for the three vegetation types in Inner Mongolia.

Discussion
Previous studies have investigated the spatial-temporal change of annual SOS and its responses to temperature and precipitation in Inner Mongolia.As one of the three most stable snow regions in China, the snow depth increases in winter and decreases in spring over the Inner Mongolia region [54,55].However, the impacts of winter snow cover dynamics on spring phenology in the Inner Mongolia grassland are still unknown.This might be because previous studies [20,22] have paid more attention to the influence of the temperature (Tmax, Tavg, and Tmin) and precipitation in different months and seasons on the vegetation phenology.Although many studies have explored the impacts of winter snow cover on alpine vegetation in the Tibetan Plateau [56,57], the impacts of winter snow cover on temperate grassland vegetation are thought to be similar to those on alpine vegetation.As the relationships between the SOS and winter snow cover might vary across different vegetation types, the possible effects of winter snow cover across temperate grassland vegetation in Inner Mongolia should be further investigated.
It is worth noting that the SOS changes were related not only to the temperature and precipitation, but also to other factors (e.g., winter snow cover, GST, SSD, and RHU).Winter snow cover could add soil moisture, arrest soil heat exchange, and have crucial effects on soil heat and moisture preservation [58,59], which would protect the over-winter survival of vegetation from wind and low air temperatures.Thus, it could influence the exchange of gas, moisture, and dissolved matter.Furthermore, it could affect the activity of the microbes in the soil as well as the transformation of soil organic matter and nutrients.In general, deeper snow depth may lead to a later SOS.The reason might be that a deeper snow depth would have a longer SCD and later SCED, which would not be beneficial for soil respiration in winter and the accumulation of heat needed to unfold leaves in spring [60 -62].Earlier snow cover could shield vegetation from low temperatures by keeping the temperature at a favorable level and protecting the activity of the soil microorganisms during the winter, which could ensure positive vegetation growth during the growing season in the following year.Moreover, the spring phenology was more strongly and negatively associated with temperature and precipitation during the SCED and SOS, which indicates that the SOS will advance with warmer temperatures or more precipitation.In contrast, the SOS showed a significant positive correlation with Tmax and Tavg for desert steppe.Increased temperatures (Tmax and Tavg) would result in the loss of copious amounts of water by evaporation [20] and exacerbate the effects of drought in dry regions, which would delay the SOS for desert steppe [63,64].Increasing Tmin could reduce frost damage to stimulate vegetation growth.A suitable sunshine environment, which plays a role equally important as precipitation and temperature, is beneficial for photosynthesis and respiration.Meanwhile, excess SSD may also contribute to the aggravation of drought through evapotranspiration.The relative humidity exhibited a negative relationship with the SOS, and its change was usually associated with the temperature, sunshine, wind, and precipitation.

Discussion
Previous studies have investigated the spatial-temporal change of annual SOS and its responses to temperature and precipitation in Inner Mongolia.As one of the three most stable snow regions in China, the snow depth increases in winter and decreases in spring over the Inner Mongolia region [54,55].However, the impacts of winter snow cover dynamics on spring phenology in the Inner Mongolia grassland are still unknown.This might be because previous studies [20,22] have paid more attention to the influence of the temperature (Tmax, Tavg, and Tmin) and precipitation in different months and seasons on the vegetation phenology.Although many studies have explored the impacts of winter snow cover on alpine vegetation in the Tibetan Plateau [56,57], the impacts of winter snow cover on temperate grassland vegetation are thought to be similar to those on alpine vegetation.As the relationships between the SOS and winter snow cover might vary across different vegetation types, the possible effects of winter snow cover across temperate grassland vegetation in Inner Mongolia should be further investigated.
It is worth noting that the SOS changes were related not only to the temperature and precipitation, but also to other factors (e.g., winter snow cover, GST, SSD, and RHU).Winter snow cover could add soil moisture, arrest soil heat exchange, and have crucial effects on soil heat and moisture preservation [58,59], which would protect the over-winter survival of vegetation from wind and low air temperatures.Thus, it could influence the exchange of gas, moisture, and dissolved matter.Furthermore, it could affect the activity of the microbes in the soil as well as the transformation of soil organic matter and nutrients.In general, deeper snow depth may lead to a later SOS.The reason might be that a deeper snow depth would have a longer SCD and later SCED, which would not be beneficial for soil respiration in winter and the accumulation of heat needed to unfold leaves in spring [60][61][62].Earlier snow cover could shield vegetation from low temperatures by keeping the temperature at a favorable level and protecting the activity of the soil microorganisms during the winter, which could ensure positive vegetation growth during the growing season in the following year.Moreover, the spring phenology was more strongly and negatively associated with temperature and precipitation during the SCED and SOS, which indicates that the SOS will advance with warmer temperatures or more precipitation.In contrast, the SOS showed a significant positive correlation with Tmax and Tavg for desert steppe.Increased temperatures (Tmax and Tavg) would result in the loss of copious amounts of water by evaporation [20] and exacerbate the effects of drought in dry regions, which would delay the SOS for desert steppe [63,64].Increasing Tmin could reduce frost damage to stimulate vegetation growth.A suitable sunshine environment, which plays a role equally important as precipitation and temperature, is beneficial for photosynthesis and respiration.Meanwhile, excess SSD may also contribute to the aggravation of drought through evapotranspiration.The relative humidity exhibited a negative relationship with the SOS, and its change was usually associated with the temperature, sunshine, wind, and precipitation.
In addition, the frequency of extreme events has been increasing under rapid climate change.Extreme events during the SCED and SOS, such as drought and frost, will influence spring vegetation phenology in the snow cover regions.Therefore, the impacts of the extreme events on the spring phenology should be investigated.

Conclusions
This paper investigated the SOS of grassland vegetation in Inner Mongolia from 1982 to 2015, and explored the responses of SOS to the winter snow cover, climate change during snow cover melt, and growing season climate factors (including average temperature, precipitation, relative humidity, and sunshine duration).The conclusions of this study can be drawn as follows: (1) During 1982-2015, 52.7% of the Inner Mongolia grassland experienced a significant advancing trend in SOS and 34.30% exhibited a delaying trend.The average SOS occurred on DOY 120 for meadow steppe and typical steppe, and on DOY 130 in the desert steppe.All three grassland vegetation types exhibited an earlier spring season at a rate of 0.3 ± 0.74 days/year across all of grassland vegetation types, and rates of 0.27 ± 0.47 days/year, 0.32 ± 0.65 days/year, and 0.47 ± 1.48 days/year for temperate meadow steppe, typical steppe, and desert steppe, respectively.(2) Winter snow cover showed a positive correlation with the SOS.By contrast, SCOD showed an opposite correlation.Focusing on the correlation between snow cover and SOS, we found that the SOS was more strongly associated with SCED, SD, and SCOD for meadow steppe, typical steppe, and desert steppe, respectively.Furthermore, the climate during snowmelt and SOS was also a significant factor contributing to the change in SOS.Higher temperatures and more precipitation advanced SOS, whereas Tmax and Tavg showed a positive correlation with SOS for desert steppe.Increasing Tmin would reduce the number of frost events and promote vegetation growth.Sunshine hours and relative humidity showed weaker correlations.
It should also be noted that this present study may have several limitations.Although the vegetation type data are fairly reliable, the possible effects of human activities, such as grazing, cannot be excluded completely.

Figure 1 .
Figure 1.(a) The geographic location of Inner Mongolia; (b) the distribution of grassland vegetation types; and (c) the elevation, weather stations and phenology observation stations of Inner Mongolia.

Figure 1 .
Figure 1.(a) The geographic location of Inner Mongolia; (b) the distribution of grassland vegetation types; and (c) the elevation, weather stations and phenology observation stations of Inner Mongolia.

Figure 2 .
Figure 2. A schematic diagram illustrates the retrieval of spring phenology using the logistic fitting method [45].The solid line indicates the fitted logistic curve and the dashed line is the rate of change in curvature of the fitted logistic curve.Onset date of spring phenology (SOS) is defined as the first local maximum of the dashed curve.The red line indicates that the vegetation index begin to increase rapidly.

Figure 2 .
Figure 2. A schematic diagram illustrates the retrieval of spring phenology using the logistic fitting method [45].The solid line indicates the fitted logistic curve and the dashed line is the rate of change in curvature of the fitted logistic curve.Onset date of spring phenology (SOS) is defined as the first local maximum of the dashed curve.The red line indicates that the vegetation index begin to increase rapidly.

13 Figure 3 .
Figure 3. (a) Spatial distribution of the mean SOS in the grassland of Inner Mongolia during 1982-2015; (b) spatial distribution of the change trend for SOS in the grassland of Inner Mongolia during 1982-2015.The SOS change trends are termed significant for pixels at p < 0.05 level.

Figure 4 .
Figure 4. 34 year average of SOS and its change trends for different vegetation types: (a) meadow steppe, (b) typical steppe, and (c) desert steppe.Each bin represents a 10 day range of SOS.The height and color of each bin indicate the number and fitting slope (i.e., the change trend) of the pixels that fall within the bin, respectively, with the color bar of the slope on the bottom of the figure.Only statistically significant trends (p < 0.05) are included.

Figure 3 .
Figure 3. (a) Spatial distribution of the mean SOS in the grassland of Inner Mongolia during 1982-2015; (b) spatial distribution of the change trend for SOS in the grassland of Inner Mongolia during 1982-2015.The SOS change trends are termed significant for pixels at p < 0.05 level.

Figure 3 .
Figure 3. (a) Spatial distribution of the mean SOS in the grassland of Inner Mongolia during 1982-2015; (b) spatial distribution of the change trend for SOS in the grassland of Inner Mongolia during 1982-2015.The SOS change trends are termed significant for pixels at p < 0.05 level.

Figure 4 .
Figure 4. 34 year average of SOS and its change trends for different vegetation types: (a) meadow steppe, (b) typical steppe, and (c) desert steppe.Each bin represents a 10 day range of SOS.The height and color of each bin indicate the number and fitting slope (i.e., the change trend) of the pixels that fall within the bin, respectively, with the color bar of the slope on the bottom of the figure.Only statistically significant trends (p < 0.05) are included.

Figure 4 .
Figure 4. 34 year average of SOS and its change trends for different vegetation types: (a) meadow steppe, (b) typical steppe, and (c) desert steppe.Each bin represents a 10 day range of SOS.The height and color of each bin indicate the number and fitting slope (i.e., the change trend) of the pixels that fall within the bin, respectively, with the color bar of the slope on the bottom of the figure.Only statistically significant trends (p < 0.05) are included.

Figure 5 .
Figure 5.The change trends of SOS for different vegetation types with different window lengths (i.e., study period): (a) meadow steppe; (b) typical steppe and (c) desert steppe.Each dot represents mean trends of each vegetation type over a single window size (1-33 years) for each start year (1982-2014).Only statistically significant trends (p < 0.05) are included in this analysis.The red (green) color indicates positive (negative) trends.

Figure 5 .
Figure 5.The change trends of SOS for different vegetation types with different window lengths (i.e., study period): (a) meadow steppe; (b) typical steppe and (c) desert steppe.Each dot represents mean trends of each vegetation type over a single window size (1-33 years) for each start year (1982-2014).Only statistically significant trends (p < 0.05) are included in this analysis.The red (green) color indicates positive (negative) trends.

Figure 6 .
Figure 6.The interannual variability of each climate factor in different grassland types.The value in each grid indicates the rate of change for each climate factor.Values are also color-coded, with blue indicating negative values and red indicating positive values.The color of each grid corresponds to the value, with the color bar on the right of the table.The asterisks (*) indicate the climate factors trends that are statistically significant at the p < 0.1 level, and the double asterisks (**) indicate significance at the p < 0.05 level.

Figure 6 .
Figure 6.The interannual variability of each climate factor in different grassland types.The value in each grid indicates the rate of change for each climate factor.Values are also color-coded, with blue indicating negative values and red indicating positive values.The color of each grid corresponds to the value, with the color bar on the right of the table.The asterisks (*) indicate the climate factors trends that are statistically significant at the p < 0.1 level, and the double asterisks (**) indicate significance at the p < 0.05 level.

9 of 13 13 Figure 7 .
Figure 7.The correlation between SOS and climate factors for three vegetation types: (a) meadow steppe; (b) typical steppe; and (c) desert steppe.Each dot represents the correlation between SOS and the corresponding climate factor, with the red color indicating a positive correlation and blue color indicating a negative correlation.The asterisks (*) indicate the trends of the climate factors trends that are statistically significant with p < 0.1 and the double asterisks (**) indicate trends with p < 0.05.

Figure 7 .
Figure 7.The correlation between SOS and climate factors for three vegetation types: (a) meadow steppe; (b) typical steppe; and (c) desert steppe.Each dot represents the correlation between SOS and the corresponding climate factor, with the red color indicating a positive correlation and blue color indicating a negative correlation.The asterisks (*) indicate the trends of the climate factors trends that are statistically significant with p < 0.1 and the double asterisks (**) indicate trends with p < 0.05.
arid northwest inland region, which is affected by Asian monsoons.The climate in Inner Mongolia can be characterized by a cold, dry winter and warm, rainy summer