Nonlinearity and Asymmetry of the ENSO Stratospheric Pathway to North Atlantic and Europe, Revisited

Nonlinearities and asymmetries of El Niño Southern Oscillation (ENSO) stratospheric pathway to the North Atlantic and Europe are examined in large ensembles conducted with fully coupled climate models during wintertime. The analysis is centered on historical experiments of the Max Planck Institute Grand Ensemble (MPI‐GE, 95 members) and expanded to six other ensembles of more limited size. In MPI‐GE, significant responses are obtained for each ENSO phase and three different intensities (weak, moderate and strong). Overall, linear relationships are found for either El Niño or La Niña key diagnostics that characterize the pathway. These relationships are generally weaker for the cold La Niña than for the warm El Niño so that asymmetries between them develop as the events intensify. Specifically for strong events, the extra‐tropical North Pacific and stratospheric responses are asymmetric, with larger responses for El Niño. In addition, the stratospheric asymmetry in strong events seems to contribute to the asymmetry in strong events in the North Atlantic—Europe response in the troposphere in late winter. The extra‐tropical North Pacific response shows general agreement between MPI‐GE and the other large ensembles. However, this agreement is not as large when other parts of the pathway are compared. Relatively high inter‐model response spread confirms the typical model uncertainty found when examining atmospheric circulation responses which include the stratosphere in state‐of‐the‐art climate models.


Introduction
Equatorial Pacific Sea Surface Temperatures (SST) are known to exhibit large-scale warming and cooling, the most distinct component of the climate phenomenon named El Niño Southern Oscillation (ENSO).The dynamics of ENSO and its teleconnections have been extensively studied because of their global impacts on ecosystems and society (McPhaden et al., 2021;extensive review).The warm and cold ENSO phases, known as El Niño (EN) and La Niña (LN) events respectively, alternate on time scales of a few years and so do their extreme global impacts (e.g., Cai et al., 2021).Intriguing questions are to what extent the EN and LN teleconnections are symmetric between each other and each of them linear with respect to the magnitude of their corresponding SST anomalies.ENSO is indeed known to be diverse in magnitude, location and duration of anomalies ("ENSO diversity," Capotondi et al., 2015).Therefore, nonlinearities and asymmetries are intrinsic to ENSO itself and can lead to nonlinear and asymmetric teleconnections.Here, we search for nonlinearities and asymmetries along the ENSO stratospheric pathway, taking advantage of Large Ensembles (LEs) of simulations with fully coupled climate models.Central to our analysis is MPI-GE, a grand ensemble of historical and scenario experiments (Maher et al., 2019).Findings from MPI-GE are supplemented by analysis of 6 LEs (30 or more members by model) performed with Coupled Model Intercomparison Project Phase 6 (CMIP6) models.To assess the nonlinearity and Abstract Nonlinearities and asymmetries of El Niño Southern Oscillation (ENSO) stratospheric pathway to the North Atlantic and Europe are examined in large ensembles conducted with fully coupled climate models during wintertime.The analysis is centered on historical experiments of the Max Planck Institute Grand Ensemble (MPI-GE, 95 members) and expanded to six other ensembles of more limited size.In MPI-GE, significant responses are obtained for each ENSO phase and three different intensities (weak, moderate and strong).Overall, linear relationships are found for either El Niño or La Niña key diagnostics that characterize the pathway.These relationships are generally weaker for the cold La Niña than for the warm El Niño so that asymmetries between them develop as the events intensify.Specifically for strong events, the extra-tropical North Pacific and stratospheric responses are asymmetric, with larger responses for El Niño.In addition, the stratospheric asymmetry in strong events seems to contribute to the asymmetry in strong events in the North Atlantic-Europe response in the troposphere in late winter.The extra-tropical North Pacific response shows general agreement between MPI-GE and the other large ensembles.However, this agreement is not as large when other parts of the pathway are compared.Relatively high inter-model response spread confirms the typical model uncertainty found when examining atmospheric circulation responses which include the stratosphere in state-of-the-art climate models.

Plain Language Summary
The alternating sea surface temperature warming (El Niño) and cooling (La Niña) events of the Equatorial Pacific Ocean are known to affect weather and climate worldwide.Yet, impacts of El Niño and La Niña may differ greatly between events, because of differences between the intensity and pattern of the events, and the confounding influence of other weather events.Here we address how the stratospheric polar vortex and the climate of the North Atlantic and Europe regions are impacted by El Niño and La Niña and we ask to what extent these impacts are equal in magnitude and opposite in sign.Disentangling the origin of the diversity in these impacts can help improve our ability in predicting European seasonal climate.By using a large number of climate model simulations, we find that overall larger responses to warm El Niño than to cold La Niña are typical in the extra-tropical troposphere and stratosphere during Northern winter.MANZINI ET AL.

