Distribution and Attribution of Terrestrial Snow Cover Phenology Changes over the Northern Hemisphere during 2001–2020

: Snow cover phenology has exhibited dramatic changes in the past decades. However, the distribution and attribution of the hemispheric scale snow cover phenology anomalies remain unclear. Using satellite-retrieved snow cover products, ground observations, and reanalysis climate variables, this study explored the distribution and attribution of snow onset date, snow end date, and snow duration days over the Northern Hemisphere from 2001 to 2020. The latitudinal and altitudinal distributions of the 20-year averaged snow onset date, snow end date, and snow duration days are well represented by satellite-retrieved snow cover phenology matrixes. The validation results by using 850 ground snow stations demonstrated that satellite-retrieved snow cover phenology matrixes capture the spatial variability of the snow onset date, snow end date, and snow duration days at the 95% signiﬁcance level during the overlapping period of 2001–2017. Moreover, a delayed snow onset date and an earlier snow end date (1.12 days decade − 1 , p < 0.05) are detected over the Northern Hemisphere during 2001–2020 based on the satellite-retrieved snow cover phenology matrixes. In addition, the attribution analysis indicated that snow end date dominates snow cover phenology changes and that an increased melting season temperature is the key driving factor of snow end date anomalies over the NH during 2001–2020. These results are helpful in understanding recent snow cover change and can contribute to climate projection studies.


Introduction
Snow cover is an integral component of the cryosphere and represents one of the three Essential Climate Variables related to snow for the Global Observing System for Climate, which plays a crucial role in the Earth's climate system through the surface energy budget [1][2][3], atmospheric circulation [4], and hydrological cycle [5][6][7], and influences freshwater resources across a large proportion of the Northern Hemisphere (NH) [6]. Snow cover phenology (SCP) variables including snow onset date (D o ), snow end date (D e ), and snow duration days (D d ) are key indicators of seasonal variation of terrestrial snow cover over the NH and becoming increasingly valuable indicators of climate change [6,8], especially in snow-dominated cold regions [6,9]. For example, SCP has considerable impact on climate variabilities, such as alpine vegetation growth dynamics on the Tibetan Plateau [10], green-up date across the NH [11], boreal springtime carbon uptake [12], permafrost degradation in sub-arctic Sweden [9], and SCP indicators are expected to provide feedback to temperature trends [13]. Moreover, abnormal snowmelt timing in spring are creating a serious threat to water resource sustainability, including agricultural how these climate change events, which exhibit high spatial heterogeneity, influence SCP. Such information can contribute to snow water management, the sustainable development of ecosystems, and predictions of catastrophic climate-related events.
To explore the latest terrestrial snow cover phenology changes and its attribution factors under the rapid climate change background, this study integrates 8-day Level 3 snow cover fraction products derived from the Moderate Resolution Imaging Spectroradiometer Satellite (MOD10C2) from September 2000 to August 2020, at a spatial resolution of 0.05 • (approximately 5 km), as well as the binary snow mask from the Interactive Multi-sensor Snow and Ice Mapping System (IMS) during 2005-2019, with a spatial resolution of 4 km. These data are used to produce a 20-year gap-free MOD10C2-based snow cover extent (GF-MOD10C2-SCE) with finer spatial resolution over the NH. Then, using the integrated GF-MOD10C2-SCE, this study retrieves SCP matrices over the NH, including D o , D e , and D d . Moreover, through reanalysis of the land surface temperature and precipitation, this study explores the SCP anomalies over the NH in the past two decades.

Study Area
To focus on changes in SCP, the study area was confined to seasonal snow cover gridcells with snow cover greater than 7 frames using 8-day MOD10C2 for 75% of the years between 2001 and 2020, similar with Choi et al. [22]. Permanent snow-covered regions with little change in SCP are beyond the scope of this study. Moreover, areas of thin or patchy snow cover were excluded from the analysis because they may have been missed in the optical satellite images [34]. The distribution of study area and excluded area calculated from GF-MOD10C2-SCE from 2001 to 2020 is shown in Figure 1. mid-latitudes of the NH over the past decade [33]. It is thus crucial to understand how these climate change events, which exhibit high spatial heterogeneity, influence SCP. Such information can contribute to snow water management, the sustainable development of ecosystems, and predictions of catastrophic climate-related events.
To explore the latest terrestrial snow cover phenology changes and its attribution factors under the rapid climate change background, this study integrates 8-day Level 3 snow cover fraction products derived from the Moderate Resolution Imaging Spectroradiometer Satellite (MOD10C2) from September 2000 to August 2020, at a spatial resolution of 0.05° (approximately 5 km), as well as the binary snow mask from the Interactive Multisensor Snow and Ice Mapping System (IMS) during 2005-2019, with a spatial resolution of 4 km. These data are used to produce a 20-year gap-free MOD10C2-based snow cover extent (GF-MOD10C2-SCE) with finer spatial resolution over the NH. Then, using the integrated GF-MOD10C2-SCE, this study retrieves SCP matrices over the NH, including Do, De, and Dd. Moreover, through reanalysis of the land surface temperature and precipitation, this study explores the SCP anomalies over the NH in the past two decades.

