Intermittency of Arctic–mid-latitude teleconnections: stratospheric pathway between autumn sea ice and the winter North Atlantic Oscillation

There is an observed relationship linking Arctic sea ice conditions in autumn to mid-latitude weather the following winter. Of interest in this study is a hypothesized stratospheric pathway whereby reduced sea ice in the Barents and Kara seas enhances upward wave activity and wave-breaking in the stratosphere, leading to a weakening of the polar vortex and a transition of the North Atlantic Oscillation (NAO) to its negative phase. The Causal Effect Networks (CEN) framework is used to explore the stratospheric pathway between late autumn Barents–Kara sea ice and the February NAO, focusing on its seasonal evolution, timescale dependence, and robustness. Results indicate that the pathway is statistically detectable and has been relatively active over the 39-year observational period used here, explaining approximately 26 % of the interannual variability in the February NAO. However, a bootstrap-based resampling test reveals that the pathway is highly intermittent: the full stratospheric pathway appears in only 16 % of the sample populations derived from observations, with individual causal linkages ranging from 46 % to 84 % in occurrence rates. The pathway's intermittency is consistent with the weak signal-to-noise ratio of the atmospheric response to Arctic sea ice variability in modelling experiments and suggests that Arctic–mid-latitude teleconnections might be favoured in certain background states. On shorter timescales, the CEN detects two-way interactions between Barents–Kara sea ice and the mid-latitude circulation that indicate a role for synoptic variability associated with blocking over the Urals region and moist air intrusions from the Euro-Atlantic sector. This synoptic variability has the potential to interfere with the stratospheric pathway, thereby contributing to its intermittency. This study helps quantify the robustness of causal linkages within the stratospheric pathway, and provides insight into which linkages are most subject to sampling issues within the relatively short observational record. Overall, the results should help guide the analysis and design of ensemble modelling experiments required to improve physical understanding of Arctic–mid-latitude teleconnections.

1 Introduction

Autumn sea ice is a potential source of skill in predicting the winter North Atlantic Oscillation (NAO) and hence the European climate (Scaife et al., 2014; Wang et al., 2017; Hall et al., 2017). One proposed mechanism for the relationship focuses on the Barents and Kara seas, a region with seasonal ice cover that has exhibited strong negative trends during the cold season over the last decades (Cavalieri and Parkinson, 2012; Serreze and Stroeve, 2015; Onarheim and Årthun, 2017). According to this mechanism, reduced Barents–Kara sea ice triggers a wave response that constructively interferes with the climatological stationary wave pattern (Peings and Magnusdottir, 2014; Kim et al., 2014; Sun et al., 2015; Nakamura et al., 2016; Wu and Smith, 2016; Hoshi et al., 2017; Zhang et al., 2018a; De and Wu, 2019), enhancing upward propagation of planetary waves, which weakens the stratospheric polar vortex (Nishii et al., 2009; Garfinkel et al., 2010; Smith et al., 2010). Downward coupling from the stratosphere to the troposphere subsequently produces circulation anomalies that resemble the negative phase of the NAO or Arctic Oscillation (AO) (Baldwin and Dunkerton, 1999; Polvani and Waugh, 2004), along with its attendant climate effects (Hurrell, 1995).

A delayed stratospheric pathway linking sea ice and the NAO is suggested by observations, but its exact nature is somewhat unclear. The observational evidence (e.g. García-Serrano et al., 2015; King et al., 2016; Koenigk et al., 2016) hinges on lagged correlations, such as the one shown in Fig. 1a (similar to Fig. 10c in García-Serrano et al., 2015, and Fig. 6b in King et al., 2016): less Barents–Kara sea ice in November is associated with higher polar cap heights in the stratosphere (i.e. polar vortex weakening) and a subsequent downward propagation of the height anomalies into the troposphere through the winter season, consistent with the appearance of negative NAO conditions several months later. However, the stationarity and statistical significance of this signal have both been questioned when using longer records that extend back to before the satellite era (Hopsch et al., 2012; Kolstad and Screen, 2019). In fact, the strength and timing of the signal can change when the observational period in Fig. 1a is extended by several additional winters, showing a statistically insignificant autumn sea ice connection to the winter NAO via the stratosphere (Fig. 1b).

Figure 1 Lead–lag correlations (shading) between November Barents–Kara sea ice index (sign reversed) and polar cap height (70^∘ N poleward) over the October-to-February cold season using ERA-Interim reanalysis for two periods: (a) 1979/80–2010/11 and (b) 1979/80–2017/18. Hatching indicates non-significant values at the 5 % level using a two-tailed t test. Linear trends and the seasonal cycle have been removed.

Evidence from modelling experiments is even more difficult to interpret because the relationship between Barents–Kara sea ice and the NAO is not robust in simulations. Some studies find a clear stratospheric signal after removing sea ice, leading to a weakening of the polar vortex and a negative NAO (Kim et al., 2014; Nakamura et al., 2015; Sun et al., 2015). A negative NAO response to sea ice loss is also possible, although much weaker, if the stratospheric pathway is not well represented or artificially suppressed (Liptak and Strong, 2014; Sun et al., 2015; Wu and Smith, 2016; Nakamura et al., 2016; Zhang et al., 2018a; De and Wu, 2019). However, other modelling studies show a weak or even positive NAO response when sea ice is reduced (Singarayer et al., 2006; Strey et al., 2010; Orsolini et al., 2012; Cassano et al., 2014; Screen et al., 2014), and we lack a comprehensive understanding of why model results are so different (Screen et al., 2018). One reason may be that the atmospheric response depends on where and when sea ice is removed; for example, some studies have shown that sea ice loss in the Pacific sector leads to a strengthening of the polar vortex (Sun et al., 2015; Screen, 2017; McKenna et al., 2018) and that winter ice loss may be more influential than autumn ice loss in weakening and shifting the jet stream (Blackport and Screen, 2019). Other possible reasons include non-linearities with respect to the amplitude of sea ice loss (Petoukhov and Semenov, 2010; Semenov and Latif, 2015; Chen et al., 2016; Overland et al., 2016) and dependence of the atmospheric response on the background state (Smith et al., 2017, 2019; Labe et al., 2019).