10.1029/2023JD039992
2 of 12 asymmetry of the ENSO stratospheric pathway and the atmospheric extra-tropical dynamical processes involved is particularly relevant for improving seasonal regional climate predictions and their impacts.
The ENSO stratospheric pathway to the North Atlantic and Europe is a teleconnection active during Northern Hemisphere winter.The link of El Niño or La Niña SST anomalies to climate variability of the North Atlantic-European regions can be mediated by the stratospheric polar vortex because the latter is sensitive to the anomalies in planetary waves triggered by ENSO (Butler et al., 2014;Cagnazzo & Manzini, 2009;Calvo et al., 2017;Ineson & Scaife, 2009;Iza et al., 2016;Manzini et al., 2006).Although several tropospheric ENSO teleconnections have been documented in the literature (see Yang et al., 2018 for a review), the role of the stratospheric pathway in the ENSO teleconnection to Europe is crucial (Domeisen et al., 2019 and references herein).Indeed, several authors have suggested that the North Atlantic-European response to EN can only be detected in winter when the stratosphere is active (Butler et al., 2014;Ineson & Scaife, 2009;Polvani et al., 2017).
The North Pacific region is a known key region for the stratospheric pathway.Early works assessing the extratropical response to ENSO typically found a nonlinearity in sea level pressure and geopotential height response in the North Pacific-North American regions (Hoerling et al., 1997(Hoerling et al., , 2001)).However, more recent analysis seems not to support the early works (Deser et al., 2017), raising the question about whether or not the early findings are simply a manifestation of internal variability.Analysis based on experiments with Atmospheric General Circulation Models (AGCM) forced by idealized SST anomalies indicate larger responses to warm than to cold SST anomalies, but overall symmetric extra-tropical response for a "realistic blend" of Central Pacific and Eastern Pacific events (Frauen et al., 2014).However, even if responses between EN and LN are symmetric overall, the responses to their extreme manifestations may not be so.
The ENSO response along the stratospheric pathway has been reported to be linear in analysis of a LE performed with an AGCM forced by observed historical SST, at least for the limited number of events considered (Weinberger et al., 2019).By performing idealized experiments with another AGCM with imposed SST anomalies (SSTAs), such that the location of the maximum magnitude of the anomalies is controlled, the two ENSO phases have instead been found to impact differently the stratospheric pathway, with the warm phase providing for substantially stronger responses than the cold one as the magnitude of the SSTAs increased (Trascasa-Castro et al., 2019).A stronger response to the warm than to the cold ENSO phase is consistent with the finding of Frauen et al. (2014), and we note that these latter two works used related AGCMs.An asymmetric stratospheric response to strong EN and LN events is also found in Hardiman et al. (2019).In this case, however, the stratospheric response to EN changes sign in mid-winter, so that the December to January stratospheric response to cold LN anomalies is stronger in magnitude.Differently from the previous studies, Hardiman et al. (2019) analyzed simulations performed with an initialized fully coupled climate prediction model (including atmosphere-ocean coupling).These divergent findings, found with different model configurations, arise the questions of how sensitive the simulated ENSO responses are to the atmospheric model used, explicitly simulating atmosphere-ocean coupling and initialization.Motivated by these questions, we focus on extending the assessment of the ENSO responses to a large number of free running climate models, the fully coupled atmosphere-ocean models which account for the air-sea interactions that are key to the development of ENSO teleconnections.By analyzing free running climate model simulations, however, we do not address the role of initialization, a topic left for future investigations.
Taking advantage of the availability of MPI-GE and CMIP6 LEs, we assess the sensitivity of the ENSO response as the magnitude of the SSTA increases, for each phase.To this end, ENSO events are grouped in three intensity categories (weak, moderate and strong) according to the magnitude of SSTAs.The use of LEs is crucial, as it allows for a relatively large number of cases in each group, especially for LN strong events.In the course of this work, the term nonlinearity is used for the deviation from a linear function of the intra-phase responses (either EN or LN), while asymmetry is reserved for the inter-phase comparison (between EN and LN).Symmetry is defined when responses to EN and LN are opposite in sign and equal in magnitude for the same magnitude of the SSTA.

