Assessing Snow Phenology and Its Environmental Driving Factors in Northeast China

Abstract

Snow cover is an important water source and even an Essential Climate Variable (ECV) as defined by the World Meteorological Organization (WMO). Assessing snow phenology and its driving factors in Northeast China will help with comprehensively understanding the role of snow cover in regional water cycle and climate change. This study presents spatiotemporal variations in snow phenology and the relative importance of potential drivers, including climate, geography, and the normalized difference vegetation index (NDVI), based on the MODIS snow products across Northeast China from 2001 to 2018. The results indicated that the snow cover days (SCD), snow cover onset dates (SCOD) and snow cover end dates (SCED) all showed obvious latitudinal distribution characteristics. As the latitude gradually increases, SCD becomes longer, SCOD advances and SCED delays. Overall, there is a growing tendency in SCD and a delayed trend in SCED across time. The variations in snow phenology were driven by mean temperature, followed by latitude, while precipitation, aspect and slope all had little effect on the SCD, SCOD and SCED. With decreasing temperature, the SCD and SCED showed upward trends. The mean temperature has negatively correlation with SCD and SCED and positively correlation with SCOD. With increasing latitude, the change rate of the SCD, SCOD and SCED in the whole Northeast China were 10.20 d/degree, −3.82 d/degree and 5.41 d/degree, respectively, and the change rate of snow phenology in forested areas was lower than that in nonforested areas. At the same latitude, the snow phenology for different underlying surfaces varied greatly. The correlations between the snow phenology and NDVI were mainly positive, but weak correlations accounted for a large proportion.

1. Introduction

Northeast China, one of the three major areas with snow cover in China [1,2], is an important agricultural production base. Snow cover can affect the scale and yield of agriculture by changing soil moisture, insultation from deep frost and prevention of microbial decomposition of organic matter. In addition, Snow accumulation and melting are important for storing and releasing water [3,4]. Snow phenology can describe seasonal snow cover variations and is also used to study the relationship between snow cover and climate change; parameters include the snow cover days (SCD), snow cover onset dates (SCOD) and snow cover end dates (SCED) [5,6]. According to the statistics, snow phenology generally has regular interannual variations [7]. Therefore, assessing snow phenology and its driving factors in Northeast China is essential for water resource management and agricultural development in this region.

However, in contrast with research on large-scale snow phenology, systematic studies of snow phenology in the region have been limited thus far. Chen et al. [8] used the MODIS snow products across Northeast China to study the spatiotemporal variations in snow cover. Ding and Gao [4] studied the SCD in Northeast China based on meteorological station data. Yang et al. [9] used MODIS data to explore only space-time distributions of SCD. Qiao et al. [10] used MODIS data to investigate the variations in snow phenology and their impact on vegetation growth in forested areas, which occupies 40% of the total area. Shi et al. [11] discussed the SCOD and SCED in the Mollisol areas across the northeastern plains of China based on snow depth data. In addition, Huang et al., Ke et al., and Ma et al. [5,6,12] used meteorological station data to analyze the snow phenology changes in snow-covered areas of China, including Northeast China.

When studying snow phenology based on meteorological observation stations, the stations give snow depth data for only the surrounding area [13]. Moreover, due to the scarcity and uneven distribution of stations, there are large discontinuities in the spatial distribution of the obtained snow phenology data, especially in forested areas and alpine regions [5,14,15,16]. The snow depth data retrieved by passive microwaves have a long time series, but the resolution is low [17], which is not appropriate for local regional research. The same is true of the Northern Hemisphere snow cover extent (NHSCE), which is more suitable for studying large-scale snow cover variations [18,19]. In contrast, MODIS data have not only higher spatial resolution but also higher temporal and spectral resolutions [20]. Therefore, MODIS is an ideal data source for studying continuous snow phenology, whether at the global or regional scale. For cloud contamination, many scholars have also studied and recovered data under clouds through a series of methods and achieved high cloud removal accuracy [21,22,23,24,25,26].

