Tracking Snow Variations in the Northern Hemisphere Using Multi-Source Remote Sensing Data ( 2000 – 2015 )

Multi-source remote sensing data were used to generate 500-m resolution cloud-free daily snow cover images for the Northern Hemisphere. Simultaneously, the spatial and temporal dynamic variations of snow in the Northern Hemisphere were evaluated from 2000 to 2015. The results indicated that (1) the maximum, minimum, and annual average snow-covered area (SCA) in the Northern Hemisphere exhibited a fluctuating downward trend; the variation of snow cover in the Northern Hemisphere had well-defined inter-annual and regional differences; (2) the average SCA in the Northern Hemisphere was the largest in January and the smallest in August; the SCA exhibited a downward trend for the monthly variations from February to April; and the seasonal variation in the SCA exhibited a downward trend in the spring, summer, and fall in the Northern Hemisphere (no pronounced variation trend in the winter was observed) during the 2000–2015 period; (3) the spatial distribution of the annual average snow-covered day (SCD) was related to the latitudinal zonality, and the areas exhibiting an upward trend were mainly at the mid to low latitudes with unstable SCA variations; and (4) the snow reduction was significant in the perennial SCA in the Northern Hemisphere, including high-latitude and high-elevation mountainous regions (between 35◦ and 50◦N), such as the Tibetan Plateau, the Tianshan Mountains, the Pamir Plateau in Asia, the Alps in Europe, the Caucasus Mountains, and the Cordillera Mountains in North America.


Introduction
Snow cover is an indicator for global climate change.Previous studies have indicated that approximately 98% of the seasonal snow cover on the earth is located in the Northern Hemisphere [1].It has relatively large intra-annual and inter-annual variations.With a maximum extent of about approximately 4.7 × 10 7 km 2 , snow cover accounts for approximately 50% of the continental area in the Northern Hemisphere [2,3].
The snow cover extent has decreased considerably and retreated to higher altitudes, and the inter-annual snow cover variability has increased [4][5][6][7].The snow-covered area (SCA) in the Northern Hemisphere from 1970 to 2010 in March and April exhibited a decreasing trend, and, on average, it decreased by approximately 8 × 10 5 km 2 each decade due to increasing temperatures [8], and the decay has accelerated during the past 40 years [9].Snow cover duration exhibits the strongest climate sensitivity [10], and temperature increases can postpone the snow cover accumulation date and the maximum SCA change from February to January, which causes the snow ablation date to arrive earlier in the spring [11].From 1967 to 2008, the stable snow cover duration during the winter in the Northern Hemisphere gradually shortened, decreasing by approximately 5.3 days every 10 years [12].The high albedo of snow can significantly reduce the capability of the ground surface to receive solar radiation and has a positive feedback effect on the climate system; in particular, positive albedo feedback in the spring period generate a stronger signal in spring snow cover duration [13].The winter SCA, duration, snow depth, and resulting variation in the surface albedo have significant influences on the atmospheric circulation at local and regional scales [14][15][16].Between the 1980s and the 21st century, the SCA decreased at a rate of (−0.55 ± 0.21) × 10 6 km 2 per decade in spring in the Northern Hemisphere [17,18].Furthermore, the Intergovernmental Panel on Climate Change (IPCC 2013) has evaluated the variation in snow cover in the Northern Hemisphere over the past 50 years, and the SCA during the spring in the Northern Hemisphere has continued to decrease (with high confidence).In the Northern Hemisphere, the average reduction rate of the SCA in March and April is 1.6% 10 −a and is 11.7% 10 −a in June, while the SCA does not exhibit a statistically significant increase in any month [19].
Remote sensing techniques have the advantages of high temporal and spatial resolution, strong objective authenticity, low cost, and high effectiveness [20], which provide effective snow cover parameters for climate models at global and regional scales.These techniques expand the traditional, limited "point" information to regional information that indicates the actual situation and plays an important role in large-scale high-precision monitoring of the earth's snow cover.Although the GlobSnow data provide information concerning the areal extent of snow (SE) on a global scale and the snow water equivalent (SWE) for the Northern Hemisphere, the SE product was only produced between the years 1995 to 2012 with huge gaps and the SWE product did not provide the areas of mountainous regions and Greenland [21,22].Currently, Moderate Resolution Imaging Spectroradiometer (MODIS) data are most widely used.MODIS is aboard the earth observation satellites Terra and Aqua, which are the main sensors of the Earth Observation System (EOS).The spatial and spectral resolution are better than the Advanced Very High Resolution Radiometer (AVHRR) and the temporal resolution is better than the Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper+ (ETM+).Thus, MODIS is the optimal data source for real-time monitoring of large-scale SCA with high accuracy [23][24][25][26][27].The Aqua satellite also carries the improved multi-band dual-polarized passive microwave radiometer AMSR-E (Advanced Microwave Scanning Radiometer-EOS).This sensor can provide more microwave band information than previous passive microwave radiometers, such as the Scanning Multichannel Microwave Radiometer (SMMR) and Special Sensor Microwave/Imager (SSM/I), and it is less affected by adverse weather conditions.Because the spatial resolution of AMSR-E is relatively coarse, the AMSR-E data are mainly used for studies on large-scale snow cover range, snow depth and snow water equivalent (SWE) at both global and hemisphere scales [28].In addition, the implementation of the interactive multi-sensor snow and ice mapping system (IMS) provides another approach for dynamic monitoring of snow cover, and scholars in China and other countries are currently using this product to monitor snow cover [29][30][31][32].It is the current representative product for the fusion of multi-source remote sensing satellites and has broad prospects in snow monitoring by remote sensing [31,33].However, the overall accuracy of the IMS is far lower than that of the MODIS [34], and the timing of IMS snow disappearance lags behind the MODIS during the spring ablation season by several weeks, which may affect our understanding of the seasonal snow cover cycle and development of parameterization schemes for climate model [35].
Brown et al. [36] used 10 data sources to carry out multi-source snow cover work firstly, and the results showed a 46% reduction in Arctic snow cover extent in June and a 14% reduction in May over the period between 1967 and 2008.The fusion of MODIS and AMSR-E data to improve the capability for remote sensing snow cover identification is a relatively mature algorithm for monitoring snow cover [37,38].Huang et al. [39] combined a digital elevation model to apply the snow line (SNOWL) algorithm [40] to different elevation regions of the Tibetan Plateau and developed a set of cloud removal algorithms based on optical and passive microwave remote sensing data to obtain a cloud-free daily snow cover product.Yu et al. [41] combined MODIS with the IMS snow cover product, and the overall classification accuracy for the generated cloud-free daily snow cover product reached 94%.
To understand the variation characteristics of snow cover in the Northern Hemisphere in the future, the temporal and spatial dynamics of snow cover in the Northern Hemisphere was analyzed during the period from 2000 to 2015, including the number of snow-covered days (SCDs), SCA, and the snow cover inter-annual variation.In addition, this work aims to establish the basis for further revealing the mechanism of interaction between temporal and spatial variations in both snow cover and climate change in the Northern Hemisphere.The objectives of this study are (1) to generate 500-m resolution daily cloud-free snow cover products with data from MODIS, AMSR-E SWE and IMS snow product for the Northern Hemisphere from 2000 to 2015; (2) to validate the accuracy of the daily cloud-free snow cover product; and (3) to analyze and discuss the characteristic trend of snow cover temporal and spatial dynamics in the Northern Hemisphere in detail.