Large Ensembles
Monthly means of atmospheric variables and SSTs from MPI-GE and 6 Large Ensembles (LEs) of historical experiments are analyzed.MPI-GE (Maher et al., 2019) is a grand ensemble of a set of experiments, including scenarios, performed with the Max Planck Institute (MPI) Earth System Model 1.1.(MPI-ESM1.1,Giorgetta et al., 2013).Here we use the 95 members of historical experiments, whose stratosphere-ocean variability linkages we have previously analyzed in Ayarzagüena et al. (2021).Given that the CMIP6 historical experiment runs from 1850 to 2014, the MPI-GE historical time series (ending at 2005) are extended to 2014 using the RCP4.5 scenario experiment, also part of the MPI-GE.The 6 LEs are from CMIP6 models for which at least 30 members of historical experiments were available at the time of this study (see Table 1 name of the CMIP6 models for which there are LEs).Each CMIP6 LE comprises ∼30 members, so that the minimum (∼25 members) suggested by Weinberger et al. (2019) for the detection of potential nonlinearities is satisfied.
Our analysis focuses on the extended Northern winter, November to March monthly mean fields.For each ensemble, field anomalies are constructed by subtracting the respective ensemble mean, grid-point by grid-point and month-by-month.For LEs, the large number of members assures that subtracting the ensemble mean is effective in removing the mean impact of external forcings (natural and anthropogenic) and the seasonal cycle.However, impacts of external factors on variability can remain in the anomalous time series.The time dependence and robustness of the statistics over the historical experiments are evaluated considering, for each year, 30-year centered sliding periods (135 in total).In doing so, we evaluate statistics for each sliding period, pooling together the 30-year time series from all available members.Therefore, for each period there are 2850 (95 × 30) winters for MPI-GE, and between 870 and 1500 winters for the CMIP6 LEs.

ENSO Indices and Categories
ENSO events are selected using the Niño3.4index by averaging SST anomalies over (170°W-120°W, 5°N-5°S) and from November to February (NDJF).NDJF is used in order to capture the Northern winter maximum in Niño3.4SST anomalies (Hurwitz et al., 2014).For each 30-year centered sliding period, the SST anomalies of each member and each year belonging to the period are thereafter normalized by the standard deviation (STD) obtained from the ensemble of this 30-year time series.Finally, the equal intensity categories key to our analysis are defined by separately grouping the LN and EN events according to the magnitude of their normalized anomalies: weak (0.5-1 STD), moderate (1-2 STD) and strong (>2 STD) events.Sample sizes for the different phases and categories of ENSO events are generally large for both MPI-GE and the CMIP6 LEs, with a few exceptions in the case of the strong LN events, CMIP6 LEs (Table S1 in Supporting Information S1).Within these few exceptions, only in the case of MIROC6 we do not report statistics for strong LN events, because these events have not been found to occur for the large majority of the 30-year periods.
For MPI-GE, along the 30-year sliding periods, the Niño3.4STD of the SSTAs are rather stable, varying by ∼6% (estimated by the 25th-75th percentile of the STD distribution) around the median, 0.83 K.The latter compares realistically with 0.89 K, the Niño3.4SSTA STD, estimated using detrended HadlSST (1870-2021, Rayner et al., 2003).A small positive trend (0.06 K/100 years) in the Niño3.4SSTA STD suggests only a weak impact of external forcing on the variability during the historical period (Maher et al., 2018).For the CMIP6 LEs as well, variations and trends of the Niño3.4SSTA STD are relatively small (Table 1).
The selection of ENSO events using the Niño3.4index does not explicitly distinguish between ENSO diversity.However, our categories might correspond to different patterns in SSTAs, if the modeled ENSO captures at least partially the well-known asymmetric properties of the Equatorial Pacific SSTAs, namely that extreme EN events have the warmest SSTAs in the East Pacific (Niño3 region) and extreme LN events have the coldest SSTAs in the Central Pacific (Niño4 region) (Cai et al., 2021;Frauen et al., 2014).To test if that is the case for MPI-GE, we have compared the distributions of Niño3 and Niño4 SSTAs, normalized by the Niño3.4STD, for the strong events occurring in the last 30-year period  including all 95 members.We find that the most extreme (warmest) SSTAs are located in the Niño3 region for the EN strong events, and, conversely, for LN strong events the coldest SSTA are located in the Niño4 region (Figure S1 in Supporting Information S1).

