The different stratospheric influence on cold-extremes in Eurasia and North America

The stratospheric polar vortex can influence the tropospheric circulation and thereby winter weather in the mid-latitudes. Weak vortex states, often associated with sudden stratospheric warmings (SSW), have been shown to increase the risk of cold-spells especially over Eurasia, but its role for North American winters is less clear. Using cluster analysis, we show that there are two dominant patterns of increased polar cap heights in the lower stratosphere. Both patterns represent a weak polar vortex but they are associated with different wave mechanisms and different regional tropospheric impacts. The first pattern is zonally symmetric and associated with absorbed upward-propagating wave activity, leading to a negative phase of the North Atlantic Oscillation (NAO) and cold-air outbreaks over northern Eurasia. This coupling mechanism is well-documented in the literature and is consistent with the downward migration of the northern annular mode (NAM). The second pattern is zonally asymmetric and linked to downward reflected planetary waves over Canada followed by a negative phase of the Western Pacific Oscillation (WPO) and cold-spells in Central Canada and the Great Lakes region. Causal effect network (CEN) analyses confirm the atmospheric pathways associated with this asymmetric pattern. Moreover, our findings suggest the reflective mechanism to be sensitive to the exact region of upward wave-activity fluxes and to be state-dependent on the strength of the vortex. Identifying the causal pathways that operate on weekly to monthly timescales can pave the way for improved sub-seasonal to seasonal forecasting of cold spells in the mid-latitudes. Cold spells in the Northern Hemisphere mid-latitudes are influenced by the stratospheric polar vortex in two different ways. Marlene Kretschmer from the Potsdam Institute for Climate Impact Research, Germany, and collaborators use cluster analysis to show there are two dominant patterns associated with a weak polar vortex. The first is zonally-symmetric, associated with absorbed upward wave activity fluxes and is linked to a negative North Atlantic Oscillation and cold spells over Eurasia. Whilst the second is zonally-asymmetric, associated with downward reflected waves over Canada and is linked to a negative Western Pacific Oscillation and cold spells over North America. This knowledge could allow better prediction of winter weather on sub- to seasonal scales.


INTRODUCTION
Variability in the stratospheric polar vortex in boreal winter can influence the tropospheric circulation and is an important source of predictability on sub-seasonal to seasonal (S2S) timescales. 1,2 Particularly, extremely weak polar vortex states, such as sudden stratospheric warmings (SSW), can be impactful to societies as they are often associated with large-scale cold-air outbreaks in the densely populated mid-latitudes. 3 Understanding the exact coupling mechanisms between the troposphere and the stratosphere is hence central to improve S2S predictions in the midlatitudes including cold spells.
SSWs have been extensively studied and have been classified by their spatial properties 4 (split versus displaced events), by the dominant type of wave forcing 5 (wave 1 versus wave 2) or by their intensity 6 (major versus minor warmings). However, these metrics describe the stratospheric extreme events itself but do not necessarily capture differences in the tropospheric response. 7,8 Recently, Kodera et al. 8 classified SSWs according to the different coupling mechanism between vertical wave activity and the stratospheric polar vortex and distinguished between socalled absorbing and reflecting SSWs. The former are characterized by a several week long pulse of troposphere-induced upward wave activity fluxes, which is absorbed in the stratosphere leading to increased temperatures over the polar cap and an overall weakening of the zonal-mean flow. Subsequently, the circulation anomalies, characterized by a negative phase of the Northern Annular Mode (NAM), 3,[9][10][11] migrate downward into the troposphere where they can persist for up to two month. 3,9,10,12 Consequently, a negative phase of the North Atlantic Oscillation (NAO) and an equatorward-shifted jet are then frequently observed, often associated with mid-latitude cold-snaps, especially over northern Eurasia. For North America this stratospheric NAM/NAO mechanism is, however, much less robust. 10,13,14 In addition to absorbing SSWs affecting the tropospheric circulation via the downward propagation of zonal-mean anomalies, 3,12,15 Kodera et al. 8 introduced so-called reflecting SSWs. These are associated with shorter pulses of vertical wave-activity fluxes, a more quickly recovering stratospheric flow and reflection of the upward-propagating planetary waves downward. More generally, irrespective of the occurrence of SSWs, various studies showed that a sufficiently strong upper stratospheric polar vortex can reflect upward-propagating waves downward, influencing the tropospheric flow. [16][17][18][19][20][21] Analyses of individual winters have indicated that wave-reflection can lead to increased ridging over the North Pacific associated with a negative phase of the Western Pacific Oscillation (WPO) and its sea-level pressure signature, the North Pacific Oscillation (NPO). 18,19 The WPO/NPO is the Pacific analog of the NAO, and projects strongly onto North American temperature and precipitation variability. 22,23 Several cases have been discussed qualitatively, but mostly with a focus on major SSWs. 8,18,19 Overall, the possible role of the stratospheric polar vortex, in particular of wave-reflection, for winter cold-spells in North America is, thus, not well understood.
Here we present cluster analysis 13,24,25 to quantitatively study spatial patterns of the polar vortex linked to mid-latitude coldspells. We focus on the role of the stratospheric polar vortex on tropospheric circulation in boreal winter and cold-extremes in northern Eurasia and North America. Moreover, we apply a novel type of time series analysis, called causal effect networks (CEN), 26 to extract the potential causal pathways associated with wavereflection.