Overall, isolating the sea ice influence on the mid-latitudes remains a challenge in part because it is a search for causal drivers in a tightly coupled system with large internal variability (Shepherd, 2016). This internal atmospheric variability itself has well-known effects on Arctic climate over a range of timescales. Synoptic weather systems carry heat and moisture poleward from the North Atlantic and are associated with moist intrusions that have been shown to warm the Arctic and melt sea ice (Woods et al., 2013; D.-S. R. Park et al., 2015; H.-S. Park et al., 2015; Gong and Luo, 2017; Kim et al., 2017; Lee et al., 2017). Feedbacks between sea ice and the NAO acting on intraseasonal timescales can yield opposite-signed relationships depending on the time lag considered: anomalously low Barents–Kara sea ice concentrations are favoured by positive NAO conditions (Fang and Wallace, 1994; Deser et al., 2000) but are also part of an ice perturbation pattern that has been found to produce negative NAO conditions (Magnusdottir et al., 2004; Deser et al., 2004; Kvamstø et al., 2004; Strong et al., 2009; Deser et al., 2010; Wu and Zhang, 2010). The causality problem with respect to sea ice extends beyond the NAO to other mid-latitude phenomena, such as Eurasian cooling, for which one finds numerous studies arguing both for (Outten and Esau, 2012; Mori et al., 2014, 2019) and against (McCusker et al., 2016; Sorokina et al., 2016; Ogawa et al., 2018; Blackport et al., 2019) sea ice loss being responsible for the recent spate of extreme winters.

In the present study, we revisit the observed relationship between autumn Barents–Kara sea ice and the winter NAO with the goal of quantifying the robustness of the stratospheric linkage. In other words, we ask how systematically the stratospheric linkage has appeared during the satellite period. While sampling issues are unavoidable when using a short observational record with large internal variability, our analysis attempts to account for this by exploring the idea that weak but statistically significant signals may arise from a teleconnection pathway that is only intermittently active.

We begin with a description of data and methods (Sect. 2), including a Causal Effect Networks (CEN) approach that provides a statistical framework for assessing causality (applied to climate problems by studies such as Ebert-Uphoff and Deng, 2012; Runge et al., 2014; Kretschmer et al., 2016, 2018). Results showing that the pathway is indeed detectable but exhibits a high level of intermittency are presented in Sect. 3, and the implications for understanding present-day Arctic–mid-latitude teleconnections are discussed in Sect. 4. We end with some concluding remarks in Sect. 5.

2 Data and methods

2.1 Reanalysis data

The Causal Effect Networks (CEN) approach requires indices (time series) of variables representing key processes in the dynamical mechanism being studied. In our study, we use sea ice area fraction, surface sensible heat flux, surface latent heat flux, sea level pressure, meridional wind, temperature, geopotential height, and downward thermal radiation at the surface. Raw daily data for the period 1979 to 2018 are from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis (Dee et al., 2011). The seasonal cycle is removed at each grid point by subtracting the climatological daily mean to obtain anomalies of each variable, and the data are then detrended. The trend is removed through all the days of the year (1 January, 2 January, etc.). The following indices are then calculated from the reanalysis data from September to March.

• Barents–Kara sea ice (ICE): sea ice area fraction averaged over 70–80^∘ N, 30–105^∘ E (Fig. 2a);

• Barents–Kara turbulent heat flux (THF): sum of surface sensible and latent heat flux averaged over 70–80^∘ N, 30–105^∘ E (Fig. 2b), with positive defined as heat flux from the ocean to the atmosphere;

• stratospheric polar vortex strength (SPV): negative of geopotential height poleward of 60^∘ N (Fig. 2c) averaged between 10 and 100 hPa, as defined by Kretschmer et al. (2016), such that positive values of the index indicate a stronger polar vortex;

• Urals sea level pressure (URALS): sea level pressure averaged over 45–70^∘ N, 40–85^∘ E (Fig. 2d);

• downward longwave radiation (IR): downward thermal radiation at the surface averaged over 70–90^∘ N (Fig. 2e);

• poleward eddy heat flux (${\mathrm{V}}^{\ast }{\mathrm{T}}^{\ast }$): product of V^∗ and T^∗ at 100 hPa averaged over 45–75^∘ N (Fig. 2f), where V and T denote the meridional wind velocity and air temperature, respectively, and the superscript ^∗ indicates deviations from the zonal mean;

• North Atlantic Oscillation index (NAO): taken from the Climate Prediction Center, based on Rotated Principal Component Analysis of 500 hPa geopotential height (see details at https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml, last access: November 2019).

Figure 2 ERA-Interim (1979–2018) DJF climatologies (shading) of key variables and regions (black boxes) for computing area-averaged indices: (a) sea ice area fraction, 70–80^∘ N, 30–105^∘ E; (b) turbulent heat flux, 70–80^∘ N, 30–105^∘ E; (c) stratospheric polar vortex, which is defined by 10-100 hPa geopotential height, 65–90^∘ N; (d) Urals sea level pressure, 45–70^∘ N, 40–85^∘ E; (e) downward longwave radiation, 70–90^∘ N; and (f) 100 hPa poleward eddy heat flux, 45–75^∘ N. For panel (b), turbulent heat flux from the ocean to the atmosphere is defined as positive. See Sect. 2.1 for details.

Finally, the daily indices are averaged up to monthly, half-monthly and pentad means for the different analyses carried out in this study.

2.2 Causal Effect Networks (CEN)

The CEN algorithm is a causal inference framework (Runge et al., 2014, 2019) aimed at identifying causal relationships between variables of interest. It was previously used to study Arctic–mid-latitude teleconnections by Kretschmer et al. (2016, 2018). Essentially, given a set of indices such as the ones described above, a CEN is constructed following three steps: (1) identifying potential causal drivers of each index (condition selection), (2) identifying the causal drivers using these potential causal drivers as a “conditioning set”, and (3) quantifying the strength of the causal relationship. We will illustrate the algorithm using January stratospheric polar vortex strength (SPV_Jan) as an example. Readers are referred to Kretschmer et al. (2018) and Runge et al. (2019) for a full description of the CEN algorithm, also known as PCMCI.

In the first step, we find all possible drivers for SPV_Jan. A preliminary list of drivers is generated by calculating the Pearson correlation r between SPV_Jan and all other indices (including SPV itself) in the preceding months, up to a maximum lag of 2 months (i.e. November and December for this example). Indices with significant correlations are retained, where an optimal significance level is determined using the Akaike information criterion (AIC). The AIC results in the selection of a 20 % significance level for the case of SPV_Jan (note that the AIC allows for these rather liberal significance levels in the first step, but more stringent levels are used later in the second step). This leaves us with the following possible drivers: ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$, URALS_Dec, SPV_Dec and URALS_Nov. This list is sorted in descending order according to the absolute value of the Pearson correlation coefficient. Next, we test for the conditional independence of all four possible drivers with SPV_Jan by calculating partial correlations, controlling for the effect of each driver one at a time, starting from the top of the sorted list. If a driver passes the partial correlation test, it is retained in the list of possible drivers; if it does not pass, it is removed from the list, meaning it is no longer in the conditioning set. For example, the partial correlation between URALS_Nov and SPV_Jan controlling for ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$ is, following the notation of Kretschmer et al. (2016):