Data
The remote sensing products used in this study mainly included (a) the daily MODIS snow product (MOD10A1 and MYD10A1, V005 version; the V006 version does not include binary snow cover products); (b) the land cover type product MCD12Q1; (c) the AMSR-E/Aqua SWE product; (d) the Landsat-TM product; and (e) the IMS product.Detailed information of the above remote sensing data is listed in Table 1.
The land cover type product MCD12Q1 derived from the classification system of the International Geosphere-Biosphere Program (IGBP) was used to study the difference in agreement between the synthesized snow cover product and Landsat-TM within different land cover types.There are eleven types of natural vegetation, three types of developed and mosaic lands, and three types of non-vegetation land.In this study, we reclassified the IGBP into seven major classes: water body, forest, shrubland, grassland, cropland, barren, and ice and snow, which statistically account for 2.2%, 17.2%, 16.9%, 24.9%, 17.8%, 13.8% and 7.2%, respectively.Because the MODIS snow cover classification product is affected by the land cover, we used the reclassified land cover product to analyze the influence of the different land cover types on the classification accuracy of the synthesized snow cover product to establish a foundation for further improving the classification accuracy of MODIS snow cover.

Cloud Removal Algorithm
The cloud removal algorithm used in this study has the following three components: (a) the composition of the daily MODIS snow cover product [46,47]-we conducted the maximum SCA fusion processing of the data for two products, MOD10A1 and MYD10A1, from the MODIS morning and afternoon tracks; (b) adjacent temporal composite-we used the snow cover image from the previous and posterior day to eliminate the cloud pixels in the snow cover image of the current day, and we ensure a minimum amount of clouds after these processes; (c) the combination with the SWE [37,38,48] or IMS product-the AMSR-E daily SWE product (filling the gap-using the AMSR-E images of the closest date to replace missing data with maximum values)/IMS product are used to completely eliminate the remaining cloud pixels, thus the accuracy of cloud pixels reclassification were fully dependent on the accuracy of AMSR-E and IMS.Prior to the procedure, some preprocessing such as reprojection and resampling (nearest neighbor interpolation) is needed.If the value of SWE exceeds 0 mm or the value of IMS is snow, then the cloud pixel is classified as snow.Using these procedures, we eliminated all cloud pixels and generated MODIS daily cloud-free binary snow cover images (MOYDCA) with a spatial resolution of 500 m, which were used to verify the accuracy for the Northern Hemisphere.The specific flow chart is shown in Figure 1.
Remote Sens. 2018, 10, 136 4 of 22 previous and posterior day to eliminate the cloud pixels in the snow cover image of the current day, and we ensure a minimum amount of clouds after these processes; (c) the combination with the SWE [37,38,48] or IMS product-the AMSR-E daily SWE product (filling the gap-using the AMSR-E images of the closest date to replace missing data with maximum values)/IMS product are used to completely eliminate the remaining cloud pixels, thus the accuracy of cloud pixels reclassification were fully dependent on the accuracy of AMSR-E and IMS.Prior to the procedure, some preprocessing such as reprojection and resampling (nearest neighbor interpolation) is needed.If the value of SWE exceeds 0 mm or the value of IMS is snow, then the cloud pixel is classified as snow.
Using these procedures, we eliminated all cloud pixels and generated MODIS daily cloud-free binary snow cover images (MOYDCA) with a spatial resolution of 500 m, which were used to verify the accuracy for the Northern Hemisphere.The specific flow chart is shown in Figure 1.

Landsat-TM Snow Mapping
Clouds and snow both have relatively high reflectivity in the visible band but different reflectivity in the shortwave infrared band.In the shortwave infrared band, clouds have very high reflectivity, and snow has a relatively strong absorptive feature.Based on this property of snow, snow cover extractions adopt the SNOWMAP algorithm [49].We first conducted the radiation calibration and the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction in ENVI 5.3 (Exelis Visual Information Solutions, Boulder, Colorado, USA).Then, we used the second and fifth wave bands of the TM images to calculate the NDSI.The calculation method is as follows: NDSI = (Band 2 − Band 5)/(Band 2 + Band 5), where Band 2 is the 0.55-μm band of the Landsat-TM image and Band 5 is the 1.64-μm band of the Landsat-TM image.We set the threshold to 0.40 [17].The reflectivity of water bodies in the visible band and shortwave infrared band is very low compared with that of snow.Thus, the water body category can be eliminated because the reflectivity of snow at the 4th band exceeds 11%.In a Landsat-TM image, if NDSI ≥0.40 and the reflectivity of the 4th waveband exceeds 11%, the pixel is

Landsat-TM Snow Mapping
Clouds and snow both have relatively high reflectivity in the visible band but different reflectivity in the shortwave infrared band.In the shortwave infrared band, clouds have very high reflectivity, and snow has a relatively strong absorptive feature.Based on this property of snow, snow cover extractions adopt the SNOWMAP algorithm [49].We first conducted the radiation calibration and the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction in ENVI 5.3 (Exelis Visual Information Solutions, Boulder, Colorado, USA).Then, we used the second and fifth wave bands of the TM images to calculate the NDSI.The calculation method is as follows: NDSI = (Band 2 − Band 5)/(Band 2 + Band 5), where Band 2 is the 0.55-µm band of the Landsat-TM image and Band 5 is the 1.64-µm band of the Landsat-TM image.We set the threshold to 0.40 [17].The reflectivity of water bodies in the visible band and shortwave infrared band is very low compared with that of snow.Thus, the water body category can be eliminated because the reflectivity of snow at the 4th band exceeds 11%.In a Landsat-TM image, if NDSI ≥0.40 and the reflectivity of the 4th waveband exceeds 11%, the pixel is identified as snow.Then, a binary snow cover image with a resolution of 30 m is generated, where 1 represents a snow-covered pixel and 0 represents a snow-free pixel.We used the Pixel Aggregate tool in ENVI to resample aggregately the Landsat binary snow cover image using a 500-m grid size.The tool averages all pixel values that contribute to the output pixel.If the output value exceeds 50%, the pixel is identified as snow.