Metrics of the Components of the Stratospheric Pathway
The stratospheric pathway is diagnosed by three indicators of circulation responses, each representing a key step in the pathway dynamical processes.The time evolution occurring along the pathway is taken into account by a 1-month time lag between the indicators: • ZG500 index, to evaluate the tropospheric response in the North Pacific, the well-known deepening and eastward enlargement of the Aleutian Low during EN (Hoerling et al., 1997).This index is defined as the 500 hPa geopotential height anomalies area-averaged over (150°E-120°W, 30°N-65°N).In this fashion, our index is broader than the NPI (North Pacific Index) to account for the EN signature.To capture the early stage of the tropospheric forcing of the stratospheric pathway, the geopotential height anomalies are averaged over December and January (DJ).• Stratospheric Polar Vortex index (SPV) index, a measure of the stratospheric response in terms of zonal mean zonal wind anomalies averaged over 55°-70°N at 10 hPa.The stratospheric response to ENSO is known to be generally stronger in mid to late winter, thus the zonal-mean zonal wind anomalies are averaged over January and February (JF).• NAE index, to represent the lagged response in sea level pressure (SLP) over the North Atlantic and Europe (NAE).The index is defined by the SLP difference between the Atlantic-Europe midlatitudes (60°W-60°E, 30°N-55°N) and Arctic (all longitudes, 65°-90°N).The restriction to the Atlantic-Europe sector of our index is consistent with the strong link of this sector to the Arctic anomalies (Deser, 2000).Considering the typical 1-month lag of the surface response to a stratospheric vortex anomaly, for this index the SLP anomalies are averaged over February and March (FM).
Unless specified, for each 30-year centered sliding period we use normalized anomalies (over the ensemble of the 30-year time series, as for the ENSO indices) to show the relative magnitude of the responses independently on the index and models and better discern normal versus unusual values.In doing so, we compare strengths of relationships along the stratospheric pathway and across models relative to their respective variability.Responses are given by the composited anomalies, stratified by ENSO phase and category.For each 30-year period, ENSO phase and category, the 95% confidence interval of the composited indicators is based on the fitted normal probability distribution to the original data values.The two-sided Student's t-test is used for calculating the statistical significace for the composited maps.