Study Area
To focus on changes in SCP, the study area was confined to seasonal snow cover gridcells with snow cover greater than 7 frames using 8-day MOD10C2 for 75% of the years between 2001 and 2020, similar with Choi, et al. [22]. Permanent snow-covered regions with little change in SCP are beyond the scope of this study. Moreover, areas of thin or patchy snow cover were excluded from the analysis because they may have been missed in the optical satellite images [34]. The distribution of study area and excluded area calculated from GF-MOD10C2-SCE from 2001 to 2020 is shown in Figure 1.

Datasets
Three satellite-observed snow cover products, 850 ground-based daily snow depth observations, and a reanalysis dataset were employed to explore the distribution and attribution of terrestrial SCP anomalies over the NH.

Satellite-observed Snow Cover Datasets
Satellite-observed snow cover datasets are preliminary input data for SCP detection over the NH. Three snow cover datasets, including snow charts derived from the Northern Hemisphere Snow Cover Extent Climate Data Record (NH SCE CDR) v01r01 [35], MOD10C2 [36], and IMS [37] were employed in this study.

Snow Cover Charts
Monthly snow charts over the NH from 1966 through 2020, calculated from the NH SCE CDR v01r01 [35] were used to investigate the long-term SCE anomaly. These snow

Datasets
Three satellite-observed snow cover products, 850 ground-based daily snow depth observations, and a reanalysis dataset were employed to explore the distribution and attribution of terrestrial SCP anomalies over the NH.

Satellite-Observed Snow Cover Datasets
Satellite-observed snow cover datasets are preliminary input data for SCP detection over the NH. Three snow cover datasets, including snow charts derived from the Northern Hemisphere Snow Cover Extent Climate Data Record (NH SCE CDR) v01r01 [35], MOD10C2 [36], and IMS [37] were employed in this study.

Snow Cover Charts
Monthly snow charts over the NH from 1966 through 2020, calculated from the NH SCE CDR v01r01 [35] were used to investigate the long-term SCE anomaly. These snow chart data are widely used in large-scale SCE anomalies over the NH due to their consistency and long time span [35].

MOD10C2 Snow Cover Fraction Dataset
The MOD10C2 provides the 8-day composite global snow cover fraction (SCF) at 0.05 • from February 2000 to the present. MOD10C2 is an aggregation of MOD10A2 products with 500-m spatial resolution, which is an 8-day composite of MOD10A1 daily SCE maps [36]. The 8-day composite is considered useful because persistent cloudiness limits the number of days available for surface observations in many regions, particularly at high latitudes [38]. In this study, the MOD10C2 was used as the primary data in SCP retrieval because MOD10C2 is the only consistent, objective snow estimate derived from optical satellite observations over the last two decades with relatively finer spatial resolution.

IMS Datasets
The IMS snow cover product provides daily binary SCE maps for the NH from February 1997 to December 2019 at three different resolutions, i.e., 1 km (since 2015), 4 km (since 2004), and 24 km (since 1997), which are created manually by a snow analyst using data from a combination of geostationary and polar orbiting satellites in visible, infrared, and microwave spectrums [37]. The IMS dataset is widely used for data integration [39,40]. The daily rate of agreement between IMS snow maps and ground-based snow observations between 2006 and 2010 was predominantly between 80% and 90% throughout the winter seasons over the continental United States [41]. To fill the gaps in MOD10C2 at the highest achievable spatial resolution, this study used the IMS at a spatial resolution of 4 km from 2005 to 2019.

Ground-Based Snow Depth Observations
Daily ground-based snow depth measurements generated by the Global Historical Climatology Network (GHCN) [42] were used to detect the ground-observed SCP and were employed as ground truth to verify the performance of the satellite-retrieved SCP in this study. The daily GHCN contains records from over 75,000 stations in 180 countries and territories [42]. By assembling and checking observations made in multiple different nations, the GHCN provides daily snow depth records over the NH, especially at high latitudes. As they are subject to the speed of data updates, GHCN data for the last three years are incomplete. Therefore, only 850 stations covering September 2000 to August 2017 with D d greater than 60 days were used for validation. The distribution of the selected GHCN snow depth observations are shown in Figure 2. sistency and long time span [35].