$$\begin{array}{}\text{(1)}& \mathit{\rho }\left({\text{URALS}}_{\text{Nov}},{\text{SPV}}_{\text{Jan}}|{\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}\right)=-\mathrm{0.274},\end{array}$$

where

$$\begin{array}{}\text{(2)}& \mathit{\rho }(x,y|z)={\displaystyle \frac{r_{xy}-r_{xz}{r}_{yz}}{\sqrt{\mathrm{1}-r_{xz}^{\mathrm{2}}}\sqrt{\mathrm{1}-r_{yz}^{\mathrm{2}}}}}.\end{array}$$

The partial correlation is significant at the 20 % level (p value =0.105); therefore, URALS_Nov is retained as a possible driver of SPV_Jan. After going through the entire list, SPV_Dec is eliminated, leaving us with three possible drivers of SPV_Jan: URALS_Nov, URALS_Dec, and ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$.

In the second step, we retest all possible links (for all indices in the preceding 2 months, including those rejected in the first step) with SPV_Jan, controlling for the combined effect of the possible drivers (conditioning set) identified in the first step. This step helps account for false positives when working with highly interdependent time series (as is often the case with climate indices) and enhances detection power (Runge et al., 2019). Specifically, the test for SPV_Jan is as follows:

$$\begin{array}{}\text{(3)}& \mathit{\rho }\left(X,{\text{SPV}}_{\text{Jan}}|{\text{URALS}}_{\text{Nov}},{\text{URALS}}_{\text{Dec}},{\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}\right),\end{array}$$

where X represents all indices of ICE, THF, URALS, ${\mathrm{V}}^{\ast }{\mathrm{T}}^{\ast }$, SPV, and NAO in both November and December. Any X producing a significant partial correlation in Eq. (3) is regarded as a causal driver of SPV_Jan. The conditioning set excludes X when X is being tested, for example:

$$\begin{array}{}\text{(4)}& \mathit{\rho }\left({\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}},{\text{SPV}}_{\text{Jan}}|{\text{URALS}}_{\text{Nov}},{\text{URALS}}_{\text{Dec}}\right)=-\mathrm{0.453},\end{array}$$

which is significant at the 5 % level (p value =0.00629). Testing all of X leaves us with three causal drivers of SPV_Jan: URALS_Nov, URALS_Dec and ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$. Note that these are the same causal drivers identified in the first step, meaning that no new drivers are reintroduced in the second step in this case. As an additional refinement, the Hochberg–Benjamini false discovery rate (FDR) control may be used to account for the multiple testing problem (Kretschmer et al., 2018; Runge et al., 2019).

In the third step, we use a multiple regression equation to quantify the influence of causal drivers and simultaneous influences on SPV_Jan:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\text{SPV}}_{\text{Jan}}^{\mathrm{s}}& ={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}\times {\text{URALS}}_{\text{Nov}}^{\mathrm{s}}+{\mathit{\beta }}_{\mathrm{2}}\times {\text{URALS}}_{\text{Dec}}^{\mathrm{s}}\\ & +{\mathit{\beta }}_{\mathrm{3}}\times {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}^{\mathrm{s}}+{\mathit{\beta }}_{\mathrm{4}}\times Y_{\text{Jan}}^{\mathrm{s}},\end{array}\end{array}$$

where the β values are regression coefficients for the standardized regressors ${\text{URALS}}_{\text{Nov}}^{\mathrm{s}}$, ${\text{URALS}}_{\text{Dec}}^{\mathrm{s}}$, ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}^{\mathrm{s}}$, $Y_{\text{Jan}}^{\mathrm{s}}$, and the superscript s indicates a standardized index. The inclusion of Y allows us to check for significant simultaneous relationships between all indices. By standardizing, the interpretation is that changing a certain regressor by one standard deviation changes SPV_Jan by β standard deviations, provided that all other variables are held fixed.

A two-tailed t test is used for significance testing. For the AIC in step one, a significance set of (5 %, 10 %, 20 %) is used. There are no substantial changes to the main messages when using other significance sets (Fig. S4 in the Supplement). A significance level of 5 % is used in the second and third steps.

The above example illustrates how the CEN algorithm identifies and evaluates causal drivers of SPV_Jan. In order to construct the complete monthly and half-monthly CENs, we identify causal drivers for all our chosen indices (ICE, THF, URALS, ${\mathrm{V}}^{\ast }{\mathrm{T}}^{\ast }$, SPV, and NAO) during the extended winter season (NDJFM). September to December (January to March) indices are taken over the period of 1979 to 2017 (1980 to 2018). All Pearson correlations and partial correlations (first and second steps) and the multiple regressions (third step) are thus based on indices with a sample size of 39 winter seasons. A similar procedure is used for the pentad CEN but with a maximum lag of two pentads to capture processes occurring on synoptic timescales.

3 Results

This section describes results from our exploration of the ICE–NAO stratospheric pathway using the CEN framework (Sect. 3.1), including an assessment of its strength (Sect. 3.2) and intermittency (Sect. 3.3) in the observational record. We also explore processes occurring on shorter timescales and discuss how these effects may reinforce or interrupt the stratospheric pathway (Sect. 3.4).

3.1 Seasonally evolving ICE–NAO pathway

We begin by examining pathways from Barents–Kara sea ice to the NAO proposed by previous studies. The CEN analysis follows the approach of Kretschmer et al. (2016) but keeps individual months separate rather than considering the DJF period as a whole. This allows us to capture the seasonal evolution of pathways through the cold season.

The CEN (Fig. 3) shows evidence for a stratospheric pathway leading from autumn sea ice perturbations in the Barents and Kara seas to a late winter NAO response. This pathway appears using both monthly (Fig. 3a) and half-monthly (Fig. 3b) averages as input to the CEN, albeit with slight differences in timing. The half-monthly CEN in Fig. 3b is displayed such that individual half-monthly linkages (shown in Fig. S2) are aggregated into full months to allow for direct comparison to Fig. 3a.

Figure 3 Seasonal evolution of the stratospheric pathway (indicated by coloured arrows) detected by the (a) monthly and (b) half-monthly CENs. Arrows indicate causal linkages, vertical lines indicate auto-correlation, horizontal bars indicate simultaneous relationships, and colours show the sign and strength of the linkages as given by the CEN beta coefficients (see Sect. 2.2). The grey background shows other significant linkages (arrows) and autocorrelations (vertical lines) but does not include simultaneous relationships. The half-monthly CEN in panel (b) has been aggregated into full months for ease of comparison with panel (a). See Fig. S2 for the unaggregated version.

