Variability in observed snow depth over China from 1960 to 2014

Snow depth is a critical variable that reflects snow variability and has an important impact on the energy and hydrological cycles. However, snow depth changes in response to regional climate warming have not been quantified in detail across China. Here, we investigated the variability in observed snow depth based on a recently released dataset collected at ~1,500 stations from 1960 to 2014. To address the non‐normal distribution of daily snow depth data, the annual cumulative amount of snow depth (CASD) is employed, and the trend is calculated and tested for the change point and field significance. The results showed that CASD experienced a slight increase of 1.75%·decade−1 in China. An increasing trend of CASD was observed in Northeastern and Northwestern China, a decreasing trend was observed in the eastern Tibetan Plateau, Loess Plateau, and Yangtze River Basin, and these changes were more remarkable after the 2000s. Since the 2000s, the increased CASD in Northwestern and Northeastern China has been accompanied by a significant decrease in winter temperature and a sharply increasing trend in snowfall events in the northeastern region. The decreased CASD in the Eastern Tibetan Plateau and Loess Plateau after the 2000s was accompanied by a decrease in snowfall events and an increase in winter temperature. Large‐scale atmospheric circulation, such as the Arctic Oscillation (AO), has a remarkable regional effect on the variability in snow depth. Statistically significant negative correlations were found in the northern part of China (i.e., positive AO phase with reduced CASD), while significant positive correlations were shown in the middle part of China (i.e., positive AO phase with enhanced CASD). Furthermore, China experienced an obviously shorter snow season during the whole period (decreased by 3.87%·decade−1), jointly resulting from later starting dates and earlier ending dates.