Accuracy Assessment of the Snow-Cover Images
Currently, verification data of satellite product mainly consists of two types, the first is point data from meteorological observatories, and the second is satellite remote sensing data with higher resolution that can reflect the true situation on the ground.
Landsat-TM images have a higher spatial resolution than MODIS data.Moreover, the Landsat-TM and MODIS products used in the validation were for the same day.Therefore, in this study, we used the snow cover classification images from the TM as the "true value" and evaluated the consistency between TM and the synthesized snow cover images from MODIS using the Kappa coefficient.The Kappa coefficient is an index that is often used to measure the agreement or accuracy between two images, and the calculation method is as follows: where P O is the actual consistency rate, P C is the theoretical consistency rate, P O = s/n, and The total number of pixels in the image is n, the number of pixels for the snow cover in a TM image of the actual ground situation is a 1 , the number of other pixels is a 0 , the number of pixels for the snow cover in the self-defined synthesized image is b 1 , the number of other pixels is b 0 , and the number of pixels for which the corresponding pixel values of the lattices in the two images is s.Klein and Barnett [24] indicated that the Kappa test can be used on five sets to express the consistency of different levels.For 0 < Kappa ≤ 0.20, the two images have extremely low consistency; for 0.20 < Kappa ≤ 0.40, the two images have normal consistency; for 0.40 < Kappa ≤ 0.60, the two images have moderate consistency; for 0.60 < Kappa ≤ 0.80, the two images have a high consistency; for 0.80 < Kappa ≤ 1, the two images are almost completely consistent.

Analysis of the Snow Cover Variation
Parametric and non-parametric tests can be used to analyze significance in trends within hydro-meteorological time series data.Parametric trend tests require data to be independent and normally distributed and are more sensitive to outliers, while non-parametric trend tests only require the data to be independent.
The Mann-Kendall method is a non-parametric test that has been frequently used to analyze long time series data [50][51][52][53][54][55][56][57], and the method can also detect the variation trend of monotonic nonlinear data without a data distribution requirement [58].In this paper, the Mann-Kendall method was used to examine the trend and significance level of the annual SCDs in the Northern Hemisphere at the pixel scale.For a series X i = (X 1 , X 2 , • • • , X n ) with n samples, the specific test process was as follows: where where n is the number of years (n = 16), m is the number of nodes (repetitive data groups) in the series, and t i is the node width (the number of repetitive data points in the ith repetitive data group).
When n ≤ 10, we directly used the statistic S for the two-sided trend test.S > 0 indicates an upward trend; S = 0 indicates no variation; and S < 0 indicates a downward trend.For a given significance level, α, if |S| ≥ S a/2 , the trend of the series is significant; otherwise, it is insignificant.
When n > 10, the statistic S approaches the standardized normal distribution, and we used the test statistic Z for the two-sided trend test.Positive values of Z indicate an upward trend while negative Z values indicate a downward trend.For the given significance level, α, we looked up the critical value Z a/2 in the normal distribution table; if |Z| > Z a/2 , the trend of the series is significant; if |Z| ≤ Z a/2 , the trend is insignificant.
In addition, we also used Sen's median method [52,[59][60][61] to analyze the slope for the variation of the annual SCD.This method calculates the median slope for n(n − 1)/2 pairs of combinations in a series of length n, and the specific calculation formula is: where β > 0 indicates the trend is increasing, and β < 0 indicates the trend is decreasing.

Validation of the Daily Cloud-Free Snow Product
The Kappa consistency test method was adopted to analyze the accuracy of the MODIS cloud-free daily snow cover product obtained in this study.The largest Kappa coefficient between the synthesized product and the TM binary snow cover product was for the grassland and barren coverage types, and the coefficient was 0.675 in both cases; the Kappa value for the cropland cover type was 0.604.The Kappa analysis showed that their agreement was good and exhibited a high consistency.This was followed by the shrubland coverage type, with a Kappa coefficient of 0.599.The agreement was relatively good with moderate consistency.The accuracy was relatively poor for the forest coverage area, with a Kappa coefficient of 0.406, and the consistency performance was normal.Figure 2 compared the SCA based on cloud-free daily snow cover product and the Landsat-TM binary snow cover product, and showed the Kappa value for different land cover types.The results indicate that, for the grassland and barren land cover conditions, the acquired SCA was essentially consistent.The percentage of the SCA differed by only 0.78% and 0.47%, respectively.For the shrubland land cover condition, the SCA differed by 2.42%.The difference for the other land cover types was relatively large, and the difference for the forest land coverage was the largest (15.88%).
Furthermore, we compared and analyzed the difference between the MODIS cloud-free daily snow cover product and the Landsat-TM snow cover product for different land cover types and the overall average Kappa coefficients in Asia, Europe, and North America were 0.575, 0.592, and 0.584, respectively, representing moderate consistency, which were almost close to a high consistency.Therefore, for such large-scale areas in the Northern Hemisphere, this product has a relatively good reliability for studying the spatial and temporal dynamic variation characteristics of snow cover.
Therefore, for such large-scale areas in the Northern Hemisphere, this product has a relatively good reliability for studying the spatial and temporal dynamic variation characteristics of snow cover.The binary snow cover images acquired by Landsat 5 were used to validate the cloud-free daily snow cover product synthesized in this study.An example on 5 May 2007 is shown in Figure 3.The MOD10A1 and MYD10A1 images in the validation region are greatly affected by clouds, and the number of cloud pixels in the two images accounted for 33.21% and 52.68% of the total number of pixels, respectively.For the latter image, more than half of the validation region was contaminated by clouds.The SCA in the MOD10A1 image, MYD10A1 image, and MOYDCA-synthesized image accounted for 2.94%, 1.19%, and 5.56%, respectively, of the total, and the SCA in the Landsat-TM binary snow cover image accounted for 6.91% of the total.The Kappa coefficients for the three images were 0.494, 0.238, and 0.726, respectively, indicating that the standard snow cover products of MOD10A1 and MYD10A1 were seriously affected by clouds and that the degree of agreement with the Landsat-TM classification images was normal.There was moderate consistency between them, although there was a high consistency between the cloud-free snow cover classification images obtained using the cloud removal algorithm and the Landsat-TM classification images.This finding indicates that the cloud removal algorithm proposed in this study has good reliability for the dynamic monitoring of large-scale snow cover.The binary snow cover images acquired by Landsat 5 were used to validate the cloud-free daily snow cover product synthesized in this study.An example on 5 May 2007 is shown in Figure 3.The MOD10A1 and MYD10A1 images in the validation region are greatly affected by clouds, and the number of cloud pixels in the two images accounted for 33.21% and 52.68% of the total number of pixels, respectively.For the latter image, more than half of the validation region was contaminated by clouds.The SCA in the MOD10A1 image, MYD10A1 image, and MOYDCA-synthesized image accounted for 2.94%, 1.19%, and 5.56%, respectively, of the total, and the SCA in the Landsat-TM binary snow cover image accounted for 6.91% of the total.The Kappa coefficients for the three images were 0.494, 0.238, and 0.726, respectively, indicating that the standard snow cover products of MOD10A1 and MYD10A1 were seriously affected by clouds and that the degree of agreement with the Landsat-TM classification images was normal.There was moderate consistency between them, although there was a high consistency between the cloud-free snow cover classification images obtained using the cloud removal algorithm and the Landsat-TM classification images.This finding indicates that the cloud removal algorithm proposed in this study has good reliability for the dynamic monitoring of large-scale snow cover.

