Historical evolution and future trend of Northern Hemisphere snow cover in CMIP5 and CMIP6 models

Historical snow cover over the Northern Hemisphere was examined in the satellite-based NOAA-CDR data for the period of 1970–2019. Observed annual snow cover fraction (SNF) has reduced over most areas by up to 2%/decade, while annual snow cover area (SCA) has reduced by 2 × 105 km2/decade. However, SCA in the October–December season has increased by about 5 × 105 km2/decade. CMIP5 and CMIP6 historical experiments were validated against the NOAA-CDR data. Snow cover was generally well simulated in both CMIPs, with CMIP6 models performing better. The biases in SCA reduction were larger and smaller during summer and winter, respectively. The observed increase of October–November–December SCA in the 2000s was not reproduced. Climate projections of future snow cover were evaluated in CMIP6. SNF is projected to decrease in the next 80 years, under all four scenarios evaluated (SSP126, SSP245, SSP370 and SSP585). The higher the greenhouse emissions in the shared socio-economic pathways, the faster the reduction. Under the SSP585 scenario, the rate of SCA reduction is projected to exceed −1.2 × 106 km2/decade. By 2081–2100, annual (January–March) SCA is projected to decrease by more than 30% (20%). Under the SSP126 scenario, annual (January–March) SCA is projected to only reduce by about 10% (5%) relative 1995–2014 values. The reduction of Greenhouse gas emissions is critical to controlling the loss of snow cover; future snow cover only stabilizes under the SSP126 scenario, but continue to decrease under the other three scenarios.