| INTRODUCTION
Snow plays an important role in the terrestrial energy budget and hydrological cycle. It can lower land surface temperatures and stabilize thermal conditions due to its high thermal resistance and low thermal conductivity (Barnett et al., 2005;Zhang, 2005). Snow can delay the transfer of precipitation to surface runoff, infiltration, and groundwater recharge (Barnett et al., 2005;Mote et al., 2005;Harpold et al., 2017). The water stored in snow pack is important for spring and early summer runoff, especially for a dry summer (Godsey et al., 2014). Under a warmer climate, a significant decreasing trend in snow cover has been found in the Northern Hemisphere (Brown and Mote, 2009;Brown and Robinson, 2011), associated with earlier snowmelt. The fast change is closely linked with streamflow, soil moisture, and vegetation. The snowmelt-driven streamflow increased in late winter and early spring (Barnett et al., 2005;Mote et al., 2005;Stewart et al., 2005;McCabe et al., 2007;Barnhart et al., 2016) and decreased in summer (Godsey et al., 2014). The soil moisture affected by melt water in spring and summer regulated vegetation growth (Peng et al., 2010). Thaw-refreeze events occurred in early and late winter (Wilson et al., 2013), which resulted in root damage and reduced plant productivity (Bokhorst et al., 2011).
Snow variability is closely related to environmental conditions such as winter air temperature (determining snow melting) and precipitation phase (determining snow accumulation) (Mote, 2003;Mote et al., 2005;Luce et al., 2014;Harpold et al., 2017). For example, an increasing trend in snow depth was detected over Western Eurasia (Kripalani et al., 2003a;Kripalani et al., 2003b;Popova, 2007), while a significant decreasing trend in snow water equivalent was found over Northern Eurasia with an increased winter air temperature of approximately 2 C (Ye et al., 2015). In China, the increasing trend of the total amount of snowfall was found in Northeast China (Sun and Wang, 2013;Wang and He, 2013), and Northern Xinjiang and Eastern Tibetan Plateau (Sun et al., 2010), but the decreasing trend in southeastern China (Sun et al., 2010;Zhou et al., 2018). Over the Western United States, snow depth showed a decreasing trend with an earlier snowmelt season (Dyer and Mote, 2006;Fyfe et al., 2017). The amount of snow and the timing of snowmelt events were affected by the warmer winter. The warming effect from the increased winter temperature frequently surmounted the effects from increased precipitation in the low-latitude area. However, in the colder areas, changes in snow were predominantly driven by precipitation, while the temperature effect was relatively small Mote, 2006). It was assumed that in cold winter regions, snow variations would be "immune" against to a warmer climate. The increased winter temperature seemingly was still low enough to have a negative effect on the snowmelt process.
The variability in snow is expected to be marked by geographical responses to climate warming in China due to its vast and complex terrain features, such as deserts, mountains, plateaus, and plains. In addition, China has also experienced obvious climate changes caused by global warming (Ren et al., 2012a;Ma et al., 2015;Zhou and Wang, 2017). The effect of climate change on snow cover has been investigated based on remote sensing data and observed data across China (Huang et al., 2016;Ke et al., 2016), but how geographical responses have impacted snow depth variability over the whole area of China has not been documented in detail. Furthermore, the warming surface temperature has slowed since approximately the 2000s (Gillett et al., 2012;Kosaka and Xie, 2013;Trenberth and Fasullo, 2013) and is especially pronounced in boreal winter Li et al., 2015). The response of snow depth variability in China to this slowdown in the warming trend remains unclear partly because the data have not been available until recently.
Several variables have been employed in previous studies to evaluate snow depth variability, such as the maximum snow depth, annual mean snow depth, and monthly mean snow depth in wintertime. These variables are more often used in snow stable regions, where the snow depth showed an appropriately normal distribution, as shown in Figure 1a. However, one concern is that the variable of annual mean value of snow depth is easily affected by extreme values, which can occur in snow unstable regions, with a non-normal distribution and are also associated with data gaps during the snow season, as shown in Figure 1d. An alternative presentation is to employ the annual cumulative amount of snow depth (CASD) with an "S"-shaped curve of observation. With a similar shape, the useful advantage of the "S"-curve is the easier comparison with the characteristics of an irregularly humped distribution (Borradaile, 2013). The application of the cumulative distribution has been successfully implemented to reflect the changing timing of streamflow from snowmelt when it is sensitive to short-term climate variability (Moore et al., 2007;Clow, 2009). Hence, the CASD could be used to help read the statistical characteristics of snow depth and reflect the variability in snow depth in both stable and unstable snow regions.
Both snowfall and snow melting are affected by the large-scale atmospheric circulations, and therefore could play a role in variability of CASD. The Arctic Oscillation (AO) pattern has been regarded as one of the dominant patterns of Northern Hemisphere climate variability (Thompson and Wallace, 1998), and played an important role in the winter weather in the mid to high latitudes. When the AO was in the negative phase, an anticyclonic circulation was found in Southwest China, and a cyclonic circulation was centred on Lake Baikal (Zhang et al., 2015). The anticyclonic circulation would favour the persistence of snow cover, as it usually provided very low temperatures mainly due to the strong nocturnal radiation cooling, which was cold enough to accumulate the snow depth. When the AO is in the positive phase in the winter, it could lead to warming conditions over East Asia (Thompson and Wallace, 2000), as the higher pressure of approximately 55 N would strengthen winds and blow ocean warmth onto continents, accompanied by negative precipitation anomalies in Northwest and Northeast China (He et al., 2017).
At decadal to multidecadal time scales, the Atlantic Multidecade Oscillation (AMO) and the Pacific Decadal Oscillation (PDO) play more important role (Li and Bates, 2007;Wang and He, 2013;Dai et al., 2015). The primary mechanism for the AMO exerting its influence on China is by intensifying midtroposphere and upper troposphere Eurasian heating. In winter, this heating tended to decrease the land-sea thermal contrast (Wang et al., 2009) and induced a surface low-pressure response over the North Atlantic extending to the Eurasian Continent, weakening the Mongolian Cold High and East Asian cold air activity (Wang et al., 2009) and consequently weakening the winter monsoon. During the negative phase of the PDO, positive pressure anomalies were located over the East Asia-Japan-North Pacific region, which weakened the East Asian trough. The weaker East Asian trough favoured southerly winds over North Central China, produced southerly advection of warm and moist air, leading to increased precipitation over North Central China (Qin et al., 2018). The East Asian winter monsoon (EAWM) governed the winter climate of China. The strong northerly and northeasterly winds over Northern and Central China brings intense cold air masses from the inner land to the lower latitude coastal regions which results in dramatic changes in air temperature (Gong et al., 2001;Wang and He, 2013). The EAWM is also modulated by atmospheric circulations in various ways. For instance, the EAWM was weakened when the AO entered the positive phase, which caused a warm winter climate in Northeastern China (Suo et al., 2009). When the AO and PDO were in the negative (positive) phases, the EAWM was strong (weak), thus leading to colder (warmer) temperatures in China . The AMO corresponded to the cold period and a strong winter monsoon; otherwise, the AMO corresponded to the warm period and weak winter monsoon, thus influencing the oscillation of the EAWM .
Here, the daily snow depth data released by the China Meteorology Administration cover the 1960-2014 period at more than 2,300 stations in China. Then, we employed a snow depth matrix to evaluate the detailed variability in snow depth across China, including the CASD, maximum snow depth, and snowfall events. Furthermore, we analysed how the CASD was affected by winter environmental conditions, such as winter temperature and snowfall events, and connected it with atmospheric circulations. Hence, the objectives in our study were (a) to summarize the spatial-temporal climatology and interannual variability in CASD and its related variables across China, (b) to analyse how climate change may adjust the relative roles of winter temperature and snowfall in driving CASD, and (c) to quantify the impact of the mentioned atmospheric circulation modes on CASD.

