Eurasian autumn snow link to winter North Atlantic Oscillation is strongest for Arctic warming periods

In recent years, many components of the connection between Eurasian autumn snow cover and wintertime North Atlantic Oscillation (NAO) have been investigated, suggesting that November snow cover distribution has strong prediction power for the upcoming Northern Hemisphere winter climate. However, the non-stationarity of this relationship could impact its use for prediction routines. Here we use snow products from long-term reanalyses to investigate interannual and interdecadal links between autumnal snow cover and atmospheric conditions in winter. We find evidence for a negative NAO-like signal after November with a strong west-to-east snow cover gradient, which is valid throughout the last 150 years. This correlation is consistently linked to a weak stratospheric polar vortex state. Nevertheless, decadal evolution of this link shows episodes of decreased correlation strength, which co-occur with episodes of low variability in the November snow index. By contrast, periods with high prediction skill for winter NAO are found in periods of high November snow variability, which co-occur with the Arctic warming periods of the 20th century, namely the early 20th-century Arctic warming between 1920 and 1940 and the ongoing anthropogenic global warming at the end of the 20th century. A strong snow dipole itself is consistently associated with reduced Barents–Kara sea ice concentration, increased Ural blocking frequency and negative temperature anomalies in eastern Eurasia.

1 Introduction

As the leading climate variability pattern affecting winter climate over Europe (Thompson and Wallace, 1998), the North Atlantic Oscillation (NAO) has been extensively studied over the last decades (Wanner et al., 2001; Hurrell and Deser, 2010; Moore and Renfrew, 2012; Pedersen et al., 2016; Deser et al., 2017). The NAO has been defined as the variability of the pressure gradient between Iceland (representing the edge of the polar front) and the Azores (representing the subtropical high ridge). The sign of the NAO is related to weather and climate patterns stretching from local to continental scales. Since its variability has severe socioeconomic, ecological and hydrological impacts for adjacent continents, seasonal to decadal predictions of the state of the winter NAO are high-priority research for many climate science centres (Jung et al., 2011; Kang et al., 2014; Scaife et al., 2014, 2016; Smith et al., 2016; Dunstone et al., 2016; Wang et al., 2017; Athanasiadis et al., 2017).

Together with the rapid warming of the Arctic and the increased frequency of severe winters over Eurasia and North America (Yao et al., 2017; Cohen et al., 2018; Kretschmer et al., 2018; Overland and Wang, 2018), recent studies have highlighted the state of the Northern Hemispheric cryosphere as a useful predictor for the boreal wintertime (December–January–February, DJF) NAO (Cohen et al., 2007, 2014; Vihma, 2014; Garcia-Serrano et al., 2015; Cohen, 2016, Orsolini et al., 2016; Crasemann et al., 2017; Warner, 2018). Although both systems seem to be connected (Cohen et al., 2014; Furtado et al., 2016; Gastineau et al., 2017), the emerging main hypothesis connects reduced autumn Barents–Kara sea ice concentration and increased Siberian snow cover with a negative NAO state in the following winter months (Cohen et al., 2014).

The proposed mechanism behind this hypothesis is a multi-step process, starting with autumn sea ice loss for the European Arctic, followed by altered tropospheric circulation due to elevated Rossby wave numbers, vertical propagation of said Rossby waves upward into the stratosphere and consequently a weakening of the polar vortex (see Cohen et al., 2014, for an in-depth discussion). With the weakening (or the reversal) of the polar vortex, a stratospheric warming signal manifests itself. This signal propagates slowly back into the troposphere, where it manifests itself as a negative NAO, connected to the concurrent cold winters for Eurasia (Kretschmer et al., 2018).

In recent years, many components of this pathway have been investigated, especially concerning the increased frequency of cold winters over Europe and the emergence of the counter-intuitive “warm Arctic–cold continent” (WACC) pattern over Eurasia (Petoukhov and Semenov, 2010; Vihma, 2014). However, there remains substantial uncertainty about the impact of Arctic sea ice in terms of location (Zhang et al., 2016; Luo et al., 2017; Screen 2017; Kelleher and Screen, 2018), timing (Honda et al., 2009; Overland et al., 2011; Inoue et al., 2012; Suo et al., 2016; Sorokina et al., 2016; King et al., 2016; Screen, 2017; Wegmann et al., 2018a; Blackport and Screen, 2019) or whether sea ice can be used as a predictor or forcing at all based on the contrasting results of model studies (McCusker et al., 2016; Collow et al., 2017; Pedersen et al., 2016; Boland et al., 2017; Crasemann et al., 2017; Ruggieri et al., 2017; Garcia-Serrano et al., 2017; Francis, 2017; Screen et al., 2018; Mori et al., 2019; Hoshi et al., 2019; Blackport et al., 2019; Romanowksy et al., 2019).