Introduction
Snow cover plays an important role in the earth's climate system and is the largest component of the terrestrial cryosphere (Lemke et al 2007), with a maximum extent of about 45 × 10 6 km 2 . Because of its high surface albedo and heat-insulation effect, Snow cover has important effects on climate for the energy exchange between the land surface and atmosphere (Barnett et al 1988, Yang et al 2001, Chapin et al 2005, Euskirchen et al 2007, Hancock et al 2013, Beniston et al 2018, Kang et al 2019. Snow cover can serve as an indicator of atmospheric circulation and climate variations because it provides part of the lower boundary condition of the atmosphere, forcing the regional and global atmosphere (Barnett et al 1989, Bamzai and Shukla 1999, Fasullo 2004, Wu and Kirtman 2007, Wu et al 2009, Allen and Zender 2011, Henderson et al 2018, Luo and Wang 2019. Snow cover over land has a cooling effect on the climate system (Flanner et al 2011), but with the reduction of snow cover over recent decades, the magnitude of this cooling influence has declined (Letterly et al 2018, Zhang et al 2019, Pulliainen et al 2020. Snow cover also plays an important role in hydrological cycle as the source of water for streamflow and groundwater recharge, and freshwater from snow melt is a volatile natural resource (Stewart et al 2004, Rauscher et al 2008, Koster et al 2010, Mahanama et al 2012, Barnhart et al 2016, Niittynen et al 2020 Due to the significant impacts of snow cover on the climate system, it is important to investigate the climatological characteristics and perform future climate projections of snow cover variations. Observational evidence shows that Northern Hemisphere (NH) snow cover extent (SCE) has declined over the past 90 years (Brown et al 2017, with the greatest reduction (∼3.3 million km 2 ) observed during the 1980s (Vaughan et al 2013). Snow ablation and accumulation are affected by many factors, such as the air temperature and surface radiation over the land surface. Most of the loss in snow occurs over 'temperature sensitive regions' , where changes in SCE are closely linked to temperature variability (Karl et al 1993, Déry andBrown 2007). In addition to temperature and precipitation, aerosol and vegetation cover are also the key drivers of snow cover variability (Bond et al 2013, Pearson et al 2013. Coupled climate models are the primary tools that are used to study snow cover evolution and project snow cover variability in future. The Coupled Model Inter-comparison Project (CMIP) through its six iterations (CMIP1 to CMIP6) provides a very useful framework for the modeling community to assess the performance of climate change simulations in state-of-the-art climate models in order to improve model fidelity. The results from CMIP5 models reproduced the reduction of snow cover extent over the recent decades, but tend to underestimate the current rate of decline compared to observations (Frei et al 2003, Derksen and Brown 2012, Zhu and Dong 2013. This underestimation, which is also found in the CMIP3 simulations (Roesch 2006), may be associated with underestimated extratropical spring warming (Brutel-Vuilmet et al 2013), and reduced sensitivity of snow to temperature . Derksen and Brown (2012) and Mudryk et al (2014) illustrate other aspects of the observed trends that are not well captured by general circulation models (GCMs).
The most recent state-of-the-art climate model experiments are now becoming available as part of the CMIP6 ensemble (Eyring et al 2016). Compared to CMIP5, the number of different models, model variants, and ensemble sizes of individual models have all increased in CMIP6. Future scenario simulations in CMIP6 use emission scenarios coordinated by the ScenarioMIP project (O'Neill et al 2016). The new set of emissions and land use scenarios are, known as shared socioeconomic pathways (SSPs) (Riahi et al 2017). However, there is a lack of studies on the performance of CMIP6 models in the simulation of NH snow cover, and its future climate projections.
To fill in this knowledge gap, we investigate three critical questions regarding climate model simulations of snow cover extent using observation data and CMIP5/CMIP6 simulations: Firstly, what are the characteristics of snow over variability over different areas in NH in recent decades? Secondly, is there a significant improvement for the simulation of snow cover extent by CMIP6 models compared with CMIP5 models? Finally, what is the evolution trend of NH snow cover under different SSPs projected by CMIP6 models in the future? In the paper, the time lag between temperature and snow cover was considered when defining the seasons for the snow cover analysis (January-February-March = JFM, April-May-June = AMJ, July-August-September = JAS, and October-November-December = OND).
The paper is organized as follows. Section 2 briefly describes the snow data sources (CMIP5 and CMIP6 models, observations) and the methodology of comparison. In section 3, an evaluation of the long-term annual and seasonal climatology of snow characteristics is presented; the simulations of snow cover in CMIP5 and CMIP6 are evaluated; projection of future trend shown by NH snow cover is also conducted in this section. Section 4 summarizes the results and suggests some possible lines of future research.

Data and method
In the paper, climate data processing and figure plotting were made using the NCAR Command Language.

Observations
For snow cover, we used the satellite-based weekly SCE dataset of Robinson and Estilow (2012) for the period of 1967-2019. This dataset is based on NOAA's National Climatic Data Center (NCDC) Climate Data Record (CDR) and available for download at http://climate.rutgers.edu/snowcover. It is the longest available satellite record of NH SCE and has been used in numerous studies Brown and Robinson 2011, Derksen and Brown 2012, Brown and Derksen 2013, Déry and Brown 2007. The SCE data is a merged product from the reanalysis of observational snow cover maps and remote sensing measurements. From October 1966 to June 1999 SCE was estimated by visually assessing satellite shortwave imagery. The manual approach was replaced by the National Ice Center daily Interactive Multisensor Snow and Ice Mapping System from June 1999 (Helfrich et al 2007), which relies on the visible satellite imagery, and utilizes station observations and passive microwave data to improve snow monitoring during cloudy or nighttime conditions. Weekly SCE maps were generated on an 89 × 89 grid overlaid on a NH stereographic map with spatial resolution ranging from 16 000 km 2 to 42 000 km 2 (Robinson et al 1993). Each grid cell was flagged as snow covered if more than 50% of its area

Method
The observed and modeled SCE was examined over the NH in terms of climatology, linear trends, spatial pattern of trends, and time change of trends. SCE was described in terms of snow cover fraction (SNF) and snow cover area (SCA). Analysis was performed for the whole year, JFM, AMJ, JAS and OND, and individual months. The 3 month seasons (JFM, AMJ, JAS and OND) were chosen considering the time lag between temperature and snow cover. The 1995-2014 period was used as the baseline against which the CMIP6 SSPs were compared. In order to evaluate the model performance, the Taylor diagram technique is used. The Taylor diagram represents the correlation (R), the centered rootmean-square-error (RMSE) and the amplitude of the standard deviations (STD). The diagram provides a way of graphically summarizing how closely a pattern matches observations (Taylor 2001). The mean, linear trend, linear correlation, spatial correlation, and other common statistical methods were used to evaluate the differences between historical simulations and observations, as well as between different models. The linear trends are calculated using least square method  (1): According to equation (1), RMSE is the root mean square error of the grid, X sim,i is the simulated value of the grid i, and X obs,i is the corresponding reanalysis data.
All data was re-gridded through bilinear interpolation onto 1.0 • × 1.0 • grids to facilitate comparison between different models or observation. Multi-model ensemble means (MMEs) refer to simple arithmetic means. All model values and fractions presented in this study are MMEs. Figure 1 shows the observed SNF over the NH during 1970-2019. Means and trends are shown for the whole year, JFM, AMJ, JAS and OND. From figures 1(a)-(d), climatological snow cover patterns had a strong dependence on latitude and seasons. Greenland was fully snow covered in all seasons. In JFM, this was also the case for northern Eurasia (EA) and northern North America (NA). Elsewhere and in the other seasons, the SCF was smaller. For JAS, the SNF was less than 50% over most area in NH except Greenland. Climatological snow cover patterns were similar between AMJ and OND with 30%-80% snow covered land in NH except Greenland.

Historical evolution
The trends of snow cover during 1970-2019 were seasonally variable (figures 1(f)-(j)). In terms of annual mean SNF, reductions occurred over the northern regions of NA, northern regions of EA, and south-western regions of EA, while increases occurred over the central-southern regions of NA and northern regions of China. In terms of JFM means, clear reductions occurred over south-western regions of EA. Over northern China, snow cover increased. There was significant decreasing trend of SNF over most parts of the NH except for the Tibetan Plateau where snow cover increased for AMJ. For JAS, there was no obvious trend of snow cover over most parts of the NH which likely reflects the virtual lack of snow in this season. For OND, increases occurred over most of the NH, with reductions occurring only over the south-western regions of EA. Figure 2 shows the time series of NH SCA during 1970-2019, as anomalies from the mean. Annual mean SCA has clearly reduced over the NH, at rates of up to −0.20 × 10 6 km 2 /decade. EA and NA contributed reductions of −0.14 × 10 6 km 2 /decade and −0.06 × 10 6 km 2 /decade, respectively. JFM mean SCA has also reduced for NH and EA with the trends of −0.09 × 10 6 km 2 /decade and −0.12 × 10 6 km 2 /decade respectively, but there was an insignificant increasing trend for NA. Significant decreasing SCA trends were found for AMJ and JAS; −0.90 × 10 6 km 2 /decade and −0.35 × 10 6 km 2 /decade, −0.60 × 10 6 km 2 /decade and −0.20 × 10 6 km 2 /decade, and −0.30 × 10 6 km 2 / decade and −0.15 × 10 6 km 2 /decade for the NH, EA, and NA, respectively. In contrast, OND mean SCA has increased over NH, with a rate of 0.54 × 10 6 km 2 /decade. The corresponding trends for EA and NA were 0.35 × 10 6 km 2 /decade and 0.19 × 10 6 km 2 /decade, respectively. The positive trend was most obvious after 1990.
From a climatological perspective (figure 3), SCA over the NH was highest in January at 47 × 10 6 km 2 , and lowest in August at 3 × 10 6 km 2 . SCA over both EA and NA was also highest in January and lowest in August. This illustrates how SCA over the NH was strongly controlled by the seasonal temperature cycle, particularly so over EA. Figure 3(b) shows the trends of SCA during 1970-2019 by month, and supports the above analysis of primarily reduction during spring-summer and some increment during autumnwinter. SCA reduction was highest in June, being −1.3 × 10 6 km 2 /decade, −0.7 × 10 6 km 2 /decade, and −0.6 × 10 6 km 2 /decade over the NH, EA, and NA, respectively. SCA increment was highest in November, being 0.7 × 10 6 km 2 /decade, 0.4 × 10 6 km 2 /decade, and 0.3 × 10 6 km 2 /decade over the NH, EA, and NA, respectively. In summary, SCA over the NH was low during summer and high during winter, while the trend was reduction during spring-summer (significant from March to August) and increment during autumn-winter (significant from November to December).

Model validation
The CMIP5 and CMIP6 ensembles were compared against observations to evaluate their ability to simulate snow cover in the NH. Spatial biases were  Figures 4(a)-(e) shows that during all five periods, SNF in the CMIP5 MME was biased low over most of the NH, except for the Tibetan Plateau and parts of north China. The AMJ and JAS biases were smaller than the annual bias, while the JFM and OND biases were larger than the annual bias. The largest bias occurred over the southwestern region of EA and over the Tibetan Plateau, with an absolute magnitude as high as 20% for JFM and OND. Figures 4(f)-(j) show that SNF in the CMIP6 MME was biased in a pattern similar to that in the CMIP5 MME. However, biases were clearly smaller in the CMIP6 MME, particularly over the northern region of EA. While both CMIP5 and CMIP6 simulated SNF better in the absolute sense during AMJ and JAS than during JFM and OND, this reflected the lower spring-summer SNF itself.
The improvement of CMIP6 over CMIP5 was examined for individual models using Taylor diagrams (figure 5). From figure 5(a), the spatial correlation between simulated and observed annual SNF was greater than 0.8 for all models. Of the few models that achieved spatial correlations greater than 0.9, more were from CMIP6. Similar results were obtained for JFM, AMJ and OND SNF. For JAS SNF, more models performed poorly, especially for CMIP5. Model performance was best for JFM SNF, and most of the better models were from CMIP6. Figure 6 shows the time series anomalies of NH SCA during 1970-2014, comparing simulated against observed. Observed annual, JFM, AMJ and JAS SCA has decreased over the NH, as well as separately over both EA and NA. SCA reached its lowest value around 1990 before recovering before 2000. Most CMIP5 and CMIP6 models reproduced the decreasing trends for annual and spring-summer SCA, but as expected did not reflect observed variability like the 1990 low. In contrast to the JFM, AMJ and JAS SCA, the observed OND SCA did not present any clear trend until around 2000, when it clearly began to increase. This increase was seen over both EA and NA, i.e. the whole NH. Neither the CMIP5 nor the CMIP6 MME reproduced the OND SCA trends, only showing decreasing trends throughout the 1970-2014 period. In summary, both CMIP5 and CMIP6 models reproduced the trends of decreasing SCA over the NH, but biases from observation worsened in the last few decades. This may be due to models simulating excessive global warming trends. , SCA climatology over the NH was well simulated in both CMIP5 and CMIP6, being highest in April and lowest in August. The bias from observation was smaller during summer in both CMIP5 and CMIP6. CMIP5 showed a large winter bias, with the MME being up to 8 × 10 6 km 2 less than observation for the whole NH in January. CMIP6 performed better than CMIP5, showing a relatively minor winter bias. From figures 7(d)-(f), observed SCA over the NH has significantly decreased in spring-autumn and increased in autumn-winter. Both CMIP5 and CMIP6 ensembles reproduced the reduction in spring-summer, although weaker than observed reduction in May-June-July. Neither CMIP5 nor CMIP6 ensembles reproduced the increment in autumn-winter. The differences between CMIP5 and CMIP6 trends were relatively minor. Although there is substantial variability in the monthly trends between the individual CMIP5 and CMIP6 models (see the error bars representing the intermodel STD), the observed increase in NH SCA in October-December was well beyond the variation between the model simulations.
In conclusion, SCA climatology and annual trend over the NH were generally well-simulated by both CMIP5 and CMIP6 ensembles. However, the observed increment of OND SCA was not reproduced, especially the trend in the 2000s, which may result from the overestimation of the warming trend during this period by CMIP models (Zhu et al 2017). Regarding the discrepancy between the observed OND increase and model-simulated decrease in SCA, a previous study  identified the NOAA-CDR as the only data set with increasing snow cover trends in OND when compared with six other observational SCE data sets in the years 1981-2010. So, the NOAA NOAA-CDR climate data record trends should be interpreted with considerable caution in these months. Analysis of multiple datasets (Mudryk et al 2015) or an ensemble product (e.g. CanESN, CESM, Mudryk et al 2017) helps to overcome the limitations of individual datasets and provides a better assessment of the climate state. CMIP6 models showed a marked improvement over CMIP5 in reproducing the winter climatology of SCA over the NH. Hence, we expect CMIP6 climate projections of SCA over the NH to be more reliable. The results are summarized in table 3.

Future projections
Twenty CMIP6 models under four SSPs were evaluated for the next 80 years (2021-2100). Figure 8 shows a comparison between SNF projected in 2081-2100 against the baseline of 1995-2014, revealing projected widespread reductions in SNF over the whole NH. Under all four scenarios, SNF reductions are projected over the western regions of EA, the Tibetan Plateau, and NA. Reductions are larger during AMJ and OND than JFM and JAS. Under the SSP126 scenario, SNF reduces by 1%-8%, mostly over the western regions of EA and the western regions of NA. Under the SSP245 scenario, SNF reductions are larger, with JFM reductions larger than 10% over the western regions of EA and the western regions of NA. Under the SSP370 scenario, SNF reductions are larger still, with JFM and OND reductions larger than 15% over the western regions of EA, the Tibetan Plateau, and the south-western regions of NA. Under the SSP585 high emissions scenario, SNF reductions are largest, with JFM and OND reductions larger than 20% over the western and northern regions of EA and central-south regions of NA.
From figure 9, SCA under all four SSPs decreased over the whole NH, as well as EA and NA independently. The four SCA trajectories are effectively indistinguishable until 2040, after which they diverge strongly. SCA is stabilized under SSP126, but continues to decrease under the other three scenarios; the higher the emissions, the larger the reduction. Like SNF, SCA reductions are larger during AMJ and OND than the other seasons.  table 4. Annual mean SCA trends over the NH under the SSP126, SSP245, SSP370 and SSP585 scenarios are −0.18 × 10 6 km 2 /decade, −0.44 × 10 6 km 2 /decade, −0.74 × 10 6 km 2 /decade, and −0.99 × 10 6 km 2 /decade, respectively. The reductions are largest for OND and smallest for JJA, but there is higher SCA climatologically during OND. The reductions are slightly larger for EA compared to NA which may be result of that the land area of EA is larger than that of NA. Figures  10(d)-(f) show the change in SCA projected at 2081-2100 relative to the 1995-2014 baseline. The results are consistent with the trends. Although the absolute reduction in OND SCA is the largest, the reduction relative to the seasonal SCA is larger for the JAS season.
Although the absolute reduction of SCA is smaller for JAS, the climate impact may still be large, because the climatological SCA is already low while solar irradiance is highest during summer. Going from SSP126 to SSP585 and relative to the baseline period, reductions of SCA over the NH for JFM, AMJ, JAS, OND are 5%-20%, 14%-39%, 30%-59%, and 10%-25%. The largest per cent reduction of NH SCA relative to the baseline period will happen in JAS. The changes are similar for the SCA between NH, EA and NA in annual, JFM, AMJ and OND, but not for JAS. Going from SSP126 to SSP585 and relative to the baseline period, the reductions of SCA over EA and NA for JAS are 50%-89%, and 21%-47%, respectively.  , showing monthly statistics for regional land SCA. Left column: monthly mean SCA in 10 6 km 2 ; right column: trends of monthly mean SCA in 10 6 km 2 /decade. Top row: NH; middle row: EA; bottom row: NA. In the left column, black: observation (OBS); blue: CMIP5 MME; red: CMIP6 MME. Colored shading represent the range of one STD for CMIP5 (blue) and CMIP6 (red). In the right column, white: NOAA-CDR (white); blue: CMIP5; red: CMIP6. Error bars indicate the 95% confidence intervals for the observation and the one STD for CMIP5 and CMIP6.

Discussion and conclusions
The satellite-based NOAA-CDR SCE data was evaluated in terms of climatology and trends over the period of 1970-2019. The whole of the NH was evaluated, as well as the EA and NA sub-regions within the NH. Statistics were calculated for the whole year, as well as the four seasons of JFM, AMJ, JAS, and OND.
The historical experiments of CMIP5 and CMIP6 then were validated against NOAA-CDR SCE using the same statistics. Finally, changes in SCE for the next 80 years (2021-2100) were evaluated under four SSPs of CMIP6. From NOAA-CDR, climatological SNF over the NH was larger with higher latitude. SCA was highest in January and lowest in August. Over the 1970-2019 period, there were significant decreasing trends of snow cover in NH in terms of annual mean and seasonal means, except for OND. The decreasing trend of snow cover was larger over EA than over NA. In contrast, mean SCA for OND has increased over NH, at rates up to 0.50 × 10 6 km 2 /decade. The positive trend was most obvious after 1990.
Verification of CMIP5 and CMIP6 historical experiments with NOAA-CDR for the 1970-2014 period confirmed that the CMIP models could reproduce the historical spatio-temporal characteristics of SCE over the NH. Spatial correlations of SNF for individual models were higher than 0.8, even exceeding 0.9 in some models. Biases were generally negative over the northern regions of EA and NA, but generally positive over the Tibetan Plateau and the northwestern regions of China. Annual SNF biases were smaller in CMIP6 compared to the CMIP5. Both CMIPs reproduced the decreasing trend of SCA over most of the validation period. However, the observed increment of OND SCA was not reproduced by CMIP5 and CMIP6 models, especially the trend in the 2000s. This may be due to overestimation of the warming trend during this period by CMIP models, or the uncertainty in observation. CMIP6 models showed a marked improvement over CMIP5 in reproducing the winter climatology of SCA over the NH. Hence, we expect CMIP6 climate projections of SCA over the NH to be more reliable.
Twenty CMIP6 models under four SSPs were evaluated for the next 80 years (2021-2100). Projections indicate future reductions in SNF over the NH, particularly over the northwestern regions of EA, the western regions of NA, and the Tibetan Plateau. The rate of reduction is strongly related to the amount of emissions in the scenario, being smaller in scenarios with lower emissions. Trends of regional SCA are not distinguishable between the four scenarios before 2040, but rapidly diverge after then. SCA stabilizes eventually in SSP126, but continues to drop in the other scenarios, and especially strongly so during JFM and OND. Under SSP585, annual SCA trend is −1.10 × 10 6 km 2 /decade and OND SCA trend is −1.50 × 10 6 km 2 /decade.
Although the largest absolute reduction of SCA is projected to occur during OND, relative to the 1995-2014 baseline, the most significant relative reduction of SCA over the NH will happen during JAS. Under the SSP126 scenario, annual and JAS SCA are reduced by about 10% and 30%, respectively; Under SSP585, these values are 30% and 59%, respectively. Both the projected trend and the relative reduction are significantly larger over EA than over NA. Since SSP585 is the most extreme scenario, and actual SCA reduction should be ameliorated by controls on greenhouse gas emissions.
This study has performed relatively simple validation on the snow cover simulation, since the simulation results were only compared against one satellitebased dataset of SCE, and only aggregated statistics were examined to evaluate model biases. We hope that more studies will continue this line of enquiry and perform more detailed evaluations of regional climate changes in snow cover. The most apparent weakness is the complete reliance on one satellite-based SCE as observation. When the historical simulations fail to reproduce the increasing trend of OND SNF in the 2000s, it is not clear whether this is due to flaws in the models, or biases in the satellite data. A detailed comparison using multiple observation sources would be a necessary follow-up study. Analysis of multiple datasets (Mudryk et al 2015) or an ensemble product (e.g. CanESM, CESM; Mudryk and Derksen 2017) helps to overcome the limitations of individual datasets and provides a better assessment of the climate state. Future continuation of this work will focus on the analysis of multiple observation datasets to further elucidate the evolution of snow cover variability over the NH. Another avenue of investigation is an attributive study of snow cover variability over the NH using large ensemble simulations and the CMIP6 Detection and Attribution Model Intercomparison Project (DAMIP; Gilett et al 2016) experiment.

Data availability statement
The data generated and/or analyzed during the current study are not publicly available for legal/ ethical reasons but are available from the corresponding author on reasonable request. their constructive comments and suggestions. Their advices will benefit the improvement of the paper and our future researches.