MOD10C2 Snow Cover Fraction Dataset
The MOD10C2 provides the 8-day composite global snow cover fraction (SCF) at 0.05° from February 2000 to the present. MOD10C2 is an aggregation of MOD10A2 products with 500-m spatial resolution, which is an 8-day composite of MOD10A1 daily SCE maps [36]. The 8-day composite is considered useful because persistent cloudiness limits the number of days available for surface observations in many regions, particularly at high latitudes [38]. In this study, the MOD10C2 was used as the primary data in SCP retrieval because MOD10C2 is the only consistent, objective snow estimate derived from optical satellite observations over the last two decades with relatively finer spatial resolution.

IMS Datasets
The IMS snow cover product provides daily binary SCE maps for the NH from February 1997 to December 2019 at three different resolutions, i.e., 1 km (since 2015), 4 km (since 2004), and 24 km (since 1997), which are created manually by a snow analyst using data from a combination of geostationary and polar orbiting satellites in visible, infrared, and microwave spectrums [37]. The IMS dataset is widely used for data integration [39,40]. The daily rate of agreement between IMS snow maps and ground-based snow observations between 2006 and 2010 was predominantly between 80% and 90% throughout the winter seasons over the continental United States [41]. To fill the gaps in MOD10C2 at the highest achievable spatial resolution, this study used the IMS at a spatial resolution of 4 km from 2005 to 2019.

Ground-based Snow Depth Observations
Daily ground-based snow depth measurements generated by the Global Historical Climatology Network (GHCN) [42] were used to detect the ground-observed SCP and were employed as ground truth to verify the performance of the satellite-retrieved SCP in this study. The daily GHCN contains records from over 75,000 stations in 180 countries and territories [42]. By assembling and checking observations made in multiple different nations, the GHCN provides daily snow depth records over the NH, especially at high latitudes. As they are subject to the speed of data updates, GHCN data for the last three years are incomplete. Therefore, only 850 stations covering September 2000 to August 2017 with Dd greater than 60 days were used for validation. The distribution of the selected GHCN snow depth observations are shown in Figure 2.

Reanalysis Temperature and Precipitation Dataset
Model resolution plays an important role in representing the orographic gradient effects on temperature and precipitation simulations [43]. ERA5-Land is a replay of the land

Reanalysis Temperature and Precipitation Dataset
Model resolution plays an important role in representing the orographic gradient effects on temperature and precipitation simulations [43]. ERA5-Land is a replay of the land component of the ERA5 climate reanalysis with a finer spatial resolution (0.10 • grid spacing) than ERA5 monthly averaged data on single levels from 1979 to the present (0.25 • spatial resolution) [44] and the Climatic Research Unit (CRU) Time-Series (TS) Version 4.04 (0.50 • spatial resolution) [45]. To address the changes of SCP over the NH, the monthly averaged 2-m temperature and total precipitation derived from the European Centre for Medium Range Weather Forecasts (ECMWF) Reanalysis v5-Land (ERA5-Land) [44] are used as surface air temperature and precipitation in the attribution analysis. Moreover, the model used in the production of ERA5-Land is the tiled ECMWF scheme for surface exchanges over land, which incorporates land surface hydrology (H-TESSEL) and has been proven suitable for modeling the ground [46,47].

