Mechanisms driving MJO teleconnection changes with warming in CMIP6

Teleconnections from the Madden-Julian Oscillation (MJO) are a key source of predictability of weather on the extended time scale of about 10-40 days. The MJO teleconnection is sensitive to a number of factors, including the mean state dry static stability, the mean flow, and the propagation and intensity characteristics of the MJO itself, which are traditionally difficult to separate across models. Each of these factors may evolve in response to increasing greenhouse gas emissions, which will impact MJO teleconnections and potentially impact potential predictability on extended time scales. Current state-of-the5 art climate models do not agree on how MJO teleconnections will change in a future climate. Here, we use results from the Coupled Model Intercomparison Project Phase 6 (CMIP6) historical and SSP585 experiments in concert with a linear baroclinic model to separate and investigate alternate mechanisms explaining why and how MJO teleconnections over the North Pacific and North America will change in a future climate, and to identify key sources of inter-model uncertainty. We find that decreases to the MJO teleconnection due to increases in tropical dry static stability alone are robust, and that uncertainty in mean state 10 winds are a key driver of uncertainty in future MJO teleconnections. We find no systematic relationship between changes in Rossby wave excitation and the MJO teleconnection. However, we find that models that predict increases (decreases) in the stationary Rossby wave number over the gulf of Alaska also predict stronger (weaker) teleconnections over North America. Uncertainty in future changes to the MJO’s intensity, eastward propagation speed, and eastward propagation extent are other important sources of uncertainty in future MJO teleconnections, although to a lesser degree than uncertainty in the mean 15 state. The overall outlook is a reduction of the boreal winter MJO teleconnection across the vast majority of CMIP6 models, especially over the North Pacific, but with some nuance over North America due to larger sensitivity to expansion of the MJO’s eastward extent.

explore the direct role that changes to the mean state winds have on changing MJO teleconnections to the broader North Pacific and North America region.
In addition to features of the mean state (its static stability and winds), the MJO teleconnection may also respond to changes in the MJO itself. The MJO teleconnection is known to be sensitive to its eastward propagation speed. In a simulation study, 60 Bladé and Hartmann (1995) showed that an eastward propagating heat source is associated with a weaker and smaller wavetrain than a stationary or westward propagating heat source. To first order, this is a linear effect: eastward (westward) propagation effectively embeds the heating in easterlies (westerlies) (see Sect. 2b of Bladé and Hartmann, 1995). Many climate models predict an increase in the MJO's propagation speed with warming (e.g., Rushley et al., 2019), which may contribute to a weakening of the MJO teleconnection with warming.
linear framework is quite simple in that almost all simulated variability is due to the addition of an external forcing, which makes the interpretation of the results very straightforward. Third, because the mean state is maintained by prescribing it at each time step, it is possible to run simulations with the temperature field out of balance with the wind field. This final point allows us to separately and cleanly quantify the uncertainty in the MJO teleconnection change from uncertainty in the thermal 95 structure of the atmosphere versus its mean winds.

Methods
We explore the sensitivity of the MJO teleconnection to changes in the MJO's propagation speed, eastward propagation extent, heating magnitude, and of the mean state dry static stability and winds using simulations with the LBM of Watanabe and Kimoto (2000). The time-integration version of the LBM that we use here is a spectral model that solves linearized equations 100 for vorticity, divergence, temperature, and surface pressure. The model takes two inputs: a mean state and a forcing (which is added to the simulated perturbation field at each time step).
The model equations are described in appendix B of Watanabe and Kimoto (2000). The model is numerically damped with biharmonic horizontal diffusion, for which we use an e-folding time of 10 minutes for the smallest resolvable waves, and with linear drag, with a damping time scale of 0.5 days applied to the lowest 3 model levels, 1 day for the top two model levels, and 105 20 days for all levels in between. We use this strong damping to inhibit the growth of baroclinic waves, particularly because in some simulations we intentionally prescribe mean state winds that are not in balance with the mean state temperature. We use a horizontal truncation of T42 (roughly 2.8 • horizontal resolution) and 20 vertical levels.
We will make the convention for composite analysis that "historical" and "future" mean states refer to January conditions during 1984-2014, and 2071-2100, respectively, of the historical and SSP585 simulations from models participating in CMIP6.