MPI-GE
The time series of composited anomalies (Figure 1, upper row) shows that the main features of the responses are stable along the historical period, but there are exceptions, most notably for the responses to LN strong events, which shows some vacillations at times relatively large (ZG index before 1960 and NAE index after 1980).Very slight overall trends for the responses of ZG500 and NAE indices to EN strong events, can also be detected.
Consequently, there can be rare occasions for which claims based even on a large ensemble size (95 members) but one 30-year period only may not be fully definitive.
Main features however emerge rather explicitly, for each metric of the components of the stratospheric pathway.
For each category, responses to EN and LN are oppositely signed, meaning that overall, there are significant responses to the opposite phases of ENSO, as expected.The ZG500 responses to EN categories are distinct, their strength increasing with the SSTA magnitude.Less so but still distinguishable by category are also the SPV and NAE responses for moderate and strong EN events.Instead, for each metric, LN responses overlap so that there is no clear distinction in the responses across the categories, although they do share the same sign.
The ZG500 and SPV responses to EN events are linear as the Niño3.4SSTA increases.A linear fit of the intra-phase EN magnitude responses can indeed be already found by only including 50% of the composite means (see the size of the bars, Figures 1d and 1e).Although Figure 1 (bottom row) shows that a linear relationship, albeit with a much smaller rate of increase, can be drawn for LN responses as well, it cannot be concluded in this case that the intra-phase LN response is linear, given that the responses to the three categories overlap (Figure 1, upper row).
The ZG500 and SPV metrics show responses to the strong events that are in magnitude more than double for EN than for LN (Figures 1a,1b,1d,and 1e).By looking at North Pacific geopotential height composite maps for the 1985-2014 period (Figure 2), we find that the difference shown by the ZG500 metric arises as the event magnitude increases, because of the anomalous pattern: larger size and eastward extension of the EN response.
Note that only the magnitude of the geopotential height maxima of the responses to strong EN and LN are comparable (∼60 m).Concerning the difference shown by the SPV metric for the strong events, its doubling for EN with respect to LN events captures the large EN versus LN difference in the zonal mean zonal wind response in the entire stratosphere (Figure 3, 1985-2014 period, right column). Figure 3 shows that the tropospheric subtropical jet responds more strongly already to weak EN compared to weak LN.For the 1985-2014 period, the stratospheric response is actually larger for weak LN compared to EN.As the SSTA magnitude increases, the stratospheric response becomes comparable (moderate events) and thereafter it is substantially larger for the EN events, so that to lead to the remarkable asymmetry for the strong events.Such large differences open the question of the role of sudden stratospheric warming (SSW) events.As noted in Ayarzagüena et al. (2021), the zonal mean zonal winds are biased low in strength in MPI-GE (Figure S2 in Supporting Information S1).Therefore, as in Ayarzagüena et al. ( 2021) we select SSW using the standard wind reversal definition with the additional condition of a 5-day duration.In doing so, we find that intra and inter-phase differences in the SSW frequency responses are generally minor in this ensemble, but for a substantial reduction (factor ∼2) between moderate and strong LN events (occurring for both JF and the entire winter season, not shown).The response to strong LN is therefore weak in spite of a decrease in SSW frequency.Overall, an SSW frequency response does not appear to be a necessary feature of the stratosphere response to ENSO, although SSWs can be facilitated or hampered.
Back to our metrics, we note that the ZG500 and SPV responses to LN remain well below those of EN also if linearly extrapolated to the larger value of the EN median, strong events (from ∼2.3 to ∼2.6 STD of the Niño3.4index, Figures 1d and 1e).In summary, we conclude that the ZG500 and SPV responses to the strong EN and LN events are asymmetric.In contrast, their responses to the weak events are classified as symmetric.For the moderate events, symmetry/asymmetry of the responses may depend on the 30-year period considered, with weak asymmetry suggested by comparing the medians of the response magnitude (Figures 1d and 1e).
Turning now to the last component of the stratospheric pathway, the NAE index, the main finding is that for each ENSO phase, the grouping of the responses by the three categories hardly leads to distinct responses (Figure 1c).Although a linear fit can be drawn between the intra-phase LN responses (Figure 1f), given the overlaps of the confidence intervals of the composite means, the intra-phase LN response is judged nonlinear, that is, largely not increasing with the SSTA intensity.The magnitudes of the responses to weak and moderate EN and LN categories are all comparable, leading to infer symmetry in this type of events.Only between the responses to moderate and strong EN events, the NAE metric suggests a distinction, as visually judged by the composite time series and indicated by the sharp increase in the median of their responses (Figures 1c-1f, respectively).In this case, a nonlinear intra-phase response is weakly suggested by the fact that the 50% of the moderate composite means do not intersect the linear fit between the three medians.The NAE response to strong events can therefore be asymmetric, with a stronger response for the warm than for the cold ENSO phase, as it is seen for instance in the SLP composite maps for the 1985-2014 period (Figure 4, right column).At the surface, however, the signatures of both the ENSO tropospheric and stratospheric pathway contribute to the responses.Given that the asymmetry of the stratospheric response extends to the eddy-driven jet in the troposphere (Figure 3, note the significant response at high latitudes), it is reasonable to understand the asymmetry of the SLP response to strong events as a downward stratospheric impact.Conversely, the overall symmetric SLP response for weak events illustrated in Figure 4 (middle column) could be due to contributions from both the stratospheric and tropospheric pathway.The tropospheric pathway has been shown to be symmetric in initialized climate prediction experiments (Hardiman et al., 2019).
Recalling that the ZG500, SPV and NAE responses to either EN or LN are independently calculated, the question now is whether the symmetry or asymmetry found in the three indicators of circulation responses is related.
Given that the ZG500 and SPV responses are linked through the response of quasi-stationary planetary waves, we have diagnosed DJ anomalous planetary waves 1 and 2 in geopotential height at 500 and 10 hPa, latitude band 20 o -90°N.Waves are diagnosed by Fourier analysis.We use the 500 hPa geopotential height anomalous wave 1 and wave 2 to reconstruct ZG500 index responses, and we compare these response reconstructions with the ZG500 index responses themselves (Figure 5a).For compactness, we show the distributions of the responses in terms of box-plots.The reconstructed ZG500 indices for either wave 1 or wave 2 show responses along the lines of those of the original ZG500 index: distinct responses by category for EN events but not for LN events and asymmetric/symmetric responses to strong/weak events.Therefore, locally, over the ZG500 index regional area, and relative to their own variability, the responses of both planetary waves are involved in the asymmetry of the geopotential height anomalous patterns, specifically the larger size and eastward extension of the EN response (Figure 2, right panels).
To estimate the global impacts of the wave response, we now turn to compare the amplitude of the non-normalized geopotential height anomalous waves at 500 hPa and averaged over 35°-50°N (Figure 5b).Notable is the sharp increase in the distribution of anomalous wave 1 amplitude from the moderate to strong EN (compare ENM and ENS blue boxes, Figure 5b).Instead, anomalous wave 2 amplitudes for moderate and strong EN events are comparable, as shown by the overlaps of the interquartile ranges of their distributions.Given the large overlaps between the interquartile ranges of anomalous wave 2 and wave 1 for weak and moderate events, both EN and LN, the only asymmetric response that emerges by this comparison is the one between the strong EN and LN events, and only for wave 1.To test if the asymmetry in the amplitude of wave 1 response can be linked to the anomalous diabatic heating associated with the atmospheric convective response to ENSO, we diagnose the deep tropical precipitation response (Figure 5c).We find that the distributions of the composited anomalous tropospheric geopotential height waves (Figure 5b) are consistent with those of the precipitation responses.Only the distributions of the responses in precipitation wave one power are distinct across the event categories.For wave 1, the responses to weak and moderate events overlap, leading to symmetry, whereas those to strong events do not, leading to asymmetry with larger responses for the warm El Niño phase.
The impact of the response of the quasi-stationary planetary waves on the stratospheric polar vortex is exerted in terms not only of amplitude but also of phase.The phases of wave 1 and 2 respond less differently across the event categories, while reproducing the well-known differences between their responses to EN and LN, examples are shown only for the strong events at 500 and 10 hPa, and only for the 1985-2014 period (Figure 6).As expected, at 500 hPa, wave 2 responses are in quadrature for both EN and LN and wave 1 responses tend to be out of phase for LN and in phase but with an eastward shift more pronounced at low latitudes, for EN.At 10 hPa, the locations of the wave 2 responses do not shift with respect to those at 500 hPa, meaning the anomalous waves are barotropic and do not contribute to the meridional heat flux response.Wave 1 responses instead tilts with elevation, so that to be fully in phase with the climatological waves for EN, and out of phase for LN at 10 hPa.The statistical significance in LN wave responses is also low.
The asymmetry in the amplitude of wave 1 response, together with the quasi-in-phase relationship of the anomalous and climatological wave 1, and the low induced by wave 2 in the central North Pacific (Figure 6, lower panel, third and fourth columns) can explain the asymmetry in the ZG500 index, strong events, by leading to the elongated lobe (South Eastward) of the EN geopotential height responses (Figure 2).The geopotential height anomaly for LN remains rather constrained, instead (Figure 2).
In turn, at 10 hPa, the wave 1 amplitude asymmetry and the in-phase relationship of the anomalous and climatological waves 1, can explain the asymmetry in the SPV index, strong events.The asymmetry of the stratospheric vortex response for the strong events can subsequently contribute to an asymmetry in the NAE response, as previously mentioned.