| Data and quality control
The daily snow depth was collected from the China National Stations' Fundamental Elements Datasets V3.0, obtained from the National Meteorological Information Center, China Meteorological Administration (CMA, http://www.cma.gov.cn/), included more than 2,300 meteorological stations (Figure 1). We chose daily observations from September 1, 1959 to August 31, 2014. According to the Specifications for Surface Meteorological Observations (China Metrological Administration, 2007), snow depth was recorded as a rounded-up centimetre integer and measured with a ruler at a daily frequency.
We excluded the non-snow sites from the following analysis. A non-snow site is defined by the annual mean value of snow depth that is less than 1 cm, referring to the criteria of (Li and Mi, 1983;He and Li, 2012). We did not consider the corresponding non-snow sites in the following analysis. There are 550 non-snow sites that are mainly located in Southern China. Furthermore, to reduce uncertainty caused by the short time period of observations, we consider stations with at least 40 years of records ( Figure 1). Quality control has been strictly implemented in the dataset, including resolving the problems of incorrect and missing data, restoring historical basic meteorological data, climatic and numerical range checks, internal consistency checks, main change range checks, and temporal and spatial consistency checks (Ren et al., 2012b). The quality control flag was provided in the dataset (i.e., flag = 0 means correct data). In our study, valid data with a flag of 0 were adopted. Records with other flags, such as missing data or failing to meet the quality control criteria, were excluded.
For each site, the missing data and flagged data should be less than 20% over 10 hydrological years (a hydrological year spanned from September 1 of the previous year to August 31 of the current year) and should be less than 20% each winter time in a hydrological year. After culling, we retained approximately 1,800 stations, among which 1,445 stations (81.45%) have 55-year records and 327 stations (18.4%) have 40-55-year records. The daily air temperature and precipitation records at the corresponding sites were also applied to the culling rule. Finally, 1,506 stations (81.8%) were retained with 55 year in air temperature records and 1,445 stations (78.2%) in precipitation, 337 stations (18.2%) with 40-55 year in temperature and 437 stations (23.1%) in precipitation.

