Increased European heat waves in recent decades in response to shrinking Arctic sea ice and Eurasian snow cover

In recent decades, unprecedented extreme summer heat waves have occurred in Europe, and they have exhibited an increasing trend since 1970s. Although previous studies have suggested that these recent hot European summers could have been instigated by the underlying surface thermal conditions, the possible influence of shrinking Arctic sea ice and Eurasian snow cover on heat waves are not well understood. Herein, we present evidence obtained via observational analyses and numerical experiments indicating that the interdecadal increase in European heat waves is closely linked to the reductions in Arctic sea ice concentration (ASIC) and Eurasian snow cover fraction (EASC) across mid–high latitudes via the excitation of the anomalous Eurasian wave train. The combined effects of declined ASIC and EASC, accompanied by the drier soil and the stronger heat flux, tend to weaken the poleward temperature gradient at mid–high latitudes and affect the midlatitude jet stream and transient eddy activities. These dynamic and thermodynamic circulations increase the likelihood of more persistent European blocking events that favor frequent and strengthened heat waves. Further projection analysis of simulations from 13 CMIP5 climate models suggests that Europe may experience more hot summers as the ASIC and EASC continue to decline over the next century.


INTRODUCTION
In recent decades, especially since 2003, Europe has experienced an exceptional number of pronounced heat wave events (HW) in summer, [1][2][3][4] and these events have had disastrous impacts on local communities and economies. [5][6][7] The causes of these European heat waves are complex, and many studies have primarily attributed them to the possible impacts of large-scale atmospheric circulation patterns, such as the North Atlantic Oscillation (NAO), 8 and the midlatitude planetary waves, [9][10][11] as well as ocean and land surface thermal conditions, such as the El Niño-Southern Oscillation (ENSO), 12 North Atlantic warming, 13 and Eurasian soil moisture variations. 14,15 Other studies have noted that anthropogenic global warming, particularly the Arctic temperature amplification, has contributed to the increasing number of summer heat wave events in Europe. 3,16 During this process, heat waves are often associated with anomalous circulation patterns characterized by the weakening of the zonal mean jets, the quasiresonant amplification of planetary waves, [9][10][11]17 and persistent blocking conditions 12,18 over the midlatitude Eurasia that produce prolonged hot surface conditions. The weakening of the zonal mean jets is generally due to the reduction of the equator-to-pole thermal gradient associated with changing surface conditions. Quasi-stationary midlatitude waves are most likely amplified by the quasi-resonance between free and forced waves in waveguides with zonal wavenumbers of 6~8. For example, the European heat wave in 2003 and the Russian heat wave in 2010 were associated with the strong amplification of midlatitude quasi-stationary waves. 9,10 The occurrence of a 2010 Russian heat wave was associated with persistent blockings in the European-Western Russian region. 12,18 These frequent and amplified extreme heat waves in Europe have been accompanied by record declines in Arctic sea ice 19 and Eurasian snow cover. 20 In particular, in spring and summer they have dramatically diminished since the 1980s. Studies have reported the underlying mechanisms as the shrinking cryosphere-induced modulation of midlatitude circulation [21][22][23][24] and adjustment of the surface energy and water balances. 22 Limited studies have also explored the impacts of shrinking summer cryosphere on European heat waves, including Arctic sea ice and highlatitude snow cover [25][26][27] and Tibetan Plateau snow cover. 23 However, whether there is robust atmospheric response to the continuous loss of spring-summer Arctic sea ice and Eurasian snow cover has yet to be examined. Furthermore, whether and how the increase in European summer heat waves is related to the combined impacts of diminishing sea ice and snow cover remain unclear. Therefore, in the present study, we conducted observational analyses and numerical experiments to address these questions. Figure 1 shows the linear trends of summer heat wave magnitude index (HWMId 4 ) in Europe (30°-75°N, 0°-45°E) between 1980 and 2015, calculated from the gridded daily maximum temperature data from Hadley Center-Global Historical Climatological Network (HadGHCN) and National Centers for Environmental Prediction/ National Center for Atmospheric Research (NCEP/NCAR), respectively. The HWMId is defined as the duration and magnitude of heat wave lasting for at least 3 days (see Methods). Generally, their spatial distributions exhibit uniform increase in heat waves across Europe and present an activity center located in Southern and Eastern Europe with a maximum value of >0.8/decade (Fig. 1a, b). Figure 1c displays the interdecadal change of the European areaaveraged HWMId during 1980-2015. Although entire period shows an increasing trend of heat wave, as reported by previous studies, 1-4 decadal changes are also evident in these indices. It can be seen that heat waves occurred more frequently during 2003-2015, particularly these intense heat wave extremes observed in the summers of 2003, 2007, 2010, and 2015, 12,28 but less heat waves are found during 1980-1997. The decadal transition around the year 1997 is also supported by the moving ttest and 9-year smoothing methods that applied to the HWMId indices. These results indicated that European heat waves occur less frequently in the early period, with small magnitudes and interannual variations, whereas occur with high frequency and strong interannual variations in the later period. To understand the interdecadal variations of European heat waves and associated causes, we conduct composite analyses separately for the high-HW period (2003-2015) and the low-HW period (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997).

