Spatiotemporal variation of snow depth in the Northern Hemisphere from 1992 to 2016

Abstract. Snow cover is an effective indicator of climate change due to its impact on regional and global surface energy and water balance, and thus also weather and climate, hydrological processes and water resources, and the ecosystem as a whole. The overall objective of this study is to investigate changes and variations of snow depth and snow mass over the Northern Hemisphere from 1992 to 2016. We developed a long term Northern Hemisphere daily snow depth and snow water equivalent product (NHSnow) by applying the support vector regression snow depth retrieval algorithm, using passive microwave remote sensing data from the period. NHSnow product was evaluated along with the other two snow cover products (GlobSnow and ERA-Interim/Land) for its accuracy across the Northern Hemisphere. The evaluation results show that NHSnow performs comparably well with relatively high accuracy (bias: −0.59 cm, mean absolute error: 15.12 cm, and root mean square error: 20.11 cm) when benchmarked against the station snow depth measurements. Further analyses were conducted across the Northern Hemisphere using snow depth, snow mass, and snow cover days as indices. Analysis results show that annual average snow mass have a significant declining trend, with a rate of about 19.72 km3 yr.−1 or a 13 % reduction in snow mass. Although spatial variation pattern of snow depth and snow cover days exhibited slight regional differences, they generally reveal the decreasing trend over the most area of the Northern Hemisphere. Our work provides evidence that rapid changes in snow depth and snow mass are occurring since the beginning of the 21st century, accompanied by dramatic climate warming.



