Quantifying the Timescale and Strength of Southern Hemisphere Intraseasonal Stratosphere‐troposphere Coupling

Abstract The Southern Hemisphere zonal circulation manifests a downward influence of the stratosphere on the troposphere from late spring to early summer. However, the strength and timescale of the connection, given the stratospheric state, have not been explicitly quantified. Here, SH zonal wind reanalysis time series are analyzed with a methodology designed to detect the minimal set of statistical predictors of multiple interacting variables via conditional independence tests. Our results confirm from data that the variability of the stratospheric polar vortex is a predictor of the tropospheric eddy‐driven jet between September and January. The vortex variability explains about 40% of monthly mean jet variability at a lead time of 1 month and can entirely account for the observed jet persistence. Our statistical model can quantitatively connect the multidecadal trends observed in the vortex and jet during the satellite era. This shows how short‐term variability can help understand statistical links in long‐term changes.

(Section 4.2). Text S4 explains the procedure needed to map the linearly forced timeseries, generated using the standardized linear coefficients of Eq. 1, back onto its physical variables (Section 4.3).
Table S1 lists all the links (significant and not) computed in the final analysis of Section 4.2.
Text S1 The algorithm used to produce the present analysis (PCMCI, Runge (2018)) is a modification of the more commonly used causal algorithm (PC, Spirtes and Glymour (1991)) with the additional benefit of providing an estimate of link strength (Runge et al., 2014).
The latter property is closely related to the computation of the conditioning set Z M IT .
Given a system of N variables X t with t = 1, . . . , T , consider each pair of lagged variables (X t , Y t+τ ) at each lag τ = 1, . . . , τ max . To determine which pairs are predictor-predictand (as defined in the main text, Section 3) the partial correlation ρ(X t , Y t+τ |S t <t+τ ) has to be computed. By definition the conditioning set is S t <t+τ = {X t+τ −1 , X t+τ −2 , . . . } \ {X t } (up to a maximum lag in the past τ max ) and it can be very high dimensional (N τ max − 1).
While S is a sufficient set to detect conditional independence, it is usually much larger than necessary. The high dimensionality not only represents a computational cost, but also reduces the detection power and thus the reliability of the partial correlation itself (Runge, 2018).
Since in many cases the true links are sparse, PCMCI operates a preliminary step to extract a low dimensional subset of S. This first step, similar to the PC algorithm, removes irrelevant conditions for each of the N variables by iterative independence testing. In this way, to each variable X is assigned a (generally) low dimensional set of potential predictorŝ P X . The conditioning subset for the pair (X t , Y t+τ ) is therefore the union of both variables' potential predictors: After this preliminary stage, the algorithm goes through the rejection/estimation step (MCI). The null hypothesis is that each pair is connected. By computing the partial correlation ρ M IT , i.e. with Z M IT as conditioning set, the null hypothesis of a connection is rejected if the computed p-value is larger than a user defined threshold α. MIT stands for Momentary Information Transfer and refers to the property of the measure being lag-specific (momentary), thus removing autocorrelations. This step allows to reliably quantify the strength of the links.
Compared to PC, this method allows to maximise the detection power of the partial correlation and to give an estimate of the strength of the (significant) links. For a detailed discussion of this result we refer to Runge (2018).

Text S2
The following two-variable linear system, studied in Runge et al. (2014), motivates the conclusion of Section 4.1 in the main text from a general perspective. The system describes a stochastic process of two auto-dependent variables (through a, b, with |a|, |b| < 1 for stationarity) that are one-way coupled, since X contributes to the state of Y (through c) but not vice-versa. The noise terms are independent and Gaussian with covariance matrix Σ and zero mean. The analytical expression of the lagged cross-correlation ρ(X t , Y t+τ ) = f (τ ; a, b, c, Σ) depends on all coefficients of the system (see Table 2 in Runge et al. (2014)). Relevant to our analysis is that the amplitude and position of its peak -often used as a diagnostic for lagged interaction -can be modulated according to the values of the auto-dependent coefficients a and b. Also, the function f does not markedly decay for a specific sign of the lag, thus potentially suggesting that Y directly influences X. This is because correlation measures the collinearity between lagged time series, which can be due to all sorts of effects. On the other hand, the partial correla- where δ τ,1 is the Kronecker delta function, Table 2 in Runge et al. (2014)). The effect of both autocorrelations has been removed (no a, b dependence) thanks to the conditioning set The difference between the correlation and partial correlation function is largely due to the auto-dependence coefficients. In the vortex-jet case, the large auto-dependence of the vortex (Fig.1b, top-left) can be seen as analogous to a large coefficient a herein. The fact that conditioning on a Z M IT = P oV m−1 (or P oV m+τ −1 ) results in a small downward and upward ρ M IT shows the vortex auto-dependence is indeed responsible for the large and persistent cross-correlations.

Text S3
The autocorrelation of vortex and jet analytically derived from Eq. 1 read:

Text S4
The coefficients computed for Eq. 1 in Section 4.2 come from standardised time series.
The internal variability of the standardised vortex is Var [PoV ] = σ 2 P /(1 − a 2 ), according to Eq. 1 ( to estimate it we use coefficients calculated in Section 4.2). To ensure that the trending synthetic vortex time-series is appropriately scaled we need to impose that it has the same variability-to-trend ratio as the observed vortex. The ratio is R = vortex variability/ strength increase due to the trend. Keeping Var [PoV ] fixed while imposing R obs = R synth , this results in the synthetic trend coefficient beingγ = K · γ P obs with the Note that K depends only on vortex parameters. In particular, the observed jet trend is not used in any way. Specifically, R obs = Var [PoV]/T obs γ P obs , with γ P obs (ms −1 /month) obtained from reanalysis as described in Section 4.3, T obs =22(years in one Epoch)· 12 (months in a year) and V ar [PoV] in the variance of the PoV monthly anomalies from reanalysis. R synth = Var [PoV ]/T synthγ , with T synth = 30 time steps (an arbitrary length chosen to be close to the small sample size from which each calendar day-specific reanalysis trend is estimated). Because the analytic prediction for the synthetic jet trend isγc/(1 − a), we can convert it to physical values by K −1 multiplication.