Spatio-Temporal Variation of SCA in the Northern Hemisphere
Figure 4 shows the variation curve of the maximum, average and minimum SCA in the Northern Hemisphere from 2000 to 2015.The maximum SCA in the Northern Hemisphere rapidly decreased in 2007, and it declined with fluctuations after reaching a peak in 2008.The variation in the minimum SCA of the Northern Hemisphere is of particular interest because the area contains perennial snow cover.Dynamic monitoring of the perennial snow area in the Northern Hemisphere plays a vital role in monitoring global climate and hydrological changes.Figure 4 shows that both the maximum and minimum SCA in the Northern Hemisphere exhibited a downward trend during the study period, and the tendency rates were −0.057% a −1 and −0.045% a −1 , respectively.Since 2000, the inter-annual variation in the SCA of the Northern Hemisphere has exhibited a downward trend with fluctuations, with an average rate of −0.066% a −1 .(5.80%).The 4th report of the IPCC noted that the largest average SCA in the Northern Hemisphere typically occurs in January, while the smallest occurs in August [62], and this study is consistent with these conclusions.In addition, based on the figure, the SCA fluctuation in the Northern Hemisphere was relatively large from February to June, and the fluctuation exhibited a negatively skewed distribution from February to April, which was mainly concentrated in the spring.This also indicated that the spring SCA had a downward trend during the study period.The 5th assessment report (AR5) of the IPCC noted that, since the middle 20th century, the SCA in the Northern Hemisphere has significantly decreased, indicating that the main cause is the reduction in the SCA in the spring [19].This figure shows that the largest and smallest SCA in the Northern Hemisphere occurred in January and August, respectively, i.e., approximately 5.5 × 10 7 km 2 (55.38%) and 5.76 × 10 6 km 2 (5.80%).The 4th report of the IPCC noted that the largest average SCA in the Northern Hemisphere typically occurs in January, while the smallest occurs in August [62], and this study is consistent with these conclusions.In addition, based on the figure, the SCA fluctuation in the Northern Hemisphere was relatively large from February to June, and the fluctuation exhibited a negatively skewed distribution from February to April, which was mainly concentrated in the spring.This also indicated that the spring SCA had a downward trend during the study period.The 5th assessment report (AR5) of the IPCC noted that, since the middle 20th century, the SCA in the Northern Hemisphere has significantly decreased, indicating that the main cause is the reduction in the SCA in the spring [19].The snow cover in the Northern Hemisphere exhibits seasonal variability.We conducted statistical analyses on the SCA for all four seasons (winter, spring, summer, and fall) in the Northern Hemisphere.The largest SCA was in the winter, which is typically the main time period for snow accumulation.Winter was followed by spring and fall, with the smallest SCA occurring in the summer.The trend in variation for the SCA during the four time periods is shown in Figure 6.Since the early 21st century, the SCA in the Northern Hemisphere during the spring, summer, and fall has decreased.The reduction in the magnitude of the SCA has been relatively large in the spring and summer, with tendencies of −0.108% a −1 and −0.110% a −1 , respectively; the SCA has slightly decreased in the fall, with a tendency of −0.035% a −1 .Considering that the global climate has become warmer, the SCA during the winter in the Northern Hemisphere did not show any obvious variation trend.The snow cover in the Northern Hemisphere exhibits seasonal variability.We conducted statistical analyses on the SCA for all four seasons (winter, spring, summer, and fall) in the Northern Hemisphere.The largest SCA was in the winter, which is typically the main time period for snow accumulation.Winter was followed by spring and fall, with the smallest SCA occurring in the summer.The trend in variation for the SCA during the four time periods is shown in Figure 6.Since the early 21st century, the SCA in the Northern Hemisphere during the spring, summer, and fall has decreased.The reduction in the magnitude of the SCA has been relatively large in the spring and summer, with tendencies of −0.108% a −1 and −0.110% a −1 , respectively; the SCA has slightly decreased in the fall, with a tendency of −0.035% a −1 .Considering that the global climate has become warmer, the SCA during the winter in the Northern Hemisphere did not show any obvious variation trend.

Spatio-Temporal Variation of SCDs in the Northern Hemisphere
SCDs represent the accumulated number of days of a particular pixel covered by snow during one year.The SCDs data from 2000 to 2015 in the Northern Hemisphere in Figure 7 showed that the area with SCDs ≤ 10 accounted for approximately 3.28 × 10 7 km 2 (33%) of the total area.The SCA of these regions exhibited an upward trend over the 16-year period, increasing from 3.18 × 10 7 km 2