The interplay between Arctic sea ice and Siberian snow is much less explored. Ghatak et al. (2010) showed that reduced autumn polar sea ice leads to the emergence of increased Siberian winter snow cover, especially so in the eastern part of Eurasia. This dipole signal was amplified in coupled climate model runs for the 21st century, where sea ice is substantially diminished. In an observational study, Yeo et al. (2017) point out that the moisture influx from the open Arctic ocean into the Eurasian continent contributes to the increase in snow cover, a mechanism described by Wegmann et al. (2015). Gastineau et al. (2017) found that reduced sea ice is connected to a distinct November snow dipole over Eurasia, both in reanalysis and model data. They further state that the snow component is a statistically more powerful predictor than sea ice for the atmosphere in the following winter. This relationship was also found in a range of climate models, albeit with weaker links. Xu et al. (2019) found the same correlation in observational and model data, however looking at winter climate only. Based on their analysis, the authors state that the enhanced snow cover in winter is a product of the negative NAO rather than a precursor. Sun et al. (2019) highlight the importance of elevated North Atlantic sea surface temperatures for the development of a Eurasian snow dipole in autumn. This warming of the North Atlantic favours reduced sea ice cover for the European part of the Arctic, which triggers a high-pressure anomaly over the northern Ural Mountains via increased ocean-to-atmosphere heat fluxes, transporting cold air masses towards the south of its eastern flank.

The possible impact of the Siberian snow on the stratosphere and eventually on the NAO is well discussed in Henderson et al. (2018). Although observational NAO prediction studies with Siberian snow showed great success in the past (Cohen and Entekhabi, 1999; Saito et al., 2001; Cohen et al., 2007, 2014; Han and Sun, 2018), links between snow and the stratosphere still seem to be missing or too weak in model studies (Furtado et al., 2015; Handorf et al., 2015; Tyrrell et al., 2018; Gastineau et al., 2017; Peings et al., 2017), whereas nudging realistic snow changes to high-resolution models seems to improve the prediction skill (Orsolini and Kvamsto, 2009; Orsolini et al., 2016; Tyrrell et al., 2019). Moreover, even though the stratosphere–surface connection is now reasonably well established (Kretschmer et al., 2018), the timing and location of the snow cover used for the prediction is, as with sea ice, still debated (Yeo et al., 2016; Gastineau et al., 2017). As an additional caveat, Peings et al. (2013) and more recently Douville et al. (2017) showed that the proposed autumn snow-to-winter NAO relationship is non-stationary for the 20th century. A possible modulator for that relationship might be the phase of the Quasi Biennial Oscillation (QBO) (Tyrrell et al., 2018; Peings et al., 2017; Douville et al., 2017). Peings (2019) argues that neither snow nor sea ice anomalies trigger the stratospheric conditions needed to produce winter extremes and that instead high tropospheric blocking frequency over northern Europe leads to the cryosphere anomalies.

Here, we follow up on the definition of a November Eurasian snow cover dipole (Ye and Wu, 2017; Gastineau et al., 2017; Han and Sun, 2018), which was identified to provide predictive power for the following winter months at the end of the 20th century. It is, however, unclear if this prediction skill is stable for time periods further back than 30 years and how it evolves in periods of high Arctic sea ice cover. In this study we address the question of (a) the non-stationarity of the Eurasian snow cover to winter European surface climate relationship in the 20th century, (b) the importance of snow versus sea ice as a predictor and (c) possible precursors/modulators of the sea-ice–snow–stratosphere chain. With this we aim to contribute to the understanding of impacts of cryosphere variability on mid-latitude circulation (Francis, 2017; Henderson et al., 2018; Cohen et al., 2020). To this end, we utilize centennial reanalyses and reconstruction data, where we focus on the transition from October to November to DJF to facilitate the idea of seasonal prediction.

This paper is organized as follows. Sect. 2 describes the data and methods used. In Sect. 3, we introduce the snow cover indices and their interannual prediction value. Section 4 investigates interdecadal shifts in the correlation between snow cover and NAO as well as possible determining factors. The results are discussed in Sect. 5 and finally summarized in Sect. 6.

2 Data and methods

2.1 Atmospheric reanalyses

To evaluate long-term reanalyses, we use snow cover, snow depth and atmospheric properties from the MERRA2 reanalysis (Gelaro et al., 2017). MERRA2 has a dedicated land surface module and was found to reproduce local in situ snow conditions over Russia very well (Wegmann et al., 2018b). For a detailed description of how MERRA2 computes snow properties, see, e.g., Orsolini et al. (2019).

To cover the 20th century and beyond, we include two long-term reanalyses in this study, namely the NOAA-CIRES 20th-century reanalysis Version 2c (20CRv2c) (Cram et al., 2015) as well as the Centre for Medium-Range Weather Forecasts (ECMWF) product ERA-20C (ERA20C; Poli et al., 2016). From the ERA20C product we use snow depth, whereas from 20CRv2c we investigate snow depth and snow cover. Both reanalyses were found to represent interannual snow variations over Eurasia remarkably well. For an in-depth discussion of their performance and their technical details concerning snow computation, see Wegmann et al. (2017b). We also performed the same analysis using the coupled ECMWF reanalysis CERA20C (Laloyaux et al., 2018) but found no added knowledge gain over ERA20C. Thus, we do not include CERA20C in any further analysis.

We use detrended anomalies of these three reanalysis products to extend the October and November index proposed by Han and Sun (2018) into the past, where the November index is in essence the snow dipole described by Gastineau et al. (2017) using maximum covariance analysis (Fig. 1). Whereas the October index is just calculated as the field average snow cover, the November index is computed as the difference between the eastern and the western field averages. It should be noted that Han and Sun (2018) found the November index to be linked to a negative NAO and colder Eurasian near-surface temperatures, whereas the October index was correlated with warmer than usual temperatures over Eurasia and a southward-shifted jet. However, since many studies focus on northern Eurasian October snow cover as the predictor for winter climate, we will include it nonetheless. MERRA2 and 20CRv2c offer snow cover as well as snow depth as a post-process output; however, ERA20C only offers snow depth. We refrained from converting it to snow cover ourselves but found the index based on snow depth to be extremely similar (also see Supplement Fig. S1) to the same index using snow cover. Moreover, comparing snow indices from reanalyses with snow indices using the NOAA Climate Data record of Northern Hemisphere snow cover extent (Robinson et al., 2012), which incorporates satellite data, does not highlight any meaningful differences (Supplement Fig. S2). All snow indices are normalized and linearly detrended with respect to their overall time period. Generally, we found the long-term reanalyses to be of comparable quality to MERRA2 during the overlapping periods.

