Strengthened Causal Connections Between the MJO and the North Atlantic With Climate Warming

While the Madden‐Julian oscillation (MJO) is known to influence the midlatitude circulation and its predictability on subseasonal‐to‐seasonal timescales, little is known how this connection may change with anthropogenic warming. This study investigates changes in the causal pathways between the MJO and the North Atlantic oscillation (NAO) within historical and SSP585 simulations of the Community Earth System Model 2‐Whole Atmosphere Community Climate Model (CESM2‐WACCM) coupled climate model. Two data‐driven approaches are employed, namely, the STRIPES index and graphical causal models. These approaches collectively indicate that the MJO's influence on the North Atlantic strengthens in the future, consistent with an extended jet‐stream. In addition, the graphical causal models allow us to distinguish the causal pathways associated with the teleconnections. While both a stratospheric and tropospheric pathway connect the MJO to the North Atlantic in CESM2‐WACCM, the strengthening of the MJO‐NAO causal connection over the 21st century is shown to be due exclusively to teleconnections via the tropospheric pathway.

source and propagation of the Rossby waves into midlatitudes. Focusing first on changes to the MJO itself, Bui and Maloney (2019) show that ratio of MJO-induced circulation to precipitation in CMIP5 climate models decreases as the climate warms, and Hsiao et al. (2020) provide evidence that these changes may already be detectable in the observations. These results are important as they suggest that MJO teleconnections, which are directly excited by the MJO divergent circulation, may weaken as the climate warms Wolding et al., 2017). Thus, the response of the MJO alone to climate warming suggests that skill provided by the MJO could disappear in the coming decades.
Changes in Rossby wave propagation with climate warming are also expected to impact MJO teleconnections. Zhou et al (2020) investigated changes in MJO teleconnections within CMIP5/CMIP6 climate models and demonstrated a robust increase in teleconnections to the west coast of the United States. Their reasoning was that the robust extension of the North Pacific jet-stream changes the Rossby wave propagation paths, acting to shift the teleconnection centers eastward to more directly impact the U.S. west coast. How Rossby wave propagation to the North Atlantic and Europe may change is perhaps more complicated. MJO teleconnections to these regions are known to be dynamically driven by two different pathways, a direct tropospheric pathway and an indirect stratospheric pathway . The tropospheric pathway is communicated via tropospheric propagation of Rossby waves from the MJO region to the midlatitudes, which takes ∼10-15 days (Cassou, 2008;L'Heureux & Higgins, 2008;Lin et al., 2009). The stratospheric pathway instead involves MJO-excited Rossby waves propagating up into the stratosphere and disturbing the stratospheric polar vortex (a time lag of ∼15-30 days) Garfinkel et al., 2014Garfinkel et al., , 2012Kang & Tziperman, 2018;Weare, 2010). It is the signal from the disrupted polar vortex which then propagates downward into the troposphere, impacting the North Atlantic circulation (e.g., Baldwin & Dunkerton, 2001;Charlton-Perez et al., 2018;Garfinkel et al., 2014;Garfinkel & Schwartz, 2017;Kidston et al., 2015;Schwartz & Garfinkel, 2020). These multiple pathways connecting the MJO to the North Atlantic suggest that changes in teleconnection strength could come about due to changes in Rossby wave propagation through the troposphere, through the stratosphere, stratosphere-troposphere coupling, or changes to the MJO itself.
Here, we investigate how MJO teleconnections to the North Atlantic may change under anthropogenic climate change using three historical and five SSP585 simulations of the Community Earth System Model 2-Whole Atmosphere Community Climate Model (CESM2-WACCM). We invoke two data-driven approaches, namely, the STRIPES index and graphical causal models to quantify MJO teleconnection pathways within CESM2-WACCM and how they may evolve with climate warming. These approaches collectively indicate that the MJO's influence on the North Atlantic strengthens in the future and that this change is due exclusively to the tropospheric pathway.

CESM2-WACCM Simulations
We analyze simulations from the latest generation of the Community Earth System Model (CESM2; Danabasoglu et al., 2020) (Hurrell et al., 2013) have resulted in improved historical simulations including major reduction in low-latitude precipitation, and better representation of the MJO and extratropical atmospheric circulation, making it an ideal model for our study (Ahn et al., 2020;Danabasoglu et al., 2020;Simpson et al., 2020). This study was also repeated for the "low top" version-CESM2-and the general conclusions remain the same as those found for CESM2-WACCM (see Figures S6 and S7).
We analyze data from ensemble members 1, 2, 3 with historical (1850-2014) forcing and ensemble members 1, 2, 3, 4, 5 with SSP585 (2015-2099) forcing. SSP585 represents the high range of possible futures and exhibits an end of century radiative forcing of 8.5 W/m 2 (O'Neill et al., 2016). We focus on daily sea-level pressure, zonal wind at 850, 250, and 50 hPa, and precipitation. All fields are re-gridded to a 2° by 2° grid prior to analysis to speed-up computation and reduce data storage. While we center our focus on winter subseasonal variability in December, January, and February (DJF), the causal inference approach uses time-shifted/ lagged variables that are allowed to extend to November and March as appropriate.

Climate Indices
When computing each of the three climate indices below (i.e. MJO, VORTEX, and NAO) every ensemble member and 40-year period are treated completely separately to ensure that the analysis is agnostic to the number of ensemble members available. Prior to any analysis, we remove the third-order polynomial trend from each calendar day for each grid point for each field, which removes both the trend over the period and the seasonal cycle. This is done to ensure that all climate indices reflect subseasonal variability, rather than a long-term trend in the mean fields. Ensemble-averages are only performed as a final step prior to plotting the results.
The North Atlantic oscillation (NAO) is defined as the leading empirical orthogonal function (EOF) of area-weighted monthly mean SLP anomalies over the region 25°N-90°N, 90°W-30°E. The principal components are defined such that a positive value refers to a positive NAO state (low pressure over the poles and high pressure over the mid-North Atlantic) and have been standardized to have a mean of zero and standard deviation of one in December-February. NAO sea-level pressure composites from CESM2-WACCM are provided in Figure S1.
The state of the stratospheric polar vortex (VORTEX) is defined as the area-weighted average of anomalous daily 50 hPa zonal wind north of 65°N (e.g., Baldwin et al., 1994;Barnes et al., 2019;Charlton-Perez et al., 2018). The index is standardized by subtracting the mean value over DJF and dividing by the standard deviation over DJF for each 40-year period.
The MJO is quantified by the real-time multivariate MJO index (RMM index; Wheeler & Hendon, 2004) defined as the two leading EOFs of meridionally averaged tropical (15°S-15°N) anomalous fields of 250 and 850 hPa zonal wind and precipitation. We use model precipitation rather than outgoing longwave radiation (OLR) as is convention since precipitation is more directly simulated by the model. Numerous climate model diagnostics studies use precipitation anomalies to assess MJO simulation skill, and model MJO skill metrics are higher when evaluated with precipitation compared to OLR (Ahn et al., 2020(Ahn et al., , 2017. In addition to removing the seasonal cycle, we subtract the previous 120-day running mean from each field for each gridpoint to remove any low-frequency modes of variability prior to EOF analysis. Each of the three fields are also standardized by removing the mean and dividing by the standard deviation over all seasons. Figure S2 shows the structure of the two leading EOFs in the historical and future periods and the historical EOFs compare well with those from observations (Wheeler & Hendon, 2004). The two leading principal components (denoted as RMM1 and RMM2) are standardized to have mean of 0 and standard deviation of 1. The different phases of the MJO are then defined by the phase space of RMM1 and RMM2 in the conventional way (Wheeler & Hendon, 2004). We define an MJO event as any day where the MJO amplitude As has been documented previously, under a warming climate the MJO EOFs shift eastward compared to the historical period in CESM2-WACCM due to the eastward extension of tropical Pacific warm sea-surface temperatures Subramanian et al., 2014;Zhou et al., 2020). Due to the shift in the EOFs with warming, we are careful when determining the sign and convention of RMM1 and RMM2 by maximizing the spatial correlation of EOF1 and EOF2 with their historical counterparts.

Sensitivity to the Remote Influence of Periodic EventS Index
The Sensitivity to the Remote Influence of Periodic EventS (STRIPES) index is designed to capture regional sensitivities to the remote influence of periodic events in a single number (Jenney et al., 2019). Here, we use the STRIPES index to quantify the impacts of the MJO (i.e., magnitude and consistency) on the midlatitude circulation at each midlatitude gridpoint. Specifically, we compute the SLP anomalies 0-35 days following every MJO event of a particular phase. If there is a strong influence of the MJO on that grid point, a tilted stripe will be present in a plot of SLP anomalies as a function of lag (0-35 days) and MJO phases 1-8. A variance calculation of these tilted anomalies produces the STRIPES value, where we focus on tilts that correspond to MJO propagation speeds of 5-7 days per phase. A detailed explanation and visualization of the STRIPES index calculation can be found in Jenney et al. (2019). A higher STRIPES value indicates a larger variance in SLP anomalies, and thus, larger remote influence of the propagating MJO. The ensemble-mean STRIPES index is computed as the average STRIPES value across ensemble members.

Graphical Causal Models
Graphical causal models based on Bayesian networks have been successfully used in climate science to gain insights into atmospheric teleconnections (Ebert-Uphoff & Deng, 2012;Kretschmer et al., 2016;Runge et al., 2019) and cross-scale interactions , as well as for climate model evaluation (Nowack et al., 2020;Vázquez-Patiño et al., 2020). These approaches provide a compact visual representation of the salient interactions between a set of random variables by representing the variables as nodes of a Directed Acyclic Graph (DAG) and the direct causal relationships between them as arrows/directed edges. In this study we use a temporal extension (Chu et al., 2005;Ebert-Uphoff & Deng, 2012) of the "PC-stable" algorithm (Colombo & Maathuis, 2014;Spirtes et al., 2000) to efficiently derive a DAG from our model data. However, the interactions identified from this data-driven approach (as opposed to targeted simulations) are only potential interactions and are not guaranteed to be true causal interactions due to reasons such as hidden common causes. Nevertheless, this method provides several advantages over traditional techniques such as correlation analysis and lagged regression by allowing one to easily distinguish direct interactions from indirect ones in a multivariate setting while still accounting for the memory and feedbacks. See Ebert-Uphoff and Deng (2012) and Samarasinghe (2020), for explanations. Barnes et al. (2019) demonstrated that this causal model approach is capable of distinguishing between the tropospheric and stratospheric pathway of MJO teleconnections to the North Atlantic within observations.
In this study, we derive a temporal DAG by representing the MJO, NAO, and the polar vortex and lagged copies of these variables as separate nodes in a graph. The PC-stable algorithm starts with a fully connected undirected graph initially assuming interactions between every pair of variables, however, instantaneous connections are not allowed. Next, a conditional independence test is used to iteratively disprove the assumed interactions. A pair of variables X and Y cannot have a direct causal connection if they are conditionally independent given any subset of nodes S (excluding X and Y) in the graph. If there is such a subset S, the edge between X and Y is eliminated. We statistically test for zero partial correlation using Fisher's Z-test as the conditional independence test with a confidence level alpha of 0.005, thus focusing on the average linear dependencies between variables. Finally, we orient the edges such that the interactions are from the past to the future. Our temporal model uses 12 lags each for the MJO, NAO, and VORTEX variables with each lag being 5 days apart. We discard one time slice to handle initialization issues following Ebert-Uphoff and Deng (2012).
To understand how the MJO-NAO teleconnection changes over time, we investigate the interactions identified by PC-stable for each 40-year long window in the historical and future periods, with a 10-year shift between adjacent windows. Prior to the causal analysis, we detrend each 40-year-long NAO and VORTEX time series by removing a first-order polynomial fit of each calendar day and then smooth the data with a 5-day, backward-looking average. We define binary variables for each MJO phase. These variables indicate an MJO event happening in the phase of interest with "1," and "0" otherwise. Unless otherwise noted, we conduct the causal analysis for pairs of phases to increase sample sizes for each group (i.e., phases 2/3, 4/5, 6/7, 8/1).
For each ensemble member, we derive a separate DAG for each 40-year window. We quantify the robustness of each potential causal connection via the "temporal consistency fraction" (tcf ) defined below.

10.1029/2020GL091168
A fraction closer to 1.0 indicates that the interaction is consistently repeating in the temporal model, and thus, robust within that time frame. Figure 3a shows a summary of the potential causal interactions learned by PC-stable that have a tcf of 0.6 or greater for a sample 40-year window of an ensemble member. We average this fraction over the ensemble members as a metric to distinguish robust causal signals in the CESM2-WACCM model while accounting for internal variability. See Figure S5 for details.

Results
Figure 1 compares composites of anomalous sea-level pressure 15 days following MJO phase 7 events between the historical (1850-1889) and future (2055-2094) period for all ensemble members. The SLP teleconnection pattern over the North Atlantic compares well with that observed following MJO phase 7 (Henderson et al., 2016; their Figure 6), and generally represents the negative phase of the NAO (e.g., Hurrell, 1995). SLP anomalies increase substantially between the historical and future period, suggesting either stronger or more consistent teleconnections between the MJO and North Atlantic under future climate warming. This strengthening is also present when the two periods are more evenly compared using three ensemble members each ( Figure S3) or when alternative lags are considered ( Figure S4).  (Figures 2a and 2b) (e.g., Cassou, 2008;Lin et al., 2009;Mori, 2008). While the location of the teleconnection hotspots appear similar between the two periods, the magnitude of the STRIPES index over these hotspots increases substantially by the end of the 21st century under SSP585 compared to the early historical period (Figures 2c and 2d). That is, consistent with Figure 1, the STRIPES analysis quantifies a strengthening of the MJO teleconnections with warming, with the largest changes occurring over the Gulf of Alaska and the North Atlantic. While the focus of this study is on the North Atlantic, we return to the North Pacific response in Section 5.
Our results support a strengthening of the teleconnection between the MJO and NAO under SSP585, however, the methodology thus far does not allow us to distinguish between the tropospheric and stratospheric pathways. To do this, we compute graphical causal models of the MJO, stratospheric polar vortex, and the NAO (see Section 3 for methodology). An example result is shown in Figure 3a for MJO phases 6/7 within a single ensemble member for the 1970-2009 period. Arrows denote potential causal connections between climate phenomena, pointing in the direction of cause to effect, and numbers denote the time lag of the connection in days. Arrows that loop back on themselves signify temporal autocorrelation. The graphical model in Figure 3a demonstrates that CESM2-WACCM simulates a direct tropospheric pathway from the MJO to the NAO with a time lag of ∼15 days following phases 6/7. Furthermore, CESM2-WACCM also simulates a stratospheric pathway, evidenced by a connection between the MJO and the VORTEX (time lag of 15 days) and then VORTEX to the NAO (time lag of 5 days), resulting in a total time lag of 20 days or so. These connections and time lags are consistent with what is observed over 1979-2016 in reanalyzes .
Graphs such as that shown in Figure 3a are computed for each 40-year period for each ensemble member for each pair of MJO phases. We summarize the results by computing the tcf (see Section 3) of each connection for each time lag, period and ensemble member, and then average across members. In doing so, we obtain an estimate of the consistency of the causal connection for lags of 5-40 days as a function of time, as shown in Figure 3b. Colored lines denote different causal connections (i.e., arrows in Figure 3a) and it is clear that the direct, tropospheric MJO→NAO connection (pink) substantially strengthens over the 20th and 21st centuries. The stratospheric pathway, however, shows no evidence of strengthening over the 21st century, with the MJO→VORTEX connection (purple) weakening mid-20th century and then returning to 1800 values, and the VORTEX→NAO connection (brown) remaining constant over the entire 250 years. That is, the causal models suggest that the changes in MJO teleconnection strength under SSP585 forcings can be attributed exclusively to the tropospheric pathway, with the stratospheric pathway remaining largely unchanged.
Both the MJO→NAO and MJO→VORTEX causal connections show pronounced increases in their tcfs in the middle of the 20th century (Figure 3b, pink and purple lines). This mid-century peak may be related to the strong impact of aerosols on the North Atlantic during this period (Booth et al., 2012). While this explanation may explain the enhancement of the MJO→NAO causal connection, it is less clear for the MJO→VORTEX. We speculate that it is possible that the North Atlantic aerosol forcing, or the multidecadal atlantic variability associated with it, also impacted the stratospheric polar vortex (e.g., Omrani et al., 2014) and thus its connection to the MJO. Additional analysis is required to determine whether this is actually the case.
The STRIPES and graphical causal model analyses provide strong evidence of a strengthening of the tropospheric MJO to NAO teleconnection under climate warming within CESM2-WACCM. At a fundamental level, this could be brought about by (1) a strengthening of the Rossby wave source in the tropics, (2) a strengthening of the wave propagation within the midlatitudes, or both. While a systematic study of the relative importance of each of these mechanisms is outside the scope of this study, Zhou et al. (2020) suggest that the extension of the North Pacific jet with warming can explain shifts in MJO teleconnections over the western United States in CMIP5 and CMIP6 simulations. On seasonal timescales, Fereday et al. (2020) also concluded that an extension of the Pacific and Atlantic jet-streams under RCP8.5 lead to a strengthening of the tropospheric pathway between the El Nino-Southern Oscillation and the NAO. We also find a robust extension of the North Pacific and North Atlantic jet-streams within CESM2-WACCM (contours in Figure 2c). This extension of the jet-stream from the North Pacific, across North America, and into Europe may enhance the waveguide, supporting Rossby wave propagation from the tropics to the North Atlantic (mechanism #2 above). In fact, changes in the STRIPES index align well with the extensions to the jet, especially over the North Atlantic basin (Figure 2c).
On the other hand, enhancement and shift of the Rossby wave source, for example, through changes in the magnitude and location of MJO heating, could also be at play. While Zhou et al. (2020)   precipitation anomalies associated with the MJO increase under SSP585 in CESM2-WACCM (not shown), so further study is required to truly separate these two mechanisms within the simulations.

Discussion and Conclusions
This study is one of the first to provide evidence of strengthened MJO teleconnections to the North Atlantic under SSP585. This result is perhaps surprising, as recent studies have argued that robust increases in tropical static stability may lead to a weakening of MJO teleconnections as the climate system warms Wolding et al., 2017). MJO-to-North Atlantic teleconnections are known to occur via both a stratospheric and tropospheric pathway, and our graphical causal model approach allows us to clearly distinguish between the two. Thus, our results documenting a strengthening of the tropospheric pathway, and theoretical arguments suggesting weakening teleconnections via increased tropical static stability, may be reconciled if the strengthening of the tropospheric pathway is driven predominantly by changes to the waveguide outside of the tropics. To truly disentangle between these competing mechanisms, however, additional model experiments are likely necessary.
The focus of this study was MJO impacts on the North Atlantic, however, our hemispheric STRIPES analysis leaves open questions regarding a strengthening of the MJO-Gulf of Alaska teleconnection in CESM2-WACCM (Figure 2). A widening of the North Pacific jet-stream on its polar flank (Figure 2c) may reduce the Rossby wave reflection there and instead be more conducive to wave propagation and breaking (e.g., Ambrizzi et al., 1995). The mechanism behind this change is still not clear to the authors, given that one might expect the extended North Pacific jet to also shift the SLP anomaly eastward, something we do not see in Figures   of Alaska region when only the MJO heating is modified to its future state in a linear baroclinic model, yet another possible reason for this response. Thus, further work is needed to understand this Gulf of Alaska response.
Our study employs two data-driven approaches for summarizing and quantifying MJO teleconnections to midlatitudes, and perhaps the methods used in this study are of as much interest as the results themselves. The STRIPES analysis allowed us to summarize many lag and phase diagrams typically used to study MJO teleconnections, while the causal model approach helped disentangle the stratospheric and tropospheric pathways. Both methods provide a straightforward way to compare various teleconnection pathways and strengths within observations and climate models, as well as how the teleconnections may change with time.