Spatio-Temporal Variation of SCDs in the Northern Hemisphere
SCDs represent the accumulated number of days of a particular pixel covered by snow during one year.The SCDs data from 2000 to 2015 in the Northern Hemisphere in Figure 7 showed that the area with SCDs ≤ 10 accounted for approximately 3.28 × 10 7 km 2 (33%) of the total area.The SCA of these regions exhibited an upward trend over the 16-year period, increasing from 3.18 × 10 7 km 2 (32.03%) in 2000 to 3.46 × 10 7 km 2 (34.86%) in 2015.The area with 120-180 SCDs also exhibited an upward trend.The region with 10 < SCDs ≤ 60 and SCDs ≥ 180 exhibited a downward trend.The regions with SCDs ≥ 350 were covered by perennial snow or glaciers throughout the year, and this part of the snow cover is an important indicator of global climatic change.There was not an obvious change in regions with other SCDs.In addition, Figure 8 shows the spatial distribution of the average annual SCDs in the Northern Hemisphere from 2000 to 2015, and Table 2 shows the percentage by region of different SCDs in Asia, Europe and North America, respectively.The results indicate that the spatial distribution of the annual average SCDs was based on latitudinal zonality.
In Asia, the instantaneous SCA with SCDs ≤ 10 was approximately 9.42 × 10 6 km 2 (26.1%), which primarily occurred near the equator, including Southeast Asia, Southern China, the Tarim Basin in Xinjiang of China, most of South Asia, Eastern Iran, and the Arabian Peninsula in the Western Asia.The unstable SCA (10 < SCDs ≤ 60) accounted for 4.98 × 10 6 km 2 (13.8%) of the total SCA, which mainly occurred in the Korean Peninsula, the Northern and most Western regions of China, the Southern islands of Japan, and Uzbekistan and its surrounding area in Central Asia.The stable SCA (60 < SCDs ≤ 350) accounted for the largest proportion, 2.17 × 10 7 km 2 (60.1%), and it was mainly distributed in the Northeastern region of Inner Mongolia of China, North of Xinjiang, the Tibetan In addition, Figure 8 shows the spatial distribution of the average annual SCDs in the Northern Hemisphere from 2000 to 2015, and Table 2 shows the percentage by region of different SCDs in Asia, Europe and North America, respectively.The results indicate that the spatial distribution of the annual average SCDs was based on latitudinal zonality.
In Asia, the instantaneous SCA with SCDs ≤ 10 was approximately 9.42 × 10 6 km 2 (26.1%), which primarily occurred near the equator, including Southeast Asia, Southern China, the Tarim Basin in Xinjiang of China, most of South Asia, Eastern Iran, and the Arabian Peninsula in the Western Asia.The unstable SCA (10 < SCDs ≤ 60) accounted for 4.98 × 10 6 km 2 (13.8%) of the total SCA, which mainly occurred in the Korean Peninsula, the Northern and most Western regions of China, the Southern islands of Japan, and Uzbekistan and its surrounding area in Central Asia.The stable SCA (60 < SCDs ≤ 350) accounted for the largest proportion, 2.17  Between 2000 and 2015, the instantaneous SCA in the European region accounted for 7.98 × 10 5 km 2 (7.9%) of the total SCA in Europe.The areas were mainly distributed in the Northern islands of the Mediterranean, the Iberian Peninsula, and the Western region of France.The unstable SCA accounted for 2.62 × 10 6 km 2 (25.9%) of the total SCA in Europe and was distributed in Britain, the Western European Plains, the Central European Plains, the Apennines and Balkan Peninsula, and the Northern region of the Caucasus Mountains.Similarly, the area of stable SCA was the largest in Europe, 6.69 × 10 6 km 2 (66.2%), and it was broadly distributed in Central and Northern Europe and the Alps and Caucasus Mountains.The perennial SCA with annual SCDs > 350 were only distributed in the Novaya Zemlya of Russia and Svalbard of Norway.
The regions with SCDs ≤ 10 in North America were mainly the Southern region of the United States and countries at lower latitudes, which accounted for 3.36 × 10 6 km 2 (13.6%).The unstable SCA Between 2000 and 2015, the instantaneous SCA in the European region accounted for 7.98 × 10 5 km 2 (7.9%) of the total SCA in Europe.The areas were mainly distributed in the Northern islands of the Mediterranean, the Iberian Peninsula, and the Western region of France.The unstable SCA accounted for 2.62 × 10 6 km 2 (25.9%) of the total SCA in Europe and was distributed in Britain, the Western European Plains, the Central European Plains, the Apennines and Balkan Peninsula, and the Northern region of the Caucasus Mountains.Similarly, the area of stable SCA was the largest in Europe, 6.69 × 10 6 km 2 (66.2%), and it was broadly distributed in Central and Northern Europe and the Alps and Caucasus Mountains.The perennial SCA with annual SCDs > 350 were only distributed in the Novaya Zemlya of Russia and Svalbard of Norway.
The regions with SCDs ≤ 10 in North America were mainly the Southern region of the United States and countries at lower latitudes, which accounted for 3.36 × 10 6 km 2 (13.6%).The unstable SCA only accounted for 2.12 × 10 6 km 2 (8.6%), and it was distributed in the mid-latitude area of the United States.The stable SCA accounted for 1.61 × 10 7 km 2 (65.1%), and it was distributed in Canada, the high-latitude region in the mainland of the United States, and Alaska.The perennial SCA with SCDs > 350 was mostly in Greenland and Canadian Arctic Archipelago, with an area covering 3.14 × 10 6 km 2 (12.7%) of the total, and this was also the main distribution region of perennial snow cover in the Northern Hemisphere.Because Greenland remains snow and ice covered during all season, we do not analyze Greenland.