Figure 1 (a) Regions for October and November snow index used in this study. (b) Linearly detrended and standardized October snow index comparison for the 20th century for snow cover (SC) and snow depth (SD) variables. Panel (c) same as (b) but for the November snow dipole.

Besides snow properties we use detrended atmospheric and near-surface anomaly fields from all three reanalyses. Moreover, as Douville et al. (2017), we use the field-averaged (60–90^∘ N) 10 hecto-pascal (hPa) geopotential height (GPH) anomalies in ERA20C as a surrogate for polar vortex (PV) strength. Although ERA20C only assimilates surface pressure, correlation between this stratospheric index in ERA20C and MERRA2 during the overlapping time periods is higher than 0.9.

The ERA20C 10 hPa November–December mean GPH shows remarkable interannual agreement with state-of-the-art reanalyses that assimilate upper-air data for the period 1958–2010 (see Supplement Fig. S3). Moreover, MERRA2 and ERA20C 10 hPa GPH anomalies agree best over the northern polar regions with correlation coefficients of >0.9 for the period 1981–2010 (see Supplement Fig. S3). This fact supports the extended value of the ERA20C polar stratosphere. Before 1958, the quality of the ERA20C stratosphere is difficult to assess, but the comparison with reconstructions of 100 hPa GPH zonal means shows very good agreement for late autumn and winter months (see Supplement Fig. S4). As the 20CRv2c ensemble mean dilutes the interannual variability signal back in time with increased variability within the ensemble members, we use the deterministic run of ERA20C for the following stratosphere analyses.

We use 6-hourly 500 hPa GPH fields (GPH500) to calculate monthly blocking frequencies according to Rohrer et al. (2018). Blockings are computed according to the approach introduced by Tibaldi and Molteni (1990) and are defined as reversals of the meridional GPH500 gradient. In accordance with Scherrer et al. (2006) the one-dimensional Tibaldi and Molteni (1990) algorithm is extended to the two dimensions by varying the latitude between 35 and 75^∘ instead of a fixed latitude.

• i.

  GPH500 gradient towards the pole:

  $$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{GPH}\mathrm{500}{\mathrm{G}}_{\mathrm{P}}& ={\displaystyle \frac{\mathrm{GPH}{\mathrm{500}}_{\mathit{\phi }+d\mathit{\phi }}-\mathrm{GPH}{\mathrm{500}}_{\mathit{\phi }}}{d\mathit{\phi }}}\\ & <-\mathrm{10}{\displaystyle \frac{\mathrm{m}}{{}^{\circ }\mathrm{lat}}}\end{array}\end{array}$$
• ii.

  GPH500 gradient towards the Equator:

  $$\begin{array}{}\text{(2)}& \mathrm{GPH}\mathrm{500}{\mathrm{G}}_{\mathrm{E}}={\displaystyle \frac{\mathrm{GPH}{\mathrm{500}}_{\mathit{\phi }}-\mathrm{GPH}{\mathrm{500}}_{\mathit{\phi }-d\mathit{\phi }}}{d\mathit{\phi }}}>\mathrm{0}{\displaystyle \frac{\mathrm{m}}{{}^{\circ }\mathrm{lat}}}\end{array}$$

Blocks by definition are persistent and quasi-stationary high-pressure systems that divert or severely slow down the prevailing westerly winds in the mid-latitudes. They influence regional temperature and precipitation patterns for an extended period. Therefore, not all blocks that fulfil the two above-mentioned conditions are retained. We only include blocks that have a minimum lifetime of 5 d and a minimum overlap of the blocked area of 70 % ($A_{t+\mathrm{1}}\cap A_t>\mathrm{0.7}\cdot A_t)$ in our blocking catalogue. This largely follows the criteria defined by Schwierz et al. (2004).

2.2 Climate reconstructions

To be as independent as possible with regard to the reanalyses we use a wide array of climate index reconstructions for the 20th century:

• Atlantic Multidecadal Oscillation (AMO): for the AMO index we take October values based on the Enfield et al. (2003) study. We choose October to allow for a certain feedback lag with the atmosphere and to have decent prediction value for the upcoming snow and NAO indices.

• El Niño–Southern Oscillation (ENSO): we choose the ENSO3.4 reconstruction based on the HadISSTv1 Rayner et al. (2003) sea surface temperatures (SSTs). As with the AMO, we select October values to allow for a reaction time in the teleconnections.

• NAO: we use the extended Jones et al. (1997) NAO index for DJF from the Climate Research Unit (CRU).

• Sea ice: we use the monthly sea ice reconstruction by Walsh et al. (2017), which covers the period 1850–2013, to create a Barents–Kara (65–85^∘ N, 30–90^∘ E) sea ice index for November.