| Matrix of snow depth quantification
The variability in observed snow depth was evaluated by a snow depth matrix, including the annual CASD, the maximum snow depth, the snowfall event, the length of snow season, the wintertime temperature, and the fraction of days with air temperatures less than 0 C during the snow season (designated F T < 0 C ). The variables described for each station were calculated in a hydrological year.
In this analysis, we employed the annual CASD and maximum snow depth in a hydrological year to jointly reflect the variability in snow depth in China. The CASD is the cumulative value of daily snow depth in one hydrological year. The variable CASD could be regarded as an indicator for the total snow water mass for a year, which could not only reflect the whole snow depth process but also include detailed information in the regions where snow depth showed large fluctuations. Specifically, the CASD can capture snowstorms or severe snow events, while the mean value of snow depth is insensitive (Gao and Yang, 2009). This method has been applied successfully to reflect the changing timing of streamflow from snowmelt (Moore et al., 2007;Clow, 2009) because the timing of streamflow is sensitive to short-term climate variability, such as fluctuations in precipitation (Stewart et al., 2005).
The snowfall event was estimated by a widely used method from Dai (2008) based on precipitation and air temperature. A snowfall event could be considered when the precipitation is larger than 0 cm and the snow occurrence frequency is more than 50%. The snow occurrence frequency was estimated based on a hyperbolic tangent function, which first considered the surface air temperature effect and then the elevation effect on the precipitation reflected in the parameters of the function (details can be found in figure 2 of Dai 2008). The snow season was defined as the duration between the first occurrence of snow depth that was larger than 1 cm and the last snow depth occurrence recorded in spring. The wintertime in this study referred to the period from December to February.
Atmospheric circulations and teleconnection modes, which have a potential influence on snow depth, were represented by several indices. The indices used in this analysis were obtained from publicly established available sources: the AO and the AMO index were obtained from the National Oceanic and Atmospheric Administration Climate Diagnostic Center (https://www.esrl.noaa.gov/psd/data/timeseries/ ), and the PDO index was obtained from the University of Washington's Joint Institute (http://jisao.washington. edu/pdo/PDO.latest). The EAWM was described by the index proposed by Wang and Chen (2014), based on the sea level pressure obtained from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis dataset 2.5 × 2.5 grid, which is available from 1948 to present.

| Trend analysis
We adopted the nonparametric Mann-Kendall (MK) test (Mann, 1945;Kendall, 1948) to detect the significance of trends in the time series of snow-related variables. The MK test is minimally affected by the unnormalized distribution of variables (Hamed, 2008). However, data should be assumed to be independent in the MK test. According to Von Storch (1999), autocorrelation would lead to a rejection of the null hypothesis of no trend when the null hypothesis is actually true. To eliminate this concern, the "trend-free pre-whitening" method proposed by Yue et al. (2002) and Yue and Wang (2004) was applied prior to the MK test to preserve the magnitude of the "real" trend. This combination can provide an accurate trend estimate of the autocorrelation process and hence is widely used in hydrological and meteorological time series (Dadaser-Celik and Cengiz, 2013; Gocic and Trajkovic, 2013;Tan et al., 2017). The slope is estimated using the Sen estimator (Sen, 1968), considering its robustness against outliers. In our study, significance levels of α = .01 and α = .05 are adopted.

| Field significance test
The field significance is important in large-scale studies to assess whether the trends occurred by chance. For a given data set, "field significance" is considered when a certain number of stations have local significance trends at a certain significance level, that is, the 5% level in our case (Livezey and Chen, 1983;Wilks, 1997). The local significance was tested by the method described in Section 3.1. However, local significance tests may suffer from the "multiplicity" problem, especially meteorological variables that are spatially interdependent (Wilks, 2011). To address this issue, we conducted a field significance test using the moving block bootstrap (Wilks, 1997). This method can incorporate the interdependence between stations and hence has been widely used for field significance tests in other regions (Kiktev et al., 2003;Alexander et al., 2006;Westra et al., 2012;Li et al., 2018) . The field significance test is calculated as follows: First, the null hypothesis assumed that the pattern of trends estimated from the actual station data was due to climate noise. Second, a matrix was formed with rows of time and columns of stations. Then, we resampled the matrix 1,000 times (with the replacement of rows) by moving block bootstrap with a block size of 5. Finally, we calculated the trend using the method described in Section 3.1 for each random matrix. We calculated the percentage of 1,000 bootstrapped field with a significance trend to form a distribution. If the percentage of stations with actual significance trends lay in the upper 5% tail of the distribution, the actual trend was considered to be field significant. Considering the importance of block length (Kiktev et al., 2003), our resampling result showed that a block length of more than 4 was robust, after checking for a length between 3 and 12, and found no significant change was found for a length of more than 4.

| Regional time series and change point test
Considering the large geographic extent and complex terrain in China, snow depth could have developed differently in different regions, especially in some featured and hot-spot regions. Hence, in our study, we selected five subregions, as shown in Figure  To describe the temporal variability in snow depth, we aggregated snow depth data and its related variables into the regional time series. First, we gridded the area of China onto 1 × 1 longitude-latitude grids and used the cosine value of the grid latitude as area weights instead of just applying a simple average to the variables due to the uneven distribution of the stations (shown in Figure 2). Any grid box with missing values was given zero weight. Then, the regional time series of variables were created by the method from Jones and Hulme (1996). For snow-related variables, the large gradient in mean and variance mainly exists among stations, so the weight standardized anomaly (also known as the "z score") (Equations (1) and (2)) was adapted to the regional averages: where S ij is the snow variable in year j at the ith station and S i and σ i are the mean and SD of the ith station, respectively. ΔS ij is the standard value in year j at the ith station; and ΔS R is the regional normalized data.
For winter temperature, the regional time series was conducted using the method in Equations (3) and (4), as the temperature was a continuous variable and the average-weight anomaly was appropriate: where T ij is the temperature in year j at the ith station, T i is the time mean, ΔT ij represents the anomalies of temperature, w i is the weight at the ith station, and ΔT R is the regional anomaly value.
Then, the regional time series of variables were transformed back to original units by denormalizing the snow depth-related time series (Equations (5) and (6)), and an absolute temperature for winter temperature was added (Equation (7)).
The change point of the regional time series was detected by a nonparametric method that is robust against outliers (Taylor, 2000). This method has been widely used to identify the timing of shifts in hydrological studies (Mavromatis and Stathis, 2011;Yang et al., 2012;Reid et al., 2016). The cumulative sum of the time series was first calculated based on Equation (8). A sudden change in the direction of the cumulative sum indicated that a change point occurred in the time series. The confidence level is determined by a 1,000 times bootstrap analysis to test if the change point happened by chance. The change point was significantly regarded with a confidence level of 95%.
where S i is the cumulative sum of the difference and S 0 = 0, X is the average value during the period of 1960-2014 (55 years).

| Climatology of snow depth
To illustrate the climatology of the observed snow depth in China, we analysed the mean value of the annual CASD, maximum snow depth value, snowfall events, and the length of snow season during the period from 1960 to 2014 ( Figure 3). Figure 3, the 55-year mean value of CASD was 229.0 cm. The CASD increased from south to north with increasing latitude. In the low-latitude and low-altitude regions, such as the Yangtze River Basin, the mean value of CASD ranged from 10 to 50 cm. Snowfall events in such regions were below 10 days per year. In the high altitudes of mountainous regions, the CASD reached above 200 cm, such as in the Northeast, Northwest, Tibetan Plateau, and Inner Mongolia Plateau. These regions experienced more frequent snowfall in winter, with a range of 20-30 days. The maximum observed snow depth in most of Northern China was approximately 15 cm ( Figure 3b). Noticeably, the CASD and maximum snow depth increased sharply in the regions above the latitude of 40 N (Figure 3a). This finding is consistent with the previous results that the snow depth and snow season north of 40 N were quite stable (Peng et al., 2010;Zhong et al., 2014).

As shown in
From Figure 3c, in accordance with the CASD pattern in the northern part of China, the Northeast, Northwest, Eastern Tibetan Plateau, and Inner Mongolia Plateau were the major snow stable regions with a length of snow season more than 90 days per year. For the North China Plain and Loess Plateau, the length of the snow season was approximately 10-60 days per year, which was unstable but showed an annual periodicity. The unstable snow regions were mainly in the southern part of China, such as the middle and lower parts of Yangtze River Basin and the southwestern part of China, with a very short duration of the snow season (1-10 days).

| Spatial patterns in the snow depth trends
During the 1960-2014 period, 647 stations showed an increasing trend in CASD, and 605 stations showed a decreasing trend (Figure 4a). The stations with increasing trends were mainly located in the northern part of China, and were especially significant in Northwestern and Northeastern China and in the Northeastern part of Inner Mongolia, while decreasing trends were found in middle China, such as the North China Plain and Loess Plateau. Slight decreasing trends were found in the CASD in the low part of the Yangtze River Basin. The histogram from Figure 4a shows that the statistically significant increasing trends in CASD were field significant.
To further interpret the spatial variation in CASD, the spatial pattern of winter climate conditions was jointly analysed, including the snowfall events, winter temperatures, and the fraction of days with temperatures below 0 C during the snow season, as shown in Figure 4. There were 261 stations with increasing trends in snowfall events, while 442 stations had decreasing trends. Increasing snowfall events were found in the northern part of China, while decreased snowfall events occurred in the Eastern Tibetan Plateau and North China Plain. It was deduced from Figure 4b that the observed trends in snowfall events were not field significant. As shown in Figure 4c, winter temperature experienced a widely increasing trend across China, which was field significant at the 5% level. The significant trends were mainly located in the northern and central parts of China. The trends in F T < 0 C are shown in Figure 4d. There were 364 stations with increasing trends and 775 stations with decreasing trends. The increasing trend in F T < 0 C in the northeast and northwest parts showed insignificance, while the trend in the central part and North China Plain showed a significant decreasing trend. The field significance test showed that the observed significantly decreasing trends in F T < 0 C were statistically significant.

| Temporal variation in regional mean snow depth
As shown in Figure 5, the CASD showed a slight increase of 1.75%Ádecade −1 across China during the 1960-2014 period. Overall, the winter temperature showed a F I G U R E 3 Spatial patterns in the climatological mean values for snow-related variables from 1960-2014: (a) annual cumulative amount of snow depth (CASD, cm), (b) duration of snow depth season (day), (c) annual maximum snow depth (cm), and (d) number of snowfall events (day). The blue lines summarize the annual CASD and length of the snow season with 1 latitude/longitude shown on the y and x-axes, respectively F I G U R E 4 The spatial distribution of snow-related indices and field significance. The left column represents trends in four variables from 1960-2014: (a) annual cumulative amount of snow depth (CASD, cm), (b) snowfall events (day/year), (c) winter temperature ( C) and (d) fraction of days with T < 0 C during the snow season from 1960 to 2014 in China. The red (blue) symbols in the left column indicate increasing (decreasing) trends, and the solid symbols are at a significance level of 10%. The middle and right columns represent the field significance test results. The histogram shows the percentage of stations with significantly increasing trends (middle column) and significantly decreasing trends (right column), which were calculated from 1,000 bootstrap resamples. Solid red (blue) circles represent significantly increasing (decreasing) results from the observed dataset. The grey dashed lines represent the position of the 5% tail of the distribution. When the circle is at the upper position of the 5% tail, the trends were regarded as the field significance significant increase (0.3 CÁdecade −1 ) and a slight decrease in F T < 0 C (−0.3%Ádecade −1 ). The mean value of winter temperature for the whole period was −2.04 C below the freezing point. The snowfall event showed little trend with large variability.
For subregions (Figure 6), the CASD showed a significantly increasing trend in the northwest and northeast at 9.81 and 7.96%Ádecade −1 , respectively. Slight decreasing trends were found in the other three subregions: −2.72%Á decade −1 in the middle and low parts of the Yangtze River Basin, −9.41%Ádecade −1 in the Eastern Tibetan Plateau, and −7.98%Ádecade −1 in the Loess Plateau, respectively. As shown in Figure 6 and Table 1, there were significant change points in the time series of CASD in the northwest (2000), northeast (1999) and Eastern Tibetan Plateau (1998). These change points occurred around the 2000s. Before the 2000s, the CASD showed an increasing trend, with increased snowfall events in these three regions. After the 2000s, the CASD increased more sharply in the northwest and northeast. There was a decreasing trend in snowfall events (−0.26 dayÁyear −1 ) and a decreasing trend in winter temperature (−0.11 CÁyear −1 ) in the northwest, whereas there was a slightly increased snowfall event (0.12 dayÁyear −1 ) and a decreased winter temperature (−0.12 CÁyear −1 ) in the northeast. For the Eastern Tibetan Plateau, the CASD showed a decreasing trend after the change point, with decreased snowfall events and increased winter temperatures.

| Relationship with the atmospheric circulation indices
To quantify the impact of large-scale atmospheric circulations on snow depth patterns in China, we examined the correlation between the CASD and major potential climate mode indices. Based on the correlation map in Figure 7, it can be considered that the large-scale patterns made the joint contribution to the snow depth variability pattern, as each of them addressed the impacts in different parts of China. Among those related atmospheric circulations, the correlation map of CASD with the AO (Figure 7a) showed the most similar pattern with the trend in CASD (Figure 4a).
A distinctive correlation pattern between the AO index and CASD was found. Specifically, in the northwest and northeast parts of China, the CASD has a negative correlation with the AO index, while in the Loess Plateau and North China Plain, the CASD has a positive correlation with the AO index. Noticeably, this pattern was consistent with the spatial trend in CASD in Figure 4. This finding indicated that when the AO was in the negative phase, the CASD increased in the northern part but decreased in the central part. In contrast, AO entered its positive phase, and snow depth decreased in the northern part but increased in the central part. The previous studies are also consistent with our results. They pointed out that the positive winter AO is associated with positive precipitation anomalies, especially in central China, and the negative precipitation anomalies in Northwest and Northeast China (Gong and Wang, 2003;Song et al., 2019). In addition, the studies on the relationship of AO and winter precipitation from Wen et al. (2009) and Matsuo and Heki (2012) showed similar results to ours.
The AO pattern was regarded as one of the dominant patterns of Northern Hemisphere climate variability (Thompson and Wallace, 1998), and played an important role in the winter weather in the mid to high latitudes. The winter climate was also significantly influenced by the large decadal to multidecadal climate variations, such as the AMO and PDO (Li and Bates, 2007;Wang and He, 2013;Dai et al., 2015). The correlation pattern between the AMO and CASD was significantly negatively correlated over the north plain of China but positively correlated in the Eastern Tibetan Plateau (Figure 7b). A decreasing trend in the CASD was obvious in the north plain of China when the AMO entered the warm phase. Meanwhile, an increased CASD in Eastern Tibetan Plateau was influenced by a warm AMO phase. Figure 7c showed that the CASD variability was significantly correlated with the PDO pattern in the central part of China, such as the Loess Plateau. A decreased CASD pattern was partly influenced by the PDO positive pattern. The EAWM is governed by the winter climate of China. The map from Figure 7d also illustrated a similar correlation pattern between the CASD and EAWM but was more significant on the northern coast area. It concluded that the pattern of snow depth was influenced by the circulations on different parts of China.

| Trends of the snow season length
In China, the length of the snow season showed an obvious decrease of 3.87%Ádecade −1 during the whole period. Approximately 61% of stations showed obvious trends in the duration of the snow season ( Figure 8) from 1960 to 2014, with field significance. A total of 1,067 stations showed a negative trend, and 30% of these stations showed a significant trend. Significant negative trends were mainly observed in Northwestern and Northeastern China, the Loess Plateau, and the Eastern Tibetan Plateau. Furthermore, the earlier ending dates and the later starting dates of snowpack resulted in a shorter period of snow depth (Figure 8b,c). 80% of all stations (971 stations) showed a negative trend on the last day of snow cover, and 16% of these trends were significant, which means that the last date of snow cover occurred earlier. Such trends were found in areas of China with frequent snowfall, such as the Loess Plateau, North China Plain, Northeastern and Northwestern China, and the eastern part of the Tibetan Plateau. Furthermore, 741 stations (60%) presented a positive trend of first day snow cover, of which 18% were significant. Most stations with significantly positive trends were in Northwestern and Northeastern China, the Loess Plateau, and the Tibetan Plateau, which meant a later starting date during the entire snow depth cover period. Hence, most parts of China experienced an obviously shorter snow season from 1960 to 2014, which was jointly determined by later starting dates and earlier end dates.

| DISCUSSION
Large variability in snow depth and the length of the snow season have been documented across China. Different patterns of snow depth variability have been shown in unstable and stable snow regions, responding to fluctuations in snow amount and temperature during the snow season. These findings were consistent with previous studies, although the variability in snow-related variables was based on multi-resource data. For example, in the Tibetan Plateau, the CASD showed an initial increase before 1980 and was subsequently followed by a decreasing trend. A similar pattern was also found based on the annual mean snow depth (Xu et al., 2016), and a small increase in snow cover was reported between 1951 and 1997 by Qin et al. (2006) using observations and remote sensing data. This small difference may be caused by the coverage of different time periods. We found that the CASD showed an evident increase in the northern part of China, especially in the northeast and northwest, which was also reported by Zhong et al. (2018) based on the variable of annual mean snow depth, with a greater increasing rate in higher-latitude regions.
The spatial pattern of snow cover days or snow season was complex in China, with an earlier end date and later The statistical information of four variable trends for the whole period and the mean value for two phases in five subregions Note: The bold numbers are trends at the 10% significance level. The term "trend" refers to the trend in each variable during the period of 1960-2014. The term "mean before" refers to the mean value of the period before the change point, and the term "mean after" refers to the mean value of the period after the change point. The four variables are cumulative amount of snow depth (CASD), snowfall event, winter temperature and F T < 0 C . The five subregions are Northwest China (NW), Northeast China (NE), the eastern Tibetan plateau (ETP), the loess plateau (LP), middle and below Yangtze River basin (YR).
start date, which was consistent with the conclusions from Ke et al. (2016). In addition, a decreasing trend in snow cover days was also found in the Loess Plateau based on MODIS data from 2003 to 2008 (Jin et al., 2015). Snow cover days increased in South China but decreased in the Tibetan Plateau from 2000 to 2014 based on remote sensing snow depth data (Huang et al., 2016). These regional characteristics were largely consistent with our findings (shown in Figure 8).
In this study, we used the snowfall event, winter temperature and F T < 0 C during the snow season, which resulted in suitable factors to better understand the fluctuation in snow depth. During the study period, the magnitude of the trend of CASD was different in subregions but was largely consistent with the trends in snowfall events. The snowfall events showed increasing trends in northeastern and northwestern regions, and also in Eastern Tibetan Plateau, in agreement with the previous studies (Liu et al., 2012;Liu et al., 2013;Wang and Zhou, 2018). The decreasing trend of snowfall events was found in Yangtze River Basin regions that the similar result was also shown in . Winter temperature showed decreasing trends after the 2000s in the northeast (−0.12 CÁyear −1 ) and northwest (−0.11 CÁyear −1 ). For the Eastern Tibetan Plateau, and even though the winter temperature increased after the 2000s (0.04 CÁyear −1 ), the CASD showed a decrease, which was consistent with a decreasing trend in snowfall events. Noticeably, the CASD showed a significant increase since the 2000s in the northwest and northeast and a decrease in the Eastern Tibetan Plateau, and the magnitude of winter temperature changes in different directions, with trends of −0.11, −0.12, and 0.04 CÁyear −1 , respectively. The mean value of temperature in those regions remained below the freezing point, which was as low as −14.06, −9.05, and −4.9 C, respectively, and the number of days with T < 0 C exhibited an increasing trend for the whole snow season, which was in favour of snow depth accumulation. Hence, we found that the increased CASD anomaly was jointly affected by winter temperature and snowfall in Northern China, but we assumed that winter precipitation could be the main F I G U R E 8 The spatial distribution of snow cover trends and the field significance test. The left column represents (a) the trend in the length of the snow season (day/a), (b) the trend in the first day of the snow season (day/a) and (c) the trend in the last day of the snow season (day/a). In the left column, red (blue) symbols indicate an increasing (decreasing) trend, and the solid symbols are at a significance level of 10%. In the middle and right columns, the histogram represents the percentage of stations with a significantly increasing trend (the middle column) and a significantly decreasing trend (the right column), based on 1,000 bootstrap resamples. The grey dashed lines represent the position of the 5% tail of the distribution. Solid red (blue) circles represent increasing (decreasing) significance test results from the observational dataset. When solid circles lay in the upper 5% tail of the distribution, the actual trend was field significant factor influencing snowpack variability at locations that were above freezing, as stated in previous studies (Mote, 2006;Brown and Robinson, 2011). In other words, when the temperature was below the freezing point, intensified precipitation would result in increasing snow depth (Ye et al., 1998;Qin et al., 2009;Peng et al., 2010). Therefore, the CASD was sensitive to the snowfall amount in those subregions where the winter temperature was below the freezing point.
Variability in the atmospheric pattern or teleconnection showed remarkable regional effects in the fluctuation in snow depth pattern (Clark et al., 2001;Ge and Gong, 2008). These atmospheric circulations caused regional variations in snow depth by altering temperature, precipitation and the dynamic and thermal effects of the atmosphere (Seager et al., 2010). From Figure 7, the trend in CASD in China was well matched with the pattern of the AO, which was consistent with previous studies (Zhang et al., 2015;Douville et al., 2017).
In general, the physical mechanisms of atmospheric circulation could force the interannual variability in snow depth. However, the impacts of circulation patterns were too complex to be easily explained by single climate atmospheric circulation or determined by relying on only observation data; rather, the circulation patterns were interacting and jointly influencing the snow depth variation and subsequently impacted the winter climate in China. In conclusion, to investigate the connection between the variability in snow depth and climate circulations in winter, numerical climate model experiments may be a more appropriate method beyond observations.

| CONCLUSIONS
In our study, the spatial pattern and temporal variability in observed snow depth and its response to regional climate change were investigated across China and its subregions from 1960 to 2014 based on a recent observation data set collected at approximately 1800 stations. A snow depth matrix was employed, including the annual CASD, maximum snow depth, snowfall event, length of snow season associated with winter temperature and fraction of days with T < 0 C during the snow season (referred to as F T < 0 C ). In addition, the relationships between snow depth and climate circulation modes were analysed.
The climatology of the CASD showed an increase at higher latitudes. The increase was much sharper at latitudes above 40 N with a maximum snow depth of approximately 15 cm. Regions located above the latitude of 40 N featured stable snow conditions, with an average snow season of more than 120 days. The periodical snow regions are mainly located in middle and low-latitude areas, such as the Loess Plateau and North China Plain. The Yangtze River basin was marked as an unstable region with a snow season shorter than 10 days.
The CASD showed a slight increase of 1.75%Ádecade −1 over China during the 1960-2014 period. It also showed spatial heterogeneity in snow-featured regions: at high latitudes, such as Northwestern and Northeastern China, the CASD increased significantly by 9.81 and 7.96%Á decade −1 , respectively, whereas in the Eastern Tibetan Plateau and Loess Plateau and Yangtze River Basin, the CASD showed the decreasing trends of −9.41, −7.98, and −2.72%Ádecade −1 , respectively. Those changes were especially evident since the 2000s: in Northwestern and Northeastern China, the increasing trends in CASD corresponded to an evidently decreased winter temperature and an increased F T < 0 C , even with different trends of snowfall events after the 2000s. In the Eastern Tibetan Plateau and Loess Plateau, the decreasing trends of in CASD were consistent with a decreased snowfall event and an increased winter temperature. The snow depth in the Yangtze River basin was influenced by the variability in snowfall, as the average winter temperature was above the freezing point. Additionally, the length of the snow season showed a significant negative trend (−3.87%Á decade −1 ) for the majority of China and was seen in conjunction with the later beginning dates of the snow depth season and the earlier ending dates of the snow depth season.
The analysis of the snow depth-climate relationship indicated that the AO pattern was correlated most with the variability in snow depth. Significant negative correlations were found in the northern part of China, while significant positive correlations were found in the middle part of China. The distribution pattern of snow depth variability was not controlled by a single climate variable but by a combination of multiple circulation variables. Our results provide a direct and detailed view of observed snow depth variation and its pattern over China and hence form a basis for the understanding of its response to a warmer climate.