Snow Cover Change in the Northern Hemisphere
We used the Mann-Kendall method to further analyze the SCD variation tendency and its significance for the different continents of the Northern Hemisphere from 2000 to 2015 at the pixel scale in Figure 9. Table 3 lists the detailed information.To validate the accuracy of the results, we also used Sen's median method to calculate and derive the SCD variation slope for different continents, and the results were shown in Figure 10.Of note, the areas where SCDs reduction was significant in the Northern Hemisphere were distributed primarily in mountainous areas at 35 • -50 • N.
From 2000 to 2015, the SCDs over 1.79 × 10 7 km 2 (49.6%) of the area in Asian exhibited an increasing trend (Z > 0), but the increase was very gentle.The areas with a significant increase only accounted for 9.03 × 10 5 km 2 (2.5%) (p < 0.05), and they are mainly distributed in the Southeastern region of Siberia and Kamchatka Peninsula in Russia, Japan, and the Northeastern region of China, the Sichuan Basin, and Southeastern region of China.The area in which SCDs decreased accounted for 1.02 × 10 7 km 2 (28.3%) of Asia (Z < 0), and the regions with a significant SCD decrease only accounted for 6.35 × 10 6 km 2 (17.6%) (p < 0.05), which was mainly distributed in the Western mountainous region of Eastern Siberia in Russia, the Central Siberian Plateau, the Central and Western regions of the Mongolian Plateau, the Altai, Tianshan, Kunlun, and the Qilian Mountains of China, the Himalayas, the Southeastern region of the Tibetan Plateau, the Pamir Plateau, the Asia Minor Peninsula of West Asia, the Great Caucasus, and the Zagros Mountains.
The slope of the SCD variation in the Northern Hemisphere calculated using Sen's median method is shown in Figure 10.In Asia, SCDs over 23.8% of the area exhibited an increasing trend (β > 0), and SCDs over 27.7% of the area exhibited a decreasing trend (β < 0).In particular, the increased SCD rate for 22.6% of the region was less than 5 day/year, and these regions were mainly distributed across the Western Siberian Plain, the Southeastern region of Russia, most of Central Asia, the Northeastern and Central region of China, and some sporadic areas East of the Tibetan Plateau and in Southwest China.The regions with an increase rate greater than 5 day/year were sporadically distributed in the West of Japan and the Tibetan Plateau region.The SCD reduction rate in 27.9% of the region was less than 5 day/year and was mainly distributed in Central Siberia in Russia, the Asia Minor Peninsula, the Mongolian Plateau, and the Western region of China.The areas where the reduction rate exceeded 5 day/year were mainly distributed in the Great Caucasus Mountains, the Pamir Plateau, the Tianshan, Kunlun, and Qilian Mountains in China, the Himalayas, and the Southeastern region of the Tibetan Plateau.There was a relatively high consistency between the results from Sen's median method and the results from the Mann-Kendall method, especially for the spatial distribution of the tendency variation.
The results from combining Figures 9 and 10 indicated that the annual SCD in most regions of Asia exhibited an increasing trend.From a spatial distribution perspective, the area exhibiting an upward trend was mainly the unstable SCA with SCDs < 60, while the SCA with SCDs > 180 exhibited a downward trend.This finding indicated that the annual SCD in the seasonal snow cover distribution area in Asia exhibited an upward trend; the annual SCD in the high-elevation mountainous regions with the presence of perennial snow exhibited a downward trend.
From 2000 to 2015, the areas in Europe where SCDs increased (Z > 0) and decreased (Z < 0) were essentially the same, i.e., 4.48 × 10 6 km 2 (44.4%) and 4.32 × 10 6 km 2 (42.8%), respectively.The area with significant SCDs increase, accounting for 3.23 × 10 5 km 2 (3.2%), was distributed in the Western European Plains, Britain, and Ireland, and the annual SCDs in these regions was <60 days.The area with a significant SCD decrease, accounting for 5.35 × 10 5 km 2 (5.3%), was mainly distributed in the Southern region of the Balkan Peninsula and the Alps, and the annual SCDs value was >180.The variation in the SCD slope in Europe calculated using Sen's median method is shown in Figure 10.The increasing trend in the annual SCDs in Northern Europe and most regions of Western Europe was <5 day/year, while the regions with a decreasing rate of >5 day/year were mainly concentrated in the Alps and the Balkan Peninsula.The combination of Figures 9 and 10 shows that the overall SCD value in Europe did not exhibit an obvious upward or downward trend.
Remote Sens. 2018, 10, 136 15 of 22 exhibited a downward trend.This finding indicated that the annual SCD in the seasonal snow cover distribution area in Asia exhibited an upward trend; the annual SCD in the high-elevation mountainous regions with the presence of perennial snow exhibited a downward trend.
From 2000 to 2015, the areas in Europe where SCDs increased (Z > 0) and decreased (Z < 0) were essentially the same, i.e., 4.48 × 10 6 km 2 (44.4%) and 4.32 × 10 6 km 2 (42.8%), respectively.The area with significant SCDs increase, accounting for 3.23 × 10 5 km 2 (3.2%), was distributed in the Western European Plains, Britain, and Ireland, and the annual SCDs in these regions was <60 days.The area with a significant SCD decrease, accounting for 5.35 × 10 5 km 2 (5.3%), was mainly distributed in the Southern region of the Balkan Peninsula and the Alps, and the annual SCDs value was >180.The variation in the SCD slope in Europe calculated using Sen's median method is shown in Figure 10.The increasing trend in the annual SCDs in Northern Europe and most regions of Western Europe was <5 day/year, while the regions with a decreasing rate of >5 day/year were mainly concentrated in the Alps and the Balkan Peninsula.The combination of Figures 9 and 10 shows that the overall SCD value in Europe did not exhibit an obvious upward or downward trend.The results from the Mann-Kendall trend analysis indicated that the regions where SCDs in North America increased and decreased during 2000-2015 were 1.39 × 10 7 km 2 (56.4%) and 7.16 × 10 6 km 2 (29.0%), respectively.The area of the regions with a significant SCD enhancement was 6.67 × 10 5 km 2 (2.7%), which was distributed in the Saskatchewan region of Canada.The regions with a significant SCD decrease accounted for 1.75 × 10 6 km 2 (7.1%), and the regions were mainly distributed in the Cordillera Mountains of North America with some relatively scattered areas in the Central and Eastern areas of Canada.The increasing trend in SCDs in 37.9% of the area was <5 day/year and was mainly distributed in Southern Canada and the Northern region of mainland United States.The decreasing rate in SCDs in 27.0% of the area was <5 day/year, and the area was mainly in the Cordillera Mountains of North America and the Northern region of Canada with some sporadic areas in the Northern islands of Canada and Greenland.The decreasing rate in SCDs in a few regions in the Western mountainous area was >5 day/year.The results in Figures 9 and 10 indicate that the annual SCDs in most regions of North America exhibited an upward trend.The annual SCDs in North America exhibited a gradual increase with increasing latitude, and the regions where the annual SCDs exhibited a downward trend were mainly the stable SCA where SCDs were >180, including Northern Canada, the Western mountainous area of the United States, and Alaska.This result indicated that the annual SCDs for seasonal snow cover in North America exhibited an The results from the Mann-Kendall trend analysis indicated that the regions where SCDs in North America increased and decreased during 2000-2015 were 1.39 × 10 7 km 2 (56.4%) and 7.16 × 10 6 km 2 (29.0%), respectively.The area of the regions with a significant SCD enhancement was 6.67 × 10 5 km 2 (2.7%), which was distributed in the Saskatchewan region of Canada.The regions with a significant SCD decrease accounted for 1.75 × 10 6 km 2 (7.1%), and the regions were mainly distributed in the Cordillera Mountains of North America with some relatively scattered areas in the Central and Eastern areas of Canada.The increasing trend in SCDs in 37.9% of the area was <5 day/year and was mainly distributed in Southern Canada and the Northern region of mainland United States.The decreasing rate in SCDs in 27.0% of the area was <5 day/year, and the area was mainly in the Cordillera Mountains of North America and the Northern region of Canada with some sporadic areas in the Northern islands of Canada and Greenland.The decreasing rate in SCDs in a few regions in the Western mountainous area was >5 day/year.The results in Figures 9 and 10 indicate that the annual SCDs in most regions of North America exhibited an upward trend.The annual SCDs in North America exhibited a gradual increase with increasing latitude, and the regions where the annual SCDs exhibited a downward trend were mainly the stable SCA where SCDs were >180, including Northern Canada, the Western mountainous area of the United States, and Alaska.This result indicated that the annual SCDs for seasonal snow cover in North America exhibited an upward trend.However, the annual SCDs in the stable SCA distributed in the high-latitude areas and high-elevation mountainous areas exhibited a downward trend.The results derived in this study are consistent with the conclusions of the IPCC 4th assessment report that the warming of the Western mountainous regions in North America has caused the SCA to decrease.