Cluster analysis
We perform hierarchical clustering on the daily geopotential height anomaly field at 100 hPa (Z100) poleward of 60°N in winter (January and February) from 1979 to 2018 (Methods). We chose such a lower stratospheric level as those are especially crucial for troposphere-stratosphere coupling, 11,12,14 and we focus on January and February because these months show the largest vortex variability. The choice of the number of clusters is usually somewhat subjective and metric dependent. However, here the goal of our clustering approach is to identify specific events of interest and use those as a starting point for a more detailed statistical analysis. We discuss the robustness of the clustering in the SI.
We detect five clusters capturing lower stratospheric variability, presented by their composites calculated over all days that were assigned to the cluster and ordered by mean polar cap height in Fig. 1. Cluster 1 represents extremely strong polar vortex states, with negative geopotential height anomalies over the entire polar cap. Clusters 2 and 3 show subsequently weaker and displaced vortex patterns. In the remainder of the manuscript, we restrict ourselves to discussing only the two weakest vortex clusters 4 and 5.
Cluster 4 is characterized by zonally asymmetric Z100 anomalies with strong positive anomalies over Siberia, but negative anomalies over Canada and the North Atlantic. It is detected on approximately 14% of all winter days with a total of 48 events (defined as consecutive days of the same cluster), giving an average duration of 7 days. Though cluster 4 events are characterized by a disturbed polar vortex (in terms of a high polar cap height anomaly), only 4 out of 48 cluster 4 events are associated with major SSWs (see SI for a discussion on the robustness of the clustering and a detailed comparison with other metrics). In contrast, cluster 5 represents a completely disturbed polar vortex with positive zonally symmetric Z100 anomalies over the entire polar cap. Approximately 16% of all winter days are assigned to the weakest polar vortex events described by cluster 5, coming from only 32 events with a mean-duration of nearly 12 days. Most of the observed extremely weak vortex states, defined as major SSWs, coincide with cluster 5 events (see SI).