CMIP6 LEs
Here we ask how general the MPI-GE results are.To this end, we calculate the responses of the three metrics, ZG500, SPV and NAE, for the three categories and the two ENSO phases, in six CMIP6 LEs.For compactness, we show the results in the form of box plots and we include the MPI-GE as well for reference (Figure 7).Time series of the composite means by model, index and category are shown in the Supporting Information S1 (Figure .Shown is the fraction of total power by wavenumber, for the first six wavenumbers, and colors denote phase and category of events, as in Figure 1.Each box-plot shows the interquartile range (bounds of the box), the median (solid line inside the box), and the maximum and minimum in the distribution that are not outliers (whiskers, 99.3%).The circles correspond to outliers. 10.1029/2023JD039992 9 of 12 S3).Given the smaller ensemble size of the CMIP6 LEs, the distributions of the responses tend to be broader.This is the case, for instance, when comparing MPI-GE with MPI-ESM1-2-LR (here denoted MPI-LE).Nevertheless, common features emerge.
The ZG500 index tends to intensify in absolute values as the EN magnitude increases, so as to suggest a linear intra-phase relationship, for all the CMIP6 LEs except for MIROC6, supporting the MPI-GE findings.The responses to LN events are more varied, although there is a tendency for the response to overlap for at least two categories, here again the exception is MIROC6.In this latter case, there are not enough strong LN events to estimate a distribution (Table S1 in Supporting Information S1), but we note that the median of their N3.4SSTA STD is the largest of the six LEs (Table 1).The IPSL-CM6A-LR (here IPSL) and MPI-LE responses for strong EN and LN appear to be asymmetric (compare the medians, below 0.5 for LN but close to −1 for EN), whereas those of CanESM5 and CNRM-CM6-1 are judged symmetric (magnitude of the respective medians comparable).
Turning to the stratosphere, the SPV diagnostic shows that three LEs (CNRM-CM6-1, MIROC6, and NorCPM1) do not report a clear ENSO response, independently on the relative size of their respective ZG500 anomalies.An  issue in MIROC6 or NorCPM1 may be a very weak upward propagation of the signal from the Aleutian low region in the troposphere (typically associated with ENSO among other phenomena) to the stratosphere, as suggested by the negligible interannual correlation coefficients between DJ Z500 and JF SPV (Table S2 in Supporting Information S1).Among the three LEs with a response in the stratosphere, CanESM5 tends to show a more symmetric response to the strong ENSO events.This agrees with its symmetric response in ZG500, supporting a direct link of asymmetry/symmetry between ZG500 and SPV, if a model is able to simulate the stratospheric response to ENSO.The IPSL response in the stratosphere is asymmetric for the strong EN and LN events, consistently with its asymmetric response for ZG500.MPI-LE responds similarly to MPI-GE, as expected.
The NAE responses are well known to include not only stratospheric linkages to ENSO but also tropospheric impacts.Among the three LEs without a response in the stratosphere, we note that only one LEs (NorCPM1) shows asymmetry in the strong events, with a stronger response for the cold events.A larger LN response is at odds with previous studies (Frauen et al., 2014), although here only the FM response is reported.The other two LEs, without a response in the stratosphere, tend to show symmetric responses, in agreement with their symmetric responses in the North Pacific region (ZG500).All the three LEs with a response in the stratosphere tend to show a different symmetry/asymmetry for the strong events from what is found for the ZG500 and SPV responses.This contrast may indicate the importance of the variety of tropospheric bridges that connect ENSO with the European weather and climate and the high internal atmospheric variability in the NAE region.