Coloured arrows in the network diagrams highlight the ICE–NAO stratospheric pathway, where red indicates positive relationships and blue indicates negative relationships (the exact values correspond to the beta coefficients in the multiple regression equation, e.g. Eq. 5). Grey arrows show other linkages that are statistically significant, including some tropospheric pathways that also contribute to the ICE–NAO relationship. A figure including all the identified causal linkages and autocorrelations appears in the Supplement (Fig. S1). For the monthly CEN, the stratospheric pathway is

$$\begin{array}{rl}\downarrow {\text{ICE}}_{\text{Oct}}& \Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec/Jan}}\Rightarrow \downarrow {\text{SPV}}_{\text{Jan/Feb}}\\ & \Rightarrow \downarrow {\text{NAO}}_{\text{Feb/Mar}},\end{array}$$

where we use the notation A⇒B to indicate index A as a “driver” of index B and ↓ and ↑ to represent a decrease or increase, respectively, of the indices. The pathway is described for the case of a negative sea ice perturbation leading to a negative NAO. For the half-monthly CEN, the pathway may be summarized as follows:

$$\begin{array}{rl}\downarrow {\text{ICE}}_{\text{Oct/Nov}}& \Rightarrow \uparrow {\text{THF}}_{\text{Nov}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec/Jan}}\\ & \Rightarrow \downarrow {\text{SPV}}_{\text{Dec/Jan/Feb}}\Rightarrow \downarrow {\text{NAO}}_{\text{Feb/Mar}}.\end{array}$$

Using the finer half-monthly resolution in the CEN prevents shorter timescale processes (such as linkages through THF) from being averaged out.

The CEN results illustrate how the stratospheric pathway unfolds through the winter season. The timing is in general agreement with previous observational studies, suggesting that the involvement of the stratosphere introduces a few months' delay in the NAO response to Barents–Kara sea ice variability (Kim et al., 2014; García-Serrano et al., 2015; Jaiser et al., 2016; King et al., 2016; Kretschmer et al., 2016; Yang et al., 2016). The causal linkages are consistent with the idea that Arctic sea ice reduction enhances upward wave activity through constructive interference between forced Rossby waves and the climatological stationary waves (Garfinkel et al., 2010; Smith et al., 2010). The resulting increase in wave-breaking in the stratosphere decelerates the polar vortex (Charney and Drazin, 1961), which in turn leads to tropospheric circulation anomalies and surface impacts via downward coupling (Baldwin and Dunkerton, 1999). Some features of the pathway, such as the relatively long lagged relationship of autumn sea ice to December Urals sea level pressure, are not well understood, an issue that will be further discussed in Sect. 4.

We will focus on the stratospheric pathway from ICE_Oct to NAO_Feb in the monthly CEN, as this timing yields the strongest negative ICE–NAO correlation (Fig. S5). The correlation between ICE_Nov and NAO_Jan is equally strong, but the causal pathway goes through the troposphere only (Fig. S1b) and is not a focus of this study. Results from the half-monthly CEN yield consistent messages and will be brought into the discussion where relevant.

3.2 Strength of the pathway

An interesting question is how to assess the strength of the ICE–NAO stratospheric pathway as a whole, and what insights may be gained by such an assessment.

The CEN analysis yields a set of beta coefficients (colours of the arrows in Fig. 3) that describe the strength of individual causal linkages in our network. Following Runge et al. (2015), the total causal effect of the stratospheric pathway from ICE_Oct to NAO_Feb may be calculated by summing over the product of beta coefficients along the two relevant chains of linkages from Fig. 3a.

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}\downarrow {\text{ICE}}_{\text{Oct}}& \stackrel{-\mathrm{0.326}}{\Rightarrow }\uparrow {\text{URALS}}_{\text{Dec}}\stackrel{\mathrm{0.390}}{\Rightarrow }\uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}\stackrel{-\mathrm{0.368}}{\Rightarrow }\downarrow {\text{SPV}}_{\text{Jan}}\\ & \stackrel{\mathrm{0.426}}{\Rightarrow }\downarrow {\text{NAO}}_{\text{Feb}}\left(\mathrm{0.0199}\right)\end{array}\\ {\displaystyle }\begin{array}{rl}\downarrow {\text{ICE}}_{\text{Oct}}& \stackrel{-\mathrm{0.326}}{\Rightarrow }\uparrow {\text{URALS}}_{\text{Dec}}\stackrel{-\mathrm{0.449}}{\Rightarrow }\downarrow {\text{SPV}}_{\text{Jan}}\\ & \stackrel{\mathrm{0.426}}{\Rightarrow }\downarrow {\text{NAO}}_{\text{Feb}}\left(\mathrm{0.0624}\right)\end{array}\end{array}$$

The total causal effect ($\mathrm{0.0199}+\mathrm{0.0624}=\mathrm{0.0823}$) tells us that a single standard deviation perturbation in ICE_Oct yields a like-signed response of 8 % of a standard deviation in February NAO (Runge et al., 2015). One might question the interpretation of the contemporaneous $\uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}$ linkage in the first chain as a causal effect, but the fact that it also shows up in the half-monthly CEN as a linkage from the first half of the December to the second half of December (Fig. S2) supports the point.

A comparison between the stratospheric and tropospheric ICE–NAO pathways shows that the latter are generally stronger in the CEN framework. Table 1 summarizes the causal effect of the three full pathways (Fig. S6). Our main stratospheric pathway of interest from ICE_Oct to NAO_Feb is comparable in strength to the pathway from ICE_Oct to NAO_Mar (0.0823 and 0.0872). The latter has both stratospheric and tropospheric chains, accounting for 30 % ($\mathrm{0.0258}/(\mathrm{0.0614}+\mathrm{0.0258})$) and 70 % ($\mathrm{0.0614}/(\mathrm{0.0614}+\mathrm{0.0258})$) of the total causal effect, respectively. The $\downarrow {\text{ICE}}_{\text{Jan}}\Rightarrow \downarrow {\text{NAO}}_{\text{Mar}}$ tropospheric pathway is the strongest in terms of the total causal effect (0.137), primarily because it involves fewer linkages. Overall, the larger causal effect of the tropospheric pathways is perhaps unsurprising, given that the stratospheric pathway may be disrupted by internal variability (noise) from both the troposphere and the stratosphere.

Table 1 A summary of the casual effect of all ICE–NAO pathways. The $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \downarrow {\text{NAO}}_{\text{Mar}}$ pathway consists of both tropospheric and stratospheric branches; n/a stands for “not applicable”.