RESULTS
To investigate the relationship between the European HWMId and Arctic sea ice concentration (ASIC)/Eurasian snow cover fraction (EASC), the composite anomalies of fields from the preceding spring to summer in high-and low-HW periods are calculated in Fig. 2. Difference between these two periods (high-HW minus low-HW) is also shown in Figure. During the high-HW period, in spring low sea ice occur in regions adjacent to the Barents Sea, in contrast, during the low-HW period high sea ice are seen in the northern Barents Sea (Fig. 2a, b). Composite difference in ASIC between the two periods reveals a reduced ASIC in the Barents Sea (Fig. 2c), similar to that during the high-HW period. The retreat of ASIC continues until summer, with spatial range expanding to the entire basin (Fig. 2d, f). When it comes to EASC, reduced snow coverage is found over regions extending from Southern and Eastern Europe to Siberia and regions around the northwestern Tibetan Plateau in spring (Fig. 2g, i). In summer, snow cover exists only in highlatitude Eurasia and the northwestern Tibetan Plateau regions, thus reduced snow coverage is  barely visible there (Fig. 2j, l). The opposite distribution features of ASIC and EASC are observed during the low-HW period (Fig. 2e, h, k). These analyses imply that the recent frequent and intensified hot European summers may be linked to the continuous loss of spring-summer ASIC and EASC.
To clarify the possible combined effects of ASIC and EASC on the European heat waves, we define a combined ASIC and EASC index (hereafter referred to as ICE/SNOW, see Methods) following a previous study. 27 Figure 1d shows the historical time series of normalized spring and summer ICE/SNOW indices from 1980 to 2015. These ICE/SNOW indices exhibit a significant decreasing trend over entire period, and the decreasing trend is accelerated after 2000, similar to the interdecadal evolution of HWMId. The spring ICE/SNOW index is correlated with the HWMIds calculated from the HadGHCN and from the NCEP/NCAR data (the Kendall Tau correlation are R = 0.32 and 0.33, respectively; p < 0.05). The correlation is slightly higher for the summer ICE/SNOW index (R = 0.38 and 0.30, respectively; p < 0.05). These correlations indicate a plausible connection between the shrinking ASIC and EASC and interdecadal variation of European heat waves, with the former in spring likely forcing the latter.
The loss of mid-high-latitude EASC in spring and summer may have a thermodynamic influence on atmospheric circulation by allowing the underlying soil to dry out and by strengthening the surface heat flux exchange from the surface to the atmosphere. 22,24 Figure 2m-p show composite differences in soil moisture and surface sensible heat flux from spring to summer between the high-and low-HW periods. As the EASC vanishes, the soils over widespread regions of highlatitude Eurasian show consistent desiccation from spring ( Fig. 2m) to the following summer (Fig. 2n). This dried soil further changes the underlying surface heat balance, 28,29 and leads to significantly stronger upward heat flux from the soil to the atmosphere (Fig. 2o, p). Meanwhile, the low ASIC leads to a stronger and persistent latent heat flux until summer over the Barents Sea region (figure not shown), which might modulate the midlatitude circulation in summer. 19,21 This persistent dry soil and/or high heat flux from spring to summer can be regarded as a "thermodynamic land/ ocean memory" 14,15 that connects the changing spring ASIC and EASC with subsequent summer atmospheric circulation.
We next explore the 200-hPa geopotential height (Z200), mean zonal wind (U200), and synoptic-scale transient eddy activity (STEA) patterns that related with the decadal changes in European heat waves. As shown in Fig. 3a, the Z200 pattern obtained from reanalysis is characterized by a classic "Omega" blocking pattern, 30 with a clear zonal " + -+ " teleconnection wave train over the mid-high-latitude Eurasia. Two positive centers are located in Southeastern Europe and the East Asia-Okhotsk Sea, and a negative center is located in the region to the west of Lake Balkhash. The height patterns at 500 hPa and 700 hPa are quite similar, although the positive anomalies are relatively weaker (Supplementary Fig. 1a-b), suggesting a barotropic tropospheric response. Areas of strong ridges depict more meridional flow and indicate a slow-moving Rossby wave. Such a blocking-like circulation pattern can guide us to further detect historical Omega blocking events (see Methods), particularly on an interdecadal timescale. Figure 4 shows the variations in Omega blocking duration (BDI) and maximum duration (BMI) and correspondent de-trended components. The BDI is defined as the number of Omega blocking event days throughout summer, and the BM is defined as the maximum number of blocking days in individual Omega blocking events throughout the summer. The BDI displays a slightly increasing trend over the period 1980-2015 (Fig. 4a) and is highly correlated with the HWMId (R = 0.32; p < 0.05). However, the increasing trend of BMI is inapparent (Fig. 4c), and BMI is weakly related to the HWMId (R = 0.23; p < 0.1). A comparison between the BDI and BMI calculated form the NCEP/NCAR reanalysis demonstrates that they are relatively high during the high-HW period (15 days for BDI and 8.7 days for BMI) and low during the low-HW period (7 days for BDI and 5.6 days for BMI, Fig.  4b, d). This suggests that the interdecadal increase in summer heat waves in Southern and Eastern Europe corresponds well to the occurrence of frequent Omega blockings.
We present evidence showing that the recent increases in European heat waves and blockings may be related to the combined influence of reductions in the ASIC and EASC via changes in large-scale circulation and synoptic-scale eddy activities. 12,31 The JJA-mean STEA 30,32 is defined as the high-pass Z200 for the 2-8-day Lanczos filter. Accompanying the occurrence of more blockings over Europe, the accelerated U200 was located along the Arctic coast, whereas the decelerated U200 extended from Southern Europe to Lake Balkhash (Fig. 3b). To some extent, the stronger zonal wind in highlatitudes and the weaker one in midlatitudes encourage the development of blockings and anticyclonic circulation 33 in Southern and Eastern Europe. In addition to the predominant role of large-scale zonal wind, the STEA plays a secondary role in the formation and maintenance of European blockings. 31 In Fig. 3c, the ASIC and EASC loss correspond to significantly stronger baroclinic eddies in the regions from the south of the Barents Sea to the Siberian Plateau, and clearly decrease over Southwestern and Southeastern Europe. Relative to the climatological position, the anomalous STEA tends to move northward. Such a decrease in the STEA over Europe may further drive positive eddy-vorticity forcing (figure not shown), which is spatially coherent with the enhanced Z200 in Fig. 3a, thus reflecting positive STEA-Z200 feedback in the upper troposphere. Under these conditions, rising daily maximum air temperature (TMX) from NCEP/NCAR reanalysis govern most regions of Europe, with magnitudes of up to 1.5°C and higher (Fig. 5a). These physical processes suggest that the frequent occurrence of summer blockings and strengthened heat waves in Europe may result from ASIC and EASC loss via changes in both large-scale Rossby waves and STEAs. These aforementioned statistical linkage is analyzed based on reanalysis data, but it does not guarantee the causation relationship. Further modeling analysis is needed to address this question.
To elucidate the causality between the declining ASIC and EASC and increased European heat waves, we successfully detected associated circulation signatures based on a large ensemble of atmospheric general circulation model simulations. Three pairs of comparative experiments 34 Table 1). These experiments are conducted by prescribing the ASIC with the observed monthly composites for low-and high-ASIC years, and by altering snowfall rate for low-and high-EASC years, as the sea surface temperature (SST) fixed as the control runs (i.e., using a monthly varying climatology). In each experiment, 50 ensemble members are integrated from 1 March to 31 August, with different initial conditions from the control run. Atmospheric responses to changed ASIC and/or EASC are ultimately compared with those in observation.
In the sea ice experiments, the tropospheric circulation anomalies associated with ASIC anomalies in observation are largely captured by the model ensemble-mean response (Supplementary Fig. 2), although the positive center on Z200 is confined to Eastern Europe (Fig. 3d) in comparison with observations (Fig. 3a). The decrease in ASIC is coupled with a meridional dipole anomaly of the U200 (Fig. 3e) and a decreased STEA across the Eastern Europe-Ural Mountain regions and the south of the Mediterranean Sea (Fig. 3f), that tend to move southward relative to the observations. The amplified high-pressure ridge in Eastern Europe is favor of frequent (Fig. 4b) and persistent (Fig. 4d) blockings, where the BDI and BMI are greater in the LICE experiment than those in the HICE experiments. We also note that these indices are lower than those in the high-HW period in observation, indicating that ASIC forcing has a limited contribution to European blockings. This suggests that the impact of melting ASIC on the summer heat waves is confined to Eastern Europe, leading to roughly 1.0°C higher TMX (Fig. 5b), although the variations in heat waves cannot be completely explained by the imposed ASIC anomaly alone.
In the snow cover experiments, to obtain a more realistic snow distribution, we roughly reduced the snowfall rate by half in the Eurasian mid-high latitudes and the northwestern Tibetan Plateau in the LSNOW experiments and doubled it in the HSNOW experiments. By using this experimental design, this model largely mimics the spatial patterns of spring-summer Eurasian snow cover anomalies ( Supplementary Fig. 3), which are consistent with the observational results. As the EASC declines, a " + − + " teleconnection wave train is evidently reproduced in Z200, with an anomalous ridge over Southern and Eastern Europe (Fig. 3g). The weakening of U200 in the Southeastern Europe is generally consistent with that in observation (Fig. 3h). Meanwhile, a zonal belt-like pattern of decreased STEA prevails across Southern Europe (Fig. 3i). Under these conditions, the TMX tends to rise as much as 1.0°C (Fig. 5c). The European BDI and BMI in the LSNOW experiments are greater than those in the HSNOW experiments (Fig. 4b, d). The fact that the BMI is also greater than that in the LICE experiment suggests that the EASC forcing can induce prolonged blockings and then explain more changes in the heat waves over Southern Europe than ASIC forcing. However, the BDI is weaker in the LSNOW experiment than that in the LICE experiment, indicating that EASC contributes less to changes in the blocking days relative to ASIC.
To explore the combined effect of the ASIC and EASC on summer heat waves in Europe, we imposed anomalous ASIC and EASC forcings in the model simultaneously. These two forcings are identical to those in the separate experiments. This experiment qualitatively reproduces the snow cover distribution in the Eurasian region ( Supplementary Fig. 4), and mimics that in observation (Fig. 2i, l) and in snow cover experiments (Supplementary Fig. 3). The prevalence of a "− + − + " zonal wave train pattern and positive height anomalies over the Southern and Eastern Europe (Fig. 3j) suggests that in response to ASIC and EASC loss, European blockings may strengthen the heat waves. The responses of the Z200, U200 and STEA (Fig. 3k, l) show patterns similar to their observational counterparts, and their amplitudes reach more than half the values of the observations, suggesting that the combination experiment can explain >50% of the circulation changes contributing to heat waves. Associated with the change in atmosphere circulation, positive TMX anomalies exceeding 1.5°C are found throughout Europe (Fig.  5d). The BDI and BMI indices are greatest among all these experiments, although the amplitudes are smaller than their observational counterparts (Fig. 4b, d). These results illustrate the importance of combined effect of ASIC and EASC on the interdecadal increase in European summer heat waves.
As noted in observations, decreasing trends of ASIC and EASC have been observed in recent years against a global warming background, and they may continue to decline in the future. 35 Contrarily, the summer blockings and heat waves in Europe may increase in the future 4 based on the mechanisms discussed above.  Table 2). Our results show significant decreasing trends in the spring and summer ICE indices, and these decreases could reach approximately 35 and 30%, respectively, by 2100 (Fig. 6a, b). The SNOW indices in spring and summer successfully reproduce the snow cover retreat for the observational period of 1980-2005 and the projected period of 2006-2100, with SNOW decreasing by approximately 26 and 80%, respectively, by 2100 (Fig. 6c, d). The above results show that the retreat of the ASIC and EASC at the beginning of the twenty-first century is not only a short-term phenomenon but will also continue or even be exacerbated in the next century. Along with that, the BDI and HWMId indices show increasing trends, and we find that the ensemble-mean simulations show good performance in simulating their variations. The  projected BDI will increase by~10 days in 2100, and the HWMId by roughly six times (Fig. 6e, f). The saturation of the BDI after 2080 may result from the parameter settings in the RCP4.5 simulation, wherein the CO 2 concentration, emissions and radiative forcing (in addition to the ASIC and EASC forcings) are set to be relatively stable after 2080. 37 The increases in BDI and associated HWMId in recent decades might be driven by thermodynamic processes, in which the drying trend of soil moisture caused by the losing EASC could enhance the surface heating, 38 resulting in the occurrence of more heat waves. Meanwhile, the surface thermodynamic conditions associated with reduced ASIC, EASC, and soil moisture may also interact with the dynamic changes through promoting weakened transient eddies 25 and quasi-resonant waves, 39 which would dynamically induce blockings and thus cause frequent heat waves. These marked linear trends provide a potential linkage between the combined ASIC and EASC and summer heat waves in Europe over the next century. Therefore, more frequent and intensified hot summers in Europe will continue in a warming scenario, although projection uncertainties remain due to the insufficient number of CMIP5 ensemble members.