We checked for autocorrelation in the time series of the snow indices, stratospheric index, Barents–Kara Sea (BKS) sea ice index (Supplement Fig. S5), AMO index and ENSO index and only found significant autocorrelation in the BKS sea ice and AMO time series. We assess the significance of a regression coefficient in a regression model by dividing the estimated coefficient over the standard deviation of this estimate. For statistical significance we expect the absolute value of the t ratio to be greater than 2 or the p value to be less than the significance level (α=0.05). The degrees of freedom are determined as (n-k), where k is the parameters of the estimated model and n the number of observations.

3 Results

3.1 Interannual links

In the following paragraphs, we investigate the year-to-year relationship between the snow indices and the following winter sea level pressure (SLP) fields. For this we use MERRA2 for a 35-year-long period ranging from 1981 to 2015, ERA20C for a 110-year-long window ranging from 1901 to 2010 and 20CRv2c for a 160-year-long window ranging from 1851 to 2010.

Figure 2 shows the linear regression fields of DJF SLP anomalies projected onto the respective snow indices in October and November. For October, we find no NAO-like pressure anomaly appears to be significantly correlated with the snow index in each of the three reanalysis products and respective time windows (Fig. 2a, b, c). Instead, negative SLP anomalies dominate northern Eurasia in MERRA2, with high-pressure anomalies towards the Himalayan Plateau. The 110-year-long regression in ERA20C shows significant negative anomalies over the Asian part of Russia, reaching as far south as Beijing. A second significant negative SLP pattern appears along the Pacific coast of Canada. Finally, SLP anomalies in 20CRv2c support the main SLP patterns shown by ERA20C but reduce the extent of negative anomalies over Eurasia and increase the extent of the negative anomalies over the North Pacific.

The DJF SLP anomaly patterns change substantially when projected onto the November snow index (Fig. 2d, e, f). All three reanalysis products show negative NAO-like pressure anomalies with significantly positive anomalies over Iceland and the northern North Atlantic and significantly negative anomalies south of ca. 45^∘ N, including Portugal and the Azores. As expected, MERRA2 shows the strongest anomalies due to the shorter regression period; however, interestingly, ERA20C, with the 110-year-long analysis period, shows less large-scale significance for positive anomalies in high latitudes compared to the 150-year-long investigation period in 20CRv2c (even though non-significant anomalies cover roughly the same area as in 20CRv2c (not shown)). This hints at decadal variations in the strength of the regression but could also be due to biases in the reanalyses.

To check for such biases we compared all reanalyses with the SLP reconstruction dataset HadSLP2r (Allen and Ansell, 2006) and found that for the regression analysis using the time period 1901–2010, 20CRv2c overestimates the polar sea level pressure response, whereas ERA20C is much closer to HadSLP2r (see Supplement Fig. S6). This would indeed support the notion of decadal variations in the strength of the relationship between predictor and predictand. However, it is worth highlighting that this overestimation for 20CRv2c is not visible for the 1851–2010 period, during which regression anomalies resemble HadSLP2r much more closely.

Figure 2 DJF sea level pressure (pascal per 1 SD) anomalies projected onto snow indices (see Fig. 1) for October (a, b, c) and November (d, e, f) for (a) and (d) MERRA2 covering 1981–2015, (b, e) ERA20C covering 1901–2010 and (c, f) 20CRv2c covering 1851–2010. Only anomalies >95 % significance level are shown.

We investigate other possible predictors for wintertime NAO via regressed anomalies onto the November BKS ice concentration, November–December mean polar GPH at 10 hPa, October AMO and October ENSO indices (Fig. 3). The periods for MERRA2 and ERA20C are identical in Fig. 2, whereas the anomaly plots for 20CRv2c use the maximum period covered in the reconstructions, namely 1851–2010 in the sea ice reconstruction, 1856–2010 in the AMO reconstruction, 1901–2010 for the polar 10 hPa GPH index taken from ERA20C and 1870–2010 for the ENSO reconstruction.

As can be seen from Fig. 3, the 35-year-long analysis in MERRA2 shows November sea ice concentration and early winter stratospheric heights regress to a similar SLP pattern as the November snow index. Positive SLP anomalies over Iceland and Greenland combined with negative anomalies over southern Europe and the adjacent North Atlantic form a negative NAO-like pattern in DJF (Fig. 3a). On the other hand, the interannual signals in the October AMO and ENSO indices do not point towards such a pressure distribution. The small interannual changes and low frequency of the AMO combined with the short sample period inhibit most of the significance, only southern Eurasia shows regions with elevated SLP. As expected, anomalies regressed on the ENSO index show significance mostly for the North Pacific and North American regions.

Looking at the regression patterns in the centennial reanalyses, the NAO-like pattern in the SLP anomalies regressed onto sea ice and stratospheric GPH can still be seen; however, the extent and strength are substantially reduced compared to MERRA2 as well as compared to the regression using November snow as a predictor. Again, ERA20C shows a decrease in the significant anomalies regressed onto sea ice compared to 20CRv2c, with possible reasons already discussed above. Elevated geopotential heights at 10 hPa consistently increase polar sea level pressure in the following winter months; however, the impact over the European and North Atlantic domain severely decreases in the centennial reanalyses.