An alternative view of the pathway strength comes from considering the amount of February NAO variance explained by the various linkages along the pathway using a multiple linear regression framework. This gives a sense of the relative importance of each linkage and how information passes through the pathway. The full pathway can be represented by the following regression equation:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\text{NAO}}_{\text{Feb}}& ={\mathit{\kappa }}_{\mathrm{0}}+{\mathit{\kappa }}_{\mathrm{1}}\cdot {\text{ICE}}_{\text{Oct}}+{\mathit{\kappa }}_{\mathrm{2}}\cdot {\text{URALS}}_{\text{Dec}}\\ & +{\mathit{\kappa }}_{\mathrm{3}}\cdot {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}+{\mathit{\kappa }}_{\mathrm{4}}\cdot {\text{SPV}}_{\text{Jan}},\end{array}\end{array}$$

where κ₀ is a constant and κ₁, κ₂, κ₃, and κ₄ are the regression coefficients for the standardized regressors ICE_Oct, URALS_Dec, ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$ and SPV_Jan, respectively. The importance of the regressors may be quantified in the following different ways.

• a.

  cumulative NAO_Feb variance explained as regressors are included, calculated by successively adding terms in Eq. (6) from left to right (orange bars in Fig. 4), e.g. for ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$:

  $$\begin{array}{}\text{(7)}& \begin{array}{rl}{\text{NAO}}_{\text{Feb}}& ={\mathit{\kappa }}_{\mathrm{0}}^{\mathrm{a}}+{\mathit{\kappa }}_{\mathrm{1}}^{\mathrm{a}}\cdot {\text{ICE}}_{\text{Oct}}+{\mathit{\kappa }}_{\mathrm{2}}^{\mathrm{a}}\cdot {\text{URALS}}_{\text{Dec}}\\ & +{\mathit{\kappa }}_{\mathrm{3}}^{\mathrm{a}}\cdot {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}};\end{array}\end{array}$$
• b.

  NAO_Feb variance explained by individual regressors, calculated via a simple bivariate regression between each regressor and NAO_Feb (blue bars), e.g. for ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$:

  $$\begin{array}{}\text{(8)}& {\text{NAO}}_{\text{Feb}}={\mathit{\kappa }}_{\mathrm{0}}^{\mathrm{b}}+{\mathit{\kappa }}_{\mathrm{3}}^{\mathrm{b}}\cdot {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}};\end{array}$$
• c.

  reduction in NAO_Feb variance explained when individual regressors are removed, calculated by removing the term from the regression equation (green bars), e.g. for ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$:

  $$\begin{array}{}\text{(9)}& \begin{array}{rl}{\text{NAO}}_{\text{Feb}}& ={\mathit{\kappa }}_{\mathrm{0}}^{\mathrm{c}}+{\mathit{\kappa }}_{\mathrm{1}}^{\mathrm{c}}\cdot {\text{ICE}}_{\text{Oct}}+{\mathit{\kappa }}_{\mathrm{2}}^{\mathrm{c}}\cdot {\text{URALS}}_{\text{Dec}}\\ & +{\mathit{\kappa }}_{\mathrm{4}}^{\mathrm{c}}\cdot {\text{SPV}}_{\text{Jan}}.\end{array}\end{array}$$

Both the blue and green bars in Fig. 4 provide a measure of the contribution of individual regressors, while comparison of these with the orange bars gives some indication of whether information from a given regressor is redundant.

Figure 4 Explanatory power of the stratospheric pathway for the February NAO assessed via multiple linear regression. Orange bars show the cumulative variance explained when including each regressor in succession from left to right. Blue bars show variance explained by the individual regressor. Green bars show the reduction in total variance explained when removing that regressor. See Sect. 3.2 for details.

The stratospheric pathway explains 26 % of the variance in the February NAO (Fig. 4). The cumulative variance explained (orange bars) increases from 11 % to 26 % as regressors are added (moving from left to right), indicating that each linkage in the pathway adds some useful information. This result is consistent with other estimates from observations but likely represents an upper limit as the Barents–Kara sea ice and NAO relationship is shown to be particularly strong during the current reanalysis period compared to the rest of the 20th century (Kolstad and Screen, 2019).

While successive linkages in the pathway add explanatory power, they are not independent. Comparing the orange and blue bars, we see that the increase in cumulative explained variance moving from left to right is much lower than the explained variance from each individual regressor. For example, while SPV_Jan explains the most NAO variance of any individual regressor (18 %), its removal from the full regression does not have much effect (3 % reduction in explained variance). However, we know that variability in upward wave activity and variability in the polar vortex are closely related, and thus in a sense it is not physically meaningful to consider one in isolation of the other. Removing both ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$ and SPV_Jan from the regression equation results in a 8 % reduction (not shown) in explained variance, which is perhaps a more representative estimate of the stratosphere's contribution. Sea ice appears to impart information that cannot be explained by the other three regressors (6 % reduction in explained NAO variance when removed), but this may also be a result of atmosphere–ice feedbacks explored in Sect. 3.4.

Overall, these analyses show a role for the stratosphere in connecting autumn ICE to late winter NAO, but it is one that accounts for a modest fraction of the total NAO variance. However, the pathway strength reported here should be considered an estimate, as there remain uncertainties associated with analysis choices, such as the time resolution of the input data and the relevant lags to include. In the next section, we will further explore reasons for this relatively weak ICE–NAO covariability.

3.3 Intermittency of the pathway

The ICE–NAO stratospheric pathway identified by the CEN comprises statistical relationships inferred from a relatively short observational record of only 39 winters. It is meaningful to ask how robust the pathway is, i.e. how systematically the relevant statistical relationships occur in the record. To assess the robustness, we perform a bootstrapping test, where bootstrap samples are created by randomly selecting 39 winters with replacement from the entire reanalysis period. The CEN of each sample is then constructed. This procedure is repeated 10 000 times.

The bootstrapping results (Fig. 5) indicate that the stratospheric pathway is intermittent. Percentages show the occurrence rate of individual segments in the pathway within the bootstrap sample population (see Fig. S7 for occurrence rates of other statistically significant linkages). By this measure, it is clear that individual segments have varying levels of intermittency, ranging from 46 % for the segment $\downarrow {\text{SPV}}_{\text{Jan}}\Rightarrow \downarrow {\text{NAO}}_{\text{Feb}}$ to 84 % for the segment $\uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}\Rightarrow \downarrow {\text{SPV}}_{\text{Jan}}$. The full stratospheric pathway (the sequence of all four segments) is detected in only 16 % of the samples, suggesting that it does not occur systematically during every winter season. An alternative three-segment pathway $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \downarrow {\text{SPV}}_{\text{Jan}}\Rightarrow \downarrow {\text{NAO}}_{\text{Feb}}$ is slightly less intermittent (22 % occurrence rate), but its physical interpretation is unclear given that there is no linkage through ${\mathrm{V}}^{\ast }{\mathrm{T}}^{\ast }$ to the polar vortex, as expected from theory. These intermittency results are a likely reason why detection of the pathway is sensitive to the choice of observational period (Fig. 1), and suggests that it may be favoured in certain background states (Overland et al., 2016; Smith et al., 2017).