DISCUSSION
The current study provides robust evidence linking the diminishing spring-summer ASIC and EASC with a greater summer HWMId in Europe in recent decades via the modulation of the European blockings using observational/reanalysis data, which is supported by the CAM large-ensemble experiments. This study provides important insights into the cause of European heat waves because previous studies mainly focused on separate summertime forcings and explained their influences through the modulation of largescale circulation. 23,27 Although the uncertainties associated with other external forcings and internal varibilities are large, in the current study the dramatic decreases in ASIC and EASC have considerable potential to modulate the decadal variations of extreme summer weather in midlatitudes. The synchronic and synergistic impacts of ASIC and EASC provoke large-scale atmospheric circulation anomalies in the southern and eastern parts of Europe, as well as blockings and heat waves. Based on both observational and modeling analyses, three mechanisms are proposed to explain such impacts: (i) a weakening of the midlatitude STEA, (ii) a "double-jet" regime in Arctic-subtropical regions due to enhanced highlatitude warming, and (iii) a midlatitude wave train pattern created by the double waveguides. These explanations of dynamical changes are consistent with those of previous observational studies 9,10 and are further confirmed by numerical experiments. An important issue should be addressed in the current study is that the diagnostics are strongly affected by the change in mean summer temperature. It is difficult, however, to distinguish whether the results are a consequence of the changing mean climate 37 or increasing extremes. 40 In Fig. 1c, we examine the variations in de-trended HWMId by removing the linear trends from the raw data to eliminate the interference of mean climate. It can be seen that the interdecadal variations and increasing trends in the raw data are largely preserved. Moreover, the original HWMId explain >50%6 0% of the local heat wave variance, while the de-trended indices account for >40%-50%. Therefore, mean climate change and increasing extreme events are both important in promoting the more frequent and intense European heat waves. Another important issue is the relative contributions of advection and turbulent heat flux (i.e., the latent and sensible heating) to the surface temperature. Supplementary Fig. 5 shows the spatial distribution of calculated relative explained variances of these contributions. The temperature advection in the mid-lower troposphere can account for~50% of the mean temperature variance in Europe and 30% of the TMX variance in Southern Europe ( Supplementary Fig. 5a, c). In contrast, the turbulent heat flux explains~30% of the mean temperature variance and 50% of the TMX variance throughout Europe ( Supplementary Fig. 5b, d).
These analyses indicate that these two contributors may explain a large proportion of temperature variations, and that they are of great importance in boosting summer hot extremes. These two questions need to be extensively discussed in future studies.