Reflective and absorbing stratospheric pathways
To test for the dynamical coupling mechanisms, we first compute the absolute and anomalous vertical wave-activity flux (WAF), here calculated as the vertical component of the Plumb fluxes 27 at 100 hPa, averaged over all cluster 4 (Fig. 2a, c) and cluster 5 events (Fig. 2b, d). During cluster 4 events, waves propagate upward over eastern Siberia with simultaneous downward propagation over Canada (Fig. 2a), which is also characterized by significant (p < 0.05, see Methods) positive and negative wave flux anomalies respectively in these regions (Fig. 2c). The enhanced upward wave activity over eastern Siberia is consistent with the associated Z100 dipole pattern of cluster 4 days, showing a shifted polar vortex towards eastern Canada (Fig. 1). During cluster 5 events, we find upward wave propagation across the hemisphere ( Fig. 2b) with strongest positive anomalies over Canada and the North Atlantic (Fig. 2d).
To assess the coupling mechanisms in more detail, we next plot the temporal evolution of different standardized stratospheric indices before, during and after cluster 4 (Fig. 3a, c, e) and cluster 5 events (Fig. 3b, d, f). Composites for the absolute values are shown in the SI. In each panel, lag 0 denotes the start day (i.e., the first day cluster 4 and cluster 5 events were respectively detected). We construct an index P 4 (P 5 ) describing the similarity between the observed polar vortex pattern and cluster 4 (cluster 5) for each day in winter: We project the daily Z100 anomaly field onto the cluster composite ( Fig. 1) and normalize the index by its multi-year standard deviation of the respective day. The upper row in Fig. 3 shows this P 4 (P 5 ) index as well as the mean polar cap height index (PCH) at 10 hPa. The difference in event-duration between cluster 4 and cluster 5 becomes evident again by the evolution of the P 4 and P 5 indices, with cluster 4 events only showing a relatively short period (14 days) of significantly high P 4 values (Fig. 3a), while the P 5 index is significantly increased in the 20 days before and after the detection of cluster 5 events (Fig. 3b). The middle row plots the vertical component of the Plumb flux at 100 hPa averaged over the latitudinal belt from 50°N to 75°N (WAF Hem , red   Fig. 3c, d). For cluster 5, however, the positive wave-activity anomalies precede the event start much earlier.
Consistently, the PCH is already significantly large before and during the onset of cluster 5 events and the vortex remains weakened afterwards (Fig. 3b). The increased hemispheric vertical wave-activity anomalies preceding cluster 5 events are preconditioned by significantly positive wave fluxes (Fig. 3f, see Fig. S4b for absolute values) over Canada~2 weeks before the event start and characterized by enhanced wave flux anomalies over the Euro-Atlantic sector with the event onset (Fig. 3d). Wave flux anomalies over Siberia are slightly positive but drop with the event start ( Fig.  3f). For cluster 4, in contrast, the hemispheric vertical wave-activity flux anomalies are significantly increased only 5 days before the event start, and become negative shortly after (Fig. 3b). Consistently, the polar vortex is briefly disrupted (i.e. the PCH rises) and recovers quickly thereafter (Fig. 3a). Also the regional Plumb flux indices show a different evolution: Vertical waveactivity flux is significantly enhanced over eastern Siberia for several days (Fig. 3e) but shows only moderate anomalies over the Euro-Atlantic sector (Fig. 3c). Moreover, wave-activity fluxes over Canada become significantly negative with the event start ( Fig. 3e, Fig. S4a for absolute values), indicating downward propagation over this region.
Overall, these findings are consistent with previous studies and support the notion of a reflecting (cluster 4) and absorbing (cluster 5) coupling mechanisms. 8,[16][17][18][19] Cluster 5 events, often coinciding with SSWs, are associated with persistent stratospheric disturbances (absorbing-type), preceded by enhanced hemisphere-wide vertical wave activity that is absorbed in the stratosphere leading to increased geopotential heights over the entire stratospheric polar cap. The composites for cluster 4 events, of which only a few are associated with major SSWs, suggest that a strong but shortlasting pulse of vertical wave-activity resulting from enhanced a Color shading shows the absolute values during cluster 4 events and the contours represent the winter climatology. b Same as in a but during cluster 5 events and with contours indicating regions over Siberia (pink), Canada (brown) and the Euro-Atlantic sector (red), which are used to calculate regional indices. c Color shading shows the anomalous values during cluster 4 events and the stippling shows significant values (p < 0.05, see Methods). d Same as c but for cluster 5 events The different stratospheric influence on cold-extremes in Eurasia. . . M Kretschmer et al. upward propagation over eastern Siberia is reflected by the stratospheric flow and descends over Canada (reflecting-type), in agreement with previous studies. 8,17 Also the precondition of the polar vortex (i.e., already weak before cluster 5 and neutral before cluster 4 events) is in agreement with earlier studies 16,17,28 and supports the finding of reflected (absorbed) waves if the stratospheric flow is sufficiently strong (weak). We find that 70% of all cluster 4 events occurred during the westerly phase of the Quasi-Biennial-Oscillation (QBO-W), which is linked to a strengthened polar vortex 29 and might thus favor the occurrence of cluster 4. Thus, the locations of the upward wave-activity flux and the initial strength of the polar vortex seem to have a central role in whether the stratospheric pattern is of reflective or of absorbing-type.
Note that several cluster 4 (and some cluster 5 events) can occur in the same winter, with the events sometimes only being separated by a few days, leading to overlapping time-ranges for the composites. Composites of all individual cluster 4 and cluster 5 events and of the associated wave flux activity are presented in the SI, partly showing substantial differences, but overall supporting the main findings. Here we decided not to apply further criteria to exclude or merge particular individual events and leave this for subsequent research.
Connection to tropospheric circulation and cold-extremes To investigate the influence of reflecting and absorbing events on the tropospheric circulation, we next compute the geopotential height anomalies at 500 hPa (Z500) during cluster 4 (reflectingtype) and cluster 5 (absorbing-type) events (Fig. 4a, b). As expected, 3,9,10,13 cluster 5 days coincide with a negative NAOtype Z500 pattern with increased geopotential heights over Iceland and the North Atlantic and anomalously low values over the Azores (Fig. 4b). Moreover, during cluster 5 days the NAO index is significantly (p < 0.01, using a two-sample Kolmogorov-Smirnov test) shifted towards more-negative values compared to all other winter days (Fig. 4d). Also the daily P 5 index is significantly (p < 0.01 according to a bootstrapping test) correlated (r = −0.5) with the NAO index. Thus, the absorbingtype cluster 5 events resemble the well-documented case of a zonally symmetric disturbed vortex including downward propagation of a negative NAM pattern. 3,9,11,12 In contrast, Z500 composites during cluster 4 days show a negative phase of the WPO with lower geopotential heights over the central Pacific and pronounced positive anomalies over the North Pacific (Fig. 4a). Moreover, Z500 anomalies over the Atlantic show a positive NAO-like dipole, but with the centers of maximum anomalies shifted southward. Comparing the WPO during cluster 4 days and all other winter days, we find again a significant negatively shifted distribution during cluster 4 days (Fig. 4c). More generally, the WPO and the daily P 4 index are significantly correlated (r = −0.5) but note that the P4 (P5) index and the NAO (WPO) only show a correlation of 0.09 (−0.02). The negative WPO pattern is consistent with the Siberian/Aleutian high in the stratosphere, descending downward due to suppression of upward wave propagation as was described in detail in Kodera et al. 8 Our results thus strongly support the hypothesis that reflecting-type events, here represented by cluster 4, can lead to increased geopotential heights over the North Pacific.
In agreement with the negative phase of the NAO during cluster 5 days, near-surface temperature composites show a pronounced pattern of significant cold anomalies over northern Eurasia (Fig.  5b). To assess the relationship between cold-extremes in this region and cluster 5 events in more detail, we count the frequency of each cluster during the 10% coldest days in northern Eurasia (averaged over 50°N-65°N, 10°E-130°E, see blue box in Fig. 5b) and normalize this number by their overall occurrence frequency (Fig. 1). A value of 1 thus indicates that the cold-extreme frequency is the same as expected by chance, implying no statistical relationship between the stratospheric cluster and coldextremes. We find that cold-extremes occur twice as often during cluster 5 days than expected by chance (Fig. 5d), indicating a strong statistical relationship between northern Eurasian coldsnaps and cluster 5 events. Yet, over North America there are no significant cold anomalies during cluster 5 events (Fig. 5b). In fact, the 10% coldest winter days in this region (averaged over 40°N -55°N to 100°W-70°W, orange box in Fig. 5a), coincide less often with cluster 5 than expected (Fig. 5c). During cluster 4 events, however, there is a pronounced pattern of negative temperature anomalies over most of Canada, the Great Lakes regions and the eastern US (Fig. 5a). The 10% coldest days in North America (orange box in Fig. 5a) occur more than twice as often during cluster 4 events than statistically expected (Fig. 5c).
Overall, the analyses thus present a strong statistical relationship between the reflecting (absorbing) cluster 4-type (cluster 5) pattern and the WPO (NAO) and associated cold-extremes over parts of North America (northern Eurasia). Moreover, in agreement with the composites, seasonal-mean cluster 4 frequency can explain a large fraction (up to 45%) of seasonal temperature variability over Canada and the US (Fig. 6a), while cluster 5 accounts for variability over large parts of Eurasia but not over North America (Fig. 6b). Thus, although in particular the reflectingtype events are only of short duration, their seasonal-mean projects strongly onto winter temperature variability. Understanding the precursors of wave-reflection might thus help to improve predictions for this region. For example, Kodera et al. 19 proposed that blocking over the North Atlantic can trigger a wave train into the stratosphere, which by suppression of upwardpropagating waves can cause blocking over the North Pacific. In subsequent research, we will assess the role of these tropospheric drivers in more detail.
Causal effect network analysis The composite analysis based on the detected cluster 4 events suggest that a stratospheric pathway contributes to the formation of North Pacific blocking, consistent with previous studies on wave-reflection. 18,19 However, clustering includes some subjective criteria which might lead to selection bias, potentially producing non-robust results when studying continuous time series. Therefore, to test the involved hypotheses and to assess the causal relationship between the considered indices in more detail, we apply causal effect network (CEN) analysis (Methods). CEN is a multi-variate statistical framework based on a causal discovery algorithm, introduced to climate science for hypothesis testing of teleconnection processes. 26 CEN detects spurious correlations (e.g., due to common drivers, indirect links or auto-correlation effects 30 ) by iteratively calculating partial correlations between different combinations of variables and at different lags. Those relationships that are found to be conditionally dependent (i.e., for which the linear relationship is significantly different from zero even when the influence of a combination of other drivers or auto-correlation effects is excluded) form the links in the CEN. These links can be interpreted as potentially causal for the set of considered processes and time-lags. Thus, CEN-analysis allows for much stronger statements towards a causal interpretation beyond simple cross-correlation analysis or the bivariate concept of Granger-causality. [30][31][32] Here we particularly want to test the hypothesis that enhanced vertical wave activity over eastern Siberia (WAF S ) leads to increased geopotential heights over the North Pacific (Pac). We also include the PCH index at 10 hPa to describe stratospheric polar vortex variability (PCH). Moreover, we include the regional WAF index over Canada (WAF C ). The regions over which the timeseries indices are calculated are displayed in Fig. 7.
Since we are interested in low-frequency sub-seasonal variability, we remove synoptic variability by calculating centered 5-day mean time-series (i.e., by calculating means over bins of 5 days) before accessing their causal links. Thus, a lag-1 relationship refers to a lag of 6-10 days which covers the timescales at which links are approximately expected based on the composites (Fig. 3) and previous studies on wave-reflection. 8,16,19 We perform CENanalysis for the winter months of January and February (JF) only. When performing conditional independence tests, we allow timelags of up to one month (lag-6). The detected causal links (i.e., the conditionally dependent relationships) are shown in Fig. 8, whereby red arrows correspond to significant (p < 0.01) positive linear relationships and blue arrows represent negative relationships. The node-color denotes the strength of the strongest autocorrelation coefficient (at lag-1). The exact values are given in Table 3 of the SI.
The CEN calculated for continuous time-series supports the findings from the event-based composites and thus confirms the proposed hypotheses. Using the terms "causes" and "leads" as conditional dependence statements, we find that an increase in vertical wave activity over Siberia (WAF S ) causes an increase in geopotential heights over the North Pacific (Pac) with a lag of 6-10 days (lag-1). Moreover, increased Siberian wave-activity flux (WAF S ) is linked to decreasing vertical wave-activity over Canada (WAF C ). These links are thus consistent with the reflecting mechanism described before. Furthermore, as expected, 15,26,33 increased upward wave-activity over Siberia (WAF S ) in winter leads to a rise in PCH, thus a weakening of the polar vortex. The other detected links, from increased North Pacific geopotential heights (Pac) to a decrease in PCH (i.e., strengthening of the polar vortex), directly or via WAF S , confirm that positive geopotential height anomalies over the North Pacific (Pac) lead to a strengthening stratospheric flow, probably through destructive interference with the climatological wave, as hypothesized by previous studies. 19,34,35 Consistently, wave activity over Canada (WAF C ), showing downward propagation in the climatological mean (Fig.  2a), decreases after Pac increases, thus also representing a weakening of the climatological state. Our results, however, show that the relationship between Pac and PCH is bi-directional, demonstrating once more the complex relationship between tropospheric blocking and stratospheric vortex variability. 34,36,37 In summary, although processes not considered or those operating at different timescales (e.g., tropical Pacific variability) might affect the results, the CEN provides robust evidence of a reflective mechanism influencing tropospheric circulation via a stratospheric pathway.