SLP anomalies regressed onto the AMO index show significant positive SLP regions for large parts of Eurasia as well as positive anomalies over the North Atlantic west of Great Britain. Interesting to note in 20CRv2c is the very strong high-pressure anomaly reaching from the BKS to the southern part of the Ural mountains, a prominent feature often found for years with positive AMO and negative sea ice concentration and frequently linked to a high frequency of Ural blockings (UBs). SLP distribution after El Niño events does not change considerably, irrespective of the dataset and time period used. A strong Pacific signal shows the northern part of the Pacific North American pattern (PNA), with negative SLP anomalies over the eastern North Pacific. Given the autocorrelation in the AMO and BKS sea ice index, the significance in Fig. 2a, b and c as well as Fig. 2g, h and i might be severely lower due to the reduced number of degrees of freedom.

Figure 3 DJF sea level pressure (Pa/SD) anomalies projected onto BKS ice concentration in November (a, b, c), polar 10 hPa GPH November–December mean (d, e, f), October AMO (g, h, i) and October ENSO indices (j, k, l) for (a, d, g, j) MERRA2 covering 1981–2015, (b, e, h, k) ERA20C covering 1901–2010 and (c, f, i, l) 20CRv2c covering 1851–2010. Regression values for BKS ice concentrations were multiplied by −1 to aid comparability. Only anomalies >95 % significance level are shown.

To investigate the vertical development of climate anomalies connected with the November snow dipole, Fig. 4 shows the zonal mean anomalies of zonal wind and temperature in ERA20C projected onto the ERA20C November snow index (for an evaluation with an upper-air climate reconstruction, see Supplement Fig. S7). The temporal evolution of the anomalies ranging from October to February shows that stratospheric warming occurs simultaneously within the same month as a positive snow cover dipole, with no stratospheric warming leading that development. Instead, significant lower-troposphere warming is shown between 60 and 90^∘ N for October. The warming signal then dominates the stratosphere and upper troposphere in December, after which the strongest anomalies subside into the lower stratosphere and tropopause in January and February. This development of atmospheric temperatures is mirrored in the evolution of the polar vortex, where a reduction in the polar vortex and strengthening of the subtropical jet is seen together with the emergence of the November snow dipole, after which the region of strongest anomalies migrates from the stratosphere to the upper troposphere.

Figure 4 Zonal mean (180^∘ E–180^∘ W, 15–90^∘ N) (left) temperature anomalies and (right) zonal mean zonal wind anomalies projected onto snow indices in November for ERA20C covering 1901–2010. Shading indicates 95 % significance level.

To address the physical reasons as to how the low sea ice and high snow indices are connected, climate anomalies are regressed onto BKS ice concentrations for November (Fig. 5). Compared to factors such as AMO and ENSO, BKS sea ice shows a distinct snow cover dipole coinciding with a high-pressure anomaly over the BKS and the northern Ural mountains, which supports a regional atmospheric blocking and cold air advection on its eastern flank. This cold air anomaly supports increased snow cover over eastern Eurasia, while relatively warm temperatures reduce the snow cover over eastern Europe. It should be noted that October BKS ice concentration shows qualitatively the same pattern for November snow cover anomalies (not shown); however, this is not statistically significant.

Figure 5 20CRv2c November anomalies projected onto BKS ice concentration in November covering 1851–2010. Regression values for BKS ice concentrations were multiplied by −1 to aid comparability. (a) November snow cover (%/SD) anomalies projected onto BKS ice concentration in November, (b) November SLP (Pa/SD) anomalies projected onto BKS ice concentration in November, (c) November atmospheric blocking (blocking per season/SD) anomalies projected onto BKS ice concentration in November and (d) November 2 m temperature (K/SD) anomalies projected onto BKS ice concentration in November. Only anomalies >95 % significance level are shown.

3.2 Interdecadal links

The interdecadal evolution of the November snow index is shown in Fig. 6. The 21-year running means of the normalized time series of AMO, BKS ice and snow hint at a multidecadal frequency, similar in wave length to the AMO and BKS ice anomalies. Even though we refrain from correlating these time series due to the 21-year filter (Trenary and DelSole, 2016), we find the possible mechanism behind the decadal co-occurrence of warm North Atlantic SSTs, reduced sea ice and increased snow cover gradient to be physically plausible (Luo et al., 2017). As Luo et al. (2017) point out, warm North Atlantic water reduces the BKS ice concentration, which decreases the meridional temperature gradient and strong westerly winds, which in turn supports high pressure over the Ural mountains and, with that, cold air advection towards eastern Eurasia. It should be noted, however, that the AMO and the November snow index are out of phase between 1880 and 1920, where uncertainties in both products are largest.

The more critical question is the interdecadal evolution of the relationship between the predictor and the predictand. Similar to Peings et al. (2013) and Douville et al. (2017), we apply a 21-year running correlation covering the period 1901–2010 to examine the stationarity of the relationship and differences between 20CRv2c and ERA20C.

Figure 7 summarizes the correlation over time for multiple pairs of climate variables. As Fig. 7b points out, the sign of the November snow to winter NAO relationship in 20CRv2c is negative throughout the whole 20th century. Periods with negative correlations can be found at the beginning and the end of the century, with relatively weak correlation during the 1930s and 1970s. The periods of strong negative correlations overlap with commonly known Arctic warming periods, the early 20th-century Arctic warming (ETCAW) and the ongoing recent Arctic warming in the context of anthropogenic global warming. Even stronger decadal variability can be seen for the running correlations between the October snow index and winter NAO-like signal (Fig. 7a), with periods of pronounced negative correlations during the early 20th-century Arctic warming and the 1980s. Emerging since the 1970s is a negative relationship shown in Fig. 7e between BKS ice reduction (multiplied by −1 to aid comparability) and the formation of a negative NAO signal in the following winter, with very weak negative correlations for the ETCAW.