METHODS Definitions
The methods used to identify the summer heat wave duration and magnitude index (HWMId), blocking duration (BDI) and magnitude (BMI) indices are described in this section. A heat wave is defined when the daily TMX is higher than the 90th percentile of TMX for the reference period of 1980-2015, 1,2 centered on a 15-day window, and lasts >3 consecutive days. Following Russo et al., 4 the HWMId is calculated as the maximum magnitude of the heat waves occurring in a summer to take into account both heat wave duration and intensity. The HWMId magnitude is the sum of the daily magnitude of the consecutive days within a heat wave. The daily magnitude is calculated as follow: with T d being the TMX on day d of the heat wave, T 25p and T 75p are the 25th and 75th percentile values, respectively, of the time series composed of TMXs from the reference period. Atmospheric blocking activities are often characteristically quasistationary and present a long duration and closed anticyclonic circulation. This activity can briefly alter the dominant westerly circulation in the troposphere. Following its definition in previous studies, 30,41 the Omega blocking index is defined as the difference in the 200-hPa geopotential height between the region with positive values (30°-60°N, 20°-50°E) and the region with negative values (35°-48°N, 60°-75°E) in the NCEP/NCAR dataset. An Omega blocking event was identified with daily blocking index exceeding one standard deviation and period >5 consecutive days. The BDI refers to the number of blocking days throughout the summer, whereas the BMI refers to the maximum number of blocking days in an individual blocking event in summer.
Following Tang et al., 27 the combined ASIC and EASC index (ICE/ SNOW) is defined as ICE/SNOW = ICE * P ICE + SNOW * P SNOW , in which the ICE index is defined as the average ASIC anomaly in Baffin Bay (65°-80°N, 70°-50°W) and the Barents Sea (74°-80°N, 30°-90°E), the spring (summer) SNOW index is defined as the average EASC anomaly in the mid-high-latitude Eurasia (35°-70°N, 40°-120°E) [the northwest Tibetan Plateau (31°-41°N, 70°-80°E)], and the P ICE and P SNOW represent the percentage slopes of the ICE and SNOW indices, respectively, from 1980 to 2015.
The statistical methods used in this paper include Kendall's tau-based slope estimator and trend analysis, composite analysis, and Kendall's tau and Pearson correlation analysis. We also use Student's t-test to assess the significance of the composite and correlation analyses.