Introduction
Seasonal snow cover is an important component of the climate system and global water cycle, and has significant impacts on the surface energy, hydrological processes and water resources, heat exchange between the ground surface and the atmosphere, and the ecosystem as a whole [1][2][3][4].
have been explored to retrieve SD [41][42][43][44][45]. However, some limitations for these retrieval algorithms that use only passive microwave brightness temperature remain, due to the diversity of the land cover types and the spatiotemporal heterogeneity of the snow properties.
Numerous studies have reported the characteristics of snow cover variability at regional or hemispheric scales [20,21,[46][47][48][49]. For example, Huang et al. [50] reported the impact of climate and elevation on snow cover variability in the Qinghai-Tibetan Plateau, including SWE, snow cover extent (SCE), and snow cover days. Hori et al. [51] developed a 38-year Northern Hemisphere daily SCE product and analyzed the seasonal Northern Hemisphere SCE variation trends. SD, which provides additional information to characterize snow cover variation, is selected as the basis of analyzing the spatiotemporal change of snow cover. Barrett et al. [52] explored intra-seasonal variability in the springtime Northern Hemisphere daily SD change in the phase of the Madden-Julian oscillation. Wegmann et al. [53] compared four long-term reanalysis datasets with Russian SD observation data; however, this study only focused on the snowfall season (October and November) and snowmelt month (April).
A previous study [24] developed an SD retrieval algorithm with a machine learning method using several source datasets over the Eurasia area, showing this algorithm to have great advantages in snow depth estimation when compared with the other four common SD retrieval methods (e.g., Chang algorithm, spectral polarization difference algorithm, artificial neural networks algorithm, and common linear regression algorithm). The primary objective of this study is to generate daily hemispherical SD datasets (hereafter referred to as the NHSnow) over 25 years (1992-2016) with a 25-km spatial resolution using the support vector regression (SVR) SD retrieval algorithm [24]. This study attempts to address the following questions: (1) How consistent are NHSnow and other sources snow cover datasets with in-situ SD measurements? (2) What is the spatiotemporal variability of SD in the Northern Hemisphere from 1992-2016? Therefore, another purpose of this paper is to provide a comprehensive description of two snow characteristics (SD and snow cover days) from 1992 to 2016.

Passive Microwave Data
SSM/I and SSMIS are the PM radiometers on board the United States Defense Meteorological Satellites Program (DMSP) satellite (data available from the National Snow and Ice Data Center, http://nsidc.org/data/NSIDC-0032). The SSM/I (F11 and F13) and SSMIS (F17) datasets from this platform have been generated into the Equal Area Scalable Earth grid (EASE-Grid) with a 25-km resolution [54][55][56][57] (Table 1). The SCE and SD derived from data of the F11 and F13 sensors have a high consistency, rendering the calibration between these two sensors for snow cover area and SD unnecessary [58]. To minimize the melt-water effect, which can change the microwave emissivity of snow, only the descending orbit (night-time) passive microwave data were used [37].

Ground-Based Data
Daily ground-based SD measurements from two sources were used to construct and verify the SD retrieval model in this study. The first dataset is the Global Surface Summary of the Day (GSOD) dataset provided by the National Oceanic and Atmospheric Administration (NOAA) (https://data. To supplement the station data that were not reported during the period of 1992 to 2016, groundbased measurements of the daily SD were gathered from an additional 635 Chinese meteorological stations available at the National Meteorological Information of the China Meteorological Administration (http://data.cma.cn/en) [24,60]. These daily records, collected since 1957, include SD in centimeters, observation time, and geographical location information.

Topographic and Land Cover Data
This study also used the topography as auxiliary information to estimate the SD. The elevation was available from ETOPO1 at a resolution of 1 arc-minute [61], found at http://www.ngdc.noaa.gov/mgg/global/. To match the resolution of the PM brightness temperature data with the 25-km spatial resolution, we resampled the ETOPO1 to a 25-km resolution ( Figure 1).
To increase the accuracy of the SD estimates for different land cover types, we used both the MODIS land cover (MCD12Q1 V051) from 2001 to 2013 [62,63], and the Advanced Very High Resolution Radiometer (AVHRR) Global Land Cover classification generated by the University of Maryland Department of Geography. The MCD12Q1 International Geosphere-Biosphere Program (IGBP) classification scheme divides the land surface into 17 types, which were reclassified into five classes, according to Xiao et al. [24].
The AVHRR imagery from the NOAA-15 satellite, acquired from 1981 to 1994 [64], was categorized into 14 land cover classes at a 1-km resolution. These data allowed us to adjust the proposed snow-depth retrieval algorithm by reclassifying the 14 native land cover classes into five classes (water, forest, shrub, prairie, and bare land) at a 25-km spatial resolution (Table A1). The MCD12Q1 is available at https://earthdata.nasa.gov/, while the AVHRR land cover data are available online (http://www.landcover.org/data/landcover/).

Snow Cover Datasets
Two kinds of snow cover datasets were utilized, based on the following two criteria: covering the Northern Hemisphere and long-term availability. We selected GlobSnow and ERA-Interim/Land, which are widely used in global and regional climate change studies [22,65,66]. These datasets were compared with the NHSnow SD data. In the data pair extraction process, the data pairs were made by choosing close (or closest) observation/measurement times.
In November 2013, the European Space Agency (ESA) released the GlobSnow Version 2.0 snow water equivalent (SWE) and snow extent (SE) data for the Northern Hemisphere [17,67]. These data include all non-mountainous areas in the Northern Hemisphere and are available online To supplement the station data that were not reported during the period of 1992 to 2016, ground-based measurements of the daily SD were gathered from an additional 635 Chinese meteorological stations available at the National Meteorological Information of the China Meteorological Administration (http://data.cma.cn/en) [24,60]. These daily records, collected since 1957, include SD in centimeters, observation time, and geographical location information.

Topographic and Land Cover Data
This study also used the topography as auxiliary information to estimate the SD. The elevation was available from ETOPO1 at a resolution of 1 arc-minute [61], found at http://www.ngdc.noaa.gov/ mgg/global/. To match the resolution of the PM brightness temperature data with the 25-km spatial resolution, we resampled the ETOPO1 to a 25-km resolution ( Figure 1).
To increase the accuracy of the SD estimates for different land cover types, we used both the MODIS land cover (MCD12Q1 V051) from 2001 to 2013 [62,63], and the Advanced Very High Resolution Radiometer (AVHRR) Global Land Cover classification generated by the University of Maryland Department of Geography. The MCD12Q1 International Geosphere-Biosphere Program (IGBP) classification scheme divides the land surface into 17 types, which were reclassified into five classes, according to Xiao et al. [24].
The AVHRR imagery from the NOAA-15 satellite, acquired from 1981 to 1994 [64], was categorized into 14 land cover classes at a 1-km resolution. These data allowed us to adjust the proposed snow-depth retrieval algorithm by reclassifying the 14 native land cover classes into five classes (water, forest, shrub, prairie, and bare land) at a 25-km spatial resolution (Table A1). The MCD12Q1 is available at https://earthdata.nasa.gov/, while the AVHRR land cover data are available online (http://www.landcover.org/data/landcover/).

Snow Cover Datasets
Two kinds of snow cover datasets were utilized, based on the following two criteria: covering the Northern Hemisphere and long-term availability. We selected GlobSnow and ERA-Interim/Land, which are widely used in global and regional climate change studies [22,65,66]. These datasets were compared with the NHSnow SD data. In the data pair extraction process, the data pairs were made by choosing close (or closest) observation/measurement times.
In November 2013, the European Space Agency (ESA) released the GlobSnow Version 2.0 snow water equivalent (SWE) and snow extent (SE) data for the Northern Hemisphere [17,67]. These data include all non-mountainous areas in the Northern Hemisphere and are available online (http://www. globsnow.info/). Processing includes data assimilation based on combining the satellite PM remote sensing data (SMMR, SSM/I, and SSMIS), spanning December 1979 to May 2016, with ground-based observation data in a data assimilation scheme to derive SWE. GlobSnow Version 2.0 (hereinafter referred as GlobSnow) provides the following three kinds of temporal aggregation level products with a 25-km spatial resolution: daily, weekly, and monthly. This dataset covers all land surface areas in a band between 35 • N~85 • N, excluding mountainous regions, glaciers, and Greenland. In order to convert between SD and SWE using GlobSnow, the snow density is held constant at 0.24 g/cm 3 [33,65,68].
ERA-Interim/Land [69] is a global land-surface reanalysis product with data from January 1979 to December 2010 based on ERA-Interim meteorological forcing. It is produced by a land-surface model simulation using the Hydrology Tiled ECMWF Scheme of Surface Exchange over Land (HTESSEL), with meteorological forcing from ERA-Interim. Dutra et al. [70] described the snow scheme and demonstrated the verification using field experiments. SWE, with units in m of water equivalent, is one of the 13 parameters provided. We converted SWE to SD using the associated snow density data. These two datasets are available online (http://apps.ecmwf.int/datasets/data/interim-land/type=an/). The spatial resolution of the reanalysis data used is 0.125 degrees.

Theoretical Basis
Snow distribution is affected by various factors, which include, but are not limited to, the following: vegetation [33,34], soil and air temperature [21,71,72], and topography, and wind [73,74]. The snow retrieval process uses a range of variables to yield snow parameters (Equation (1)) [24].
where S is the snowpack properties (e.g., snow grain size (SD)), g (·) denotes the retrieval function, a is the atmosphere factors (e.g., wind speed and air temperature), T is the topography factors (e.g., elevation, slope, and aspect), G is the ground surface environment factors (e.g., surface temperature and vegetation type), L is the location factors (latitude and longitude), and DS is the digital signal from the remote sensing sensor (e.g., PM and optical remote sensing). D is the time factor and ε is the residual error or uncertainty (the difference between the sensor observation and the ground measurement). This retrieval formula can be exemplified by general research, in which the SD variable (S) is usually derived from the PM data (DS) using a linear regression method (g(·)) [15].

Processing Flow Overview
A previous study of the machine learning SVR method for SD retrieval models performed well with reduced uncertainties compared to the measured data from ground stations across the Eurasia region, based on analyzing the correlation coefficient (R), mean absolute error (MAE), and root mean squared error (RMSE) [24]. Here, we intended to use daily observational datasets over the Northern Hemisphere by applying the SVR method (except for July and August). The SVR SD retrieval algorithm uses a non-linear regression (support vector regression, SVR) method as the retrieval function (g (·) in Equation (1)). This machine learning (SVR) method was applied to establish the complex relationship between the target variable and input variables. Following Equation (1), we used ten variables as the inputs, including PM brightness temperature (19 GHz, 37 GHz, and 85 or 91 GHz) with vertical and horizontal polarizations, geographical location (latitude and longitude), elevation, and the measured SD. This algorithm also indirectly considered seasonal variation (day of the year; D) and vegetation (land cover; L) influence factors to improve the accuracy of the estimated SD. The output parameter is Remote Sens. 2020, 12, 2728 6 of 25 the estimated SD. This SVR SD retrieval algorithm mainly involves five steps (Steps 1-5). Steps 6-8 are the SD dataset generation processes ( Figure 2). Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 28 Figure 2. Process flowchart diagram for developing the Northern Hemisphere's daily snow depth and snow water equivalent data.
Step 1. Involves pre-processing meteorological station SD measurements, PM brightness temperature data, and other auxiliary datasets (DEM, land cover, and spatiotemporal information). Before estimating the SD from the PM data, it is necessary to identify snow cover and dry snow by a set of criteria [75]; in the training stage, it should extract and match the station measurements, PM brightness temperature data (classified as snow cover), and other auxiliary information; and in the predicting stage, the PM brightness temperature and auxiliary information were preprocessed to the datasets (Step 2). Step 1. Involves pre-processing meteorological station SD measurements, PM brightness temperature data, and other auxiliary datasets (DEM, land cover, and spatiotemporal information). Before estimating the SD from the PM data, it is necessary to identify snow cover and dry snow by a set of criteria [75]; in the training stage, it should extract and match the station measurements, PM brightness temperature data (classified as snow cover), and other auxiliary information; and in the predicting stage, the PM brightness temperature and auxiliary information were preprocessed to the datasets (Step 2).
Step 3. we developed SD retrieval models based on four land cover types to segregate the land cover effects. For the years over which our study period pre-dated the MODIS data, we used AVHRR land cover as the supplement data. We reclassified MODIS and AVHRR land cover into four classes (forest, shrub, prairie, and bare-land), which were the bases for constructing the SD retrieval sub-model. Table A1 (in the Appendix A) describes the reclassification scheme for AVHRR land cover. The MODIS land cover reclassification schemes were documented in Xiao et al. [24]. Because of the relatively stable change in land cover, the MODIS land cover in 2013 was used for each year from 2013 to 2016. Similarly, the MODIS land cover in 2001 was used for each year from 1998 to 2001, and the AVHRR land cover data was used for each year from 1992 to 1997. The retrieval sub-models were established on different land cover types (forest, shrub, prairie, and bare-land).
Step 4. we converted the day of the year (D) effects into the following three snow cover stages: snow accumulation stage (September to November), stabilization stage (December to February), and ablation stage (March to June). To consider the evolution effect of the snow properties, the SD retrieval model was established based on these three snow cover stages.
Step 5. we chose SVR as the retrieval function (Equation (1)), with specific kernel functions and parameters [24]; the optimization of parameters for the selected method were performed in the training stage to improve the model performance in this study.
Step 5.1. construction of a sub-continental model. In this study, we separately constructed the SD retrieval models for Eurasia and North America, based on the information that the snow properties show a discrepancy between Eurasia and North America. Taking snow density as an example, Bilello [76] pointed out that the mean snow density in the former Soviet Union (0.21~0.31 g/cm 3 ) was lower than that in North America (0.24~0.31 g/cm 3 ). Zhong et al. [77] explained the reasons for such a difference between Eurasia and North America.
Step 5.2. selection of training sample. The accuracy of the estimated SD primarily depended on the quality of the training samples [24]. More data than needed in the training stage to train the SD retrieval model may lead to over-fitting and yield an estimated SD with a high error. In this study, we collected numerous daily records for over 25 years. A changed sample selection rule was used in this study to avoid data information redundancy, and was divided into two steps. First, the number of samples of the three groups split by snow depth values were solidly quantified, i.e., group 1 (0 cm ≤ SD < 50 cm; shallow snow), group 2 (50 cm ≤ SD < 100 cm; intermediate depth), and group 3 (SD ≥ 100 cm; deep snow). To avoid an inflated training sample in group 2 and group 3, we set a threshold value number of used datasets (3000) determined by several tests (not shown). The threshold of the number of samples (12,000) for group 1 was adopted. Table 2 details the selection rules for the training sample for each group. Second, the quality of the training sample in each group was controlled using stratified random sampling. Stratification was performed at 1-cm SD intervals, and in each interval value of the datasets, we operated random selection to choose the needed numbers of samples. Step 6. Through the above steps, we generally created 24 SD retrieval sub-models (two (continents) * four (land classes) * and three (stages)), and then the preprocessed PM brightness temperature datasets were used to predict daily SD across the Northern Hemisphere from January 1992 to December 2016 (excluding July and August).
Step 7. The generated raw daily SD data were then processed to the daily SD map. If one pixel was detected as snow cover by the detection decision tree [75], but was likely to be shallow snow with an estimated value of equal or less than 1 cm in the raw predicted SD data, then the SD value was set as 5 cm [33,78] (Figure 2). Owing to the radiometer observations, the generated SD data were only reliable in areas with seasonal dry snow cover. Areas with sporadic wet or thin snow were not reliably detected, and areas marked as snow-free may have included areas with wet snow.
Step 8. This study excluded the Greenland and Iceland regions from the generation and analysis of the NHSnow datasets (NH_SD), because of the difficulty in discriminating snow from ice ( Figure 2) [19]. Missing data and zero-data gaps occurred when generating the daily SD grid. We applied the following filter: the daily estimated SD was defined as the midpoint of a sliding seven-day average window in order to reduce noise and compensate for missing data in the daily time series. For example, the SD estimated for 4 January was an average of the assimilated scheme output for 1 to 7 January [17,33]. When finished, the sliding SD method generated full coverage daily SD datasets for the entire Northern Hemisphere ( Figure 2).

Analysis Index Description
We used the seasonal maximum SD, seasonal average SD, and annual average SD as the SD analysis indexes to characterize the SD variability across the Northern Hemisphere. The seasonal maximum SD of a pixel refers to the maximum SD for a season (autumn: September-November; winter: December-February; spring: March-May). The sasonal average SD is defined as the cumulative SD divided by the days in one season (Equation (2)). The calculation of the annual average SD is similar to the seasonal average SD (Equation (2)).
where n is the number of days of one specific period, such as, a snow cover year (September-June) or a season (autumn, winter, or spring), and i is ith day of this specific period. SD is snow depth. The snow cover detection of the NHSnow SD data was based on Grody's snow cover identification method [75]. This means that SD was derived when the ground was detected as snow cover. Therefore, snow cover days (SCDs) were defined as the number of days that the SD was over 0 cm in a snow cover year (from September 1 to June 30). When analyzing the SD and SCD variation from 1992 to 2016, we divided the variation trends into five grades (highly significant decrease, significant decrease, non-significant change, significant increase, and highly significant increase; Table 3), and the significance of the trends was based on a t-test.

Validation of Snow Depth Dataset
To provide insight into the relative performance of the SD datasets, we compared three sources of snow cover datasets (NHSnow, GlobSnow, and ERA-Interim/Land) with ground SD measurements (Figures 3-8), using three indices (bias, mean absolute error (MAE), and root mean squire error (RMSE)). We collected daily SD measurements for the common period (1992-2010) of the three datasets as the validation data. This primarily focused on the snow stabilization stage (December to February). As the snow density changes slowly over a smaller range in the snow cover stabilization stage, we used a constant value (0.24 g/cm 3 ) for GlobSnow. Note that the data used for validation in this section are completely independent of the dataset used for the SD retrieval models.

Validation of Snow Depth Dataset
To provide insight into the relative performance of the SD datasets, we compared three sources of snow cover datasets (NHSnow, GlobSnow, and ERA-Interim/Land) with ground SD measurements (Figures 3-8), using three indices (bias, mean absolute error (MAE), and root mean squire error (RMSE)). We collected daily SD measurements for the common period (1992-2010) of the three datasets as the validation data. This primarily focused on the snow stabilization stage (December to February). As the snow density changes slowly over a smaller range in the snow cover stabilization stage, we used a constant value (0.24 g/cm 3 ) for GlobSnow. Note that the data used for validation in this section are completely independent of the dataset used for the SD retrieval models.        A small bias (blue and green dots in Figure 3) between the estimated SD against the measured SD was found for the mid and low latitude regions (<60° N) for the three SD datasets. However, a large bias was found for the Arctic region and along the coast, such as the Russian coastal regions, the Russian Far East, the Korean peninsula region, and Northeast Canada. For NHSnow and GlobSnow, bias was mostly distributed near the mean value μ = 0 line with a high frequency, although some bias was greater than 100 cm (or less than −100 cm; Figure 3d,e). Positive (negative) biases indicate that the mean estimated values are greater (less) than the corresponding measured SD values. ERA-Interim/Land overestimated the snow depth in Western Siberian Plains and Eastern European Plains (around 60° N; red dots, Figure 3c).
For the analysis indexes, MAE, and RMSE, the distribution of error points of NHSnow and GlobSnow were almost the same as the distribution of its bias (Figures 3-5). We used all of the valid daily records for each evaluation station to calculate three precision indexes for the three snow cover datasets, and found that the bias, MAE, and RMSE were −0.59 cm, 15.98 cm, and 20.11 cm, respectively, for the NHSnow grid datasets. For GlobSnow, however, there was more bias (−1.19 cm), and a lower MAE (15.12 cm) and RMSE (15.48 cm; Table 4). This comparison (NHSnow vs. GlobSnow) showed a relatively good agreement, although NHSnow over-or under-estimated the SD with a larger RMSE. Overall, the performance of GlobSnow was better than the NHSnow grid data. However, part of the validation data was also applied for the GlobSnow assimilation; it is highly possible that, in this case, the GlobSnow validation may not have been completely independent. The grain size was more important than the snow density and temperature in the microwave emission A small bias (blue and green dots in Figure 3) between the estimated SD against the measured SD was found for the mid and low latitude regions (<60 • N) for the three SD datasets. However, a large bias was found for the Arctic region and along the coast, such as the Russian coastal regions, the Russian Far East, the Korean peninsula region, and Northeast Canada. For NHSnow and GlobSnow, bias was mostly distributed near the mean value µ = 0 line with a high frequency, although some bias was greater than 100 cm (or less than −100 cm; Figure 3d,e). Positive (negative) biases indicate that the mean estimated values are greater (less) than the corresponding measured SD values. ERA-Interim/Land overestimated the snow depth in Western Siberian Plains and Eastern European Plains (around 60 • N; red dots, Figure 3c).
For the analysis indexes, MAE, and RMSE, the distribution of error points of NHSnow and GlobSnow were almost the same as the distribution of its bias (Figures 3-5). We used all of the valid daily records for each evaluation station to calculate three precision indexes for the three snow cover datasets, and found that the bias, MAE, and RMSE were −0.59 cm, 15.98 cm, and 20.11 cm, respectively, for the NHSnow grid datasets. For GlobSnow, however, there was more bias (−1.19 cm), and a lower MAE (15.12 cm) and RMSE (15.48 cm; Table 4). This comparison (NHSnow vs. GlobSnow) showed a relatively good agreement, although NHSnow over-or under-estimated the SD with a larger RMSE. Overall, the performance of GlobSnow was better than the NHSnow grid data. However, part of the validation data was also applied for the GlobSnow assimilation; it is highly possible that, in this case, the GlobSnow validation may not have been completely independent. The grain size was more important than the snow density and temperature in the microwave emission modeling. The different performances for these two datasets might have been mainly caused by the variation of the snow grain size, which was used to generate the SWE of GlobSnow [17]. Che et al. [33] reported that the grain size is more important than snow density and temperature. Furthermore, ERA-Interim/Land had the worst performance of all of the three snow cover datasets, with the highest bias (5.60 cm), MAE (18.72 cm), and RMSE (37.77 cm). The smallest bias was found for the mid-latitude regions (<50 • N), and much of the bias lay at 0-100 cm for ERA-Interim/Land (Figure 3c,f). We found large MAE and RMSE in the high latitude and coastal regions (Figure 3c). Unlike NHSnow and GlobSnow, ERA-Interim/Land was more likely to overestimate SD, and appeared to be less consistent with the in situ observations across the Northern Hemisphere (Figure 3f).

Variation of Snow Depth
To understand and interpret snow cover variations better in the Northern Hemisphere, we analyzed the SD variation using the seasonal maximum SD of the NHSnow dataset from 1992 to 2016 with variation trend grading (Table 3; Figure 6). The variation trend of the seasonal maximum SD exhibited distinctly different patterns in the three seasons during 1992-2016. Spatially, we found that the areas where SD exhibited a highly significant decrease trend in autumn (Figure 6a) were mainly located in the Russian Far East, the Qinghai-Tibet Plateau, the south of the central Siberian Plateau, and the northeastern regions of Canada. Over Russia's Taimyr Peninsula and the Alaska region, SD exhibited a highly significant increase trend (0-1 cm yr −1 ). Moreover, the seasonal maximum SD in winter and spring also underwent a highly significant decrease trend in the Qinghai-Tibet Plateau and the northeastern region of Canada (Figure 6b,c). The areas with a highly significant SD decrease extended to the Western Siberian plain region. Wang and Li [79] used nearly 50 years of daily station SD observation data to analyze the trend of maximum SD in China, and found that variation trend of the seasonal maximum SD in the Qinghai-Tibet Plateau was consistent with the observations from this study. The seasonal maximum SD variations in autumn and winter in the Russian Far East area exhibited a highly significant decreasing trend, while in spring, it showed a highly significant increasing trend. This variation trend pattern of the NHSnow seasonal maximum SD in spring was consistent with the analysis results obtained using the GlobSnow datasets [80], showing increasing trends in the Scandinavian peninsula and Alaska regions, and decreasing trends in the Kamchatka Peninsula region. We quantified the proportion of each of the variation trend grades in each season, because the snow cover area varied between seasons. From Figure 7, non-significant change dominated the seasonal maximum SD variation trends of the Northern Hemisphere (more than 50%) in three seasons. There were more regions (green regions) with a (highly) significant decrease trend of seasonal maximum SD in winter and spring, accounting for 31% and 34%, respectively, while the areas showing (highly) significant increasing trends (19% and 8%) were less. As for spring, the proportion of area with an increasing trend (about 25%) in the seasonal maximum SD was slightly larger than the decreasing trend (about 24%).
Finally, we analyzed the variation pattern of the average SD across the Northern Hemisphere using the average SD for each season. The variation rate of each season's average SD exhibited great fluctuations in different regions and seasons. The winter average SD decreased with the highest rate of −0.11 cm yr −1 compared with the other two seasons during 1992-2016 (Figure 8b; Table 5). In particular, the seasonal average SD variation rate in the Arctic region fluctuates greatly from region to region far more than average (Figures 8 and 9). Additionally, the absolute values of the seasonal average SD variation rate in the Arctic regions were apparently larger than in the middle-low latitude regions (Figures 8 and 9; Figures A1 and A2 in Appendix A). These results are indeed consistent with station observations over Northern Russia [81], Northern Canada [82], and Alaska. The average seasonal SD in the Northern Hemisphere is subject to a large standard deviation. A similar conclusion can be found in the two periods of 1992-2001 and 2002-2016 ( Figures A1 and A2 in Appendix A; Table 5). The majority regions in the Northern Hemisphere showed increasing trends in seasonal average SD during 1992-2001 at the rate of 0.02 cm yr −1 ( Figure A1 and Table 5). However, the seasonal average SD variation trend during 2002-2016 was reversed, presenting decreasing trends (autumn: 0.15 cm yr −1 ; winter: 0.22 cm yr −1 ; and spring 0.07 cm yr −1 ). The absolute variation rate during 2002-2016 was apparently greater than the rate during 1992-2001 ( Figure A2 and Table 5). We found that a high absolute value of the SD variation rate (>1 cm yr −1 ) usually occurred in the Arctic and the Qinghai-Tibetan plateau for all three seasons. The snowpack inevitably contained liquid water, as snow melting events occurred especially in spring. During the spring periods, consequently, the SD and snow cover days derived from the NHSnow datasets were partly underestimated and bias might have occurred in the trend change analysis of snow cover. The breakpoint time (2001/2002) was selected because it was the beginning of the new century.  Table 5). The majority regions in the Northern Hemisphere showed increasing trends in seasonal average SD during 1992-2001 at the rate of 0.02 cm yr −1 ( Figure A1 and  Figure A2 and Table 5). We found that a high absolute value of the SD variation rate (>1 cm yr −1 ) usually occurred in the Arctic and the Qinghai-Tibetan plateau for all three seasons. The snowpack inevitably contained liquid water, as snow melting events occurred especially in spring. During the spring periods, consequently, the SD and snow cover days derived from the NHSnow datasets were partly underestimated and bias might have occurred in the trend change analysis of snow cover. The breakpoint time (2001/2002) was selected because it was the beginning of the new century.

Snow Cover Days
Investigating the variation of SCD showed a rate of −0.99 ± 2.23 day yr −1 during 1992-2016. Most areas across the Northern Hemisphere presented a prominent decreasing trend, with the decreasing rate ranging from 0 to 5 day yr −1 (Figure 10a). The decreasing trend regions are mainly distributed in Eurasia, covering a vast area of the northern Russia and central Asia. The regions with a decreasing rate greater than 5 day yr −1 were almost all located in China, for example the Qilian Mountain, central Tibetan Plateau, and Tianshan regions. The regions with increasing trends were found in the central North America, Western Europe, Northwestern Mongolia, and some areas of China. Throughout the Northern Hemisphere (Figure 10a), the areas with decreasing trends covered the most regions, spanning from 25 • N through 85 • N, with a mean decreasing rate of approximately 1 day yr −1 . The most notable variation (either decrease or increase) occurred in the Arctic regions (Figure 10b).

Snow Cover Days
Investigating the variation of SCD showed a rate of −0.99 ± 2.23 day yr −1 during 1992-2016. Most areas across the Northern Hemisphere presented a prominent decreasing trend, with the decreasing rate ranging from 0 to 5 day yr −1 (Figure 10a). The decreasing trend regions are mainly distributed in Eurasia, covering a vast area of the northern Russia and central Asia. The regions with a decreasing rate greater than 5 day yr −1 were almost all located in China, for example the Qilian Mountain, central Tibetan Plateau, and Tianshan regions. The regions with increasing trends were found in the central North America, Western Europe, Northwestern Mongolia, and some areas of China. Throughout the Northern Hemisphere (Figure 10a), the areas with decreasing trends covered the most regions, spanning from 25° N through 85° N, with a mean decreasing rate of approximately 1 day yr −1 . The most notable variation (either decrease or increase) occurred in the Arctic regions (Figure 10b). Unlike the SCD variation rate pattern, the SCD variation trend showed that the majority of regions (approximately 64%) demonstrated a non-significant change over the Northern Hemisphere (Figure 10b; Figure A3 in the Appendix A). The highly significant decrease and the significant decrease trend, which appeared in the northwest of Hudson Bay in Canada, the Kamchatka peninsula, the Eastern European plains, Northern Russia, the Iranian plateau, and several regions in China (the Tibet Plateau, Mount Tianshan, and Northeast China Plain), accounted for 33% of the total snow cover area. In addition, areas with a highly significant decrease and significant increase trend (3%) were only found in a limited area of North America and the eastern Qinghai-Tibetan Plateau regions.
In terms of the changes in SCD and the annual average SD, they were either in the same or opposite directions over the past few decades. We found that the area where the annual average SD changed significantly (increase or decrease) was apparently larger than that of SCD (54% vs. 36%; Figure A3 in the Appendix A). Generally, more than 60% of a snow cover area has the same variation trends for the SCD and annual average SD (Figure 11), and the remaining area (40%) had different variation trends (non-significant change or a significant increase or decrease). Interestingly, the areas Unlike the SCD variation rate pattern, the SCD variation trend showed that the majority of regions (approximately 64%) demonstrated a non-significant change over the Northern Hemisphere (Figure 10b; Figure A3 in the Appendix A). The highly significant decrease and the significant decrease trend, which appeared in the northwest of Hudson Bay in Canada, the Kamchatka peninsula, the Eastern European plains, Northern Russia, the Iranian plateau, and several regions in China (the Tibet Plateau, Mount Tianshan, and Northeast China Plain), accounted for 33% of the total snow cover area. In addition, areas with a highly significant decrease and significant increase trend (3%) were only found in a limited area of North America and the eastern Qinghai-Tibetan Plateau regions.
In terms of the changes in SCD and the annual average SD, they were either in the same or opposite directions over the past few decades. We found that the area where the annual average SD changed significantly (increase or decrease) was apparently larger than that of SCD (54% vs. 36%; Figure A3 in the Appendix A). Generally, more than 60% of a snow cover area has the same variation trends for the SCD and annual average SD (Figure 11), and the remaining area (40%) had different variation trends (non-significant change or a significant increase or decrease). Interestingly, the areas with the opposite trend of variation in SCD and annual average SD were approximately 20%, i.e., areas with a significant increasing trend in the annual average SD while showing the SCD with a significant decreasing trend in the corresponding areas (pink area in Figure 11; 4%), and areas with a significant decreasing trend of annual average SD, while showing the SCD with a significant increasing trend (blue area in Figure 11; 16%). These regions with the opposite trend (increasing for annual average SD but decreasing for SCD) mostly appeared in the northern Central Siberian Plateau, the eastern Scandinavian Peninsula, Western Alaska, and the northern Greater Khingan Mountains of China. Zhong et al. [83] reported similar results across the northern Eurasian continent using ground-based measurements. As a reference, we also analyzed the variation trends of SCD and annual maximum SD ( Figure A4 in the Appendix A), showing that the areas with the opposite trend of variation in SCD and the annual maximum SD accounted for 21%. The pattern of these variation trend areas was also the same as the pattern for the SCD and annual average SD, except for a few areas (Figure 11 vs. Figure A4). The increase in SD may be explained by the increasing frequency of extreme snowfall [83]. Additionally, a recent study found that the greater the SWE, the faster the melting rate, leading to a shortened SCD in the Northern Hemisphere [80], along with an increasing winter air temperature [81,82]. Furthermore, the areas showing an annual average SD significantly decreased, but SCD significantly increased, accounting for 16% of the total snow cover area, mainly distributed in regions around 60 • N and the margin of the Qinghai-Tibetan Plateau.
with the opposite trend of variation in SCD and annual average SD were approximately 20%, i.e., areas with a significant increasing trend in the annual average SD while showing the SCD with a significant decreasing trend in the corresponding areas (pink area in Figure 11; 4%), and areas with a significant decreasing trend of annual average SD, while showing the SCD with a significant increasing trend (blue area in Figure 11; 16%). These regions with the opposite trend (increasing for annual average SD but decreasing for SCD) mostly appeared in the northern Central Siberian Plateau, the eastern Scandinavian Peninsula, Western Alaska, and the northern Greater Khingan Mountains of China. Zhong et al. [83] reported similar results across the northern Eurasian continent using ground-based measurements. As a reference, we also analyzed the variation trends of SCD and annual maximum SD ( Figure A4 in the Appendix A), showing that the areas with the opposite trend of variation in SCD and the annual maximum SD accounted for 21%. The pattern of these variation trend areas was also the same as the pattern for the SCD and annual average SD, except for a few areas (Figure 11 vs. Figure A4). The increase in SD may be explained by the increasing frequency of extreme snowfall [83]. Additionally, a recent study found that the greater the SWE, the faster the melting rate, leading to a shortened SCD in the Northern Hemisphere [80], along with an increasing winter air temperature [81,82]. Furthermore, the areas showing an annual average SD significantly decreased, but SCD significantly increased, accounting for 16% of the total snow cover area, mainly distributed in regions around 60° N and the margin of the Qinghai-Tibetan Plateau. Figure 11. Analysis results of the variation trends for SCD and the annual average SD over the Northern Hemisphere in 1992-2016. Grey indicates that the areas are undergoing a significant increase (or decrease) in SCD (or annual average SD) and non-significant change for annual average SD (or SCD); green indicate that the areas are undergoing a significant decrease in the annual average SD and SCD; blue indicates that the areas are undergoing a significant decrease in the annual average SD and a significant increase in SCD; cyan indicates that the areas are undergoing a non-significant change in the annual average SD and SCD; violet indicates that the areas are undergoing a significant increase in annual average SD and a significant decrease for SCD; orange indicates that the areas are undergoing significant increase in annual average SD and SCD. The significant increase (or decrease) here means that this variation trend passes the 95% confidence level test. The numbers in the brackets represent their proportion of the total snow cover area.

Discussion
While these three gridded data (NHSnow, GlobSnow, and ERA-Interim/Land) do a fairly good job in small SD accumulations (shallow and mid-deep snow cover) regions, they all struggle to capture SDs with a low bias, MAE, and RMSE in high SD accumulation (deep snow) regions ( Figures  3-5). There are insurmountable saturation effects in deep snowpacks when using passive microwave remote sensing data [24,84]. All of the SD retrieval algorithms based on the attenuation of radiation have insensitivity to deep snow. As a result, the snow cover variation could not be adequately captured in areas with deep snow; thus, we should be cautious when interpreting the validation result in deep snow regions.
It should be noted that snow-covered area of NHSnow was dry snow, as determined by Grody's series of rules [75]. SCD is determined on the identified snow-covered area [24,75]. The threshold of SD may have an effect on the specific value of the generated SCD. The sensitivity of SCD to the SD Figure 11. Analysis results of the variation trends for SCD and the annual average SD over the Northern Hemisphere in 1992-2016. Grey indicates that the areas are undergoing a significant increase (or decrease) in SCD (or annual average SD) and non-significant change for annual average SD (or SCD); green indicate that the areas are undergoing a significant decrease in the annual average SD and SCD; blue indicates that the areas are undergoing a significant decrease in the annual average SD and a significant increase in SCD; cyan indicates that the areas are undergoing a non-significant change in the annual average SD and SCD; violet indicates that the areas are undergoing a significant increase in annual average SD and a significant decrease for SCD; orange indicates that the areas are undergoing significant increase in annual average SD and SCD. The significant increase (or decrease) here means that this variation trend passes the 95% confidence level test. The numbers in the brackets represent their proportion of the total snow cover area.

Discussion
While these three gridded data (NHSnow, GlobSnow, and ERA-Interim/Land) do a fairly good job in small SD accumulations (shallow and mid-deep snow cover) regions, they all struggle to capture SDs with a low bias, MAE, and RMSE in high SD accumulation (deep snow) regions (Figures 3-5). There are insurmountable saturation effects in deep snowpacks when using passive microwave remote sensing data [24,84]. All of the SD retrieval algorithms based on the attenuation of radiation have insensitivity to deep snow. As a result, the snow cover variation could not be adequately captured in areas with deep snow; thus, we should be cautious when interpreting the validation result in deep snow regions.
It should be noted that snow-covered area of NHSnow was dry snow, as determined by Grody's series of rules [75]. SCD is determined on the identified snow-covered area [24,75]. The threshold of SD may have an effect on the specific value of the generated SCD. The sensitivity of SCD to the SD threshold changes were analyzed using two datasets (September 2014 to June 2015; September 2015 to June 2016; Figure 12). We found that the low threshold of SD (<3 cm) used to define the SCD did not significantly change the value of SCD, while a high value for SD could have significant decrease for SCD and snow cover extent, as reported by [49]. When compared to the SCD derived from optical sensor data, the specific quantity of SCD and the SCD variation rate derived from the NHSnow data was overestimated. As the optical (MODIS or AVHRR) and microwave sensors (SSM/I or AMSR-E) respond to different parts of the electromagnetic spectrum, the estimated snow cover area and SCD from these two sensors would vary somewhat [51,85]. Notwithstanding, another explanation for these discrepancies may be the snow cover identification algorithm [25,86,87]. The microwave radiation characteristics of snow cover are similar to that of precipitation, cold desert, and frozen ground [75]. Commission and omission errors in the NHSnow datasets may result from coarse spatial resolution, snow characteristics and topography [21], and precipitation [75,86], especially over frozen ground [88]. Several rules for NHSnow algorithm development were applied to distinguish snow from precipitation, cold desert, and frozen ground [24]; it was impossible to entirely remove interference factors in each image at the current studies. Moreover, the precondition of NHSnow is dry snow, which means almost no wet snow was considered in the SCD variation analysis [89]. The poorer performance of the microwave derived datasets was anticipated because of the documented difficulties in monitoring snow cover over forested and mountainous terrain [34,73].
Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 28 threshold changes were analyzed using two datasets (September 2014 to June 2015; September 2015 to June 2016; Figure 12). We found that the low threshold of SD (<3 cm) used to define the SCD did not significantly change the value of SCD, while a high value for SD could have significant decrease for SCD and snow cover extent, as reported by [49]. When compared to the SCD derived from optical sensor data, the specific quantity of SCD and the SCD variation rate derived from the NHSnow data was overestimated. As the optical (MODIS or AVHRR) and microwave sensors (SSM/I or AMSR-E) respond to different parts of the electromagnetic spectrum, the estimated snow cover area and SCD from these two sensors would vary somewhat [51,85]. Notwithstanding, another explanation for these discrepancies may be the snow cover identification algorithm [25,86,87]. The microwave radiation characteristics of snow cover are similar to that of precipitation, cold desert, and frozen ground [75]. Commission and omission errors in the NHSnow datasets may result from coarse spatial resolution, snow characteristics and topography [21], and precipitation [75,86], especially over frozen ground [88]. Several rules for NHSnow algorithm development were applied to distinguish snow from precipitation, cold desert, and frozen ground [24]; it was impossible to entirely remove interference factors in each image at the current studies. Moreover, the precondition of NHSnow is dry snow, which means almost no wet snow was considered in the SCD variation analysis [89]. The poorer performance of the microwave derived datasets was anticipated because of the documented difficulties in monitoring snow cover over forested and mountainous terrain [34,73].

Conclusions
This study applied the revised support vector regression snow-depth retrieval algorithm, which uses passive microwave remote sensing and other auxiliary data to generate long term (1 January 1992, to 31 December 2016) Northern Hemisphere daily snow depth (SD) datasets (NHSnow) with a 25-km spatial resolution. When compared to the other snow cover datasets (GlobSnow and ERA-Interim/Land) benchmarked against the ground SD measurements, the NHSnow SD datasets had a relatively smaller bias (−0.59 cm), mean absolute error (15.12 cm), and root mean square error (20.11 cm).
NHSnow datasets with full coverage and high accuracy can provide useful information about the spatiotemporal distribution of snow cover. We quantificationally investigated the spatial and temporal change characteristics in snow cover (SD and SCD) across the Northern Hemisphere. The analysis results show that the SD variation pattern in the variation rate or variation trend (increase or decrease) varies greatly among the different areas. For the seasonal maximum SD, areas with a

Conclusions
This study applied the revised support vector regression snow-depth retrieval algorithm, which uses passive microwave remote sensing and other auxiliary data to generate long term (1 January 1992, to 31 December 2016) Northern Hemisphere daily snow depth (SD) datasets (NHSnow) with a 25-km spatial resolution. When compared to the other snow cover datasets (GlobSnow and ERA-Interim/Land) benchmarked against the ground SD measurements, the NHSnow SD datasets had a relatively smaller bias (−0.59 cm), mean absolute error (15.12 cm), and root mean square error (20.11 cm).
NHSnow datasets with full coverage and high accuracy can provide useful information about the spatiotemporal distribution of snow cover. We quantificationally investigated the spatial and temporal change characteristics in snow cover (SD and SCD) across the Northern Hemisphere. The analysis results show that the SD variation pattern in the variation rate or variation trend (increase or decrease) varies greatly among the different areas. For the seasonal maximum SD, areas with a significant decreasing trend are larger than those with a significant increasing trend in three seasons, namely, autumn (34% vs. 8%), winter (31% vs. 19%), and spring (25% vs. 24%). The results show that the annual average SD decreased at a rate of −0.06 cm yr −1 from 1992 to 2016 across the Northern Hemisphere, and the variation rate in winter was the highest, reaching up to −0.11 cm yr −1 ; meanwhile, the absolute values of the variation rate of the seasonal average SD in the most Arctic regions are apparently greater than that in the middle-low latitude regions. In terms of the spatiotemporal variation trend of SCD over the Northern Hemisphere, we found that most areas show non-significant variation trends (more than 64%), while 33% of the area shows a significant decreasing trend. Interestingly, we found that only 60% of areas show the same variation trends for the SCD and annual average SD, but the area presenting a completely opposite variation trend for the SCD and annual average SD accounts for 20% of the total snow cover area, occurring in part of the Arctic regions and in the third polar region (Tibetan Plateau).
While this study shed light on the spatiotemporal variability trends of snow cover across the Northern Hemisphere using the 25-year NHSnow datasets, we cannot claim that the NHSnow dataset could completely capture the climate change signal in each region and season. Compared to the current relative best snow cover datasets, there are some deficiencies and limitations (e.g., overestimation and underestimation) for the NHSnow datasets, and further efforts should be made to improve the estimation accuracy and robustness of the SD inversion algorithm coupling the snow physical model/radiative transfer model, which appears to be a promising approach [27,90,91]. Meanwhile, a comprehensive validation analysis should also be carried out in complex terrain regions and different land cover types [21,22,92]. Further analyses and study are still needed to help understand the characteristics of SD changes in the Northern Hemisphere, e.g., analyzing and investigating the difference of snow cover variation and its response to climate change in regional and continental areas [17,[93][94][95]. It is recommended that future works focus on the climatic effects and climatological causes in snow cover changes in order to comprehensively understand the associated snow cover change mechanisms against a climate change background [5,7,50,66].  Remote Sens. 2020, 12, x FOR PEER REVIEW 22 of 28 Figure A3. The statistical analysis of the variation trend grades for snow cover days and annual average snow depth during 1992-2016. Figure A4. The analysis results of the variation trends for SCD and the annual maximum SD over the Northern Hemisphere during 1992-2016. Grey indicates that the areas are undergoing a significant increase (or decrease) for SCD (or annual maximum SD) and non-significant change for annual maximum SD (or SCD); green indicates that the areas are undergoing a significant decrease for annual maximum SD and SCD; blue indicate that the areas are undergoing a significant decrease for annual maximum SD and significant increase for SCD; cyan indicates that the areas are undergoing a nonsignificant change for annual maximum SD and SCD; violet indicates that the areas are undergoing a significant increase for annual average SD and a significant decrease for SCD; orange indicates that the areas are undergoing a significant increase for annual maximum SD and SCD. The significant increase (or decrease) here means that this variation trend passes the 95% confidence level test. The numbers in brackets represent the proportion of the total snow cover area. Figure A3. The statistical analysis of the variation trend grades for snow cover days and annual average snow depth during 1992-2016.
Remote Sens. 2020, 12, x FOR PEER REVIEW 22 of 28 Figure A3. The statistical analysis of the variation trend grades for snow cover days and annual average snow depth during 1992-2016. Figure A4. The analysis results of the variation trends for SCD and the annual maximum SD over the Northern Hemisphere during 1992-2016. Grey indicates that the areas are undergoing a significant increase (or decrease) for SCD (or annual maximum SD) and non-significant change for annual maximum SD (or SCD); green indicates that the areas are undergoing a significant decrease for annual maximum SD and SCD; blue indicate that the areas are undergoing a significant decrease for annual maximum SD and significant increase for SCD; cyan indicates that the areas are undergoing a nonsignificant change for annual maximum SD and SCD; violet indicates that the areas are undergoing a significant increase for annual average SD and a significant decrease for SCD; orange indicates that the areas are undergoing a significant increase for annual maximum SD and SCD. The significant increase (or decrease) here means that this variation trend passes the 95% confidence level test. The numbers in brackets represent the proportion of the total snow cover area. Figure A4. The analysis results of the variation trends for SCD and the annual maximum SD over the Northern Hemisphere during 1992-2016. Grey indicates that the areas are undergoing a significant increase (or decrease) for SCD (or annual maximum SD) and non-significant change for annual maximum SD (or SCD); green indicates that the areas are undergoing a significant decrease for annual maximum SD and SCD; blue indicate that the areas are undergoing a significant decrease for annual maximum SD and significant increase for SCD; cyan indicates that the areas are undergoing a non-significant change for annual maximum SD and SCD; violet indicates that the areas are undergoing a significant increase for annual average SD and a significant decrease for SCD; orange indicates that the areas are undergoing a significant increase for annual maximum SD and SCD. The significant increase (or decrease) here means that this variation trend passes the 95% confidence level test. The numbers in brackets represent the proportion of the total snow cover area. Wooded grassland Prairie (grassland) 10 Grassland 8 Closed shrub land Shrub 9 Open shrub land 11 Cropland Bare-land 12 Bare ground 13 Urban and built