Figure 5 Results of a bootstrapping test to assess the robustness of causal linkages within the stratospheric pathway. The percentages above the arrows show the occurrence rate of each linkage out of 10 000 bootstrap samples. The colours of the arrows (identical to Fig. 3) and the black lines show observed beta coefficients for each linkage for the reanalysis period. The histograms above show the corresponding distribution of beta coefficients (absolute value) in the bootstrap samples. The histogram for the $\uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \downarrow {\text{SPV}}_{\text{Jan}}$ linkage is not shown. Note that the distributions are composed only of samples in which the linkage is detected.

The existence of intermittency in the stratospheric pathway is consistent with previous suggestions that internal variability modulates the influence of Arctic sea ice on the mid-latitude circulation (Screen et al., 2014; Overland et al., 2016; Shepherd, 2016). An examination of where in the pathway the intermittency is strongest provides clues to its origins. For example, the upward coupling from sea ice to the stratosphere includes the segments $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}$, and $\uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}$, whose occurrence rates are 50 % and 74 %, respectively. The occurrence of these two linkages together is seen in about 41 % out of 10 000 bootstrap samples, meaning that most of the time when the $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}$ linkage is detected, the subsequent linkage to ${\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\mathrm{Dec}}$ follows. Conversely, when the $\uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}$ linkage is detected, it is preceded by the $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}$ linkage in only about half the cases. An obvious source of the intermittency in both segments (individually and in terms of their “combined” occurrence rate) is regional sea level pressure (SLP) variability over the Urals related to atmospheric internal variability. Similarly, the downward coupling from SPV to NAO is vulnerable to both stratospheric and tropospheric internal variability, leading to a relatively low occurrence rate of 46 %. This is consistent with the idea that not all polar vortex strengthening and weakening events affect the tropospheric circulation (Karpechko et al., 2017). Most robust is the $\uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec}}\Rightarrow \downarrow {\text{SPV}}_{\text{Jan}}$ linkage (84 %), which arises from well-known physical processes related to upward planetary wave flux and polar vortex weakening. Sea ice variability can also contribute to intermittency in the pathway through higher-frequency synoptic processes, a topic we will explore in Sect. 3.4.

The strength of the segments in the pathway also exhibits large variability among the bootstrap samples. This can be seen in histograms of the beta coefficients for all segments in the pathway (Fig. 5). While the beta coefficients exhibit ranges of up to 0.5 for any given segment, the sign is always the same, indicating that the sign of the relationship between variables is robust. The observed beta coefficients (black lines) for the reanalysis period itself fall within the spread of the distributions. Note that the distributions are composed only of samples in which the linkage of interest is detected by the CEN algorithm (i.e. a beta coefficient can be calculated from Eq. 5), which is why some of the distributions appear skewed. This is particularly true for the linkages that are least robust (the first and last segments, for which the observed beta coefficients are towards the weaker end of the distributions). Overall, these results indicate that even when the stratospheric pathway is active, there is substantial interannual variability in how it manifests.

3.4 Synoptic linkages and interactions across timescales

In the monthly CEN analysis, there are simultaneous relationships between Barents–Kara sea ice, Urals sea level pressure, and the NAO (Fig. 6) that point to linkages through shorter timescale synoptic processes. For example, the NAO shows significant negative simultaneous relationships with Barents–Kara sea ice (positive NAO with reduced ice) in December and March, reflecting a well-known pattern of atmospheric forcing on sea ice via anomalies in surface heat fluxes driven by wind and temperature variability (Fang and Wallace, 1994; Deser et al., 2000). Additional simultaneous relationships between sea ice, turbulent heat flux, and Urals sea level pressure are consistent with synoptic features related to cyclones (Boisvert et al., 2016; Wickström et al., 2019) and moist intrusions (Woods et al., 2013; H.-S. Park et al., 2015) entering the Arctic. Moist intrusions in particular appear to occur preferentially during the positive phase of the NAO (Luo et al., 2017) and have been shown to lead to enhanced downward longwave radiation, surface warming, and sea ice reductions (Gong and Luo, 2017; Chen et al., 2018). We explore the possible influences of such events within the CEN framework by using higher-frequency data to capture the relevant synoptic processes. The input data are pentad (5 d) means of Barents–Kara sea ice (ICE), Urals sea level pressure (URALS), and downward longwave radiation (IR). The maximum lag is set to two pentads (10 d) to isolate the synoptic timescale. The results are summarized in Fig. 7 by summing the number of times a linkage appears in each month from Fig. S8. The maximum count for a given linkage in a month is 12 (6 pentads in a month and up to 2-pentad lag considered). Autocorrelation is strong on these short timescales and thus is not used to reject causal linkages in the partial correlation tests.

Figure 6 Simultaneous relationships between monthly indices in November (N), December (D), January (J), February (F), and March (M). Colours indicate the sign and strength of the relationship as given by the CEN beta coefficients (see Sect. 2.2).

Figure 7 Results of the pentad CEN analysis assessing relationships between downward longwave radiation (IR), Barents–Kara sea ice (ICE), and Urals sea level pressure (URALS) aggregated into months (October, November, December, January, February and March, from left to right). The height of each bar is the number of counts. (a) Linkages from the Arctic to the mid-latitudes. (b) Linkages from the mid-latitudes to the Arctic. Red (blue) colours denote positive (negative) relationships. See Fig. S8 for the unaggregated version.

The CEN detects synoptic-scale influences from the Arctic to the mid-latitudes that reinforce linkages found in the monthly analysis. A linkage from ICE to URALS appears regularly throughout the winter season (Fig. 7a), both indirectly through IR and as a direct connection and in the correct sense to contribute to the $\downarrow {\text{ICE}}_{\text{Oct/Nov}}\Rightarrow \uparrow {\text{URALS}}_{\mathrm{Dec}}$ linkage shown in the monthly and half-monthly CENs (Fig. 3). The $\downarrow \text{ICE}\Rightarrow \uparrow \text{IR}$ linkage (blue bars, first histogram in Fig. 7a) follows from the idea that sea ice retreat exposes open ocean, which is a local evaporative source for water vapour, leading to a moister, optically thicker atmosphere (Kim and Kim, 2017; Zhong et al., 2018). The linkage $\uparrow \text{IR}\Rightarrow \uparrow \text{URALS}$ (red bars, second histogram in Fig. 7a) is consistent with a suggested mechanism, whereby the resulting surface warming weakens zonal wind locally and promotes blocking over the Urals (Luo et al., 2016). These synoptic processes, if habitually occurring, can imprint onto longer timescales but may also produce interference effects, as seen by the appearance of opposite-signed causal relationships from those described above from time to time through the winter season.