Based on the above situation, we analyzed the snow phenology of Northeast China from 2001 to 2018 in this study using the MODIS snow product. First, the daily cloud-free snow products were obtained through the conditional probability interpolation method based on a space-time cube, and accuracy was verified by ground observations. Based on this work, we explored the spatiotemporal variations in snow phenology and the relative importance of potential drivers, including climate, geography, and the NDVI, and we then discussed the roles of major factors in driving snow phenology.

2. Study Area

Northeast China is located at 115°30′~135°60′E, 38°42′~53°36′N, and the area is approximately 1.24 million km². The whole region is surrounded by the Daxingan Mountains, Xiaoxingan Mountains, and Chang-bai Mountains, with the Northeast Plain in the middle. It is a relatively independent and complete natural geographical area (Figure 1). The annual average temperature, precipitation, and elevation are −5~11 °C, 300~1000 mm, and 443 m, respectively [27,28,29]. Due to the existence of large-scale forests and climate factors, this region in China has a long snow period and thick snow cover [30]. The Daxingan Mountains, Xiaoxingan Mountains, and Changbai Mountains are the most important snowfall areas, with an average annual snowfall of more than 60 mm. The average snowfall is more than 90 mm in the northern Daxingan Mountains and Changbai Mountains, while in some plain areas, the average snowfall is less than 30 mm [31]. Here, the snowfall means the snow water equivalent.

3. Methodology

3.1. Data Sources

3.1.1. MODIS Data