DISCUSSION
Using cluster analysis, we identified (1) a pattern of extremely weak polar vortex states (cluster 5) linked to cold-spells over northern Eurasia and (2) a pattern with increased geopotential height over eastern Siberia (cluster 4) associated with coldextremes over the Northeastern US. The former is in agreement with the well-documented relationship between a strongly disrupted polar vortex (absorbing-type), a downward propagating negative NAM 3,9 and Eurasian cold spells. Combining lagged composites and CEN, we demonstrate that cluster 4-type events (reflecting-type) are linked to a negative WPO and North American cold-spells via reflected upward-propagating waves over eastern Siberia, as was hypothesized in previous studies. 8,[16][17][18][19] Our results show that the stratospheric influence on winter circulation should not exclusively be analyzed in terms of major SSWs and an associated zonally symmetric, downward propagating NAM signal. When studying North American cold-spells, it seems crucial to also consider zonally asymmetric vortex states associated with wave-reflection events (cluster 4). Although these persist on average only for approximately a week, seasonal (JF) frequencies of this cluster can explain a much larger part of North America's seasonal temperature variability than cluster 5 (Fig. 6). A better understanding of the relevant low-frequency processes and boundary conditions favorable for wave-reflection might therefore help to improve S2S predictions for this region. 38,39