At the same time, causal effects from the mid-latitudes to the Arctic are also detected, consistent with an influence from moisture transport by cyclones or synoptic moist intrusions (Fig. 7b). This is represented by the $\uparrow \text{URALS}\Rightarrow \uparrow \text{IR}$ linkage (most frequently observed in October, January and February) and the $\uparrow \text{IR}\Rightarrow \downarrow \text{ICE}$ linkage (most frequently observed in November, January and February), which reflect the transport of moist air into the dry Arctic atmosphere by the large-scale flow or by cyclones tracking into the Barents Sea. These mid-latitude-to-Arctic linkages have a uniform sign (all red bars in first histogram, all blue bars in second histogram), suggesting that the effect of the relevant processes is rather systematic despite exhibiting month-to-month variability. We also detect a direct linkage from the Urals to Barents–Kara sea ice that can be of either sign. In the slightly more frequent negative sense ($\uparrow \text{URALS}\Rightarrow \downarrow \text{ICE}$), it can be interpreted as a direct effect of warm air advection and mechanical forcing of the ice cover from enhanced southerlies over the Barents–Kara region (Sorokina et al., 2016; McCusker et al., 2016). Together, these synoptic linkages show how Urals SLP variability, which has a large internally generated component, can reinforce or interrupt the ICE–NAO stratospheric pathway.

Given that our understanding of Arctic–mid-latitude teleconnections must account for the combined influences of such linkages across regions and timescales, it is no surprise that we have yet to identify a definitive set of mechanisms. Implications of such scale interactions and how they relate to viewpoints presented in previous studies are further discussed in Sect. 4.

4 Discussion

This study quantifies the robustness of atmospheric teleconnections between the Arctic and mid-latitudes, documenting their high level of intermittency in the observational record. In a bootstrapping test, the full stratospheric pathway emerges in only 16 % of the sample populations derived from the observations (Fig. 5). The existence of intermittency is likely why studies using various analytical approaches and time periods find teleconnections that differ in pattern, timing, robustness, and apparent mechanisms (Overland et al., 2016; Francis, 2017; Cohen et al., 2018; Overland and Wang, 2018; Cohen et al., 2020). In this section we discuss some of the factors that may contribute to the intermittency. Of course, anything that influences polar vortex strength is a potential source of intermittency (including internal variability, anthropogenic forcing, tropical variability, etc.), but we focus the discussion on factors that are most directly related to our CEN results.

To be more concrete, the intermittency of the stratospheric pathway stems from the fact that it can be reinforced or interrupted by other processes. For example, reinforcement can come from tropospheric pathways also detected by the CEN algorithm (see Fig. S1a and b):

1. $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{THF}}_{\text{Nov}}\Rightarrow \uparrow {\text{URALS}}_{\text{Jan}}\Rightarrow \uparrow {\text{URALS}}_{\text{Feb}}\Rightarrow \downarrow {\text{NAO}}_{\text{Mar}}$;

2. $\downarrow {\text{ICE}}_{\text{Jan}}\Rightarrow \uparrow {\text{URALS}}_{\text{Feb}}\Rightarrow \downarrow {\text{NAO}}_{\text{Mar}}$;

3. $\downarrow {\text{ICE}}_{\text{Nov}}\Rightarrow \downarrow {\text{NAO}}_{\text{Jan}}$;

4. $\downarrow {\text{ICE}}_{\text{Jan}}\Rightarrow \downarrow {\text{NAO}}_{\text{Mar}}$.

All these tropospheric and stratospheric pathways lead from the reduction of sea ice to a negative NAO, although they differ slightly in timing. The existence of the tropospheric pathway is supported by sea ice and surface heating perturbation experiments, where negative NAO and AO responses are simulated even when the stratospheric pathway is suppressed (Wu and Smith, 2016) or not well represented (Sun et al., 2015). However, the NAO/AO response is stronger when the stratospheric pathway is active than when it emerges through the tropospheric pathway alone (Nakamura et al., 2016; Zhang et al., 2018a, b).

Another example of a factor that may contribute to intermittency is the El Niño–Southern Oscillation (ENSO). El Niño winters are associated with a deepened Aleutian low, which enhances upward propagating waves, weakens the polar vortex, and favours negative NAO conditions (Domeisen et al., 2019). As such, the stratospheric pathway may be reinforced if an El Niño develops following a low autumn ice season (both are associated with a weakened polar vortex, e.g. winter 1986/87 or 2009/10, Fig. 8); if a La Niña develops instead, the stratospheric pathway may be weakened (e.g. winter 2007/08 or 2010/11, Fig. 8). The relationship between wintertime ENSO and the NAO is rather weak (Brönnimann, 2007; Domeisen et al., 2019), consistent with Fig. 8, which shows high and low NINO3.4 values in both the lower (negative NAO) and upper (positive NAO) quadrants of the scatter plot. Given that we find no systematic phasing of ENSO with Barents–Kara sea ice variability during the reanalysis period, it is likely that ENSO contributes to intermittency in the ICE–NAO pathway.

Figure 8 Scatter plots between February NAO and the late fall (mean of October and November) Barents–Kara sea ice index for the reanalysis period. Shading indicates the DJF NINO3.4 index. Red (blue) denotes El Niño (La Niña) events.

In terms of reinforcing the stratospheric pathway, blocking over the Urals region seems to play a particularly important but not fully understood role. Enhanced Urals sea level pressure is closely linked to the Scandinavian pattern in Euro-Atlantic climate variability and is related but not directly equivalent to the occurrence of atmospheric blocking. The Urals linkage appears in the monthly CEN ($\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}\Rightarrow \uparrow {\mathrm{V}}^{\ast }{{\mathrm{T}}^{\ast }}_{\text{Dec/Jan}}$ in Fig. 3a). The latter segment from Urals sea level pressure to poleward eddy heat flux is fairly systematic (appears in 74 % of the bootstrap samples in Fig. 5) and is grounded in the idea that tropospheric precursors over the Urals lead polar vortex weakening (Cohen and Jones, 2011; Cohen et al., 2014). However, the first segment from Barents–Kara sea ice to Urals blocking is more intermittent (appears in 50 % of the bootstrap samples), and whether it is in fact a causal linkage has been questioned by a recent modelling study using ensemble nudging experiments (Peings, 2019). Interestingly, not only Barents–Kara sea ice (Fig. 5e in King et al., 2016) but also ENSO (Fig. 5e and f in King et al., 2018) has been linked to the Scandinavian pattern, which suggests another avenue for ENSO to contribute to intermittency.