Data Preparation
Details of datasets used in this study are provided in Table 1. Abbreviations used in this study are listed in Table A1. To match the spatial resolution of the satellite-retrieved SCP matrixes, the monthly averaged 2-m temperature and total precipitation derived from ERA5-Land were regridded at a spatial resolution of 0.05 • by using the resampling method of "cubic-spline" with the help of gdalwarp (http://www.gdal.org/gdalwarp.html).

Generation of Gap-Free MOD10C2-Based SCE Dataset
To overcome the shortage of optical MOD10C2 in snow discrimination caused by invalid observations and cloud contamination, we produced GF-MOD10C2-SCE dataset before SCP retrieving. First, using spatial complete IMS from 2005 to 2019, we produced an 8-day snow cover probability (SProb) map in each 8-day period of the year. For gridcell i in a given 8-day period of the year, the SProb i was calculated by averaging non-zero values in this 8-day period from 2005 to 2019. Second, the SProb map was used to fill the gaps in MOD10C2. If SCF i = 0 and SProb i > 50%, the SProb i will replace SCF i in original MOD10C2. Finally, the gridcells with SCF > 0 were used as the GF-MOD10C2-SCE and employed to retrieve the SCP matrixes over the NH. Using homogeneous GF-MOD10C2-SCE in SCP retrieval will reduce the error caused by missing observations in the original MOD10C2. The flowchart of GF-MOD10C2-SCE generation is presented in Figure 3.

Snow Cover Phenology Retrieval
We used the hydrological year in this study to explore interannual anomalies in the SCP. Based on the seasonal cycles of snow cover for northern Eurasia and North America [48], this study defined the hydrological year (t) as the period from September in year (t-1) to August in year (t). To further attribute changes in SCP, we defined the accumulation season of snow cover from September to December and the melting season from March to June in the attribution analysis.

Snow Cover Phenology Retrieval
We used the hydrological year in this study to explore interannual anomalies in the SCP. Based on the seasonal cycles of snow cover for northern Eurasia and North America [48], this study defined the hydrological year (t) as the period from September in year (t-1) to August in year (t). To further attribute changes in SCP, we defined the accumulation season of snow cover from September to December and the melting season from March to June in the attribution analysis.

Definition of Satellite-Retrieved Snow Cover Phenology
Published study have used full snow seasons and core snow seasons in SCP identification [22], in which full snow seasons is the interval between the first appearance and the last disappearance of snow cover and core snow seasons is the longest interval of the season with an unbroken string of weeks for which a grid cell is snow covered. In this study, we mainly focus on SCP changes in core snow seasons.
The definition of D o , D e , and D d in a hydrological year using original 8-day MOD10C2 SCF dataset is shown in Figure 4. To avoid the impact of missing observations (D-G) on SCP retrieval, we used GF-MOD10C2-SCE in SCP retrieve. For the 8-day GF-MOD10C2-SCE dataset, we first identified the 8-day interval frames (day index i to i + 7) between the frames when the SCE was superior to zero and when it equaled zero in core snow seasons, then defined D o and D e at i+3.5 days of the first frame when the SCE was superior to zero and when it equaled zero in core snow seasons, respectively.

Snow Cover Phenology Retrieval
We used the hydrological year in this study to explore interannual anomalies in the SCP. Based on the seasonal cycles of snow cover for northern Eurasia and North America [48], this study defined the hydrological year (t) as the period from September in year (t-1) to August in year (t). To further attribute changes in SCP, we defined the accumulation season of snow cover from September to December and the melting season from March to June in the attribution analysis.

Definition of Satellite-retrieved Snow Cover Phenology
Published study have used full snow seasons and core snow seasons in SCP identification [22], in which full snow seasons is the interval between the first appearance and the last disappearance of snow cover and core snow seasons is the longest interval of the season with an unbroken string of weeks for which a grid cell is snow covered. In this study, we mainly focus on SCP changes in core snow seasons.
The definition of Do, De, and Dd in a hydrological year using original 8-day MOD10C2 SCF dataset is shown in Figure 4. To avoid the impact of missing observations (D-G) on SCP retrieval, we used GF-MOD10C2-SCE in SCP retrieve. For the 8-day GF-MOD10C2-SCE dataset, we first identified the 8-day interval frames (day index i to i+7) between the frames when the SCE was superior to zero and when it equaled zero in core snow seasons, then defined Do and De at i+3.5 days of the first frame when the SCE was superior to zero and when it equaled zero in core snow seasons, respectively.

Definition of Ground-based Snow Cover Phenology
To avoid the impact of ephemeral snow on SCP retrieval and match the temporal resolution of 8-day GF-MOD10C2-SCE, Do was defined as the first four consecutive days

Definition of Ground-Based Snow Cover Phenology
To avoid the impact of ephemeral snow on SCP retrieval and match the temporal resolution of 8-day GF-MOD10C2-SCE, D o was defined as the first four consecutive days in the accumulation season and D e was defined as the last four consecutive days of persistent snow cover in the melting season. D d was defined as the number of days between D o and D e .

Validation of Satellite-Retrieved Snow Cover Phenology Using Ground Observations
Although there are large differences in the spatial representation of satellite-retrieved SCP and ground-based SCP, the ground observations still provide the most convincing results when compared with satellite-retrieved SCP. In this study, the root-mean-square error (RMSE), mean relative error (MRE), and bias were used as criteria to evaluate the relative accuracy of the satellite-retrieved SCP (X sat ) relative to ground-based SCP (X obs ) during 2001-2017. The RMSE, MRE, and bias of X sat relative to X obs are expressed as follows: where X sat,i and X obs,i are the satellite-retrieved SCP and ground-based SCP in a given gridcell i, respectively.

Attribution Analysis
As anomalies in D d are directly correlated to variabilities in D o and D e , we mainly focus on the attribution analysis of D o and D e in this study. To compare the contributions of D o and D e from each climate variable consistently, all variables were converted to standardized anomalies (z-score) using the mean and standard deviation in the attribution analysis.
Published studies have shown that changes in SCP are largely determined by surface air temperature and precipitation variability [6,13]. Therefore, we assumed that the interannual variability of D o was driven by the competing effects of surface air temperature (T a ) and precipitation (P a ) in the snow accumulation seasons. To attribute changes in D o , we regressed D o as the dependent variable with T a and P s as the independent variables, following Chen et al. [24].
Here, β 1 and β 2 are regression coefficients for T a and P a , and ε 1 is the residual. The contributions of T a and P a to D o anomalies are reflected by the terms β 1 × T a and β 2 × P a , respectively. The contributions of T a and P a to D o were computed by regressing the annual time series of D o z-scores against the time series of z-scores of T a and P a . The resulting regression coefficients β 1 and β 2 were multiplied by T a and P a z-scores to derive the contributions of T a and P a to the D o z-scores.
Moreover, D e mainly depends on the maximum spring snow depth to melt and the surface air temperature in the snow melting season (T m ), in which the maximum snow depth in spring is determined using T a and P a . Therefore, we assumed that the interannual variability of D e was driven by T a , P a , and T m . We regressed D e as the dependent variable, with T a , P a , and T m as the independent variables, as performed in Chen et al. [24].
Here, β 3 , β 4 , and β 5 are regression coefficients for T a , P a , and T m , respectively, and ε 2 is the residual. The contributions of T a , P a , and T m to the D e anomalies are reflected by the terms β 3 × T a , β 4 × P a , and β 5 × T m , respectively. Similarly, the contributions of T a , P a , and T m to D e were computed by regressing the annual and zonal time series of D e z-scores against the time series of z-scores of T a , P a , and T m . The resulting regression coefficients were multiplied by the time series of T a , P a , and T m z-scores to derive the contributions of T a , P a , and T m to the D e z-scores. A similar formula was employed by Peng et al. [13] and Chen et al. [24] to identify the sensitivity of SCP to climate variabilities over the NH.

Results
To explore changes in SCP over the NH, we first analyzed the monthly SCE changes over the NH from 1966 to 2020. Then, we mapped the latest 20-year averaged SCP over the NH from 2001 to 2020 and explored its changes and attribution during the corresponding period.

Observed Long-Term Anomalies in Snow Cover Extent over the NH
The monthly SCE anomalies from 1966 to 2020 are shown in Figure 5, and detailed changes in the monthly SCE over the NH are presented in Table 2 August (Table 2), with a maximum SCE reduction of up to 30% in June (−3.22 10 6 km 2 ). In contrast to the negative changes of SCE in the melting seasons, SCEs largely increased during the accumulation season from September to November; the maximum SCE increase occurred in October (2.29 10 6 km 2 ), with an amplitude of up to 12.96% greater than the monthly averaged October SCE in the period of 1966-2000. Increased SCE in the snow accumulation season coincides with increasing snow cover and widespread boreal winter cooling since 1990 [33]. Moreover, the contrasting seasonal anomalies in SCE directly corresponded to SCP anomalies during 2001-2020.
the NH from 2001 to 2020 and explored its changes and attribution during the correspond-ing period.

Observed Long-Term Anomalies in Snow Cover Extent over the NH
The monthly SCE anomalies from 1966 to 2020 are shown in Figure 5, and detailed changes in the monthly SCE over the NH are presented in Table 2. Compared with the monthly averaged SCE in the earlier period of 1966-2000, the later period of 2001-2020 exhibited large differences; the SCE decreased during the melting season from May to August (Table 2), with a maximum SCE reduction of up to 30% in June (-3.22 10 6 km 2 ). In contrast to the negative changes of SCE in the melting seasons, SCEs largely increased during the accumulation season from September to November; the maximum SCE increase occurred in October (2.29 10 6 km 2 ), with an amplitude of up to 12.96% greater than the monthly averaged October SCE in the period of 1966-2000. Increased SCE in the snow accumulation season coincides with increasing snow cover and widespread boreal winter cooling since 1990 [33]. Moreover, the contrasting seasonal anomalies in SCE directly corresponded to SCP anomalies during 2001-2020.

Climatology of Snow Cover Phenology over the NH from 2001 to 2020
To explore the distribution of SCP over the NH, we first generated 20-year averaged D o , D e , and D e values over the NH from 2001 to 2020 from the implemented GF-MOD10C2-SCE dataset.

Climatology of Snow Cover Phenology over the NH from 2001 to 2020
To explore the distribution of SCP over the NH, we first generated 20-year averaged Do, De, and De values over the NH from 2001 to 2020 from the implemented GF-MOD10C2-SCE dataset.

Climatology of Snow Cover Phenology from 2001 to 2020
The spatial distributions of the 20-year averaged Do, De, and Dd over the NH from 2001 to 2020 are shown in Figure 6. The 20-year averaged Do, De, and Dd over the study area were day of year (DOY) 279.69 (±85.37), DOY 109.13 (±48.75), and 166.29 (±75.29) days, respectively, based on satellite-retrieved SCP matrixes.  As shown in Figure 6, satellite-derived SCP clearly describes the spatial patterns of D o , D e , and D d (Figure 6c) from low latitudes to high latitudes over the NH. In addition, the impacts of elevation on SCP attributes were also well represented in the satelliteretrieved SCP maps. As a result of the temperature gradient over the NH, the snow cover appears earlier in high-latitude areas around the Arctic and high-altitude regions including the Rocky Mountains, Tibetan Plateau, Sayan Mountains, and Yablonovy Mountains. Comparably, snow cover appears later in the central United States and low latitudes of Eurasia ( Figure 6a). As shown in Figure 6b, the spatial distribution of D e is generally the opposite to that of D o , i.e., an earlier D o typically corresponds to a later D e . The combined spatial distributions of D o and D e result in longer D d in grid cells located in high latitudes and high-altitude regions (Figure 6c), but shorter D d in areas including lower-latitude regions of the United States, Eurasia, and China.

Validating Snow Cover Phenology Using In Situ Observations
Subject to the data availability of GHCN daily snow depth observations, the accuracy of satellite-retrieved SCP using ground observations was estimated for the period of  Previous research has proved that ground observations yield results that are highly dependent on the specific location (latitude and elevation). Therefore, results from ground observations mostly reflect the local conditions instead of meaningful information about the climate [49]. However, as shown in Figure 7  Previous research has proved that ground observations yield results that are highly dependent on the specific location (latitude and elevation). Therefore, results from ground observations mostly reflect the local conditions instead of meaningful information about the climate [49]. However, as shown in Figure 7, the satellite-retrieved Do, De, and Dd values generally captured the spatial variability of Do, De, and Dd over the NH, which proves the suitability of satellite-retrieved SCP for current climate change studies, also validated by the latitudinal distributions of Do, De, and Dd in Figure 7d-f.   cover (-2.63 ± 7.26 days) over the NH from 2001 to 2020 (Figure 8c and Figure 9c). In addition to the central United States, Tibetan Plateau, and Scandinavian Mountains, most grid cells within the study presented negative changes in Dd over the last two decades. Moreover, the observed Dd increase in Central North America is consistent with the increased snow cover days reported by Wang, et al. [50] from 2001 to 2015 and Allchin and Déry [51] from 1971 to 2014.

Attribution of Snow Cover Phenology Changes over the NH
The snow accumulation season temperature, Ta, precipitation, Pa, and snow melting season temperature, Tm, were employed to attribute the anomalies in the SCP over the NH during 2001-2020. The spatial distributions of changes in Ta, Pa, and Tm are shown in Figure 10a-c. The attribution analysis of the SCP is presented in Figure 10d-f. The Ta ( Figure  10a) and Tm (Figure 10b) values displayed a general warming trend over the NH from 2001 to 2020. This is reasonable according to the latest WMO Statement on the State of the Global Climate [32], as the past five years (2015-2019) are the five warmest on record. Moreover, both Ta and Tm exhibited significant latitude differences under the influence of decreasing snow cover and sea-ice extent [25][26][27][28], in which the temperature increased at a higher rate in high latitudes and at a much lower rate in middle latitudes. The anomalies in Ta and Tm over the NH during 2001-2020 agree with the land surface temperature anomalies mapped from the Goddard Institute for Space Studies (GISS) (https://data.giss.nasa.gov/gistemp/maps/) based on the latest GHCN Monthly version 4 [52]. In addition, changes in atmospheric circulation [29] and human activity [30,31] resulted in a notable increase in precipitation during the snow accumulation season, Pa (Figure 10c), in the Central United States, Tibetan Plateau, Northern Europe, and Northeastern Russia, but led to a significant decline in Pa in Central Russia, Northeast Canada, and Southern Europe.
The interannual variability analysis revealed a shortening trend in Dd over the NH from 2001 through 2020 at a rate of -1.64 days decade -1 (p < 0.05). Owing to an exceptionally strong El Niño, 2016 was the warmest year on record in all major global surface temperature datasets [53], which resulted in the shortest Dd in the past two decades. The attribution analysis of Dd anomalies over the NH during 2001-2020 is displayed in Figure 10d (Figure 8b), which led to contrasting D e anomalies over the NH from low to high latitudes (Figure 9b). The contrasting D e anomalies for the period of 2001-2020 agrees with the findings of Chen et al. [24], who observed an advance in D e at northern high latitudes (52-75 • N) but a delay in D e at northern mid-latitudes  • N) at the 90% confidence level. The combination of anomalies in D o and D e led to a shortened duration of snow cover (−2.63 ± 7.26 days) over the NH from 2001 to 2020 (Figures 8c and 9c). In addition to the central United States, Tibetan Plateau, and Scandinavian Mountains, most grid cells within the study presented negative changes in D d over the last two decades. Moreover, the observed D d increase in Central North America is consistent with the increased snow cover days reported by Wang et al. [50] from 2001 to 2015 and Allchin and Déry [51] from 1971 to 2014.

Attribution of Snow Cover Phenology Changes over the NH
The snow accumulation season temperature, T a , precipitation, P a , and snow melting season temperature, T m , were employed to attribute the anomalies in the SCP over the NH during 2001-2020. The spatial distributions of changes in T a , P a , and T m are shown in Figure 10a-c. The attribution analysis of the SCP is presented in Figure 10d-f. The T a (Figure 10a) and T m (Figure 10b) values displayed a general warming trend over the NH from 2001 to 2020. This is reasonable according to the latest WMO Statement on the State of the Global Climate [32], as the past five years (2015-2019) are the five warmest on record. Moreover, both T a and T m exhibited significant latitude differences under the influence of decreasing snow cover and sea-ice extent [25][26][27][28], in which the temperature increased at a higher rate in high latitudes and at a much lower rate in middle latitudes. The anomalies in T a and T m over the NH during 2001-2020 agree with the land surface temperature anomalies mapped from the Goddard Institute for Space Studies (GISS) (https://data.giss.nasa.gov/gistemp/maps/) based on the latest GHCN Monthly version 4 [52]. In addition, changes in atmospheric circulation [29] and human activity [30,31] resulted in a notable increase in precipitation during the snow accumulation season, P a (Figure 10c), in the Central United States, Tibetan Plateau, Northern Europe, and Northeastern Russia, but led to a significant decline in P a in Central Russia, Northeast Canada, and Southern Europe.

Discussion
Accurate estimates of SCP are highly relevant for improving atmospheric reanalysis [54], climate predictions [55], and climate projections. To overcome the shortage of individual datasets in the SCP retrieve, in situ measurements are highly dependent on the location (latitude and elevation) and limited in spatial coverage, visible and near-infrared satellite data are largely influenced by cloud coverage, and it is difficult to distinguish wet and shallow snow in PM snow maps, this study employed the combination of MOD10C2 and IMS prior to the retrieval of SCP matrices. The combination of MOD10C2 and IMS data largely reduced the error in retrieved SCP caused by missing observations in optical SCE products.
Validating satellite-retrieved SCP matrixes using field measurements involves uncertainties because there is a large spatial bias between satellite-retrieved SCP at a spatial resolution of 0.05° and ground measurements at specific locations. Moreover, the different temporal resolution between satellite-retrieved SCP and ground measurements may result in systematic bias in the cross validation process. However, the validation of satellite-  Figure 10e,f. The contributions of each variable to D o and D e variability were computed annually to display the changing influence of climate variables on underlying SCP variability. An analysis of linear correlations of D o with T a and P a over the snow accumulation season revealed that D o was largely determined by T a anomalies (Figure 10e). Meanwhile, a similar analysis revealed that changes in D e were highly correlated with T a , T m , and P a at the 95% confidence level. However, compared with contributions from T a and P a , T m dominated D e anomalies over the NH from 2001 to 2020 (Figure 10f).

Discussion
Accurate estimates of SCP are highly relevant for improving atmospheric reanalysis [54], climate predictions [55], and climate projections. To overcome the shortage of individual datasets in the SCP retrieve, in situ measurements are highly dependent on the location (latitude and elevation) and limited in spatial coverage, visible and near-infrared satellite data are largely influenced by cloud coverage, and it is difficult to distinguish wet and shallow snow in PM snow maps, this study employed the combination of MOD10C2 and IMS prior to the retrieval of SCP matrices. The combination of MOD10C2 and IMS data largely reduced the error in retrieved SCP caused by missing observations in optical SCE products.
Validating satellite-retrieved SCP matrixes using field measurements involves uncertainties because there is a large spatial bias between satellite-retrieved SCP at a spatial resolution of 0.05 • and ground measurements at specific locations. Moreover, the different temporal resolution between satellite-retrieved SCP and ground measurements may result in systematic bias in the cross validation process. However, the validation of satelliteretrieved SCP matrixes using GHCN ground-based snow depth observations indicated that satellite-retrieved SCP matrixes can generally capture the "real" distributions of D o (r = 0.57), D e (r = 0.54), and D d (r = 0.66) at the 95% significance level, which indicate the fitness of using satellite-retrieved SCP in climate change studies over the NH.
Based on satellite-retrieved SCP matrixes, this study reported a delayed D o , associated with an advanced D e over the NH from 2001 to 2020. The above findings are consistent with the published results of Choi et al. [22] and Chen et al. [24], who reported that changes in SCP are largely determined by D e anomalies for the periods of 1967-2008 and 2001-2014. However, the results are in contrast with the findings reported by Allchin and Déry [56], who observed shifting spatial and temporal patterns in the onset of seasonally snow-dominated conditions in the NH from 1972 to 2017. Because the attribution analysis is significantly influenced by temporal coverage, the different time interval may be the possible reason for the inconsistency of the analysis results.
According to the latest climate projections, the reduction of NH SCE will continue in the future [20,57], along with hemispheric land surface warming and a positive feedback on the Earth's climate system through snow-albedo mechanisms [1,18,19,58]. Therefore, the response of SCP to SCE reduction should be investigated in future studies. Moreover, reductions in terrestrial snow cover accounted for one-third of the Arctic albedo decline from 1982 to 2014 [59]. Thus, a coupled study between terrestrial snow cover and Arctic sea ice should be conducted in the future.

Conclusions
With the help of satellite-retrieved snow cover observations, ground-based snow depth records, and reanalysis climate variables, this study explored the SCP distribution and its causes over the NH in the past two decades from 2001 to 2020, which are important data for current cryosphere change studies and future climate projections.
The observed longterm anomalies in SCE over the NH demonstrated that, compared with the monthly averaged SCE in the earlier period of 1966-2000, the later period of 2001-2020 exhibited an increased SCE in the snow accumulation season from September to November, but a decreased SCE in the snow melting season from May to August in the later period of 2001-2020. The enhanced seasonal SCE differences between snow accumulation season and melting season led to corresponding SCP anomalies. The climatology of 20-year averaged D o , D e , and D d over the NH from 2001 to 2020 are well represented in spatial and temporal by satellite-retrieved SCP matrixes. Moreover, the validation results by using 850 GHCN snow depth observations revealed that satellite-retrieved SCP matrixes generally capture the spatial variability of D o , D e , and D d at the 95% significance level during the overlapping period of 2001-2017, which makes it reasonable to employ satellite-retrieved SCP in climate change studies.
Based on satellite-retrieved SCP matrixes, this study reported a delay in D o (1.69 ± 5.93 days) and an advanced D e (−0.94 ± 5.90 days) over the NH for the period of 2016-2020, compared with the period of 2001-2005. The combination of anomalies in D o and D e result in a shortened D d (−2.63 ± 7.26 days) over the NH at the corresponding period. Moreover, Pearson correlation analysis revealed that changes in D d are highly correlated with variabilities in both D o and D e . However, linear analysis displayed unnoticed changes in D o . Therefore, D d changes over the NH from 2001 to 2020 are mostly attributed to anomalies in D e . Furthermore, attribution analysis revealed that the increased melting season temperature is a key driving factor of D e anomalies over the NH from 2001 to 2020. Compared with total precipitation, surface air temperature should be of high concern in SCP analysis.
Although the distribution and attribution of NH SCP has been discussed in previous research, this study provides better results with finer spatial resolution and present the latest SCP changes, which is necessary in the rapid climate change background. Hence, we believe that our study benefits understanding of terrestrial SCP changes and related climate projection studies in response to surface warming in the past decades.   [35], and Muñoz Sabater at 10.24381/cds.68d2bb30 [44]. Moreover, the GHCN-Daily can be found at https://www.ncdc.noaa. gov/ghcnd-data-access.

Acknowledgments:
We acknowledge the three anonymous reviewers for their valuable comments and suggestions. Moreover, we would like to thank the National Snow and Ice Data Center of United States (https://nsidc.org/) for the online available MOD10C2 and IMS datasets, the Rutgers University Global Snow Lab for the online available NH snow cover extent chart product, as well as the Copernicus Climate Change Service for the online available ECMWF ERA5-Land reanalysis 2-m temperature and total precipitation datasets. Furthermore, we would like to thank all people fighting with COVID-19.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Abbreviations used in this study.