Discussion and Conclusions
We have diagnosed nonlinearity and asymmetry along the ENSO stratospheric pathway in LEs of historical experiments performed with climate models.For one LEs (MPI-GE), for which we have 95-members and that we have analyzed in detail, a clear picture emerges.However, only a few of the features found in MPI-GE are shared by the other LEs.
To diagnose the ENSO stratospheric pathway, we have evaluated three indicators of circulation responses, as functions of the SSTA intensity.Specifically, we have assessed the North Pacific extra tropical mid-troposphere response (ZG500 index), the stratospheric polar vortex response (SPV index) and the North Atlantic -Europe (NAE index) surface response.SSTAs have been diagnosed using the Niño3.4index and have been stratified in three intensity categories according to their magnitude (weak, moderate, strong) to allow for testing for linearity of the response within each phase (intra-phase comparison) and for asymmetry (inter-phase comparison) between warm El Niño and cold La Niña phases.
Based on MPI-GE results, we conclude: 1. Linear relationships as a function of the SSTA magnitude can easily be drawn for either EN or LN responses, although only in the case of EN the moderate and strong responses are well separated for all three indicators of circulation responses.Consequently, the intra-phase EN responses can be judged to be overall linear while the LN responses are not.2. Only El Niño and La Niña strong events show an asymmetric response that can be traced along the stratospheric pathway, from the North Pacific troposphere, via the stratospheric polar vortex, to the North Atlantic-Europe surface response.Larger responses are found for the strong El Niño than for strong La Niña.The origin of this asymmetry in the strong events can be traced to the asymmetry in wave 1 amplitude responses of both the extratropical geopotential height at 500 hPa and the precipitation in the deep tropics.3. Symmetry between EN and LN responses instead characterizes the significant responses of weak events, along the stratospheric pathway considered.4. Responses to moderate events are more varied.A certain degree of asymmetry is found for moderate events in the North Pacific troposphere and for the polar vortex responses, but not for the North Atlantic-European region, at the surface.Clearly, the NAE region is affected by other factors, such as the tropospheric pathway, which also contribute to the ENSO response therein.
The ZG500 and SPV asymmetry in strong events found in MPI-GE is rooted in the larger rate of increase of the EN than LN responses with SSTA intensity.The tendency for a larger response for El Niño than La Niña is consistent with and supports the results of Frauen et al. (2014).In addition, the asymmetry along the stratospheric pathway that follows is broadly consistent with the results from the numerical experiments of Trascasa-Casto et al. ( 2019).However, in contrast to Trascasa-Casto et al. ( 2019), our LN responses are well developed for weak and moderate events as well.Consequently, the asymmetry we have found might be judged to be weaker.This discrepancy may be explained by the fact that we have analyzed experiments from a fully coupled climate model and used the Niño3.4index to identify events.In coupled experiments, the Niño3.4index can capture the location asymmetry in the SSTA maxima between EN and LN events, which is enabled by the explicit simulations of the air-sea feedbacks crucial for the growth of a response.The significant responses for moderate LN events found in MPI-GE are in fact in agreement with reanalysis results (Iza et al., 2016).
The relatively large responses to weak and moderate LN events found in MPI-GE let extrapolate that by merging the responses to the three categories, the overall EN versus LN asymmetry of responses may be reduced.Thus, not distinguishing by intensity category may favor agreement with recent observation-based assessments (Deser et al., 2017) and deductions from AMIP experiments (e.g., Garfinkel et al., 2019;Weinberger et al., 2019).Seasonality of the responses may be a key factor responsible for the different asymmetry in the stratospheric response found here and in Hardiman et al. (2019), because the response to EN may change sign in mid-winter in the stratosphere as shown indeed by Hardiman et al. (2019).In their work, Hardiman et al. (2019) focused on the December to February stratospheric response, however, their Figure 3a shows that by considering JF, as we did in our analysis, their stratospheric response is as well stronger (in magnitude) for strong El Niño than for strong La Niña, reconciling their results with the current work.Related to the seasonality of the stratospheric response, at the surface, climate models have been found not to reproduce the intra-seasonal change in the pattern of the response that observations or seasonal forecasts do show (Ayarzagüena et al., 2018).So far, it is unclear if the seasonality of the response is affected by model biases, initialization or unrelated to ENSO, given the low frequency variability of global SST in the observed record.
We have touched on extending the analysis to six CMIP6 LEs.However inter-model spread of the results is relatively high.Only the extra-tropical North Pacific responses show general agreement between MPI-GE and the other large ensembles, with a linear response with increasing signal intensity for El Niño but not so for La Niña, and generally smaller responses for La Niña.Only three out of the six examined models reported a response in the stratosphere.Of these three, two tend to show the same asymmetry discussed for MPI-GE, with stronger responses for El Niño.Of interest, the symmetry/asymmetry of ZG500 is carried over to the stratosphere, supporting the notion that the ZG500 response is a good constraint for the stratospheric response.In addition, from the analysis here reported, the relationship between the asymmetry in the North Pacific response and in the stratospheric polar vortex does not appear to be necessarily affected by the presence of the Quasi-Biennial Oscillation (QBO).The MPI-GE does not simulate a QBO and only one of the three CMIP6 LEs with a stratospheric vortex response to ENSO (IPSL) simulates a QBO.Further analysis, beyond this first exploratory reporting and larger ensembles by single climate models, would be needed to assess and understand inter-model spread and the modulating effects of additional factors.
Finally, it is important to note that even though extreme La Niña events are infrequent, they do exist and, more importantly, an increase in their frequency has been observed in the last decades and is also projected to occur during the 21st century (Cai et al., 2015).Although we have analyzed historical experiments, understanding the teleconnection of extreme La Niña events to the Euro-Atlantic sector is therefore becoming more relevant under the ongoing global warming.