The ICE–URALS relationship highlights the complexity of interactions between atmospheric internal variability and Barents–Kara sea ice over a range of timescales. On synoptic scales, the pentad CEN (Fig. 7) shows linkages from reduced sea ice to enhanced Urals sea level pressure but also linkages in the opposite direction ($\uparrow \text{URALS}\Rightarrow \uparrow \text{IR}\Rightarrow \downarrow \text{ICE}$), with Urals sea level pressure altering ice cover via changes in poleward moisture transport that have been tied to synoptic moist intrusions (Woods et al., 2013; Luo et al., 2016; Gong and Luo, 2017; Lee et al., 2017). This chain of linkages can act as a positive feedback on sea ice perturbations but also provides a pathway by which blocking variability (internal to the atmosphere) may interrupt the expected troposphere–stratosphere coupling in response to autumn sea ice (for example, imagine a case where atmospheric conditions inhibit Urals blocking after a low-ice autumn). Furthermore, enhanced Urals blocking and moist intrusions can lead to highly transient perturbations in turbulent heat flux over the Barents and Kara seas. Initially, turbulent heat loss from the ocean is suppressed near the sea ice edge where moist intrusions act to weaken temperature and moisture contrasts between the atmosphere and ocean (Woods et al., 2013; Gong and Luo, 2017). But the heat flux anomaly can become positive (enhanced heat loss from the ocean) after the sea ice melts back in response to the moist intrusion, 1 to 2 weeks later (Woods and Caballero, 2016; Lee et al., 2017). On longer (monthly to seasonal) timescales, there is evidence that atmospheric variability is the main driver of heat flux variability over the Barents and Kara seas both in observations and models (Sorokina et al., 2016; Blackport et al., 2019). This perhaps explains why turbulent heat flux does not show up in the monthly CEN (Fig. 3a) but does show up in the half-monthly CEN (Fig. 3b). Across synoptic to seasonal timescales, it appears that sea ice is best thought of as an intermediary rather than a true boundary forcing, as is implied by prescribed sea ice (e.g. AGCM) experiments.

One outstanding issue involves the mechanisms that have been proposed to explain the $\downarrow {\text{ICE}}_{\text{Oct}}\Rightarrow \uparrow {\text{URALS}}_{\text{Dec}}$ linkage, which act on timescales that are inconsistent with the 2-month delay found in observations. For example, reduced sea ice may allow more heating of the atmosphere by the ocean to produce a Rossby wave train with an anomalous high over the Urals region (Honda et al., 2009), but this would be expected to manifest within a matter of days to a week. Alternatively, reduced ice may reduce local baroclinicity, which discourages cyclones from tracking into the Barents and Kara seas and produces an anomalous high due to the relative absence of low-pressure systems (Inoue et al., 2012). This mechanism could introduce some delay between the ice perturbation and sea level pressure perturbation, but a 2-month persistence of such a pattern is unlikely. Finally, reduced ice may increase atmospheric moisture content, leading to increased Eurasian snow cover, diabatic cooling, and anomalously high sea level pressure over the continent (Liu et al., 2012; Cohen et al., 2014; Garcia-Serrano and Frankignoul, 2014). Though this would plausibly lead to persistence on the required timescale, recent observational and modelling studies do not support a role for Eurasian snow in this teleconnection pathway (Kretschmer et al., 2016; Peings et al., 2017; Henderson et al., 2018), and we chose not to include it in our main analyses. Note that these mechanisms may still be responsible for contemporaneous forcing of the winter atmospheric circulation by winter sea ice variability, which has been suggested to be a stronger influence than the forcing by autumn sea ice variability (Blackport and Screen, 2019).

Lastly, our experience with the CEN offers some cautionary notes about its application to climate problems. The CEN approach was designed for hypothesis testing, i.e. to test causal pathways that are thought or known to exist, either from theory or existing evidence. It should not be used as an exploratory data analysis tool to search for causal pathways because the statistics behind the CEN do not know whether relationships are physically meaningful. One specific problem we encountered is that the algorithm may drop an existing causal linkage if a new variable is added. For example, when we introduce downward longwave radiation into the monthly CEN, its strong correlation with sea ice overrides the $\downarrow {\text{ICE}}_{\text{Jan}}\Rightarrow \uparrow {\text{URALS}}_{\text{Feb}}$ linkage (see Fig. S9 compared to Fig. S1a). Since many climate variables are highly correlated but not necessarily directly related via specific processes, the CEN’s ability to identify physically meaningful linkages depends critically on the careful selection of input variables.

5 Concluding remarks

This study uses the Causal Effect Networks (CEN) framework to quantify the robustness of the stratospheric pathway between late autumn Barents–Kara sea ice and the February NAO, documenting its high level of intermittency in the observational record. The pathway has been relatively active over the satellite period, explaining approximately 26 % of the interannual variability in the February NAO. However, this result is highly sensitive to which winters are included in the analysis. Results from a bootstrapping test show that the full stratospheric pathway appears in only 16 % of the sample populations derived from the observations. The result reflects the strong internal variability of the mid-latitude atmosphere and the likelihood that Arctic–mid-latitude teleconnections may require certain background flow conditions. On synoptic timescales, we identify two-way interactions between Barents–Kara sea ice and the mid-latitude circulation, suggesting a role for atmospheric blocking over the Urals region and moist intrusions, both of which can reduce Barents–Kara sea ice. These synoptic processes can reinforce or interrupt the stratospheric pathway, contributing to intermittency. Finally, we cannot rule out that the causal linkages found on longer timescales may be artefacts of averaging over the synoptic processes or even the result of entirely different mechanisms (Smith et al., 2017; Hell et al., 2020).

Coupled interactions between sea ice and the mid-latitude circulation involve complicated lead–lag feedbacks over a range of timescales. Applying causal inference frameworks such as the CEN can help clarify some of the important physical processes at play, but, in the end, models are required to improve our understanding. A complication is that the fidelity of climate models in representing the relevant processes is difficult to ascertain (King et al., 2016; Smith et al., 2017; Mori et al., 2019), especially those processes at fine spatial and temporal scales and their interactions across scales. However, ways forward are indicated by this study, as well as others (McCusker et al., 2016; Sun et al., 2016; Peings, 2019), that provide insight into which linkages are most robust and which are subject to sampling issues within the relatively short observational record.