110
SSP585 is a relatively high atmospheric greenhouse gas emissions scenario following shared socioeconomic pathway 5 (SSP5; O'Neill et al., 2017) and representative concentration pathway 8.5 (Moss et al., 2010). For each CMIP6 model mean state combination, we construct an ensemble of 30 separate LBM simulations, each using a different January. We use interannual ensembles in this way to isolate the climate change signal apart from noise due to internal variability.
We restrict our analysis to those 29 CMIP6 models for which there was monthly mean data for the variables required as input 115 by the LBM on the data node at https://esgf-node.llnl.gov/search/cmip6/ at the time of access. We use the r1i1p1f1 variant for all models except for GISS-E2-1-G (r1i1p1f2), HADGEM3-GC31-LL (r1i1p1f3), UKESM1-0-LL (r1i1p1f2), and MIROC-ES2L (r1i1p1f2). Lastly, the SSP585 simulation of CAMS-CSM1-0 was only carried out until 2099, and so for this model we use only 29 ensemble members.
We force the LBM with a propagating idealized horizontal dipole heating, which is preferable to a realistic heat source (for 120 example, one obtained from composites of observed MJO events), since it allows easily manipulation of the MJO's intensity and propagation characteristics. The heating is designed to be MJO-like, with the rate of change of perturbation temperature given by where T is the perturbation temperature of the forcing, φ and λ are the latitude and longitude, respectively, z is the model 125 vertical coordinate with subscripts bot and top referring to the lowest and highest model level, respectively, A is the heating amplitude, k is the zonal wavenumber, ω is the temporal frequency of the heating, and t is the time. Unless otherwise noted, we set A = 0.1 K day −1 , σ = 0.1, k = 1.9, ω corresponding to a temporal period of 40 days, b = 130 • , and c = 45 • . Figure 1 shows the heating in the mid-troposphere where it is a maximum for t = 10 days, 20 days, 30 days, and 40 days. the other with the mean state given by future values (see the "sub-experiment" column in Table 1). For example, in testing the sensitivity of the teleconnection to the mean state winds, we run two sets of experiments: one in which the dry static stability is held constant at historical values and the other with future dry static stability. The use of the many CMIP6 mean states permits the quantification of uncertainty in the future MJO teleconnection due to "model uncertainty" (e.g., Lehner et al., 2020)  Specifically, the gaussian multiplier of the perturbation forcing is the control gaussian with the maximum value sustained for an additional 20 • to the east before decreasing. Lastly, we are informed by Bui and Maloney (2018), who find that the subset of CMIP5 models that produce MJOs validating best against present-day observations tend to simulate changes of the MJO's precipitation amplitude between -10% to + 20% between the historical and future climates under RCP8.5 forcing. In 155 the experiment testing sensitivity to MJO amplitude, we use perturbed MJO heating amplitudes of 0.9 and 0.12 K day −1 , corresponding to a 10% reduction and 20% amplification of the control forcing.
With the described setup, some of the simulations became unstable. We readily admit this is an inescapable trade-off of having assumed a forced linear model framework; despite its advantages for rapidly sampling and intercomparing first-order effects, the neglect of secondary nonlinear buffering mechanisms can corrupt some integrations. Pragmatically, we thus omit 160 all LBM simulations in which, at any location, the maximum value during the second half of the integration is more than ten times the maximum value during the first half. Additionally, we omit all LBM simulations using a specific CMIP6 model's mean state if 10 or more of the 30 constituent ensemble simulations for any of the mean state combinations meets the instability criterion described above, leaving 24 of the 29 CMIP6 models for analysis. Figure A1 in Appendix A shows, for each basic state combination, the fraction of ensemble members that became unstable. The last two rows of Table 1 list the number of 165 total simulations run for this study (minus those from the five CMIP6 models we omit) and the number of simulations we omit from our analysis for each experiment (which is reassuringly less than 5% of the total simulations for each experiment). Note that the "number of unique simulations" field is blank for the mean state winds feature (second row) because the simulations are shared with that of mean state dry static stability (first row).

Mean state
We conduct experiments with the experimental setup described by Table 1.
Reassuringly, the propagating MJO-like forcing in the model excites a plausible time-varying response in the extratropics. relative to a typical observed midlatitude value is an expected consequence of the smallness of the amplitude used to force the model, which we kept small to try to minimize the number of unstable simulations.
We next introduce a scalar metric of MJO teleconnection strength that is appropriate to visualize in map form, recognizing that the magnitude of peaks and troughs of the ensemble mean response to the propagating forcing in Fig. 2a is one way to 180 quantify the magnitude of the consistent (ensemble mean signal that stands apart from the noise due to internal variability) teleconnection strength. Given that the amplitude of an infinite sine or cosine wave is the square root of twice its variance, we calculate the teleconnection strength at each point as the square root of twice the variance of the ensemble mean meridional wind at 850 hPa, ignoring the first 10 days of each simulation to avoid sampling conditions prior to the establishment of representative spread between interannually varying ensemble members. We will refer to this value as the "amplitude" of the  surface temperature variability across North America via low-level temperature advection (Seo et al., 2016), and may thus be 190 more relevant for quantifying changes to the MJO's impact on near-surface midlatitude weather in a future climate.
The baseline multi-model MJO teleconnection looks reassuringly plausible. Figure 2b shows the multi-model mean MJO the southeastern and northwestern region, with a minimum over the central region, is consistent with the pattern of the observed boreal winter MJO influence on North American temperature (e.g., Jenney et al., 2019;Zhou et al., 2012).
The remainder of Fig. 2 reveals an overall reduction of MJO teleconnectivity due to future mean state effects, decomposed into separate contributions from stability vs. winds. Figure 2c-e shows the amplitude of the modeled MJO teleconnection 200 but for the set of simulations performed by perturbing the mean state. The forcing used in all simulations is the same. The difference between the top and bottom rows is the sensitivity of the MJO teleconnection to the change in the thermal structure of the atmosphere (i.e., its dry static stability), alone. The pattern of the MJO teleconnection is the same between each top panel and the panel below it, but the teleconnection is weaker in the simulations with future dry static energy because the higher static stability weakens the MJO circulation and subsequent Rossby wave source (Sardeshmukh and Hoskins, 1988) 205 that results from the prescribed thermal forcing. The difference between the left column and the right column is the sensitivity of the MJO teleconnection amplitude to the simulated change in mean winds between the historical and future climates of CMIP6. For most regions, the change to the mean state winds between the historical and future climates leads to a reduction of the MJO teleconnection amplitude. This reduction is particularly apparent over the central Pacific. Other studies have noted an increase in the MJO's boreal winter influence near California (Zhou et al., 2020), attributing this to an eastward extension 210 of the Pacific jet. Here, we find a similar strengthening over this region. In Fig. 2c-e, we include a white line highlighting the 5 ×10 − 3 m s −1 contour to help show this strengthening, which is apparent as an eastward extension of this contour line towards the California coast.
The difference between panels (b) and (e) of Fig. 2 is the sensitivity of the MJO teleconnection to a change in the full mean state, neglecting any changes to the MJO. In isolation (apart from changes to the MJO), changes to the mean state weaken 215 the multi-model mean teleconnection amplitude almost everywhere due to both an intensification of the dry static stability and changes to the winds.
We now explore the inter-model spread of changes to the MJO teleconnection due to changes in the mean state. Figure  Our first main finding is strong corroboration of the hypothesis that the increase of the tropical dry static stability is a robust thermodynamic response of the climate to warming that tends to reduce the MJO's teleconnections in almost all models. This 230 is the first multi-model confirmation of this hypothesis, which first appears in the single-model study of Wolding et al. (2017). There is relatively good agreement in the teleconnection amplitude change between model mean states, particularly over the North Pacific. Comparison between panels a and c of Fig. 3 shows that there is greater spread downstream of the tropical heat source over North America, although reasons for this are unclear.
Despite agreeing on effects of stability, models disagree wildly on the effects of future wind changes. Figure 3b,d show the 235 change in the teleconnection amplitude due to the change in the wind between the historical and future period. The multi-model mean response of the MJO teleconnection amplitude over the North Pacific is a modest weakening, with a majority of the region also experiencing a weakening. However, the inter-model spread is relatively large, with models showing changes between −15 to 7 % K −1 multi-model mean warming. Over North America, the multi-model mean change in the teleconnection due to the wind is a small strengthening, with most of the region showing weaker teleconnections. But again, the inter-model spread is 240 large (−9 to 12 % K −1 ) with the multi-model mean is indistinguishable from zero as was found in Zhou et al. (2020).
In summary, three key results emerge concerning how mean state changes in a future warmer climate will contribute to changing MJO teleconnections to the North Pacific and North America. First, we find robust support for the expectation that increases in tropical dry static stability will contribute to weaker MJO teleconnections. Second, we find that changes to the mean state winds can contribute to very large changes in the MJO teleconnection; that is, large enough to as much as double or 245 negate the weakening that is expected due to increasing tropical dry static stability. Third, there is a large amount of inter-model spread in LBM-simulated changes to the MJO teleconnection that result from changes to mean state winds, which suggests that uncertainty in the response of the mean state winds to warming dominates uncertainty in modeled changes to the MJO teleconnection.
It is thus important to better understand the cause of teleconnection changes associated with mean state wind changes in the 250 CMIP6 models. Mean state winds can increase or decrease the strength of MJO teleconnections via an intensification of Rossby wave excitation and/or through mean state changes that permit more Rossby wave energy propagation into these regions. The following analysis will argue that the latter-wind-induced changes to Rossby wave ducting-are more important than changes to wave excitation in producing the inter-model teleconnection spread.
To show this, we begin by quantifying the strength of Rossby wave excitation using the linearized equation for the "Rossby 255 wave source" (Sardeshmukh and Hoskins, 1988), S , a term that represents the generation of large-scale vorticity, where ζ is the absolute vorticity, v χ is the horizontal divergent wind, the overbar represents mean state quantities, and the prime represents anomalies. In the subsequent analysis, we use quantities at 200 hPa to calculate S .
We do not find a systematic relationship between the difference in S and the difference in the mean teleconnection amplitude teleconnections. If changes to Rossby wave excitation were to play a key role in the difference in the teleconnection amplitude between simulations with historical and future mean state winds, then it would be reasonable to expect either a weaker or similar-sized link over North America than over the North Pacific, since most Rossby waves excited by the MJO propagate first over the Pacific before reaching North America. However, the relationship between the difference in |S | and the teleconnection amplitude difference appears stronger over North America (r = 0.43) than over the North Pacific (r = 0.27). Thus we conclude 270 that the change in |S | is not a primary mechanism that explains how mean state winds affect the teleconnection amplitude difference in our LBM simulations.
In addition to playing a key role in Rossby wave excitation, mean state winds also determine the path that Rossby waves take as they propagate (e.g., Karoly, 1983;Hoskins and Ambrizzi, 1993). Analyses of the stationary Rossby wave number have been useful for assessing how the mean state winds affect Rossby wave propagation (e.g., Henderson et al., 2017;Karoly, 275 1983;Tseng et al., 2020b;Wang et al., 2020). The stationary Rossby wave number on a Mercator projection (to account for spherical geometry) is defined as where β M is the meridional gradient of absolute vorticity on a Mercator projection, 280 and u M = u/ cos θ is the mean zonal wind divided by the cosine of latitude (θ). In Eqs. (3) and (4), a is the earth's radius and Ω is the earth's rotation rate. Hoskins and Ambrizzi (1993) showed that K s can be used to understand Rossby wave propagation: waves turn towards regions of higher K s , are generally reflected away from regions where the zonal Rossby wavelength is equal to K s or where β M ≤ 0, and are either dissipated beyond or reflected from regions where u ≤ 0. Regions of K s maxima, which tend to occur in strong westerly jets, act as Rossby waveguides. Figure 5a shows K s computed from the multi-model  results. Regions of strong positive correlation suggest a systematic (common across models) positive relationship between the local change in K s and the regional (i.e., across North Pacific or North America) teleconnection amplitude change. Similarly, regions of strong negative correlation suggest a systematic inverse relationship between the local change in K s and the regional teleconnection change. We include boxes outlining the North Pacific and North America regions for reference. Correlations linking regional changes in the local K s with the regional mean teleconnection amplitude change are generally weaker for the 295 North Pacific than for North America. We thus focus on North America, but point out that the two maps have similar patterns, suggesting that changes to the winds leading to stronger teleconnections over North America may be in part due to changes in the winds that also support stronger teleconnections over the North Pacific. Figure 5c shows a region of high positive correlation over the Gulf of Alaska. Here, models exhibiting an increase in K s tend to have an increase in the teleconnection over North America, and vice versa. Additionally, there is a region of negative correlation just to the south of this area. We interpret this 300 as a slight northward shift and broadening along the northward flank of the Pacific waveguide, and subsequent increase in the number of Rossby waves traveling northward towards North America, in the models that have stronger MJO teleconnections.
This view is consistent with the work of Tseng et al. (2020b), who demonstrate how increased K s over the Gulf of Alaska and resultant increases in northward Rossby wave propagation out of the Pacific waveguide during La Niña years leads to an enhanced MJO teleconnection over North America. 305 We thus identify a systematic mechanism that may explain how mean state wind changes lead to changes to teleconnectivity over North America across CMIP6 models. This suggests that increased confidence in how winds over this region will change 13 https://doi.org/10.5194/wcd-2021-9 Preprint. Discussion started: 10 February 2021 c Author(s) 2021. CC BY 4.0 License.
in a future climate may help reduce uncertainty in future estimates of the change in the strength of MJO teleconnections to North America. We would like to point out, however, that changes in the winds over this region alone are likely not the sole important change. That is, for some models, it is possible that changes to S may be important, as well.

MJO intensity and propagation characteristics
We now investigate how changes to MJO intensity and propagation characteristics affect the MJO teleconnection by conducting simulations in which the mean state is held constant, and either the propagation speed, eastward propagation extent, or intensity of the idealized thermal forcing are perturbed to upper and lower bounds of expected end-of-century changes to the MJO.
Instead of using mean states from all CMIP6 models, we use a subset of mean states from 10 models to minimize utilization 315 of computational resources. For these experiments, we chose CMIP6 models that produced a minimum number of unstable ensemble members in experiments used in Sect. 3.1 (see Fig. A1). linearly with the magnitude of the forcing. We verified this with simulations using a 20% increase in the magnitude of the forcing for two different model basic states. For the remainder of the 8 models used in these experiments, we thus did not actually run simulations (see Table 1), because the magnitude of the extratropical response to the perturbed forcing amplitude is independent of the mean state. If full non-linearity had been considered, the variations across the models' mean state could have also been important in the response of the teleconnection.

325
Unlike perturbations to the heating intensity, the responses of the MJO teleconnection to perturbations to the eastward propagation speed (Fig. 6b,e) and extent (Fig. 6c,f) are linearly sensitive to the mean state. In general, increases to the propagation speed produce modest decreases in the teleconnection strength over the North Pacific, with changes over North America centered on and straddling zero change. Increasing the eastward extent of the propagating MJO thermal forcing in the LBM simulations increases the teleconnection amplitude for all mean states used and has the potential to produce the largest increases 330 in the MJO teleconnection amplitude, more so over the North America region. Figure 6 highlights two sources of uncertainty that lead to uncertainty in the response of MJO teleconnections to changes to the MJO itself: uncertainty in changes to the MJO and uncertainty in the mean state. Uncertainty in the MJO's response to warming (here, between a 10% overall decrease or 20% overall increase between the present day and end of the 21st century), leads to a range of about 8% K −1 uncertainty in the change to the MJO teleconnection to the North Pacific and North America.

335
For the propagation characteristics, both uncertainty in the mean state and uncertainty in the MJO's response to warming leads to uncertainty in the teleconnection response to warming, with a range of about 5% K −1 from changes to the MJO's propagation speed and up to about 10% K −1 from the MJO's eastward propagation extent. Additionally, for each individual MJO propagation or intensity characteristic, the uncertainty in the MJO teleconnection resulting from a perturbed MJO is on the order of the model uncertainty that results from perturbations to the mean state alone (Fig. 3). This further challenges the 340 notion that it may currently be possible to have confidence in future projections of MJO teleconnections. As expected, given that the increase of the tropical dry static stability is a robust thermodynamic response to surface warming, we verify that the increase in the tropical dry static stability with warming produces robust decreases in the LBM-simulated MJO teleconnection amplitude. On the other hand, we find relatively large spread in the LBM-simulated teleconnection response to changes in the mean state wind. Nonetheless, when the mean state changes are considered together (Fig. 7), they lead to robust decreases in the MJO teleconnection to the North Pacific. Over North America, the sum of mean state changes leads 385 to a weaker MJO teleconnection for most CMIP6 model mean states, although with many of these models hovering close to zero change.
Investigation into a potential systematic link between changes to the wind and changes to MJO teleconnection revealed a possible connection related to Rossby wave propagation over the eastern Pacific, where models that display increases (decreases) in the stationary wavenumber over the Gulf of Alaska also show increases (decreases) in the teleconnection amplitude 390 over North America.
We find that for perturbations to MJO characteristics (i.e., its intensity, eastward extent, or propagation speed), the uncertainty in the MJO teleconnection that results is on the order of the uncertainty from the mean state changes alone. In particular, increases in the MJO's heating intensity and eastward extent have the potential to lead to large increases in the MJO teleconnection, with the sensitivity to extent being especially strong for teleconnection over North America for some mean states.

395
Overall, over the North Pacific, despite a relatively large range of −21 to 0.6 % K −1 , which depends on the CMIP6 model used as the mean state in the LBM simulations, in general it appears as though the MJO teleconnection to the North Pacific will likely weaken with warming. Over North America, despite a modest multi-model mean weakening result from the influence of the mean state alone, including MJO changes results in many CMIP6 models showing a range of possible teleconnection changes that is either close to or near-zero.
We have thus identified sources of uncertainty leading to uncertainty in projections of MJO teleconnections in a future warmer climate. To first order, model spread in projections of future winds is the biggest source of uncertainty in projecting the response of MJO teleconnections to warming. For teleconnections to North America in particular, we have identified a region over the eastern Pacific where changes to the winds appear to explain much of the inter-model spread. We speculate that increased confidence in changes to the mean flow over the Gulf of Alaska would help reduce uncertainty in teleconnection 405 changes over North America. Uncertainty in the mean state also affects how MJO teleconnections respond to increases in the eastward extent of MJO convection. Lastly, uncertainty in how the MJO itself will respond to climate change is nontrivialleading to a range of about 5% K −1 for any given mean state (see Fig. 7)-but interestingly secondary to inter-model uncertainty in the wave ducting properties of the basic flow.
Over North America, a possible future outcome is that of the median, near-zero area-mean total teleconnection change.

410
In this case, the pattern of the teleconnection may change, despite zero average change over the region as a whole. In this case, a reduction of the uncertainty in future projections of the mean state winds will still help reduce uncertainty in the MJO teleconnection change, due to the tight dependence of the stationary wave pattern on the mean winds.
We readily admit limitations of our choice of a linear model to explore how the MJO teleconnection responds to changes in the mean state and to the MJO itself. While this allows us to easily diagnose causal relationships, and span a previously sparsely 415 sampled causative parameter space, the linear framework is missing nonlinear relationships between the mean state and the MJO forcing. For example, the mean state exerts a strong control on MJO propagation characteristics (Jiang et al., 2020). The MJO may also play an active role in shaping the midlatitude mean state winds, for example through the modulation of the local Hadley cell (Lyu et al., 2019). Such feedbacks are not simulated in the LBM. Nonetheless, previous work has also shown that the MJO teleconnection is, to first order, linear (Mori and Watanabe, 2008). For now, we suggest this work as a compelling first 420 step towards understanding the intriguing model spread in future MJO teleconnection strength in the extratropics, and another reason to elevate the ongoing community quest to reduce uncertainty of the mean midlatitude circulation's response to climate change.
Code and data availability. Appendix A Figure A1 shows the number of simulations that became unstable (and hence omitted from our analysis; see Sect. 2), and indicates the CMIP6 models for which all simulations were excluded from our analysis for having 10 or more unstable ensemble members for a particular basic state combination.