Daily MOD10A1 (V006) data with a 500 m resolution and containing 6574 scenes during 2001–2018 were downloaded with the help of the Google Earth Engine (https://code.earthengine.google.com (30 October 2021)). The product has four data types: NDSI data, quality assessment, albedo, and orbital information. In this research, only NDSI_Snow_Cover data were used as the main data source to extract snow products, and the NDSI of snow cover always greater than 0, but not all surface features with NDSI > 0 is snow cover. Detailed attribute information is shown in

Table 1. MODIS Snow Products (V006).

Value	Attributes	Value	Attributes
0–100	NDSI_Snow_Cover	239	ocean
200	missing data	250	cloud
201	no decision	254	detector saturated
211	night	255	fill
237	inland water

[32].

3.1.2. Snow Depth Records

The snow depth records in Northeast China from 1 January 2013 to 31 December 2018 were collected from the National Cryosphere Desert Data Center (http://www.ncdc.ac.cn (30 October 2021)), and these records in snow season (October to April of the following year) were used to validate the cloud removal accuracy. The station distribution is shown in Figure 1. Generally, SD ≥ 1 cm indicates snow cover; otherwise, the station is regarded as snow-free [6,12,33].

3.1.3. Meteorological Data

Meteorological data with a 1 km resolution were acquired from the National Earth System Science Data Center, National Science and Technology Infrastructure of China (http://www.geodata.cn (30 October 2021)). The dataset is monthly data, and in this study, the temperature/precipitation in different months, the annual mean temperature (Temp-mean) and annual total precipitation (Pre-total) in hydrological years were used to explore the roles in snow phenology.

3.1.4. DEM

SRTM digital elevation data with a spatial resolution of 90 m was obtained from the Google Earth Engine (GEE). In addition, the data was resampled to 500 m in order to keep the resolution consistent with MODIS snow products. Based on this, the altitude, slope and aspect were derived.

3.1.5. Vegetation Data

To analyze the variations in snow phenology for different vegetation types in Northeast China, the MODIS vegetation product MOD13A1 normalized difference vegetation index (NDVI) from 2001 to 2018 was derived through the GEE cloud platform with spatial resolution and temporal resolution of 500 m and 16 d, respectively. In addition, the monthly/yearly NDVI maximum of all grids was then calculated by the maximum synthesis method to eliminate the interference of clouds and shadows.

3.1.6. Land Cover Data

MCD12Q1 is an annual land cover product with five classification standards. In this paper, the International Geosphere-Biosphere Program (IGBP) classification standard was chosen, and land cover types in this region were reclassified into water bodies, forested areas (including broadleaf forest, coniferous forest, mixed forest and shrub), and nonforested areas (except water and forested areas), and then we explore the differences in snow cover phenology for various underlying surfaces at the same latitude.

3.2. Methodology

3.2.1. MODIS Snow Product Cloud Removal

Due to cloud contamination, we used a noval cloud removal algorithm to obtain daily cloud-free MODIS snow products. First, a threshold of 10 (expand 100 times) was used to delineate between snow and snow-free conditions in the NDSI_Snow_Cover data layer [34], and the rest attributes were classified as clouds, except for inland water and ocean. Second, the conditional probability interpolation method based on a space-time cube was used to remove the clouds. In addition, snow probability of the cloud pixels was calculated by using the conditional probability of the central pixel and every neighboring pixel in a space-time cube of 5 × 5 × 5 under the same snow condition as the weight. Finally, the snow condition of pixels covered with clouds was recovered according to the snow probability [26]. The formulas are as follows:

$$P\left(C_{x,y}|C_n\right)=\frac{\sum 1-ABS\left(C_{x,y,t}-C_{n,t^{\prime }}\right)}{N_{x,y}}$$

(1)

$$P\left(x_0,y_0,t_0\right)=\frac{\sum P\left(C_{x_0,y_0}|C_n\right)\times S_n}{\sum P\left(C_{x_0,y_0}|C_n\right)\times V_{n^{\prime }}}$$

(2)

$$C\left(x_0,y_0,t_0\right)=\left\{\begin{array}{c}snow, P\left(x_0,y_0,t_0\right)\ge 0.5\\ snowfree, P\left(x_0,y_0,t_0\right)<0.5\end{array}\right.$$

(3)

Here, $P\left(C_{x,y}|C_n\right)$ is the conditional probability having the same snow condition for the central and n-th adjacent pixels in the space-time cube. $C_{x,y,t}$ and $C_{n,t^{\prime }}$ represent snow (C = 1) or snow-free (C = 0) conditions for days t and t′, respectively. $N_{x,y}$ are cloud-free days for the central pixel and n-th neighboring pixels within the study time; $P\left(x_0,y_0,t_0\right)$ is the snow probability of the cloud gaps. $S_n$ means that the n-th pixel has snow ($S_n$= 1) or is snow-free ($S_n$= 0); $V_{n^{\prime }}$ indicates whether the n-th pixel is cloudless ($V_{n^{\prime }}$= 1) or covered by clouds ($V_{n^{\prime }}$= 0); and $C\left(x_0,y_0,t_0\right)$ is the snow condition.

3.2.2. Snow Phenology Calculation

Snow phenology mainly includes snow cover days (SCD), snow cover onset dates (SCOD) and snow cover end dates (SCED). In this study, daily cloud-free snow products were obtained through cloud removal algorithm above, and snow phenological parameters in hydrological year, which was defined from 1 September to 31 August of the following year, were calculated pixel by pixel. The SCD was the total days when a pixel is snow in a hydrological year. The SCOD was the first day when pixel was covered with snow lasting at least five days for the first time, and SCED was the last day of at least 5 days of continuous snow. This avoided the influence of instantaneous snowfall [10,16,35].

3.2.3. Cloud-Free Snow Product Accuracy Assessment

Currently, station data from meteorological observatories are usually regarded as “truth” data and used to evaluate the cloud removal accuracy. The accuracy assessment metrics include overall accuracy (OA), underestimation error (UE) and overestimation error (OE) based on the confusion matrix (

Table 2. Confusion matrix.

		MODIS
		Snow	Snow-Free
Truth	snow	a	b
	snow-free	c	d

), which are defined as follows:

$$OA=\frac{a+d}{a+b+c+d}$$

(4)

$$UE=\frac{b}{a+b+c+d}$$

(5)

$$OE=\frac{c}{a+b+c+d}$$

(6)

The definitions of a, b, c and d are given in Table 2. OA represents the proportion that pixels are consistent with the truth and MODIS classification; UE is the proportion that pixels are snow-free in MODIS, but the corresponding pixels in the truth are covered with snow; and OE refers to the proportion that pixels are covered with snow in MODIS, but the corresponding pixels in the truth are snow-free.

3.2.4. Trend Analysis

The Mann–Kendall test and Theil–Sen median analysis were used to explore the variation trend in snow phenology from hydrological years 2001–2018. In the Mann–Kendall test, to calculate whether the trend of snow phenology was increasing or decreasing at the 0.05 confidence level, the Z value was divided into five types: significant decrease (Z < 1.96), slow decrease (−1.96 ≤ Z < 0), nonsignificant change (Z = 0), slow increase (0 > Z ≥ 1.96) and significant increase (Z > 1.96). The formulas are as follows:

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{var\left(S\right)}}，S>0\\ 0 ，S=0\\ \frac{S+1}{\sqrt{var\left(S\right)}}，S<0\end{array}\right.$$

(7)

where,

$$var\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}$$

(8)

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sgn\left(S_j-S_i\right)$$

(9)

$$sgn(S_j-S_i)=\left\{\begin{array}{c}1，S_j-S_i>0\\ 0，S_j-S_i=0\\ -1，S_j-S_i<0\end{array}\right.$$

(10)

When $Z$ > 0, the trend is upward, and when $Z$< 0, it is downward. $S_i$/$S_j$ represent the value in years i/j, n is the length of the time series. When |Z| > Z_{1 − α/2} (α is the significance level), the trend is significant in the time series. In this paper, α = 0.05 was used.

In the Theil–Sen median analysis, to explore the details of the variations in snow phenology, S_snow was divided into seven types: <−4 d/a, -<−4–−2 d/a, −2–0 d/a, 0 d/a, 0–2 d/a, 2–4 d/a and >4 d/a. The formula is as follows:

$$S_{snow}=Median\left(\frac{S_j-S_i}{j-i}\right)\text{，}\forall j>i$$

(11)

where $S_{snow}$ > 0 and $S_{snow}$ < 0 represent upward and downward trends, respectively.

3.2.5. Relative Importance of Multiple Factors to Snow Phenology

The use of geodetector is a common statistical approach that can analyze spatial variations and reveal the driving factors behind them [36]. A geodetector contains four subdetectors: factor detector, risk detector, ecological detector and interaction detector [37]. In this research, we employed a factor detector to quantify the relative contributions of vegetation and geographical (altitude, slope, aspect, latitude and longitude) and meteorological (temperature and precipitation) factors to snow phenology variations, and the dominant driving factor was then the highest contributor to snow phenology variations.

The factor detector is calculated by the following q-statistic:

$$q=1-\frac{{{\displaystyle \sum }}_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}=1-\frac{SSW}{SST}$$

(12)

$$SSW ={\displaystyle \sum }_{h=1}^L{N}_h{\sigma }_h^2,SST = N{\sigma }^2$$

(13)

where 0 ≤$q$≤ 1, and the larger the value, the greater the influence of the factor. h is the number of strata for variables or factors, N represents the number of units in stratum h, and ${\sigma }^2$ and ${\sigma }_h^2$ denote the variance in the entire study area and stratum h, respectively. SW and SST are the sum of squares within the data and the total sum of squares, respectively.

3.2.6. Correlation Analysis

To investigate the influence of major driving factors on snow phenology, we calculated the correlation coefficients between these factors and snow phenology:

$$r_{xy} =\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(14)

where x_i and y_i represent the values in the i-th year and $\overline{x}$ and $\overline{y}$ are the average values for all years. If r > 0, two variables are positively correlated, and if r < 0, it is negatively correlated. When |r| ≤ 0.3, two variables are weakly correlated or have no correlation; 0.3 < |r| ≤ 0.5 indicates that there is a moderate correlation between two variables; 0.5 < |r| ≤ 1 indicates strong correlation [38].

4. Results

4.1. Validation of the Daily Cloud-Free MODIS Snow Products

Based on above accuracy assessment metrics, the available station data in Northeast China were screened to quantitatively evaluate the cloud removal results. The accuracies of the three evaluation indicators were obtained (

Table 3. Accuracy evaluation results for ground observation stations.

Time	OA	UE	OE
1 October 2013–30 April 2014	0.95	0.03	0.02
1 October 2014–30 April 2015	0.94	0.05	0.01
1 October 2015–30 April 2016	0.92	0.06	0.03
1 October 2016–30 April 2017	0.94	0.04	0.02
1 October 2017–30 April 2018	0.93	0.05	0.02

), and Figure 2 summarizes the monthly accuracy evaluation results for ground observation stations. Overall, the OA of the MODIS cloud-free product was more than 0.9, the UE and OE were less than 0.1, and the UE values were greater than the OE values. Generally, the accuracies in the snow stabilization season were higher than that in the snow accumulation and melting season. However, in the October and April, the OA was the highest, and the UE and OE were the lowest, which may be attributed to the lack of snow during this period. There seems odd relative to others in the February, which is due to the lack of 9 days MODIS data (19 February 2016–27 February 2016). In the cloud removal processing, the missing data are thought to be cloud. Therefore, the cloud removal accuracy is very low for these data because there is no spatial and temporal neighbor data available. The analysis reveals that the cloud removal product had good reliability and met the needs of the snow phenology analysis.

4.2. Spatiotemporal Variations and Trends in Snow Phenology

4.2.1. SCD

The distribution of SCD in Northeast China was visually consistent with those of topography and latitude, showing the characteristics of “high SCD in the mountains and low SCD in the plains” and “high SCD at high latitudes and low SCD at low latitudes”. The SCD values in the northern Daxingan Mountains were more than 150 d, those of the Xiaoxingan Mountains were mainly within 120–150 d, and the SCD values of the plain area were significantly less than those of the mountainous area. The SCD values of the Sanjiang Plain were within 90–150 d, and the Liaohe Plain had the shortest SCD values of less than 60 d (Figure 3).

Figure 4 showed the histogram of annual average SCD computed for a 10-day interval from HY2001 to HY2017. The average SCD in Northeast China ranged from 0 to 230 d, with an obvious bimodal distribution, and the annual average SCD was 93 d. The frequency of SCD in the range from 10–20 d accounted for 9.95%, and that in the range from 120–130 d accounted for 13.92%.

Figure 5a indicates that the SCD in Northeast China mainly increased, accounting for 62.27% of the total, among which only 4.73% increased significantly. The regions where the SCD increased were mainly in the southern Daxingan Mountains, Xiaoxingan Mountains and Changbai Mountains. SCD decreased mainly in the northern Daxingan Mountains and the southwestern Changbai Mountains, accounting for 20.25% of the total area, and only 0.24% of the total area recorded significant decreases. Figure 5b showed the area with an SCD trend >0 d/a accounted for 60.61%, which was similar to the spatial distribution of SCD with an increasing trend; the area—revealing no trend was concentrated in the Liaohe Plain and Songnen Plain, accounting for 16.80%; the areas with trends of less than 0 d/a accounted for 18.65% and were mainly located in the northern Daxingan Mountains and the southwestern Changbai Mountains. The results of the Theil–Sen trend and Mann–Kendall significance test were consistent, which further verified the accuracy of these trends.

4.2.2. SCOD

Figure 6 shows that the SCOD was delayed from north to south, which meant an obvious increase with latitude. At high latitudes, the SCOD were mainly concentrated within 60–90 d, that is, in November of the current year, and in the northern part of the Daxingan Mountains, the SCOD were at the end of September or early October of the current year. In low-latitude areas, the SCOD were concentrated within 90–120 d, that is, December of the current year. The SCOD were later in some areas of the southern Songnen Plain, and the latest was in January of the next year. The annual average SCOD in Northeast China were mainly concentrated within 60–100 d, from November to the middle of January of the next year, accounting for 84.78% of the total area. Among them, the peak was 70–80 d, which was approximately the middle of November of the current year, with a proportional area of 30.22% (Figure 7).

Figure 8 shows the significance tests and the trend of SCOD from HY2001 to HY2017. The results of the Theil–Sen trend and Mann–Kendall significance test were consistent. The area in which the SCOD showed an early trend accounted for 31.93%, while the delayed trend accounted for 32.80%. The trends of significant advance and delay were relatively small, only 0.24% and 0.12%, respectively (Figure 8a). Overall, the area with an increase in Figure 8a corresponded to an SCOD trend > 0 d/a in Figure 8b, and that the area with a decrease corresponded to an SCOD trend < 0 d/a.

4.2.3. SCED

The distribution of SCED was consistent with that of SCD, presenting “high SCED in the mountains and low SCED in the plains” and “high SCED at high latitudes and low SCED at low latitudes” trends. From the perspective of spatial distribution, the SCED in some areas of the Daxingan Mountains, Xiaoxingan Mountains and Changbai Mountains were mainly in the range of 210 d to 240 d in approximately April of the next year. The SCED of the Songnen Plain and Liaohe Plain were relatively early, probably between January and February of the next year (Figure 9). The annual average SCED distribution in Northeast China also had two peaks, which were 130–140 d and 200–210 d, corresponding to the middle of January and late April of the next year, accounting for 5.54% and 20.50% of the area, respectively. The area with an SCED of less than 240 d accounted for 99.59%; that is, at the end of April of the next year, the snow in Northeast China almost disappeared (Figure 10).

For the SCED, 29.44% of the total area showed an early trend, while 36.70% had a delayed trend. In addition, only 1.49% and 2.53% of the area were significantly early and delayed, respectively (Figure 11a). The delayed areas were mainly distributed in the Sanjiang Plain, Daxingan Mountains, Xiaoxingan Mountains and Changbai Mountains; the areas with an advancing trend of SCED were mainly concentrated in the middle of the region and the transition area between plain and mountainous areas. Figure 11b shows the trend of SCED in Northeast China. Compared with the significance test, the overall spatial distribution characteristics and the increasing/decreasing trends of SCED were highly consistent. The area with an SCED trend > 0 d/a accounted for 36.36%; that with an SCED trend < 0 d/a accounted for 29.11%.

In general, the correlation analysis of the SCD, SCOD, and SCED revealed significant relationships between the SCED and SCD (r = 0.70). The correlation between the SCOD and SCD was −0.14, and the increasing trend of the SCD was determined by the advancement of the SCOD and the delay of the SCED. Considering that there was no significant change in the SCOD across time, the snow phenology variations in Northeast China from HY2001 to HY2017 were attributed mostly to the changes in SCED.

4.3. Roles of Multiple Factors in Snow Phenology

Figure 12 shows that the geographical and meteorological factors and the NDVI all affected the SCD, SCOD and SCED. Annual mean temperature had the greatest impact on the SCD, SCOD and SCED, followed by latitude. Precipitation, aspect and slope all had little effect on the SCD, SCOD and SCED, and all these q values were less than 0.1. Compared with the SCOD, the NDVI and longitude both had a greater impact on the SCED and SCD, with q values of 0.35 (0.15) and 0.30 (0.13), respectively. However, altitude was an important factor affecting the SCOD compared with the SCD and SCED.

We further analyzed the relative importance of monthly temperature and precipitation to snow phenology, and the results are shown in Figure 13. Obviously, the influence of monthly temperature on snow phenology is much greater than that of monthly precipitation. In terms of the roles of different months, the temperature from January to May has a stronger impact on snow phenology, while the impact from June to August is relatively weak. Precipitation has a greater impact on snow phenology in May, September and snow stable period, which the snow completely melted in May, and snowfall occurred in some areas in September.

5. Discussions

5.1. Response of Snow Phenology to Climate

Figure 14 shows the annual variations in the mean temperature and snow phenology from HY2001 to HY2017. Except for the SCOD, there were clear increases in both the SCD and SCED and a decrease in temperature. According to statistics, SCD and SCED were both strongly negatively correlated with the mean temperature; at the 95% confidence level, the correlation coefficients were −0.73 and −0.57, respectively. The correlation with the SCOD was moderate (r = 0.41), which can be explained that besides mean temperature, the effect of latitude on SCOD was also great (q = 0.55).

The spatial pattern of correlation between the mean temperature and snow phenology is presented in Figure 15. For 99.91% of the pixels, the mean temperature was negatively correlated with the SCD, which means the lower the mean temperature was, the longer the SCD in the whole study area; 59.40% of the pixels showed strong negative correlations, and areas with weak negative correlation were mainly distributed in the northern Daxingan Mountains and Xiaoxingan Mountains (Figure 15a). The correlations between the SCOD and mean temperature were mainly positive, and the weak positive correlations accounted for a large proportion (40.81%). Areas with strong positive correlations accounted for only 8.19% and were distributed mainly in the Xiaoxingan Mountains and Changbai Mountains (Figure 15b). The correlations between the mean temperature and SCED were similar to those of the SCD, and most regions had negative correlations (93.86%). The difference was that the proportion of strong negative correlations was relatively small (33.06%), and the average correlation of the SCED was −0.38, while that of the SCD was −0.52 (Figure 15c).

5.2. Geographical and Vegetation Controls on Snow Phenology

To quantitatively investigate the latitudinal zonation of snow phenology further, statistics were analyzed in combination with the underlying surface conditions. The slope and the regularity between snow phenology and latitude in the nonforested and forested areas further quantitatively proved the latitudinal zonation of snow phenology, as shown in

Table 4. Linear slope and R² of the mean SCD, SCOD, and SCED with latitudes from 39.22°N to 53.22°N in the nonforested and forested areas in Northeast China.

	SCD		SCOD		SCED
	Slope	R2	Slope	R2	Slope	R2
Northeast China	10.2	0.97	−3.82	0.89	5.41	0.96
Nonforested areas	7.8	0.76	−3.41	0.74	4.12	0.71
Forested areas	5.41	0.98	−2.02	0.93	2.91	0.97

Note: The slope represents the change rate in the days when the latitude increased by 1° (d/degree).

and Figure 16. In Northeast China, as the latitude increased by 1 degree, the average SCD increased by 10.2 d, the SCOD advanced by 3.82 d and the SCED was delayed by 5.41 d. Generally, the change rate in forested areas was slower. The change rate in the SCD, SCOD and SCED were 5.41 d/degree, −2.02 d/degree and 2.91 d/degree, respectively, which were closely related to the closed environments of the forested areas themselves.

The higher the latitude was, the longer the SCD, the earlier the SCOD and the later the SCED, in both nonforested and forested areas. At the same latitude, the SCD of forested areas was higher than that of nonforested areas. At high latitudes, the SCOD of forested areas was later than that of nonforested areas. For the SCED, at low latitudes, the snow melting time in forested areas was later than that in nonforested areas, and the SCED was basically the same at high latitudes; these changes can be explained by the observations that snow melting in high-latitude areas was mainly affected by temperature and that the SCED in this region was in approximately April, so the rapidly increasing temperature led to the melting of snow in the forested and nonforested areas.

In contrast to the mean temperature, the NDVI denoting vegetation greenness was positively correlated with the SCD, SCOD and SCED, accounting for 68.69%, 59.22% and 58.68%, respectively. However, the weak correlation (−0.3< r <0.3) between snow phenology and the NDVI was dominant, with 69.85%, 74.99% and 72.65% of pixels (Figure 17a–c), which illustrates that the NDVI was not the main factor affecting snow phenology.

5.3. Comparison with Previous Results

Shi et al. [11] reported strong spatial heterogeneity in snow phenology in the Mollisol region of Northeast China. The SCD increased from southwest to northeast gradually, and snow cover began to accumulate in mid-November and completely melted in late March during 1978–2016. In the western Changbai Mountains and northern Daxingan Mountains, the SCOD always began in early October, and the SCED in the western Daxingan Mountains, Xiaoxingan Mountains, and Changbai Mountains always occurred in May [8]; which was similar to the results of this study. In the forested areas of Northeast China, the SCOD were later, and the SCED were earlier in the plains. The SCOD occurred between late November and mid-November in the regions with high altitudes, and the SCED occurred later with increasing latitude and altitude. The earlier SCOD and later SCED led to an increase in the SCD, especially in high mountain areas. In most areas, the SCOD experienced an advancing trend, the SCED exhibited an obvious delaying trend, and the SCD showed an opposing trend from south to north in Northeast China from 2004 to 2018 [10].

This paper found that SCD and SCED increased, and SCOD basically did not change in Northeast China from 2001 to 2018. However, previous studies on snow phenology in China found that due to the increase in temperature, the SCOD in most areas were delayed, and SCED were advanced [5,6,12]. These contradictions may have been caused by the inconsistency of the time span. Although snow cover will decrease under the background of global warming, the decrease of mean temperature in a short time has led to the increase of snow cover. Studies have shown that there are strong correlations between meteorological and geographic factors and snow cover [39,40,41]. Vegetation change was also closely related to snow phenology [10,42,43]. Temperature and precipitation could affect snow cover variations [44,45]. However, in this research, the results showed that temperature was the main factor affecting the variations in snow phenology, and precipitation had little effect. In the three basins of Songhua River in Northeast China, temperature, precipitation and altitude were considered to be the three most important factors [9]. Huang et al. [12] studied the snow cover variations across China from 1951 to 2018 based on snow depth dataset and used model to analyze the driving effect of multiple factors on snow cover phenology and found that the most important factors influencing the SCD, SCOD and SCED were annual coldest monthly minimum temperature, altitude and annual mean temperature, respectively. This difference may have been caused by inconsistencies in time scales, study areas, data and research methods. In addition, we both found that precipitation had little effect on snow phenology.

6. Conclusions

Snow cover is one of the most active features on the land surface [46], and it is very important to quantitatively explore and understand snow phenology variations for local climate change. The snow phenology variations and their driving factors in Northeast China from 2001 to 2018 were assessed based on daily cloud-free snow products generated in this research. The conclusions were as follows.

(1)

The SCD, SCOD and SCED all showed the characteristics of latitudinal zonal distribution, and the SCED and SCD distributions had obvious consistency. With increasing latitude, the SCD was longer, the SCOD began earlier, and the SCED appeared later. Overall, the SCD showed mainly an increasing trend, which was mostly distributed in the southern Daxing’an Mountains, Xiaoxing’an Mountains and Changbai Mountains. The SCOD showed advanced and delayed trends that accounted for 31.93% and 32.80%, respectively. The corresponding proportions of the SCED accounted for 29.44% and 36.70%, respectively, which meant that the SCED showed a delayed trend overall. On the Liaohe Plain and Songnen Plain, the snow phenology basically did not change.

(2)

For snow phenology, the mean temperature was identified as the most important driver, followed by latitude. In terms of the roles of temperature in different months, the snow phenology is mainly affected by the temperature in winter of current year and spring of the next year. The decrease in temperature directly led to the extension of SCD, the advancement of SCOD and the delay in SCED. Precipitation, aspect and the slope all had little effect on snow phenology. Compared with the SCOD, the NDVI and longitude both had a greater impact on the SCED and SCD, while SCOD showed a greater impact from altitude.

(3)

The mean temperature was mainly negatively correlated with the SCD and SCED and mostly positively correlated with the SCOD. As the latitude increased, the snow phenology changed gradually, and the change rate in the SCD, SCOD and SCED in the whole Northeast China were 10.20 d/degree, −3.82 d/degree and 5.41 d/degree, respectively. The change rate in the snow phenology in forested areas and nonforested areas were inconsistent, and it was slower in forested areas than nonforested areas. Snow phenology was mainly positively correlated with the NDVI, but weak correlations with the NDVI accounted for a large proportion.