Attributing snow cover extent changes over the Northern Hemisphere for the past 65 years

Snow is a crucial component of the cryosphere and its has been experiencing a significant decline for the Northern Hemisphere (NH) (35–90°N) snow cover extent (SCE) in March-April over the 1948–2012 period. However, the causes of this historical snow cover decline are still uncertain. Here, based on the Global Land Data Assimilation System version 2.0 (GLDAS) reanalysis reanalyzed snow cover extent (SCE) and climate model simulations that separate different anthropogenic forcings, we attribute changes of spring SCE over the NH to anthropogenic and natural climate forcings using an optimal fingerprint technique. We find that fingerprints from greenhouse-gases (GHG) and stratospheric aerosols can be clearly detected in the historical SCE records, whereas natural forcing has not contributed to the long-term SCE trend in a discernible way. The GHG-induced warming is primarily responsible for the NH SCE decline, which, however, has been partly offset (by ∼16%) by aerosol-induced climatic cooling. Our findings confirm the negative effect of warming on SCE changes yet highlight the positive role of aerosols in slowing SCE decline over the past 65 years.


Introduction
Snow, as a vital component of the cryosphere, covers an area of about 46 million square kilometers, 98% of which is located in the Northern Hemisphere (NH) (Armstrong and Brodzik 2001). Due to its unique radiative (high reflectance) and thermal (low thermal conductivity) properties, a small variation of the snow cover can substantially affect the radiated energy budget of the land-atmosphere system and regional horizontal thermal difference (Groisman et al 1994a, Groisman et al 1994b, Déry and Brown 2007, Brown and Mote 2009, Kang et al 2010, Liang et al 2010, with consequent impacts on the local, regional and global climate systems (Fasullo 2004, Seol and Hong 2009, Wu et al 2009, Cohen et al 2012, Furtado et al 2015, Yang et al 2019, You et al 2020. Over the last few decades, the global snow cover extent (SCE) has been experiencing a significant decline trend, especially since the 1980s and during the springtime (Robinson and Frei 2000, Déry and Brown 2007, Brown et al 2010, Brown and Robinson 2011, Derksen and Brown 2012, Brutel-Vuilmet et al 2013, Rupp et al 2013. This decline in SCE has been related to the local mean climate, e. g. the IPCC Fifth Assessment Report (Stocker et al 2013) noted a significant relationship between the global SCE decline and warming. However, the causes of the historical SCE decline still remain unclear, raising uncertainties in predicting future SCE changes and relevant consequences.
Anthropogenic climate warming has been widely considered as the dominant factor that drives the recent SCE decline (Derksen and Brown 2012, Najafi et al 2016, Willibald et al 2020. The radiative response to elevated greenhouse gases (GHG) concentration causes an increased air temperature in the lower atmosphere that expedites the melt of snow cover. Besides that, natural variability (NAT) of the Earth system itself, such as volcanic eruption and solar variability, also plays an essential role in driving global climate change (Mantua et al 1997, Schlesinger andRamankutty 1994) and affects SCE. However, the impact of natural variability of the Earth system on SCE changes is not well documented in existing studies. In addition, it is important to note that both human activities and volcanic eruptions are associated with the injection of SO 2 gas into the stratosphere that favors the creation of sulfate aerosols, which weakens solar radiation reaching the Earth's surface and may lead to temperature drop (Rotstayn et al 2013, Najafi et al 2015. Studies have found that other anthropogenic forcings (OANT), mainly stratospheric aerosols, could offset about 60% of the GHG-induced climate warming in the Arctic region over the past century, which substantially slows down the melt of the Arctic sea ice (Najafi et al 2015, Mueller et al 2018. Therefore, it is reasonable to hypothesize that the GHG-induced NH SCE decline is also partly masked by the aerosol effect. Nevertheless, we are not aware of a formal detection and attribution analysis that considers the combined effects of GHG, NAT and OANT on the change of NH SCE up to now. To fill this knowledge gap, in this study, we use a long-term SCE dataset (covering a period of 1948-2012) from the GLDAS (Global Land Data Assimilation System version 2.0) reanalysis and the state-of-the-art fully coupled climate models that participated in the fifth phase of the Coupled Model Intercomparison Project (CMIP5) to explicitly attribute changes in SCE to GHG, NAT and OANT forcings for the past 65 years. Our specific objectives were (i) to examine the temporal pattern of SCE over the NH during 1948NH during -2012 to investigate the forced response of the NH SCE to anthropogenic and external natural forcings, and (iii) to explore whether and to what extent the GHG effect on SCE changes is masked by the OANT effect.

Historical snow cover datasets
The analysis of the SCE trend and its attribution are based on the reanalyzed SCE from the Global Land Data Assimilation System Version 2.0 (GLDAS_CLSM025-D) over 1948-2014 . The GLDAS SCE (SCE GLDAS ) has a spatial resolution of 0.25°and is estimated based on daily snow water equivalent (SWE) directly simulated by GLDAS and updated by assimilating the MODIS (Moderate-resolution Imaging Spectroradiometer) snow cover records from NASA'S Terra satellite after 2000 (Rodell and Houser 2004). Gridcells with SWE higher than 5 kg m −2 (∼5 mm) are considered to be covered by snow; otherwise are considered as snow-free pixels (Brutel-Vuilmet et al 2013). The monthly SCE GLDAS is derived by averaging daily SCE values within a month.
To evaluate the performance of SCE GLDAS in capturing the temporal SCE trend, two additional SCE datasets are used in this study. The first dataset is the EASE-Grid 2.0 Weekly Snow Cover and Sea Ice Extent, which is available at a 25 km spatial resolution and weekly temporal resolution from 1968 through 2012 and is produced by the National Ocean and Atmospheric Administration (NOAA) and National Climatic Data Center Climate Data Record (referred to as SCE NOAA hereafter) (Robinson et al 1993). The second dataset is the Reconstructed North American, Eurasian and Northern Hemisphere Snow Cover Extent from 1915 through 1997 produced by Brown (2002) (referred to SCE Brown hereafter). The SCE Brown is reconstructed by combining SCE NOAA with in situ observations of snow depth and climate variables from the U.S., Canada, China and the Former Soviet Union (FSU) using statistical methods (Brown 1997, Frei et al 1999. Here, we use the spring (March-April) SCE for the comparison of the three SCE datasets as SCE Brown only contains SCE for March and April.

CMIP5 climate model simulations
Climate models participated in CMIP5 were used to examine the forced response of the climate system to anthropogenic and external natural forcings (Taylor et al 2012). Three independent modeling experiments are considered, including (i) the historical simulation forced by natural (solar variations and volcanic eruptions) and anthropogenic (well-mixed greenhouse gases, aerosols, and land-use change) changes (ALLF), (ii) the GHG experiment that only considers the effect of well-mixed greenhouse gases, and (iii) the NAT experiment that merely accounts for the effect of natural forcing changes (volcanic aerosols and solar variability). As a result, the impacts of other factors (OANT) including aerosols, ozone and land use changes on SCE changes were determined as the residuals of ALLF experiment minus the sum of GHG and NAT experiments. Only seven CMIP5 models conducted all three experiments (table S2) so that we only used these seven models to avoid the uncertainties caused by varying model structures and/or parameterization schemes. In addition, the unforced control simulations of pre-industrial control simulations (CTL) for 15860 years from 33 CMIP5 models (244 non-overlapping 65-years segments) (table S2) were also used to determine the internal climate variability. Finally, all CMIP5 outputs were resampled to a 0.25°spatial resolution using the first-order conservative method (Jones 1999, https://code.zmaw.de/projects/cdo).
2.3. Attribution of SCE changes using optimal fingerprinting analysis The total least squares (TLS)-based optimal fingerprinting method (Allen and Tett 1999, Allen and Stott 2003) is used to detect the relative contribution of each external forcing factor (GHG, OANT and NAT) to SCE GLDAS changes. This approach expresses the time series of SCE GLDAS (y) as a linear combination of scaled (β i ) responses to external driving factors (x i ) and internal variability (ε) (equation (1)): Specifically, the SCE GLDAS anomalies are first regressed onto the ALLF, GHG and NAT simulations to obtain the scaling factors for ALLF, GHG and NAT simulations (equation (2)), Then the SCE anomaly from the ALLF experiments are decomposed into the SCE anomalies caused by GHG, NAT and OANT forcings (equation (3)) by considering linear additivity (Najafi et al 2015), where SCE GLDAS is a spatial-temporal vector (135-yr periods×3 sub-areas) of SCE anomalies (relative to1960-1990 average) obtained by 13 consecutive 5-year periods for 1948-2012 and 3 distinguished sub-areas to summarize the SCE spatial variability over the NH (Eurasia, Asian, North America over the NH, table S1 (available online at stacks.iop.org/ERC/3/061001/mmedia)). The fingerprint vectors contained in the external and internal forcings are constructed using the same technique. A positive scaling factor that is significantly different from zero at the 90% confidence interval (CI) indicates that the signal of the corresponding response can be detected in the SCE record, with the scaling factor closer to 1 indicating a stronger control of the corresponding factor on SCE changes. Finally, a residual consistency test (RCT) was conducted to check if the regression residuals of SCE and bestestimate combinations of signals are consistent with the internal SCE variability by two covariance matrices of internal variability (Allen and Tett 1999). The two covariance matrix estimates were divided into two groups from the full set of chunks of the CTL (one is for optimization, another is for test). The F-test was involved in the comparison of the variances of the regression residuals with the internal variability from CTL. The p-value represents whether the null hypothesis is rejected or not at a 10% significance level, so that if the p-value falls between 0.05 and 0.95 then the test is passed (Allen andStott 2003, Rupp et al 2013). If the estimated internal variability from models is different from that in observations or if model simulations had incorrect response patterns, then the test is failed. A detailed description of the residual consistency test can be found in Allen and Stott (2003) and Ribes et al (2013).

Reliability of the SCE GLDAS dataset
To check the reliability of SCE GLDAS in capturing changes of the NH SCE, we compared SCE GLDAS with the other two SCE datasets (i.e., SCE NOAA and SCE Brown ) during their respective overlapping periods. Averaged over the NH, significant positive correlations are found between SCE GLDAS and SCE NOAA (r=0.72; p<0.01) during 1968-2012 and between SCE GLDAS and SCE Brown (r=0.82; p<0.01) during 1948-1996 (figure S1). Additionally, we also compared the spatial pattern of the SCE trend between SCE GLDAS and SCE NOAA for 1968-2012 (note that SCE Brown does not provide SCE spatial patterns) and find a generally consistent pattern between the two datasets (figure S2). The above results confirm the overall reliability of SCE GLDAS in capturing the long-term SCE changes in the NH.

Changes in SCE from SCE GLDAS and CMIP5 model simulations
The spring season (March-April) SCE exhibited a significant decreasing trend over the past 65 years across the majority of the snow-covered areas in the NH, which is captured by both SCE GLDAS (in 80% of snow-covered areas) and the multi-model ensemble-mean SCE under ALLF in CMIP5 simulations (in 75% of snow-covered areas) (figures 1(a) and (b)). In both datasets, the largest decreases are found in western North America, Europe and central Eurasia, whereas a slight SCE increase is observed across a narrow latitudinal band beyond 60°N in northeastern Asia. Averaged over the NH, SCE GLDAS showed a more rapid decrease (−0.5×10 6 km 2 per decade, p<0.01, figure S3) than the ALLF simulations (−0.2×10 6 km 2 per decade, p<0.01, figure S3), which was primarily caused by an evident larger decreasing trend of SCE GLDAS before the 1960s. For more recent decades, the two datasets showed a very close SCE decreasing trend (−0.5×10 6 Km 2 per decade, p<0.01 for SCE GLDAS versus −0.4×10 6 Km 2 per decade, p<0.01 for ALLF simulation; figure S4). Moreover, it is noted that the temporary SCE increases in ALLF forcing correspond well with increased aerosol (stratospheric sulfate, SO 4 ) burden after the three most recent major volcanic eruptions, e.g. Mt. Agung (1963), El Chichón (1982) and Mt. Pinatubo (1991) (figures 2 and S3), indicating that the CMIP5 models have reasonably captured the temporary weakening of the SCE decline due to aerosol cooling.   [1948][1949][1950][1951][1952][1953]. The red, green and blue curves represent simulated SCE under ALLF, GHG and NAT forcings, respectively. The orange curve represents the OANT forcing, which is calculated as ALLF minus the sum of GHG and NAT. The solid lines are the ensemble mean of the CMIP5 models and the shaded area indicates the 5%-95% confidence interval for ALLF simulation. The black curve is the SCE anomaly from GLDAS, which is considered as observation in this study. The vertical grey dashed lines indicate the timing of the three big volcanic eruptions (Agung eruption, El Chichón eruption and Pinatubo eruption).

Detection and attribution of SCE changes
Three factors are considered in this study as potential causes of SCE changes, including (i) GHG: humaninduced well-mixed GHG, (ii) NAT: natural variability, and (iii) OANT: other factors but mainly the aerosol effect. According to the CMIP5 model simulations, as expected, GHG-induced warming has led to a widespread spring SCE decline in the NH over the 65 years ( figure 1(c)). In the ensemble-mean of the GHG-only experiments, about 83% of snow-covered regions show a reduction in spring SCE, with an average decreasing trend of −0.3×10 6 Km 2 per decade (p<0.01) (figure 2).
By contrast, spring SCE in the OANT-only simulation exhibits an upward trend across the majority of the NH (in about 66% of the snow-covered areas; figure 1(d)) for the same period, with an average increasing trend of 0.1×10 6 km 2 per decade (p<0.01) (figure 2). This suggests that the aerosol effect has partly offset the GHGwarming effect on spring SCE decline over the past six decades. Compared with GHG and OANT, NAT has not exerted an evident spring SCE trend across most of the NH. In the NAT-only simulations, increases and decreases in spring SCE are found in 48% and 52% of the snow-covered areas, respectively, leading to a nonsignificant overall NAT-induced SCE trend for the past 65 years (−0.04×10 6 Km 2 per decade; p=0.06) (figure 2). Similar results are obtained with a different spring definition (i.e., March-May) (figures S5 and S6).
We then perform the optimal fingerprinting method to detect the relative contribution of each external forcing to SCE changes. The 90% confidence intervals (CI) of the resulting scaling factors, which scale the simulated responses to best reproduce the spring SCE changes, are shown in figure 3 (represented by black circles). It is found that scaling factors of GHG (β=1.53) and OANT (β=0.60) are both higher than zero and their corresponding 90% CIs encompass one, suggesting that the effects of GHG and OANT on spring SCE changes can be clearly detected by the optimal fingerprinting method. By contrast, the scaling factor of NAT is close to zero and its 90% CI encompasses zero, indicating that the effect of NAT alone is generally undetectable or inconsistent with the observed SCE changes. Similar results are obtained with a different definition of spring (March-May) (figure S7). Moreover, we estimated that the variance of the regression residuals (in equation (2)) is not significantly different from the internal variability from CTL, as indicated by a p-value of 0.21 using an Ftest. This suggests there is no inconsistency between the regression residual and the internal variability in the climate model simulations Trends in SCE induced by each influence factor are obtained by multiplying the trends in the multi-model forced responses with the estimated scaling factors (figure 3, represented by red circles). By doing so, we estimated the GHG-induced warming effect alone has resulted in a spring SCE trend of −0.51 (ranging between −0.23 and −0.79) × 10 6 km 2 per decade over the last 65 years across the entire snow-covered NH. On the contrary, the OANT impact has led to a cooling effect that increases the spring SCE at a rate of 0.08 (ranging between 0.02 and 0.15)×10 6 km 2 per decade. This means that the OANT effect has offset approximately 16% of the GHG-induced SCE reduction. At the regional scale, this OANT-induced offsetting effect is mostly manifested in northern Europe (56%), the Tibetan Plateau (75%) and eastern Asia (83%), and also clearly detected in central North America (27%), central Europe (40%) and central Asia (43%) (figure 4). These areas are mainly distributed in the mid-latitude zone, where human activities (except for the Tibetan Plateau) are relatively intensive. However, in higher latitudes where the population density is relatively small (e. g., Alaska, Canada and northern Asia), the offsetting effect is typically lower than 5% (figure 4). As for NAT, its induced SCE trend is generally very small and not discernible at both hemisphere and regional scales (figures 2-4). Attributing SCE trend to GHG, OANT, and NAT forcings during 1948-2012 over the NH. Black dots are the best estimates of scaling factors (black) and the red dots depict the corresponding attributable trend. Error bars quantify the 10th (lower end) and 90th (upper end) confidence intervals (CI). Positive values of scaling factors with the 10th CI above zero indicates that the signal is detectable, and the scaling factor closer to 1 with smaller uncertainty ranges indicate a good agreement between CMIP5 model estimations and SCE GLDAS .

Discussion and concluding remarks
GHG-induced anthropogenic climate warming is widely considered to be responsible for the cryosphere changes including the widespread degradation of permafrost and mass losses of glaciers over the past few decades (Hirabayashi et al 2016, Guo et al 2020. Here, by analyzing observed and modeled long-term SCE using a fingerprint technique, we show that the anthropogenic warming is also the key driver of the SCE decline over the NH for the past 65 years (figures 1(c) and 2). This result is consistent with previous findings using temperature-index snowpack model sensitivity analysis on the cause of NH SCE changes (Brown et al 2003, Brown and Mote 2009, Mudryk et al 2017 and again, highlights the predominant role of temperature rise on the control the long-term SCE regimes.
However, previous studies on the attribution of SCE changes ignored an important fact that aerosol particles are concurrently emitted into the atmosphere along with GHG by human activities. Aerosols are small particles (typically with a radius<20 μm) suspended in the atmosphere, lasting for days to months, which weaken the solar radiation reaching the Earth's surface that may induce a cooling effect in the lower atmosphere (Charlson et al 1992). Additionally, aerosol particles serve as condensation nuclei, which favor the formation of precipitating clouds and consequently result in increased precipitation. Both cooling and precipitation (when fall as snow) enhancing effects of aerosols favor the growth of snow cover, which may partly offset the GHGinduced SCE decline. By using an optimal fingerprint technique, our results clearly show that the OANT forcing (mainly aerosols) has offset approximately 16% of the GHG-induced SCE decline in the NH over the past 65 years. This result aligns with the findings of Arctic sea ice changes by Mueller et al (2018), who found that OANT could offset about 23% of Arctic sea ice decline caused by GHG-induced climate warming. At the regional scale, higher offsetting effects of the OANT forcing on SCE decline are generally found in areas with more intensive human activities, such as eastern Asia, North America and Europe. In eastern Asia and North America, the high aerosol concentration is mainly constituted of sulfate aerosols that originate from fossil fuels and/or biomass burning (Huang et al 2015, Xie et al 2018, whereas in Europe, the aerosols are mainly secondary organic aerosols indirectly produced by mineral dust and vegetation emissions, and secondary inorganic aerosols indirectly produced by human emissions (e.g. ammonium salts, nitrates and sulfates) (Daellenbach et al 2020). Different from these three human activity intensive regions, the high offsetting effect of the OANT forcing detected over the Qinghai-Tibet Plateau region is mainly caused by aerosols and other atmospheric pollutants discharged from South and Central Asia that are then transported to the Qinghai-Tibet Plateau via atmospheric circulation and consequently affect its climate environment (Li et al 2016. In comparison to GHG and OANT forcings, the NAT forcing (or natural variability of the climate system) does not contribute to the long-term SCE trend in a detectable manner. This finding holds at the hemispheric scale as well as for the analyzed 13 sub-regions (figures 1(e) and 4). This is somewhat expected, as the variability Figure 4. Comparison of the magnitude of spring SCE trend under ALLF, GHG, NAT and OANT forcings over 13 regions across the NH. The regional boundaries were derived from the IPCC Fifth Assessment Report (AR5). Boxes indicate the 10th and 90th percentiles and bars represent the minimum and maximum value of all seven models. The median value is shown by the horizontal bar inside each box. of solar forcing and/or large-scale atmospheric circulations (e.g., ENSO) have a relatively short return period (e.g., ENSO has a return period of 2 to 7 years) and therefore may not exert any long-term impact. Nevertheless, both SCE GLDAS and CMIP-modelled SCE clearly demonstrate the impact of natural variability of the climate system on SCE change over a relatively short period. For example, the responses of SCE to the three big volcano eruption events during the study period are well captured by the CMIP models (the OANT line in figure 2) and reflected in the overall signal of SCE changes (ALLF and OBS lines in figure 2).
Finally, it is worthwhile mentioning that the current attribution analysis also has a few uncertainties. First, we use the GLDAS reanalysis as an observational SCE data since long-term remote sensing-based SCE is still lacking. Although comprehensively validated, SCE GLDAS still contains certain uncertainties (Broxton et al 2016). This uncertainty may be potentially overcome with the continuous of relevant satellite missions and the development of retrieval algorithms. Second, because the current CMIP models do not perform an aerosol-only simulation, the aerosol effect on SCE changes is determined by subtracting the warming and natural forcing effects from the ALLF simulation. Besides aerosols, other environmental changes (e.g., ozone and land cover changes) may also contribute to the OANT signal. However, previous studies have demonstrated that the OANT signal is dominated by aerosol changes, as the cooling effect of aerosols is much larger than the small warming effect induced by ozone changes and the effect of land cover changes is even smaller (Fyfe et al 2013, Najafi et al 2015. Nevertheless, this issue calls for more separate simulation scenarios in upcoming CMIP missions such that to more robustly attribute SCE changes in the future.