Figure 6 The 21-year running means of (a) November snow index from 20CRv2c, (b) November BKS ice concentration and (c) October AMO.

Together with the emergence of the sea-ice-to-NAO relationship, negative correlations between BKS sea ice and November snow index (Fig. 7d) as well as between stratospheric warming and winter NAO strengthen towards the end of the 20th century (Fig. 7f). This strengthening is also found in ERA20C for the correlation between November snow and a following stratospheric warming, wheres 20CRv2c shows consistently positive correlation values throughout the 20th century (Fig. 7c).

Overall, the 20CRv2c November snow index shows a more stationary relationship with tropospheric and stratospheric winter circulation than ERA20C. Possible explanations for this behaviour will be discussed in the following section.

For all of the linear relationships shown in Fig. 7 we performed a Durbin–Watson test to check for serial correlation between two variables and did not find any compelling indication of co-dependence in any case (see Supplement Table S1). Moreover, we investigated different running correlation windows (11, 21, 25 and 31 years) and find that the main outcome of the analysis is not dependent on the choice of the correlation window (see Supplement Fig. S8).

Figure 7 The 21-year centred running correlation time series between (a) October snow index and DJF NAO, (b) November snow index and DJF NAO, (c) November snow index and mean November–December polar 10 hPa GPH index, (d) November snow index and November BKS ice concentration, (e) November BKS ice concentration multiplied by −1 to aid comparability and DJF NAO, and (f) mean November–December polar 10 hPa GPH and DJF NAO index. Black dashed line indicating the 95 % confidence level for a two-sided Student's t test assuming independence and normal distribution.

Based on the results from Fig. 7 (and the overall significance of linear relationships; see Supplement Fig. S9) we investigate very basic linear multiple and simple regression models to predict the upcoming DJF NAO index sign and assess the contributions to the prediction skill by November sea ice, November snow cover and November–December mean stratospheric conditions. For the period 1901–2010 we investigate three different multiple regression models with

• a.

  $\text{DJF NAO}\left(t\right)=a\mathrm{1}\times \text{Nov. snow cover}\left(t\right)+b\mathrm{1}\times \text{Nov. BKS sea ice}\left(t\right)+c\mathrm{1}\times \text{ND}\mathrm{10}\text{hPa GPH}\left(t\right)$,

• b.

  $\text{DJF NAO}\left(t\right)=a\mathrm{1}\times \text{Nov. snow cover}\left(t\right)+b\mathrm{1}\times \text{Nov. BKS sea ice}\left(t\right)$,

• c.

  $\text{DJF NAO}\left(t\right)=a\mathrm{1}\times \text{Nov. snow cover}\left(t\right)+b\mathrm{1}\times \text{ND}\mathrm{10}\text{hPa GPH}\left(t\right)$ and one simple linear regression model,

• d.

  $\text{DJF NAO}\left(t\right)=a\mathrm{1}\times \text{Nov. snow cover}\left(t\right)$,

where DJF NAO is the standardized NAO index calculated by EOF analysis of 20CRv2c SLP data, Nov. snow cover is the November 20CRv2c snow cover index, Nov. BKS sea ice is the Walsh et al. (2017) November BKS sea ice index and ND 10 hPa GPH is the ERA20C November–December mean 10 hPa GPH index, with a1, b1 and c1 being the constants determined by the least-squares calculations. Moreover, we perform (b) and (d) also for the period 1851–2010.

Figure 8 shows original and predicted normalized DJF NAO values together with the 21-year running correlation of both indices. Overall correlation values are low but significant for the 110-year time period (ranging from 0.41 to 0.38), but specific periods of high correlation emerge for both Arctic warm periods, the first one being centred around 1925 and the second one being centred around the year 2000 with both periods reaching correlation coefficients above 0.6. The multiple regression prediction model with three different predictors performs best, with a significant correlation to the original NAO variability of 0.41 for 110 years (Fig. 8a). Nevertheless, November snow cover seems to add most of the prediction skill, since the decrease in correlation coefficient between the multiple regression model with three predictors and the simple linear regression model with just November snow cover as a predictor is 0.03. Moreover, periods of high correlation coefficients align with periods of strong negative relationships in Fig. 7b.

For the same empirical prediction model using 160 years, the overall correlation coefficients decrease to around 0.3. As expected, the same periods of increased prediction skill emerge (Fig. 8e and f) and the added prediction skill of sea ice is low. It should be noted, however, that sea ice increases prediction skill during the current Arctic warming period, as well as the end of the 19th century, with the second-highest correlation coefficients centred around 1890 (not shown).

Figure 8 Comparison of 1901–2010 20CRv2c DJF standardized NAO values based on EOF analysis with predicted values from multiple and simple linear regression models showing (a) multiple linear regression model with November snow cover index, November BKS sea ice index and ND 10 hPa geopotential height index with an overall correlation of 0.41, (b) multiple linear regression model with November snow cover index and ND 10 hPa geopotential height index with an overall correlation of 0.4, (c) multiple linear regression model with November snow cover index and November BKS sea ice index with an overall correlation of 0.39, and (d) simple linear regression model with November snow cover index and November BKS sea ice index with an overall correlation of 0.38. Panels (e) and (f) same as (c) and (d) but for the period 1851–2010. Left y axis indicates standard deviation; right y axis indicates correlation coefficient. Red dashed line indicates 95 % significance level for a 21-year period.