CAM model experiments
We used the NCAR CAM3.1 model for the numerical experiments. The resolution was T42 spectral truncation (with a horizontal spatial resolution of 128 × 64) and 26 vertical levels (in which the top layer is~3.7 hPa). 42 Studies have shown that CAM3.1 can capture the summer climate and circulation features over Eurasia. 22,43 Three pairs of comparative experiments, i.e., LICE and HICE experiments, LSNOW and HSNOW experiments, and LICESNOW and HICESNOW experiments (Supplementary Table 1), were driven by the identical anomalous boundary conditions of the ASIC and EASC forcings. The experimental protocol was as follows.
Forcings: atmosphere-only experiments require ASIC and SST data to be prescribed. In the LICE and HICE experiments, the annually repeating monthly ASICs from March to August were derived from HadISST for the periods of 2003-2015 and 1982-1997, respectively. In the LSNOW and HSNOW experiments, the snowfall rates in northern Eurasia (40°-70°N, 40°-120°E) were halved, and in the northwestern Tibetan Plateau (31°-41°N, 70°-80°E) were doubled. In the LICESNOW and HICESNOW experiments, the combined ASICs and snowfall rate anomalies were imposed simultaneously. All other external forcings (e.g., SST) were set as annually repeating cycles of climatology. Initial field and integration: to obtain relatively reliable and robust simulation results, 34 in each experiment, 50-member ensembles were integrated forward for 60 years, with the first 10 years discarded as spin-up. The initial conditions were from a control run that was forced with the climatological ASIC and SST for 1980-2015. Presentation of results: The simulation results are presented as the composite differences between respective pairs of experiments, i.e., LICE-HICE, LSNOW-HSNOW, and LICESNOW-HICESNOW.

DATA AVAILABILITY
The datasets analyzed in this study are available on-line. The atmospheric circulation and daily maximum surface air temperature data are from the NCEP/NCAR reanalysis dataset 44 . HadGHCN daily maximum surface air temperature data  are from the Met Office Hadley Center 45 . ASIC and SST data with a spatial resolution of 1.0°× 1.0°are from the Hadley Center in the United Kingdom (HadISST). 46 The version v01r01 of the weekly Northern Hemisphere Snow Cover Extent data are from the NOAA National Centers for Environmental Information (NCEI). 47 Snow water equivalent data for areas north of 35°N in the Northern Hemisphere with a spatial resolution of 25 ×25 km are from the Finnish Institute of Meteorology 48 (http://www.globsnow.info/). Soil moisture and surface heat flux data are from the ERA-Interim reanalysis dataset 49 (https://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=sfc/). The soil moisture data includes four soil layers (7, 21, 100, and 289 cm). The data used in the present study are primarily from the surface layer (0-7 cm). CMIP5 data are from https://esgfnode.llnl.gov/projects/cmip5/. The modeling data that support the findings of this study are available from the corresponding author upon request.