Discussion
The synthesized snow product is a fusion of multisource data in order to improve the existing standard snow cover and snow water equivalent products.Optical remote sensing products have high spatial resolution, but they are seriously affected by clouds.Passive microwave data are available that penetrate clouds, but these data have coarse resolution.On this basis, the biggest innovation in this study is that the optical product MODIS with the assistance of the passive microwave remote sensing products AMSR-E and IMS were used to integrate the advantages of different cloud removal algorithms and generate a large-scale cloud-free daily snow cover product at high resolution (500 m) for the Northern Hemisphere.
This study of the variation characteristics of the temporal and spatial dynamic in snow cover covering 2000-2015 provides further evidence supporting significant snow declines in the Northern Hemisphere.This conclusion is also consistent with the previous studies, which indicated that the variation SCA in the Northern Hemisphere was well-defined during the 20th century [4,5,9,11,13,63].Snow cover trends are significantly dependent on latitude, elevation, and the increasing spring solar radiation that can reinforce snow melting processes [12].Average stable snow cover duration has decreased at a rate of 5.3 days decade −1 in the Northern hemisphere, and variations in this region are attributed primarily to a progressively earlier offset, which has advanced poleward at a rate of 5.5 days decade −1 [12].The decrease of −0.58 × 10 6 km 2 decade −1 in spring snow cover across the period 1967-2014 provides further evidence that the Northern Hemisphere is experiencing large reductions in springtime snow cover that are particularly strong over the past decade [64].Moreover, other researchers also reported that the earlier loss of spring snow cover in the Northern Hemisphere is linked to the rise in surface air temperatures [65].These spring air temperature increases will continue to advance snow cover declines and thus may further affect the spring snow offset dates in the Northern Hemisphere.
The data sources adopted in the above studies were mostly the National Oceanic and Atmospheric Administration's the National Environmental Satellite, Data and Information Service (NOAA's NESDIS) digitized weekly data of the SCA or re-analysis data with low spatial and temporal resolution.Although these data can reflect the trend in snow cover to some extent, the accuracy is not high, and there has not been an analysis on the latest variations in snow cover in the Northern Hemisphere for the 21st century.The findings of this study are consistent with previous studies, but, more importantly, this study contains more detailed information at a higher resolution.It should be noted that the conclusion of this study is the trends in places with shorter snow cover days are positive and those with longer are negative; the reason may be due to the inverse temperature trends in Northern Hemisphere continents [66].Nonetheless, further work is needed to analyze the relationship between variations in snow cover and temperatures under conditions of climate-warming in detail.
Utilizing the MODIS V005 snow product to produce a cloudless snow cover products is a relatively mature algorithms [37,38,[46][47][48].Still, the MODIS V005 snow product is a bit outdated and the V006 product offers the NDSI instead.For future work, we plan to introduce the new V006 product in the later time series.Meanwhile, a longer time series of data is needed to more clearly understand further definitive conclusions concerning trend characteristics in the Northern Hemisphere.
The accuracy of integrated products is consistent with the standard products in the areas without the cloud cover.On this basis, cloud pixels were reclassified for the cloud-covered areas.Although the combination of MODIS and AMSR-E/IMS can take advantage of both high spatial resolution of optical data and cloud transparency of passive microwave data, the combination might reduce the spatial resolution of MODIS products because cloud pixels in MODIS are replaced by AMSR-E/IMS images, which have coarse spatial resolution; therefore, the accuracy of the cloud pixel reclassification fully depends on the of AMSR-E and IMS.In previous studies, AMSR-E and IMS data were widely used for monitoring snow, and both of them have relatively high accuracy [28][29][30][31].On the other hand, because of the poor performance in the forest area, the Kappa test of integrated snow products shows moderate consistency with Landsat, and snow monitoring in this regions also remains a major problem for virtually all types of global remote sensing snow products, for which overall accuracy were seriously affected due to the forest canopy [37].The overall result is still not satisfying, although the accuracy of snow monitoring can be improved by correcting for the vegetation index [67].For future improvement of the snow mapping algorithm, the influence of snow physical characteristics, forest density and solar radiation in forest regions needs to be further discussed in order to improve the accuracy of large-scale snow monitoring.

Conclusions
Along with the constant threat of global climate warming, studies related to cryosphere change have become hot topics with increasing interest from researchers.Large-scale variations in snow cover are important indicators of climate change, and snow cover can also affect global-and regional-scale radiation and energy equilibrium because of its high reflectivity.In this study, we combined the MODIS standard daily snow cover products MOD10A1 and MYD10A1 with the passive microwave remote sensing data of the AMSR-E SWE product and the data of the multi-source remote sensing product IMS to develop a new snow cover mapping algorithm to obtain a cloud-free daily snow cover product for the Northern Hemisphere.Moreover, we used Landsat-TM images to validate the accuracy of the cloud-free daily snow cover product synthesized in this study, and the average Kappa was 0.581 with moderate overall consistency.By completely removing cloud interference, we effectively improved the capability to correctly monitor snow cover range over a large scale, although different land cover types can affect the accuracy of the product.This new product was used to analyze the characteristics of the temporal and spatial dynamic variations in snow cover across the Northern Hemisphere, including SCA, SCDs, and snow cover inter-annual, monthly, and seasonal variations.The conclusions of this study include the following aspects: (1) The maximum, minimum, and annual average SCA in the Northern Hemisphere all exhibited a downward trend with fluctuations.The variation trend in snow cover in the Northern Hemisphere exhibited significant inter-annual and regional differences.(2) The largest average SCA in the Northern Hemisphere occurred in January, while the smallest was in August.For the monthly variations, between 2000 and 2015, the SCA in the Northern Hemisphere for January, July, and October exhibited an upward trend, while the SCA in February, March, April, June, August, and December exhibited a downward trend.As for the seasonal variation, the SCA in the Northern Hemisphere during the spring, summer, and fall seasons exhibited a downward trend.The reduction magnitude of the SCA was relatively large in the spring and summer, and the SCA in winter did not exhibit an obvious trend.(3) The spatial distribution of the annual average SCDs was based on latitudinal zonality.The upward trend region was mainly confined to the mid-to low-latitude seasonal SCA with unstable snow cover.
(4) A decreasing trend in the SCA was observed in regions with perennial snow cover, including the high-latitude or high-elevation mountainous regions in the Northern Hemisphere (between 35 • and 50 • N), such as the Tibetan Plateau, the Tianshan Mountains, the Pamir Plateau in Asia, the Alps in Europe, the Caucasus Mountains, and the Cordillera Mountains in North America.Overall, the perennial snow in the Northern Hemisphere is transitioning to seasonal snow cover.

Figure 1 .
Figure 1.Moderate Resolution Imaging Spectroradiometer (MODIS) daily cloud-free snow cover products flow chart.Note: Tij, Yij, and Nij, represent the pixel values of the ith column and the jth row of the current day image, previous day image, and posterior day image, respectively.

Figure 1 .
Figure 1.Moderate Resolution Imaging Spectroradiometer (MODIS) daily cloud-free snow cover products flow chart.Note: T ij , Y ij , and N ij, represent the pixel values of the ith column and the jth row of the current day image, previous day image, and posterior day image, respectively.

Figure 2 .
Figure 2. Comparison of the cloud-free daily snow product and Landsat Thematic Mapper (TM) for different land cover types.

Figure 2 .
Figure 2. Comparison of the cloud-free daily snow product and Landsat Thematic Mapper (TM) for different land cover types.

Figure 3 .
Figure 3. Examples of snow-cover maps from MODIS composites and Landsat-TM on 5 May 2007.

Figure 3 .
Figure 3. Examples of snow-cover maps from MODIS composites and Landsat-TM on 5 May 2007.

4. 2 .
Figure4shows the variation curve of the maximum, average and minimum SCA in the Northern Hemisphere from 2000 to 2015.The maximum SCA in the Northern Hemisphere rapidly decreased in 2007, and it declined with fluctuations after reaching a peak in 2008.The variation in the minimum SCA of the Northern Hemisphere is of particular interest because the area contains perennial snow cover.Dynamic monitoring of the perennial snow area in the Northern Hemisphere plays a vital role in monitoring global climate and hydrological changes.Figure4shows that both the maximum and minimum SCA in the Northern Hemisphere exhibited a downward trend during the study period, and the tendency rates were −0.057% a −1 and −0.045% a −1 , respectively.Since 2000, the inter-annual variation in the SCA of the Northern Hemisphere has exhibited a downward trend with fluctuations, with an average rate of −0.066% a −1 .

Figure 4 .
Figure 4. Annual maximum, average and minimum snow-covered area in the Northern Hemisphere from 2000 to 2015.

Figure 5
Figure5presents a box plot for the monthly SCA in the Northern Hemisphere from 2000 to 2015.This figure shows that the largest and smallest SCA in the Northern Hemisphere occurred in January and August, respectively, i.e., approximately 5.5 × 10 7 km 2 (55.38%) and 5.76 × 10 6 km 2 (5.80%).The 4th report of the IPCC noted that the largest average SCA in the Northern Hemisphere typically occurs in January, while the smallest occurs in August[62], and this study is consistent with these conclusions.In addition, based on the figure, the SCA fluctuation in the Northern Hemisphere was relatively large from February to June, and the fluctuation exhibited a negatively skewed distribution from February to April, which was mainly concentrated in the spring.This also indicated that the spring SCA had a downward trend during the study period.The 5th assessment report (AR5) of the IPCC noted that, since the middle 20th century, the SCA in the Northern Hemisphere has significantly decreased, indicating that the main cause is the reduction in the SCA in the spring[19].

Figure 4 .
Figure 4. Annual maximum, average and minimum snow-covered area in the Northern Hemisphere from 2000 to 2015.

Figure 5
Figure 5 presents a box plot for the monthly SCA in the Northern Hemisphere from 2000 to 2015.This figure shows that the largest and smallest SCA in the Northern Hemisphere occurred in January and August, respectively, i.e., approximately 5.5 × 10 7 km 2 (55.38%) and 5.76 × 10 6 km 2 (5.80%).The 4th report of the IPCC noted that the largest average SCA in the Northern Hemisphere typically occurs in January, while the smallest occurs in August[62], and this study is consistent with these conclusions.In addition, based on the figure, the SCA fluctuation in the Northern Hemisphere was relatively large from February to June, and the fluctuation exhibited a negatively skewed distribution from February to April, which was mainly concentrated in the spring.This also indicated that the spring SCA had a downward trend during the study period.The 5th assessment report (AR5) of the IPCC noted that, since the middle 20th century, the SCA in the Northern Hemisphere has significantly decreased, indicating that the main cause is the reduction in the SCA in the spring[19].

Figure 5 .
Figure 5. Box plots of the monthly averages of snow-covered areas between 2000 and 2015.

Figure 5 .
Figure 5. Box plots of the monthly averages of snow-covered areas between 2000 and 2015.

Figure 6 .
Figure 6.The average area proportion of snow cover in each season in the Northern Hemisphere from 2000 to 2015.

Figure 6 .
Figure 6.The average area proportion of snow cover in each season in the Northern Hemisphere from 2000 to 2015.

22 ( 32 .
03%) in 2000 to 3.46 × 10 7 km 2 (34.86%) in 2015.The area with 120-180 SCDs also exhibited an upward trend.The region with 10 < SCDs ≤ 60 and SCDs ≥ 180 exhibited a downward trend.The regions with SCDs ≥ 350 were covered by perennial snow or glaciers throughout the year, and this part of the snow cover is an important indicator of global climatic change.There was not an obvious change in regions with other SCDs.

Figure 7 .
Figure 7. Area percentage of the different snow-covered days during the 16-year period from 2000 to 2015 in the Northern Hemisphere.

Figure 7 .
Figure 7. Area percentage of the different snow-covered days during the 16-year period from 2000 to 2015 in the Northern Hemisphere.
× 10 7 km 2 (60.1%), and it was mainly distributed in the Northeastern region of Inner Mongolia of China, North of Xinjiang, the Tibetan Plateau, Mongolia, most of Russia, Northern Kazakhstan, and the Eastern region of Turkey.The perennial SCA with SCDs > 350 only accounted for 7.22 × 10 3 km 2 (0.02%) of the total SCA, and the areas were mainly distributed in the regions of Severnaya Zemlya in Russia, West of the Tianshan Mountains in Xinjiang (China), the Kunlun and Himalayan Mountains on the Tibetan Plateau, and the Pamir Plateau.Remote Sens. 2018, 10, 136 13 of 22 Plateau, Mongolia, most of Russia, Northern Kazakhstan, and the Eastern region of Turkey.The perennial SCA with SCDs > 350 only accounted for 7.22 × 10 3 km 2 (0.02%) of the total SCA, and the areas were mainly distributed in the regions of Severnaya Zemlya in Russia, West of the Tianshan Mountains in Xinjiang (China), the Kunlun and Himalayan Mountains on the Tibetan Plateau, and the Pamir Plateau.

Figure 8 .
Figure 8. Spatial distribution of the average annual snow-covered days (SCD) during the period 2000-2015 in the Northern Hemisphere.

Figure 8 .
Figure 8. Spatial distribution of the average annual snow-covered days (SCD) during the period 2000-2015 in the Northern Hemisphere.

Figure 9 .
Figure 9. Significance of the variation in the annual SCDs in the Northern Hemisphere based on the Mann-Kendall method.

Figure 9 .
Figure 9. Significance of the variation in the annual SCDs in the Northern Hemisphere based on the Mann-Kendall method.

Figure 10 .
Figure 10.Variation slope of the average annual SCDs in the Northern Hemisphere based on Sen's median method from 2000 to 2015.

Figure 10 .
Figure 10.Variation slope of the average annual SCDs in the Northern Hemisphere based on Sen's median method from 2000 to 2015.

Table 1 .
The information of remote sensing data.

Table 2 .
Percentage of different snow-covered days by area in Asia, Europe and North America.

Table 3 .
Percent variation by area in Asia, Europe and North America based on the Mann-Kendall method.