4 Discussion

We used a variety of reanalyses and reconstructions to address some of the open questions regarding the relationship between Eurasian snow cover and the state of the NAO in the following winter.

Given the highly discussed research topic of Northern Hemisphere sea ice cover and snow cover impact on mid-latitude circulation (Cohen et al., 2020), as well as the highlighted need to investigate relationships over several decades (Kolstad and Screen, 2019), we investigated a promising November west–east snow cover dipole over Eurasia (Gastineau et al., 2017; Han and Sun, 2018) and its relationship to the DJF NAO state up to the middle of the 19th century to cover 150 years of internal and external climate forcings. Given the importance for seasonal prediction, we addressed the question of the stationarity of said relationship as well as its context within other common Northern Hemispheric predictors.

Compared to Gastineau et al. (2017) and Han and Sun (2018), we could extend the reanalysis study period from 35 to 150 years and highlighted the consistently negative sign of the snow–NAO relationship in the 20CRv2c dataset. Partial correlations for 110 years show that reduced BKS sea ice shows a similar response in DJF SLP anomalies; however, its statistical importance, and therefore its quality as being the prime predictor, is less than the November snow index (see Supplement Table S2 for partial correlations). This is also found in simple multiple regression prediction models, whereas the November snow cover index incorporated the major share of the prediction power. Extending the analysis of Gastineau et al. (2017) to 150 years further underlines the lack of snow–atmosphere feedback in most of the CMIP5 models and reduces the probability that the snow–NAO link is due to random internal variability at the end of the 20th century.

Moreover, the monthly evolution of vertical temperature anomalies related to a high snow cover supports the theoretical framework (Cohen et al., 2014; Henderson et al., 2018) for a Eurasian snow-cover-to-stratosphere link in reanalyses for at least the 20th century and probably before. We found surface cooling and snow cover expansion east of the sea ice anomaly, where cold air is advected on the eastern side of a Ural blocking anomaly (Fig. 5). The increased geopotential heights and the related Rossby wave energy reach the stratosphere (Supplement Fig. S7), where a stratospheric warming and a slowdown of the polar vortex manifests itself (Fig. 4). These anomalies reach the troposphere in January and February, where they express themselves as a negative NAO signal (Fig. 2). It is noteworthy that all of these features are significantly correlated with the November snow cover index for more than 100 years.

Peings et al. (2013) and the follow-up study by Douville et al. (2017) found that the October and October–November mean snow cover over a broader region of northern Eurasia and its relationship to the wintertime NAO is indeed not stationary over time. We found a strong relationship between the reduced variance of the snow index time series with the reduction in correlation strength of snow cover and the wintertime NAO (Fig. 9). The reduction in variance is even stronger in ERA20C than in 20CRv2c, which would explain the less stationary correlations in ERA20C. Furthermore, such periods of low snow variability coincide with a reduction in polar vortex variability, hinting even more so at possible links between November snow and stratospheric temperatures in the following month. Together with the snow cover index, the November BKS sea ice index shows increased variability with strengthened negative correlation to DJF NAO at the end of the 20th century (see Supplement Fig. S11).

These periods of increased variability in the November snow cover index co-occur arguably with the common Arctic warming periods of the 20th century, the ETCAW (Wegmann et al., 2016; Hegerl et al., 2018) and the recent ongoing Arctic warming, with peak variance and correlation values centred around the years 1920 and 2000. Interestingly, October snow cover index and BKS sea ice index variability peaks slightly after the ETCAW around the year 1945. Analysing temperature anomalies (not shown) for all three periods reveals more continental warming over Russia for the period 1911–1930, whereas warming between 1936 and 1955 is located very much on the Kara Sea coast of Russia. Both the October snow index and the BKS sea ice index are thus impacted by the locally increased near-surface temperatures during the latter period. Generally, Arctic warming periods appear to increase the variability of cryospheric predictors considerably and thus strengthen their value in seasonal prediction frameworks. Given the importance of stratospheric variability for seasonal prediction and the apparent relationship between snow cover variability and stratospheric variability (Fig. 9), it can be expected that the cryosphere–stratosphere pathway is also considerably stronger in Arctic warm periods than for cold periods. Moreover, in our statistical analysis, we found no indication of a stratospheric precursor of November snow cover anomalies.

Figure 9 The 21-year running standard deviation time series of (a) October snow index and (b) November snow index in ERA20C and 20CRv2c (snow cover and snow depth). Dashed black line shows running standard deviation of 10 hPa November–December mean GPH over the polar regions.

In accordance with the shorter time frame analysis of Sun et al. (2019), decadal variability of the November snow cover index seems mostly dominated by low-frequency variability in the AMO and subsequently reduced or increased polar sea ice concentration. This mechanism is also supported by the results of Luo et al. (2017), who highlighted the decadal relationship between a positive AMO, reduced sea ice and increased Ural blocking for the second half of the 20th century. Looking at this mechanism on an interannual basis, we showed a robust strengthening of the November snow dipole with decreasing BKS ice concentration, circulation changes over the BKS region and consequently cold air advection towards the eastern part of the snow dipole region for a period of 150 years. With this, our results support recent studies, which point to the counter-intuitive mechanism of Arctic warming and increased continental snow cover via sea ice reduction and circulation changes (Cohen et al., 2014; Wegmann et al., 2015; Yeo et al., 2016; Gastineau et al., 2017).