Figure 1 .
Figure 1.MPI-GE results.Upper row: Time series of the responses to EN or LN by category, for (a) DJ ZG500 index, (b) JF SPV index, and (c) FM NAE index.The responses (circles) are given by compositing normalized anomalies.Each composite mean is based on a 30-year period and 95 members and the bar denotes the 95% confidence interval of the mean values based on the fitted normal probability distribution to the original data.Composites for weak, moderate and strong events are respectively shown in pink, magenta and red for EN; and azure, blue, and black for LN.Lower row: Intra-phase linear fits in function of Niño3.4SSTA magnitude, for the magnitudes of the (d) DJ ZG500 index; (e) JF SPV index; and (f) FM NAE index.The fits are calculated over the three medians of the distributions of the composites (i.e., the distributions made of the circles shown in the upper panels) by category, EN (red) and LN (blue).The horizontal and vertical bars show the interquartile ranges around the respective medians.

Figure 2 .
Figure 2. MPI-GE results.DJ 500 hPa geopotential height (m) composited anomalies for the 1985-2014 period.(Upper row) LN and (lower row) EN composites by category: (left) weak, (middle) moderate and (right) strong events.Shading shows statistically significant (p < 0.05) values.Black box shows the region considered in the calculation of the ZG500 index.

Figure 3 .
Figure 3.As Figure 2 but for JF zonal mean zonal winds (ms −1 ), latitude (degrees North) and pressure (hPa) sections.Thick black line highlights the latitudes and the level considered for the SPV index.

Figure 4 .
Figure 4.As Figure 2 but for FM mean sea level pressure (hPa).Thick black lines enclose the two regions considered for the computation of the NAE index.

Figure 5 .
Figure5.MPI-GE results.(a) Box-plot showing the distribution of composites, by ENSO phase and category, for the DJ ZG500 index and its wave 1 and wave 2 reconstructions, normalized anomalies.Colors denote phase and category of events, as in Figure1.(b) Box-plot showing the distribution of the amplitude of anomalous DJ 500 hPa geopotential height (m) wave 1 (blue) and wave 2 (magenta), averaged over 35-50ºN, by ENSO phase and event category.LN amplitudes shifted to the right, and the separation in the x axis is based on the Niño3.4mean value of each category.(c) Box-plot showing the distribution of composited anomalous DJ precipitation power averaged over 5°S-5°N.Shown is the fraction of total power by wavenumber, for the first six wavenumbers, and colors denote phase and category of events, as in Figure1.Each box-plot shows the interquartile range (bounds of the box), the median (solid line inside the box), and the maximum and minimum in the distribution that are not outliers (whiskers, 99.3%).The circles correspond to outliers.

Figure 6 .
Figure 6.MPI-GE results.DJ geopotential height (m) planetary wave 1 and 2 components, (upper row) at 10 hPa and (lower row) at 500 hPa.Wave 1 responses are shown in the first and third columns, wave 2 in the second and fourth columns.Strong LN at left and strong EN at right.In each panel, composited anomalies for the 1985-2014 period are shown in shading when statistically significant (p < 0.05).Contours show the climatological waves.

Figure 7 .
Figure 7. Box plot showing the distributions of composites by phase and category, (a) DJ ZG500 index; (b) JF SPV index; and (c) FM NAE index.Results for MPI-GE and the 6 CMIP6 LEs.Colors denote phase and category of events, as in Figure1.Each box-plot shows the 25th and 75th percentiles (bounds of the box), the median (solid line inside the box), and the maximum and minimum in the distribution that are not outliers (whiskers, 99.3%).The crosses correspond to outliers.

Table 1 CMIP6
Models, Number of Members and Niño3.4SSTA STD Statistics by LE