METHODS Clustering
We use daily-mean ERA-Interim 40 data from 1979 to 2018, leap days excluded. We perform hierarchical clustering on the daily geopotential heights anomalies field at 100 hPa (Z100) poleward of 60°N in winter (January and February). Climatological anomalies for each day are calculated by subtracting their multi-year mean. To account for the denser grid towards the pole, we apply area-weighting by multiplying each value with the cosine of its latitudinal location. The cluster algorithm starts with n clusters (the starting vectors) and then iteratively merges two clusters until only one cluster (the mean over all vectors) exists. In each step, the clusters with minimal distance are merged and their mean is calculated. Here we use Ward's metric criteria, meaning that the two clusters to be merged at each step are those which result in the minimal increase in variance in the merged cluster, over all possible unions of clusters. Hierarchical clustering has the advantage that the number of clusters does not need to be chosen a priori but can be decided based on the dendogram (Fig. S2). Nevertheless, the choice of the 'optimal' number of clusters usually requires some subjective criteria.

Composites
Before calculating composites, we detrend the respective variables over the considered time-period to prevent selection biases. Significance is tested creating 1000 artificial time-series at each grid-point by randomly selecting and shuffling blocks of the original time-series (with a blocklength of five days, i.e., beyond synoptic variability). For each newly created time-series, we pick as many days as were used to form the composite but we also keep the start days and length of the identified events from the original time-series to account for a potential increase in auto-correlation during long-lasting cluster events. We apply Benjamini-Hochberg false discovery rate (FDR) correction to the spatial composites to account for the increase of false positives due to multiple testing. 41 Given the spatial correlation of the fields, we use α FDR = 2*0.05, to get a global α level of 0.05.