Peings (2019) performed model experiments with nudged November Ural blocking fields, BKS ice and snow anomalies. The author found that UB events are not triggered by reduced sea ice, but in fact lead sea ice decrease. Moreover, more November snow alone did not lead to an increase in blocking frequency, nor to a stratospheric warming. The study highlights the UB events as a primary predictor for a negative NAO and the WACC pattern. On the other hand, Luo et al. (2019) established a causal chain via a stratospheric pathway from reduced sea ice to a reduced potential vorticity gradient and increased blocking events leading to cold extremes over Eurasia. We computed the field average of the blocking frequency within the domain of Peings (2019) (10^∘ W–80^∘ E, 45–80^∘ N) and could find a strong correlation with the WACC pattern over time, but only for DJF blocking events (not shown).

We found a correlation of November UB events with wintertime NAO, which, however, is still weaker than the relationship with the November snow dipole, as well as our BKS ice index (see Supplement Fig. S10). Moreover, blockings within the domain of Peings (2019) (10^∘ W–80^∘ E, 45–80^∘ N) are not related to a snow dipole whatsoever, neither in October nor in November (see Supplement Fig. S10). That said, we want to highlight the fact that the blocking pattern emerging in Fig. 5 is mostly outside of the boundaries of this UB index (10^∘ W–80^∘ E, 45–80^∘ N) and thus might not be caught by our study. Furthermore, Peings (2019) applies a very general snow cover increase in his nudging experiment, rather than a snow dipole with a west-to-east gradient. Finally, although we focused here on the connection to the NAO, we did not find strong significant correlations between autumn snow and winter WACC. As pointed out by Peings (2019), the most important driver for the WACC signal is the Ural blocking, for which we found strong correlations throughout the 20th century (not shown).

Overall, we advocate the importance of the signal-to-noise ratio rather than mean states for the evolution of the November snow to winter NAO relationship. In our statistical analysis, we did not find any indication of a centennial relationship between the autumn ENSO or autumn QBO sign and the variability of the relationship between November snow cover and DJF NAO (not shown). As mentioned above we found the strongest influence to be the increased variability of the system due to energy uptake.

That said, a source of uncertainty is the disagreement between ERA20C and 20CRv2c when it comes to the stationarity of the relationship. 20CRv2c shows negative correlation throughout the whole 20th century, whereas ERA20C flips the sign of the correlation in the late 1930s and late 1970s. The same relationship but using October snow shows high agreement between the two datasets, and the same is the case for the correlations between snow and stratospheric GPH. We therefore conclude that the information stored in the November snow cover in 20CRv2c is slightly different to the information stored in the ERA20C snow depth. Wegmann et al. (2017b) found that Eurasian November snow depth shows much larger disagreement between 20CRv2c and ERA20C than the same snow depth in October. In the same study, the authors found decadal trends (although linear trend subtraction for all predictor time series was done for this study) in ERA20C snow depth, which might impact the running correlations. Finally, since snow depths are relatively low in October, differences between using snow cover and snow depth might be less important from an energy transfer point of view.

The disagreement between ERA20C and 20CRv2c may also be related to uncertainties and inhomogeneities in both reanalyses. Many studies showed that both ERA20C and 20CRv2c are not suitable for trend analysis and may include radical shifts in atmospheric circulation, particularly over the Arctic (e.g. Dell'Aquila et al., 2016; Rohrer et al., 2019). However, Rohrer et al. (2019) showed that although trends in centennial reanalyses may be spurious, at least in the Northern Hemisphere year-to-year variability of mid-tropospheric circulation is in agreement even in the early 20th century.

5 Conclusions

Several reconstruction and reanalysis datasets were used to examine the link between autumn snow cover, ocean surface conditions and the NAO pattern in winter for the whole 20th century and into the 19th century. We found evidence for a manifestation of a negative NAO signal after November with a strong west-to-east snow cover gradient, with this relationship being significant for the last 150 years. Interdecadal variability for this relationship seems to be linked to Arctic warm periods, which increase the variability of the cryospheric predictors considerably. As a result, increased variability in the predictors helps to generate a better seasonal prediction estimation.

Furthermore, our analysis of centennial time series supports studies pointing out the link of autumn snow to stratospheric circulation as well as the co-occurrence between reduced BKS ice concentration and increased snow cover in eastern Eurasia. The latter mechanism is triggered via the development of an atmospheric high-pressure anomaly adjacent to the BKS sea ice anomaly, which transports moisture and cold air along its eastern flank into the continent. The interdecadal evolution of the November snow index also points towards co-dependence with high North Atlantic SSTs and subsequently reduced sea ice.

Extending the investigation period from 35 to 110 and up to 150 years increases the confidence in recently proposed physical mechanisms behind cryospheric drivers of atmospheric variability and decreases the probability of random co-variability between the Arctic cryosphere changes and mid-latitude climate.

For future studies regarding seasonal prediction, we emphasize the use of the November snow dipole concerning a prediction of the winter NAO state. Nevertheless, periods of weak correlation might occur again, especially since it is uncertain how the sea-ice-to-snow relationship will change with stronger anthropogenic global warming, once the Arctic is ice-free in summer or the local warming is strong enough to override the counter-intuitive snow cover increase. Thus, further studies are needed to investigate the interplay between Arctic sea ice and continental snow distribution. Future experiments should take into account year-to-year variability and realistic distribution of snow cover if links to the stratosphere are to be examined.