Climate indices
The Western Pacific Oscillation index (WPO) is constructed by subtracting the (area-weighted) Z500 anomalies over the region 50°N-70°N and 200°E -235°E (Pac) from the region 25°N-40°N and 140°E-210°E. The North Fig. 6 Explained variance of winter temperature by cluster 4 and cluster 5. a R 2 values for regression of winter (JF) mean temperature at each grid-point on seasonal-mean cluster 4 frequency and, b on seasonal-mean cluster 5 frequency. Before calculation the regression models, the linear trends of the regressors and the temperature were removed. Significant (p < 0.05) models according to two-sided Student's t-test are indicated in hatches The different stratospheric influence on cold-extremes in Eurasia. . . M Kretschmer et al.
Atlantic Oscillation index (NAO) is calculated by subtracting the (areaweighted) sea-level pressure anomalies over the region 55°N-90°N and 90°W -60°E from the region 20°N-55°N and 90°W-60°E. 42 We calculate the vertical wave-activity fluxes (WAF) as described in Plumb et al. 27 and compute an index (WAF Hem ) as the area-weighted zonal-mean WAF at 100 hPa from 50°N-75°N. Next to this hemispheric index, we calculate regional components over the Euro-Atlantic sector (WAF Euro-Atl ; 40°W-35°E), eastern Siberia (W S ; 120°E-185°E) and Canada (WAF C ; 225°E-300°E). The polar cap height index (PCH) is computed as the area-weighted polar cap mean geopotential height anomaly northward of 60°N at 10 hPa. Indices are standardized, by dividing each value by its multi-year standard deviation. The QBO index at 50 hPa is taken from Freie Universität Berlin (http://www. geo.fu-berlin.de/met/ag/strat/produkte/qbo/qbo.dat).

Causal effect networks (CEN)
CEN is a multi-variate approach aiming to detect causal relationships amongst a set of time-series. Using partial correlations, it iteratively tests for conditional independence of two indices conditioned on a combination of other time-series at different lags (see Kretschmer et al. 26 for a detailed introduction). Thus, links in the CEN are those for which the linear relationship between two time-series cannot be explained by the combined influence of other included indices or by auto-correlation effects.
Here, CEN is constructed using the 2-step causal discovery algorithm PCMCI. 31 The first step is a condition-selection based on the PCalgorithm 30,43 which is run for each variable: It starts with all variables at all lags as potential causal drivers (=parents) and iteratively removes parents whose partial correlation with the target variable, conditional on iteratively chosen subsets of the other predictors, vanishes. Compared to the PC-algorithm implemented in Tigramite 1.0 (used in Kretschmer et al. 26 ), it utilizes a recent improvement of the stability of the PCalgorithm. 44 The significance level at which partial correlations are deemed non-significant is here chosen by comparing the parents obtained with the algorithm for different levels (0.05, 0.1, 0.2, 0.3, 0.4, 0.5) using the Akaike information criterion (AIC). In the second step, the momentary conditional independence (MCI) test, the parents with the highest score are then used as conditions when calculating partial correlations. 31 The MCI test is conducted for each pair of variables and at all time lags up to a maximum lag and yields a p-value for each link. To account for multiple testing, we adjust the p-values using the Benjamini-Hochberg FDR controlling procedure. Finally, to assess the strength of a link, standardized regression models are fitted to the significant parents of each target variable. The regression coefficients corresponding to each parent then yields the strength of the causal links in the CEN. Causal effect network (CEN) constructed for winter 5-days-mean indices of the polar cap height mean index at 10 mb (PCH), regional indices of vertical wave activity over eastern Siberia (WAF s ) and Canada (WAF C ), and Z500 over the North Pacific (Pac). The node-color denotes the strength of the auto-correlation coefficient at lag-1. Red (blue) arrows correspond to significant (p < 0.01) positive (negative) causal (i.e., conditionally dependent) relationship and the lag is 6-10 days (i.e., lag-1). Exact values of the links are given in the SI