# Towards understanding the robust strengthening of ENSO and more frequent extreme El Niño events in CMIP6 global warming simulations

Ulla K. Heede, Alexey V. Fedorov

► To cite this version:

Ulla K. Heede, Alexey V. Fedorov. Towards understanding the robust strengthening of ENSO and more frequent extreme El Niño events in CMIP6 global warming simulations. Climate Dynamics, 2023, 61, pp.3047-3060. 10.1007/s00382-023-06856-x. insu-04188257

HAL Id: insu-04188257

https://insu.hal.science/insu-04188257

Submitted on 25 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Towards understanding the robust strengthening of ENSO and more frequent extreme El Niño events in CMIP6 global warming simulations

Ulla K. Heede<fnref n="1" />,<fnref n="3" /> · Alexey V. Fedorov<fnref n="1" />,<fnref n="2" />

Received: 6 June 2022 / Accepted: 4 June 2023 / Published online: 25 June 2023

© The Author(s) 2023, corrected publication 2023

## Abstract

The El Niño Southern Oscillation (ENSO) has profound impacts on weather patterns across the globe, yet there is no consensus on its response to global warming. Several modelling studies find a stronger ENSO in global warming scenarios, while other studies suggest ENSO weakening. Using a broad range of models from the Coupled Model Intercomparison Project phase 6 (CMIP6) and four types of warming experiments, here we show that the majority of the models predict a stronger ENSO by century-end in Shared Social Pathway (SSP) experiments, and in idealized 1pctCO₂ and abrupt 4xCO₂ experiments. Several models, however, do predict no change or ENSO weakening, especially in the idealized experiments. Critically, the strongest forcing (abrupt-4xCO₂) does not induce the strongest ENSO response, while differences between the models are much greater than those between warming scenarios. For the long-term response (over 1000 years) the models disagree even on the sign of change. Furthermore, changes in ENSO sea surface temperature (SST) variability are only modestly correlated with the tropical Pacific mean state change. The highest correlation for ENSO SST amplitude is found with the mean zonal SST gradient in the SSP5-8.5 experiment (R = -0.58). In contrast, changes in ENSO rainfall variability correlate well with changes in the mean state, as well as with changes in ENSO SST variability. When evaluating the Bjerknes Stability Index for a subset of models, we find that it is not a reliable predictor of ENSO strengthening, as this index tends to predict greater stability with warming. We argue that the enhanced ENSO stability is offset by increases in atmospheric noise or/and potential nonlinear effects. However, a robust inter-model mechanism that can explain a stronger ENSO simulated with global warming is still lacking. Therefore, caution should be exercised when considering ENSO changes based on a single model or warming scenario.

**Keywords** El Niño · ENSO · Tropical Pacific · CMIP6 · Climate change · Extreme events

## 1 Introduction

The El Niño Southern Oscillation (ENSO) is the strongest interannual oscillation in Earth's climate, which has dramatic impacts on weather patterns and extreme events in the tropics and beyond (Cai et al. 2021; McPhaden et al. 2020). El Niño conditions are characterized by warm central and/or eastern Pacific sea surface temperature (SST)

anomalies, reduced easterly trade winds and the corresponding east-west slope of the ocean thermocline, and eastward shift of atmospheric convection that modifies precipitation patterns on both sides of the tropical Pacific. By exciting atmospheric planetary waves, ENSO also affects a range of whether phenomena outside the tropical Pacific via teleconnections, and the signature of ENSO has been observed on all seven continents (Yeh et al. 2018).

Whether ENSO is already changing or will change with future warming, and whether this change will lead to stronger and/or more frequent El Nino events is of great importance for timely implementation of societal adaptation and mitigation strategies. However, despite the importance of ENSO in Earth's climate, much uncertainty remains concerning these questions.

It has long been known that ENSO can be affected by changes in the mean state of the tropical Pacific (Fedorov

<figure><img src="image_2.png" /></figure>

<footnote>Ulla K. Heede
ulla.heede@colorado.edu</footnote>

<footnote n="1">¹ Department of Earth and Planetary Science, Yale University, New Haven, CT, USA</footnote>

<footnote n="2">² LOCEAN/IPSL, Sorbonne University, Paris, France</footnote>

<footnote n="3">³ Cooperative Institute for Research in Environmental Sciences (CIRES), University of Colorado, 4001 Discovery Dr, Boulder, CO 80303, USA</footnote>
and Philander 2000; Jin et al. 2006; DiNezio et al. 2012; Fedorov et al. 2020; An and Jin 2000; Zhao and Fedorov 2020; Guilyardi 2006; Guilyardi et al. 2009, 2012). Therefore, it is expected that anthropogenic climate change will influence ENSO by modifying the mean state. Specifically, a vast majority of models participating in the Coupled Model Intercomparison Project, phase 6 (CMIP6) predict that the Walker cell will eventually slow down with global warming, leading to enhanced warming of the eastern and central equatorial Pacific by century-end, dubbed the eastern equatorial Pacific warming pattern (e.g. Heede and Fedorov 2021; Xie et al. 2010). This pattern can be understood in terms of competition between a number of atmospheric and oceanic mechanisms, which acts to reduce equatorial trade winds and eastern Pacific upwelling (Heede et al. 2020; Heede and Fedorov 2021; Vecchi et al. 2006), even though on transient timescales the Walker cell may actually strengthen (Seager et al. 2019; Heede and Fedorov 2023).

Despite the model consensus on the Walker circulation weakening by the late twenty-first century and the established links between the mean state and ENSO characteristics, there has been no consensus among climate models regarding the response of ENSO to global warming. A number of studies noted a stronger ENSO in warming climates in CMIP5 models, particularly when utilizing precipitation metrics (Huang and Xie 2015; Cai et al. 2015) or metrics that isolate Eastern Pacific (EP) or Central Pacific (CP) El Niño events (Cai et al. 2018). The tendency for a stronger ENSO in future warming scenarios has also been reported for subsets of CMIP6 models (Fredriksen et al. 2020; Cai et al. 2021; Brown et al. 2020). Cai et al. (2018) pointed out two mechanisms that could result in stronger EP El Niño events: a stronger coupling between wind stress and SST in the eastern Pacific, and a stronger coupling between ocean thermocline and wind stress due to enhanced upper ocean stratification.

Contrary to Cai's et al. (2018) findings, several recent studies have argued that ENSO amplitude may actually decrease with warming (Wengel et al. 2021; Callahan et al. 2021; Kohyama et al. 2018). These studies typically rely on abrupt CO₂ quadrupling experiments and argue that increased thermodynamic damping and decreased positive ocean feedbacks, among other mechanisms, lead to a more stable ENSO in warmer conditions. Notably, Kohyama et al. (2018) specifically argues that enhanced ocean stratification should result in a stiffer thermocline and a decrease in ENSO nonlinearity, contradicting the findings of Cai et al. (2018). Finally, Maher et al. (2023) find from large ensemble CMIP6 SSP experiments that some models exhibit a time dependent response where the ENSO amplitude initially increases, then decreases. These disagreements highlight the fact that changes in different feedbacks resulting from changes in the

mean state and ocean-atmosphere coupling can partially compensate one another, leading to unforeseen outcomes.

The goal of this study is to try to reconcile these divergent views of how ENSO will respond to global warming and to offer a more comprehensive analysis of how ENSO changes within the CMIP6 models, compared to previous studies (i.e. Fredriksen et al. 2020; Cai et al. 2022), by using a broader range of models and climate scenarios. As studies finding a stronger ENSO primarily rely on the twenty-first century SSP scenarios, whereas studies suggesting a weaker ENSO use the abrupt CO₂ quadrupling experiments, here we present an overview of how models respond across these idealized and more realistic forcing experiments. This approach enables us to address uncertainties associated with the types of experiments used.

Further, we aim to explore how the various mechanisms proposed for a stronger or weaker ENSO may manifest across different models and document how ENSO changes might be linked to the tropical Pacific mean state. Finally, to examine whether the transient model response differs from the long-term response in the abrupt-4xCO₂ experiments as proposed by Callahan et al. (2021), we analyze three 1000-year simulations. To our knowledge, no study has compared ENSO changes across such a broad range of idealized and realistic forcing scenarios in CMIP6 models, as well as attempted to link those changes to tropical mean state changes in a comprehensive manner.

# 2 Methods

## 2.1 CMIP6 archive

We analyze two types of experiments from the Climate Model Intercomparison Project Phase 6 (CMIP6) archive (Eyring et al. 2016). The first type consists of two hypothetical or idealized CO₂-only experiments: abrupt-4xCO₂ rise and 1pctCO₂ gradual CO₂ rise (1% per year) where CO₂-increases are relative to the pre-industrial level of 280 ppm. The second type includes more realistic, full-forcing experiments following two future scenarios: SSP5-8.5 which is a high-emission scenario; and SSP1-2.6, in which emissions peak and decline in the twenty-first century (O'Neill et al. 2016).

We use a total of 20 CMIP6 models in our overview analysis (Table 1), utilizing surface temperature, surface precipitation, and zonal wind stress. The criterion for including a given model into the analysis is whether it has surface temperature data available for at least three ensemble members for the SSP5-8.5 scenario. Some models do not have all datasets available for all experiments, and models are excluded when variables are not available. For the idealized warming experiments (4xCO₂ and 1pctCO₂), where ensemble members are

<figure><img src="image_3.png" /></figure>
<table><caption>**Table 1** Overview of models used and their associated marker</caption><thead><tr><th>Model</th><th>Marker</th><th>Model</th><th>Marker</th></tr></thead><tbody><tr><td>ACCESS-CM2</td><td>●</td><td>FGOALS-g3</td><td>▲</td></tr><tr><td>ACCESS-ESM1-5</td><td>▼</td><td>FIO-ESM-2-0</td><td>+</td></tr><tr><td>CanESM5-CanOE</td><td>▲</td><td>HadGEM3-GC31-LL</td><td>★</td></tr><tr><td>CanESM5</td><td>▲</td><td>HadGEM3-GC31-MM</td><td>●</td></tr><tr><td>CESM2-WACCM</td><td>►</td><td>IPSL-CM6A-LR</td><td>+</td></tr><tr><td>CESM2</td><td>Y</td><td>KACE-1-0-G</td><td>✗</td></tr><tr><td>CNRM-CM6-1</td><td>Y</td><td>MIROC6</td><td>✗</td></tr><tr><td>CNRM-ESM2-1</td><td>Y</td><td>MIROC-ES2L</td><td>✗</td></tr><tr><td>EC-Earth3</td><td>Y</td><td>MPI-ESM1-2-LR</td><td>|</td></tr><tr><td>EC-Earth3-Veg</td><td>■</td><td>UKESM1-0-LL</td><td>—</td></tr></tbody></table>

<table><tr><td>ACCESS-CM2</td><td>●</td></tr><tr><td>ACCESS-ESM1-5</td><td>▼</td></tr><tr><td>CanESM5-CanOE</td><td>▲</td></tr><tr><td>CanESM5</td><td>▲</td></tr><tr><td>CESM2-WACCM</td><td>►</td></tr><tr><td>CESM2</td><td>Y</td></tr><tr><td>CNRM-CM6-1</td><td>Y</td></tr><tr><td>CNRM-ESM2-1</td><td>Y</td></tr><tr><td>EC-Earth3</td><td>Y</td></tr><tr><td>EC-Earth3-Veg</td><td>■</td></tr></table>

generally not available, we use a single ensemble member, and for the SSP scenarios we use the ensemble-mean results based on three members. We analyze 85 years of simulation for the SSP experiments, 150 years for the idealized experiments, and 200 years for the piControl experiment respectively.

## 2.2 Metrics

We define the amplitude of ENSO SST variability as the standard deviation of temperature anomalies band-passed in the frequency range between 1.5 and 7 years and averaged for the central-eastern equatorial Pacific (180–280 °E). This region is chosen to capture ENSO changes across both the central and eastern Pacific following Ferrett and Collins (2019). Unless specified differently, the equatorial Pacific is defined as a band between 5 °S and 5 °N. Similarly, we define the strength of ENSO rainfall variability as the standard deviation of precipitation anomalies band-passed in the same frequency range and averaged for the same region. The Butterworth IIR 4-th order filter with second-order-series is used in the xr-sipy frequency filter wrapper. Before applying the filter, the climatology was subtracted from the data.

While the issue of differences between Central Pacific (CP) versus Eastern Pacific (EP) El Niño events has been raised as a potential culprit when using set boxes to characterize ENSO (Cai et al. 2018), we find that the estimated ENSO changes are not very sensitive to the boundaries of the chosen box. In a complementary analysis we also use the standard deviation of the E-index and C-index as defined in Cai et al. (2018) as

$$E = (PC1 - PC2)/\sqrt{2}$$

$$C = (PC1 + PC2)/\sqrt{2}$$

where PC1 and PC2 are the first two principal components computed from EOF analysis for the region 15 °S to 15 °N, 140 °E to 80 °W. In the warming scenarios, we compute the two indices using EOFs from the piControl following Frederiksen et al. (2020).

Further, we define extreme El Niño events as events characterized by warm temperature anomalies in the central-eastern Pacific (180–280 °E) that exceed 2 standard deviations. While 2.5 standard deviations are occasionally used to define an extreme event (e.g. Yu and Fedorov 2020), we chose 2 to get a more robust statistics as most models have few events exceeding 2.5 standard deviation in the timeframe between 85 and 200 years.

All changes are computed relative to the piControl experiment. While some studies compare SSP scenarios to the historical period (i.e. Cai et al. 2022), we compare all types of experiments to the piControl to ensure a fair comparison between idealized experiments and future warming scenarios. We have conducted an additional analysis indicating that, for the majority of the models, the strength of ENSO in the piControl experiment is comparable to that of the historical experiment over the period 1850–2015.

Finally, we define the mean state east-west Pacific SST gradient along the equator as the difference between time mean anomalies of the western Pacific (120–180 °E) and the central-eastern Pacific (180–280 °E). In addition to the mean SST gradient, we also investigate links between changes in ENSO SST variability and ocean mean stratification. We define ocean stratification following Cai et al. (2018) as the as the difference in potential temperature between mean
temperature of the upper 75 m and the temperature at 100 m in a box of 5 °S to 5 °N and 150 °E to 140 °W.

## 2.3 Bjerknes Stability Index

To evaluate possible links between mean state changes and the ENSO response to global warming, we utilize the Bjerknes Stability Index (Jin et al. 2006). The Bjerknes Index (BJ index, in units of yr⁻¹) assumes that ENSO dynamics are controlled by the recharge oscillator physics (Jin 1997) where SST perturbations in the eastern Pacific are either amplified or dampened through a series of feedbacks:

$$BJ = CD + TD + ZA + ED + TC$$

Here, CD is damping by the mean currents given by:

$$CD = - \left( \frac{\langle \bar{u} \rangle}{L_x} + \frac{\langle -2y\bar{v} \rangle}{L_y^2} + \frac{\langle H(\bar{w})\bar{w} \rangle}{H_m} \right)$$

Angle brackets indicate a spatial average for the central-eastern Pacific (180–280 °E) and bars indicate a time-mean over the chosen period of the experiment. $\bar{u}$, $\bar{v}$ and $\bar{w}$ are the mean zonal, meridional and vertical velocities in the mixed layer. $L_x$ and $L_y$ are the zonal and meridional scale lengths. $H_m$ is the depth of the ocean mixed layer (50 m), and $H(\bar{w})$ is the Heaviside step function, which ensures that only positive values of vertical velocity are taken into account.

TD is thermodynamic damping given by:

$$TD = -\alpha$$

Here, $\alpha$ is the linear regression coefficient from the mixed layer temperature anomalies to net anomalous heat flux $E_{net}$ at the ocean surface. $E_{net}$ is defined as:

$$E_{net} = E_{\text{incoming\_shortwave}} + E_{\text{incoming\_longewave}} - E_{\text{outgoing\_shortwave}} \\ - E_{\text{outgoing\_longewave}} - E_{\text{latent\_heat}} - E_{\text{sensible\_heat}}$$

To have units of yr⁻¹, $\alpha$ is further multiplied by $\frac{1}{\rho C_p H_m}$, where $\rho$ is sea water density and $C_p$ is specific heat capacity of seawater.

ZA, EK and TC correspond to the zonal advection, Ekman and thermocline feedbacks respectively, given by:

$$ZA = \mu_a \beta_u \langle -\frac{\partial \bar{T}}{\partial x} \rangle$$

$$EK = \mu_a \beta_w \langle -\frac{\partial \bar{T}}{\partial z} \rangle$$

$$TC = \mu_a^* \beta_h \langle -\frac{\bar{w}}{H_m} \rangle$$

where $\mu_a$ and $\mu_a^*$ are the coefficients of linear regression of wind stress onto temperature anomalies. $\beta_u$, $\beta_w$, and $\beta_h$ are the coefficients of linear regression of zonal current, vertical current and thermocline anomalies, respectively, onto wind stress anomalies. $\frac{\partial \bar{T}}{\partial x}$ is the mean zonal temperature gradient and $\frac{\partial \bar{T}}{\partial z}$ is the mean vertical stratification. $\bar{w}$ is the mean vertical velocity. All regression coefficients are for the central-east equatorial Pacific (180–280 °E), except $\mu_a^*$ which is calculated for the full east-west extent of the Pacific equatorial band (130–280 °E). Following Kim et al. (2014) and Ferrett and Collins (2019), we use temperature at 50 m depth as a proxy for thermocline anomalies. In a supplementary analysis, we repeat our calculation of the Bjerknes Index, but for the Niño3 region (210–280 °E) to illustrate how the relative roles of the damping and feedback terms change when a different averaging region is used. Again, the equatorial region is defined between 5 °S and 5 °N.

## 2.4 Atmospheric noise

Finally, an important contribution to ENSO dynamics comes from atmospheric noise, which can sustain a damped oscillation even when the BJ-index is negative. Here we define atmospheric noise in the tropical Pacific following Philip and Oldenborgh (2009) as:

$$\epsilon(x, y, t) = \tau_x(x, y, t) - \sum_{i=1}^{2} A_i(x, y)_i T'_i(t)$$

where $\epsilon(x, y, t)$ denotes stochastic forcing by random wind stress variations, $\tau_x(x, y, t)$ denotes zonal wind stress anomalies, $T'_i(t)$ are SST anomalies averaged over selected regions, and $A_i(x, y)_i$ are the domain-wide wind stress patterns corresponding to those SST anomalies. Here we use $i=1$ for the Niño3 region, and $i=2$ for the Niño4 region. The patterns $A_i(x, y)$ are obtained by regressing monthly wind stress anomalies onto the corresponding monthly SST indices. To estimate the noise magnitude, we then average $\epsilon$ over the equatorial Pacific region (130–280 °E, 5 °S–5 °N) and compute the standard deviation over time.

# 3 Results

## 3.1 Changes in SST and rainfall variability, and extreme events, in a broad range of warming experiments

We find that all models show an enhanced ENSO SST variability in the SSP5-8.5 scenario relative to the piControl experiment (in the ensemble-mean sense). Additionally, most models exhibit a similar or slightly smaller increase of
ENSO SST amplitude in the SSP1-2.6 scenario, although the differences between the two scenarios are small compared to variations among the models. The ENSO response in the 1pctCO₂ scenario is also well correlated with the SSP5-8.5 response, but it is more muted.

On average, the abrupt-4xCO₂ scenario yields a lower correlation with the SSP5-8.5 response and also has the largest spread compared to other warming scenarios. In fact, six models in the 4xCO₂ scenario exhibit a reduction in ENSO SST variability, with CESM2 having the most dramatic reduction (Fig. 1). In contrast, two models (ACCESS-ESM5-1 and FIO-2) show an increase in ENSO SST variability in the abrupt-4xCO₂ experiment greater than that observed in any ensemble member of the SSP scenarios.

As expected, changes in extreme El Niño events, defined through their Niño3 SST index, are linked to changes in ENSO SST amplitude. Indeed, the relationship between ENSO amplitude and extreme events demonstrates a sharp slope, such that a 25% increase in ENSO SST amplitude results in an approximately 100% increase in the frequency of extreme El Niño events (Fig. 2). The strongest correlation between the occurrence of extreme events and ENSO

SST variability is found in the abrupt-4xCO₂ experiments (Fig. 2).

While the response of ENSO SST variability is highly dependent on both the experiment and model used, there is a universal increase in ENSO rainfall variability across all models and all experiments. In contrast to ENSO SST amplitude, the increase in ENSO rainfall variability is largest, on average, in the abrupt-4xCO₂ experiment (Fig. 3).

The next question is whether and how these ENSO changes are related to changes in the mean state of the tropical Pacific and other factors. We find that the correlations between changes in ENSO SST amplitude and in the mean state zonal SST gradient or ocean stratification are relatively weak (Fig. 4). The strongest correlation between the mean state zonal gradient and ENSO SST amplitude changes is observed in the SSP5-8.5 experiment (R = -0.58, p = 0.01), while the strongest correlation between ocean stratification and ENSO SST amplitude changes is found in the abrupt-4xCO₂ experiment (R = -0.56, p = 0.05). Other correlations are much smaller and not significant. Using the Niño3.4 index, E-index or C-index did not improve these correlations (Supplementary Figs. 1 and 2).

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Overview of ENSO changes in CMIP6 models across different experiments. **a** A bar chart showing the time-mean amplitude of ENSO SST variability (Methods) in individual models for the piControl, abrupt-4xCO₂, 1pctCO₂, SSP5-8.5 and SSP1-2.6 experiments. Thin errorbars for SSP experiments represent the maximum and minimum values across three ensemble members. **b–d** Changes in ENSO SST amplitude (in units of °C) in the 1pctCO₂, abrupt4xCO₂ and</figcaption></figure>

SSP1-2.6 scenarios, respectively, versus the SSP5-8.5 scenario. These changes are calculated relative to the piControl time-mean ENSO SST amplitude. Each marker + color combination denotes the same model across all plots (Table 1). The corresponding correlations and p-values are shown in each bottom panel. The high correlation values imply a general consistency of ENSO response across different types of warming experiments
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Overview of changes in extreme El Niño events. **a** A bar chart showing the time-mean frequency of extreme ENSO SST events (Methods) for the piControl, abrupt-4xCO₂, 1pctCO₂, SSP5-8.5 and SSP1-2.6 experiments, respectively. Thin errorbars for SSP experiments represent the maximum and minimum values across three ensemble members. **b–d** Relative changes (%) in the frequency</figcaption></figure>

of extreme events versus changes in ENSO SST amplitude in the 1pctCO₂, abrupt4xCO₂ and SSP5-8.5 scenarios, respectively. Note that the changes are calculated relative to the piControl time-mean ENSO SST variability. Each marker+color combination denotes the same model across all plots (Table 1). Note the robust increase in the frequency of extreme El Niño events in most models

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Overview of ENSO precipitation changes. A bar chart showing the time-mean amplitude of ENSO rainfall variability (Methods) in individual models for the piControl, abrupt-4xCO₂, 1pctCO₂, SSP5-8.5 and SSP1-2.6 experiments, respectively. Thin errorbars for SSP experiments represent the maximum and minimum values across three ensemble members</figcaption></figure>

Correlations between changes in ENSO rainfall variability and the mean zonal SST gradient in the gradual experiments (1pctCO₂ and SSP5-8.5) are generally higher (Fig. 5a–c) than those for ENSO SST variability, except for

the abrupt-4xCO₂ experiment. However, the highest correlations are observed between changes in ENSO rainfall variability and ENSO SST variability (Fig. 5d–f). In other words, a greater ENSO amplitude in terms of SST corresponds to

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Linking ENSO SST variability changes to changes in the mean east-west SST gradient and ocean stratification. **a–c** Changes in ENSO SST amplitude versus changes in the mean zonal SST gradient for the 1pctCO₂, abrupt-4xCO₂ and SSP5-8.5 scenarios, respectively. Each marker+color combination denotes the same model (Table 1). The corresponding correlations and p-values are shown at the top of each panel. Negative values on the horizontal axes indicate the weakening of the SST gradient. The link between ENSO SST variability and the mean zonal gradient is generally weak, even though there is a tendency towards stronger ENSO when this gradi-</figcaption></figure>

ent relaxes in the 1pctCO₂ and SSP5-8.5 scenarios. The strongest correlation is observed in the SSP5-8.5 experiment (R = -0.58). **e–g** As in the top panels but with the horizontal axes showing changes in mean ocean stratification. Stratification is defined as the potential temperature difference between the upper 75 m and at 100 m in the box between 5 °S to 5 °N and 150 °E to 220 °W. The link between changes in ocean stratification and changes in ENSO SST variability is weak. The strongest correlation is seen in the abrupt-4xCO₂ experiment (R = -0.56)

higher rainfall variability. We also investigated whether changes in ENSO rainfall variability might be linked to the eastern equatorial Pacific warming, mean tropical warming, and mean rainfall, but no statistically significant relationships were found.

## 3.2 Changes in the Bjerknes Stability Index, thermodynamic damping, and atmospheric noise

To investigate the connection between mean state changes and ENSO SST variability response in a more quantitative manner, we compute the Bjerknes Stability Index (BJ Index) for a subset of models that have ocean subsurface data available. First, we compute this index using the central-eastern equatorial Pacific (180–280 °E) as the averaging region, following Kim et al. (2014). This analysis reveals that all models in the subset show an increase in mean state stability with warming. For five models, out of seven used, increase

in thermodynamical damping α dominates the changes, reducing the Bjerknes Index and hence increasing ENSO stability (Fig. 6). For two models (MIROC6 and MIROC-ES2L), the increase in α is more modest and reduction in the thermocline feedback dominates the change of the Bjerknes Index with warming.

For all models, the decrease in the Bjerknes Index is higher for abrupt-4xCO₂ than for SSP5-8.5, suggesting that a rapid strong increase in radiative forcing stabilizes the system in terms of linear stability. While the strengthening of the damping terms in response to the warming is robust across the models, there is no consensus on changes to the feedback terms. Five models show an increase in the Ekman feedback, four models show an increase in the thermocline feedback, while all seven models show only small changes in the zonal advection feedback.

We next recompute the Bjerknes Index but for the smaller Niño3 region (Supplementary Fig. 3). In this case, the computed Bjerknes Index is generally more negative, due to weaker positive feedbacks compared with those computed
<figure><img src="image_9.png" /><figcaption>**Fig. 5** Linking ENSO rainfall variability changes to changes in SST mean state gradient and ENSO SST variability. **a–c** Changes in ENSO rainfall variability versus changes in the zonal SST mean gradient for the 1pctCO₂, abrupt-4xCO₂ and SSP5-8.5 scenarios, respectively.</figcaption></figure>

**d–f** As in the top panels but with the horizontal axes showing changes in ENSO SST amplitude. Each marker + color combination denotes the same model (Table 1). The corresponding correlations and p-values are shown at the top of each panel

<figure><img src="image_10.png" /><figcaption>**Fig. 6** The Bjerknes Stability Index and its contributing terms for a subset of CMIP6 models. **a** ENSO SST amplitude for a subset of models for the piControl, abrupt-4xCO₂ and the SSP5-8.5 experiments. **b–d** and **g–i** Negative and positive feedback terms for the same experiments. **e** Atmospheric noise amplitude. **f** The total magnitude</figcaption></figure>

of the Bjerknes Index. All feedback terms and the BJ index are in units of year⁻¹. Noise is in N m⁻² and ENSO SST amplitude is in °C. The Bjerknes index is computed for the region 180–280 °E (c.f. Supplementary Fig. 3 using the Niño3 region for averaging)
for the full central-eastern Pacific region. For the global warming simulations, 5 out of 7 models still show a more stable Niño3 Bjerknes Index, which is similar to the central-eastern Pacific Bjerknes Index. However, two models (MIROC6 and MIROC-ES2L) now show a more unstable Niño3 Bjerknes Index with global warming, driven by an increase in the zonal advection and thermocline feedbacks. This could indicate that changes in the mean state of these two models lead to a more unstable ENSO in the Niño3 region, which might be able to explain a stronger ENSO SST variability in these two particular models. Nevertheless, for these two models the use of the Bjerknes Index for explaining ENSO strengthening becomes problematic, given the strong sensitivity of the results to the choice of the averaging region.

Next, to examine the role of potential changes in non-linearities that could counteract increased linear stability, we calculate changes to wind stress-SST coupling ($\mu_a$) and thermodynamic damping ($\alpha$) for three different ranges of temperature anomalies (corresponding to large positive, moderate, and large negative values, respectively). We find that the thermodynamic damping $\alpha$, while increasing in the central range of temperature anomalies (–1 to 1 °C), actually decreases for temperature anomalies above 1 °C in CanESM5, CESM2, CNRM-CM6-1, HadGEM3-CG31-LL (Supplementary Fig. 4). This suggests that changes in $\alpha$, which dominate the Bjerknes Index response, may be overestimated in some models because the thermodynamic damping could become less efficient for larger temperature anomalies. The coupling coefficient $\mu_a$ becomes stronger in the warming experiments in the majority of models in the central temperature anomaly range (–1 to 1 °C), while in IPSL-CM6A-LR, MIROC6 and MIROC-ES2L, $\mu_a$ reduces for temperature anomalies below –1 °C (Supplementary Fig. 5). One however has to be cautious in interpreting these results since the correlations between the terms used when computing regressions may become very small for the large positive and large negative ranges of temperature anomalies.

To investigate the causes of discrepancy between changes in the Bjerknes Index and ENSO response, we also analyze changes to atmospheric noise within the models for which wind stress data is available. Atmospheric noise, including westerly wind bursts that occur frequently in the tropical Pacific, is believed to play an important role in sustaining ENSO and especially extreme El Niño events (e.g. Puy et al. 2019; Yu and Fedorov 2020; Fedorov 2002; Fedorov et al. 2015; Hu and Fedorov 2019; Capotondi et al. 2018).

We find that on average, correlations between changes in atmospheric noise and ENSO SST amplitude (Fig. 7) are significant and higher than those between changes in the mean state and ENSO SST amplitude. This indicates that changes in atmospheric noise could play an important role in ENSO response to global warming and potentially

explain the general strengthening of ENSO, even though there are a few outliers in the abrupt-4xCO₂ experiment. However, quantifying the noise/ENSO amplitude relationship is difficult because (1) the computation of the Bjerknes index and noise amplitude is sensitive to the choice of the averaging regions and (2) the question of causality is difficult to resolve as the noise amplitude could increase simply as a result of strengthening of ENSO, and not vice versa. Nevertheless, the recent work of Liang and Fedorov (2022) does suggest that westerly wind bursts could become stronger in a warming climate regardless of changes to ENSO amplitude.

## 3.3 Long-term warming experiments

Lastly, we investigate three abrupt-4xCO₂ experiments, each lasting 1000 years, for three models (ACCESS-ESM1-5, CESM2 and IPSL-CM6A-LR). Surprisingly, even after such a long time, ENSO response does not converge in these models. In particular, CESM2 shows a drastic reduction in ENSO strength, with its amplitude reaching only about 1/3rd of that in the control by the end of the simulation. In contrast, in ACCESS-ESM1-5, ENSO remains stronger than in the control throughout the simulation period (Fig. 8), while in IPSL-CM6A-LR, ENSO shows minimal changes over the first 600 years, followed by a slight strengthening.

# 4 Discussion and conclusions

We have investigated the ENSO response to global warming across 20 models and four different scenarios, encompassing both realistic and idealized radiative forcings. Our findings reveal that most of the models exhibit enhanced ENSO SST variability in all warming scenarios considered. Moreover, in the high-forcing scenario SSP5-8.5, all models demonstrate a stronger ensemble-mean ENSO. However, significant variations exist among the models, surpassing the differences observed between the warming scenarios. Most models do show a generally consistent response across different types of warming experiments indicating that the simulated robust ENSO strengthening is indeed a result of radiative forcing, and the spread among models is primarily caused by factors other than internal variability within the models.

It is noteworthy that for the majority of models, the SSP scenarios with a gradual CO₂ increase exhibit a stronger ENSO response to warming on average than the abrupt-4xCO₂ scenario. In the last scenario, roughly 40% of the models analyzed show a weakening or no change in ENSO SST variability, which implies that caution should be exercised when using such experiments to make projections of ENSO future changes in the twenty-first century.
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Overview of changes in atmospheric noise in the tropical Pacific. **a** A bar chart showing the time-mean amplitude of model atmospheric noise (Methods) for the piControl, abrupt-4xCO₂, 1pctCO₂, SSP5-8.5 and SSP1-2.6 experiments. Thin errorbars for the SSP-experiments represent the maximum and minimum values across</figcaption></figure>

three ensemble members. **b–d** Changes in the noise amplitude versus changes in ENSO SST amplitude in the abrupt-4xCO₂, 1pctCO₂, and SSP5-8.5 experiments, respectively. The corresponding correlations and p-values are displayed at the top of each panel

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Simulated long-term changes in ENSO SST variability. The plot shows ENSO SST amplitude estimated for consecutive 200-year time intervals of 1000-year abrupt-4xCO₂ simulations using three different climate models</figcaption></figure>

We have attempted to relate the projected changes in ENSO SST variability to the simulated changes in the mean state characteristics, specifically to the mean east-west equa- torial SST gradient and ocean vertical stratification. We find that in the SSP5-8.5 scenario, ENSO SST amplitude tends to increase when the mean SST gradient is reduced, but the correlation is modest (R = -0.58, p = 0.01) and disap- pears in other types of experiments considered. Likewise, we find some correlation between ENSO SST amplitude and ocean stratification in the abrupt 4xCO₂ experiment (R = -0.56, p = 0.05), but it decreases in the SSP5-8.5 sce- nario (R = -0.26, p = 0.05).

Despite this uncertainty in ENSO SST response, ENSO rainfall variability in the tropics increases universally across all experiments and models, which has important conse- quences for adaptation and mitigation because changes in drought or flooding conditions may cause more damage than changes in SST and because ENSO remote telecon- nections depend of tropical latent heat release. The fact that rainfall variability increases robustly with warming is expected given that a warmer atmosphere can hold more water following the Clausius-Clapeyron relation and since the eastern Pacific warming can increase sensitivity of pre- cipitation to SST anomalies (Yun et al. 2021). However, we find that changes in the background state alone cannot fully

<figure><img src="image_13.png" /></figure>
account for the modeled change in ENSO rainfall variability as changes in ENSO SST variability play a critical role as well. This highlights the need to better understand the response of ENSO SST variability to improve predictions for ENSO rainfall response.

Another important factor to consider is the occurrence of extreme events, which can change radically even for a small increase in ENSO SST amplitude. In fact, a relatively modest change in ENSO SST variability can easily lead to the doubling of the number of extreme El Niño events (e.g. Yu and Fedorov 2020). For the SSP5-8.5 experiment, which shows the largest increase in extreme events, the frequency of extreme events increases by 100% from 1 extreme event on average every other decade to one event per decade. For the other experiments, the increase is between 50 and 100%. The strongest correlation between ENSO SST and extreme events is observed in the abrupt-4xCO₂ experiments. This close connection between extreme events and ENSO SST amplitude is evident even in a relatively short timeframe of 150 years. The fact that the occurrence of extreme events is closely connected to the change in ENSO SST variability has important consequences for adaptation and mitigation since extreme events often cause more damage for society than a change in mean conditions (Trenberth 2012). Given the sensitivity of extreme events to small changes in ENSO SST, it is therefore crucial to improve projections of ENSO SST amplitude.

In an attempt to explain the robust strengthening of ENSO in warming scenarios we compute changes in the Bjerknes Stability Index but find it to be a relatively poor predictor for changes in ENSO in the small subset of models used (7 models), which is similar to findings in other studies using the Bjerknes Index for individual models (Manucharyan and Fedorov 2014; Ferrett and Collins 2019). Like Callahan et al. (2021), we find increased thermodynamic damping to be the most important term in the Bjerknes Index for 5 out of 7 models, yielding a more stable Bjerknes index in those models. However, the fact that all models show a stronger ENSO SST variability in the SSP5-8.5 scenario points to other effects counteracting this stability increase. For example, here we show that the thermodynamic damping may be overestimated in the models due to a nonlinear relationship between SST and surface heat fluxes. We show that the coupling between SST and surface heat fluxes decreases in warming experiments when estimated for temperature anomalies above a 1 °C threshold in 4 out of 7 models, which could explain why ENSO SST amplitude can increase despite a stronger thermodynamic damping estimated for arbitrary temperature anomalies.

Two models (MIROC6 and MIROC-ES2L) show a decrease in the Bjerknes Stability Index for the Niño3 region

(but not for broader equatorial regions) as opposed to other models for which the Bjerknes Index was calculated. These models, together with EC-Earth3 and EC-Earth3-Veg, are outliers among the 20 models analyzed in that they exhibit a very strong increase in ENSO amplitude in the warming experiments. Thus, such decrease in stability for the Niño3 region, driven by stronger positive feedbacks in the eastern Pacific, may explain why some models have a large increase in ENSO amplitude; however, this factor cannot explain the robust increase in ENSO across the vast majority of models and experiments.

We suggest here that atmospheric noise, including westerly wind bursts, may have a crucial role in driving projected changes in ENSO, which is supported by a strong correlation between changes in noise and ENSO SST variability, as found in this and other studies (Wengel et al. 2018; Lopez et al. 2022). However, the key question remains whether the increase in noise is in fact driving a stronger ENSO, or if ENSO itself generates a stronger noise in a warming climate (Kug et al. 2008; Eisenman et al. 2005). The former possibility is supported by the recent study of Liang and Fedorov (2022).

This result highlights a gap in our understanding of what drives ENSO changes in response to global warming in GCMs. For example, we cannot *consistently* explain those changes using a simple linear framework of the Bjerknes Stability Index that would link those changes to changes in the tropical mean state. This is exemplified by the fact that mean state changes are largest in the abrupt-4xCO₂ scenario, and yet by century-end ENSO amplitude increases more in the SSP5-8.5 scenario. In fact, 6 models in the former scenario actually show a weaker ENSO, but all models in the latter scenario suggest a stronger ENSO in a warming climate.

This shortcoming in our understanding of ENSO changes with global warming calls for the development of new tools to analyze projections of the ENSO response to warming in the future. One such method includes a more comprehensive linear stability analysis that would compute the full leading eigen modes of the system, making it possible to consider the sensitivity of unstable ENSO related modes to changes in the background state very accurately. Previously, this has been done for Zebiak-Cane (Zebiak and Cane 1987) type models (e.g. Fedorov and Philander 2000, 2001; An and Jin 2000), and ocean GCMs (Sévellec and Fedorov 2013, 2015), but yet to be completed for coupled GCMs. Another potential approach could be the application of machine learning techniques. In particular, neural networks could be used to disentangle nonlinear relationships between changes in certain background climatic variables and ENSO and to highlight which variables are most important (e.g. Mamalakis et al. 2022).
Overall, our results point towards a robust increase in ENSO activity in the SSP5-8.5 and SSP1-2.6 scenarios, but a common overarching mechanism to explain these changes is still lacking. Models that show a strong increase in ENSO amplitude, such as MIROC6 and MIROC-ES2L, may be driven by stability changes over the Niño3 region in combination with a minimal or no change in thermodynamic damping. In contrast, in models with a moderate change in ENSO, the linear Bjerknes Index decreases, primarily because of increased thermodynamic damping, resulting in a more stable system. This suggests that atmospheric noise and/or nonlinear changes may contribute to a stronger ENSO in these models. On the other hand, in CESM2, which exhibits a weaker ENSO in the 4xCO₂ experiment, the Bjerknes Index does show the largest reduction, driven by increased thermodynamic damping and a decrease in the thermocline feedback. Yet, the same model exhibits an ENSO strengthening, albeit small, in SSP5-8.5 and SSP1-2.6 relative to the piControl.

On average, ENSO amplitude also increases in idealized CO₂ scenarios, such as the 1pctCO₂ or abrupt-4xCO₂ scenario, although the change is less consistent across models. Furthermore, contrary to the findings of Callahan et al. (2021), we do not find evidence that models converge to a weaker ENSO in the abrupt 4xCO₂ experiments over longer time scales. In fact, ENSO remains stronger than the control in several models we considered even after 1000 years of computation. On the whole, the abrupt 4xCO₂ scenario produces the widest spread in ENSO projections, and as such may not be the most reliable indicator of changes to come. These findings underscore that, despite the robust strengthening of ENSO by century-end in a broad range of models and warming scenarios of CMIP6, there is still a large uncertainty in ENSO future response to global warming. This uncertainty should be addressed by evaluating ENSO drivers across multiple warming experiments and multiple models.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00382-023-06856-x.

**Author contributions** UKH and AVF contributed equally to designing the research. UKH performed the data analysis and, together with AVF, interpreted the results. UKH wrote the manuscript and edited it together with AVF.

**Funding** U.K.H. is supported by a NASA FINEST Fellowship (80NSSC20K1634). A.V.F. is supported by grants from NASA (80NSSC21K0558, 80NSSC18K1396), NOAA (NA20OAR4310377) and DOE (DE-SC0023134). Additional funding is provided by the ARCHANGE project (ANR-18-MPGA-0001, France). We also acknowledge a generous gift to Yale University from T. Sandoz. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**Data availability** CMIP6 data are available at https://esgf-node.llnl.gov/search/cmip6/. Code used for the analysis is available at https://github.com/ubbu36/ENSO_overview_in_CMIP6.

## Declarations

**Conflict of interest** The authors declare no competing interests.

## References

* An S-I, Jin F-F (2000) An eigen analysis of the interdecadal changes in the structure and frequency of ENSO mode. Geophys Res Lett 27(16):2573–2576. https://doi.org/10.1029/1999GL011090

* Brown JR, Brierley CM, An S-I, Guarino M-V, Stevenson S, Williams CJR, Zhang Q et al (2020) Comparison of past and future simulations of ENSO in CMIP5/PMIP3 and CMIP6/PMIP4 Models. Preprint. Climate Modelling/Modelling only/Pleistocene. https://doi.org/10.5194/cp-2019-155.

* Cai W, Santoso A, Wang G, Yeh S-W, An S-I, Cobb KM, Collins M, Guilyardi E, Jin F-F, Kug J-S (2015) ENSO and greenhouse warming. Nat Clim Chang 5(9):849–859

* Cai W, Wang G, Dewitte B, Lixin Wu, Santoso A, Takahashi K, Yang Y, Carréric A, McPhaden MJ (2018) Increased variability of Eastern Pacific El Niño under greenhouse warming. Nature 564(7735):201–206. https://doi.org/10.1038/s41586-018-0776-9

* Cai W, Santoso A, Collins M, Dewitte B, Karamperidou C, Kug J-S, Lengaigne M et al (2021) Changing El Niño-Southern Oscillation in a warming climate. Nat Rev Earth Environ 2(9):628–644. https://doi.org/10.1038/s43017-021-00199-z

* Cai W, Ng B, Wang G, Santoso A, Lixin Wu, Yang K (2022) Increased ENSO sea surface temperature variability under four IPCC emission scenarios. Nat Clim Chang 12(3):228–231. https://doi.org/10.1038/s41558-022-01282-z

* Callahan CW, Chen C, Rugenstein M, Bloch-Johnson J, Yang S, Moyer EJ (2021) Robust decrease in El Niño/Southern oscillation amplitude under long-term warming. Nat Clim Chang 11(9):752–757. https://doi.org/10.1038/s41558-021-01099-z

* Capotondi A, Sardeshmukh PD, Ricciardulli L (2018) The nature of the stochastic wind forcing of ENSO. J Clim 31(19):8081–8099

* DiNezio PN, Kirtman BP, Clement AC, Lee S-K, Vecchi GA, Wittenberg A (2012) Mean climate controls on the simulated response of ENSO to increasing greenhouse gases. J Clim 25(21):7399–7420

* Eisenman I, Lisan Yu, Tziperman E (2005) Westerly wind bursts: ENSO's tail rather than the dog? J Clim 18(24):5224–5238

* Eyring V, Bony S, Meehl GA, Senior CA, Stevens B, Stouffer RJ, Taylor KE (2016) Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. Geosci Model Dev 9(5):1937–1958. https://doi.org/10.5194/gmd-9-1937-2016
Fedorov AV (2002) The response of the coupled tropical ocean-atmosphere to westerly wind bursts. Q J R Meteorol Soc 128(579):1-23.
https://doi.org/10.1002/qj.200212857901

Fedorov AV, George Philander S (2000) Is El Niño changing? Science 288(5473):1997-2002

Fedorov AV, George Philander S (2001) A stability analysis of tropical ocean-atmosphere interactions: bridging measurements and theory for El Niño. J Clim 14(14):3086-3101. https://doi.org/10.1175/1520-0442(2001)014%3c3086:ASAOTO%3e2.0.CO;2

Fedorov AV, Shineng Hu, Lengaigne M, Guilyardi E (2015) The impact of westerly wind bursts and ocean initial state on the development, and diversity of El Niño events. Clim Dyn 44(5):1381-1401.
https://doi.org/10.1007/s00382-014-2126-4

Fedorov AV, Hu S, Wittenberg AT, Levine AFZ and Deser C (2020) ENSO low-frequency modulation and mean state interactions. El Niño Southern Oscillation in a Changing Climate, 173-98

Ferrett S, Collins M (2019) ENSO feedbacks and their relationships with the mean state in a flux adjusted ensemble. Clim Dyn 52(12):7189-7208. https://doi.org/10.1007/s00382-016-3270-9

Fredriksen H-B, Berner J, Subramanian AC, Capotondi A (2020) How does El Niño-southern oscillation change under global warming—a first look at CMIP6'. Geophys Res Lett 47(22):e2020GL090640.
https://doi.org/10.1029/2020GL090640

Guilyardi E (2006) El Niño-Mean state-seasonal cycle interactions in a multi-model ensemble. Clim Dyn 26:329-348

Guilyardi E, Wittenberg A, Fedorov A, Collins M, Wang C, Capotondi A, Oldenborgh GJV, Stockdale T (2009) Understanding El Niño in ocean-atmosphere general circulation models: progress and challenges. Bull Am Meteor Soc 90(3):325-340

Guilyardi E, Cai W, Collins M, Fedorov A, Jin F-F, Kumar A, Sun D-Z, Wittenberg A (2012) New strategies for evaluating ENSO processes in climate models. Bull Am Meteor Soc 93(2):235-238

Heede UK, Fedorov AV (2023) Colder eastern equatorial pacific and stronger walker circulation in the early 21st century: separating the forced response to global warming from natural variability. Geophys Res Lett 50(3):e2022GL101020. https://doi.org/10.1029/2022GL101020

Heede UK, Fedorov AV, Burls NJ (2020) Timescales and mechanisms for the tropical pacific response to global warming: a tug of war between the ocean thermostat and weaker walker. J Clim. https://doi.org/10.1175/JCLI-D-19-0690.1

Heede UK, Fedorov AV (2021) Eastern equatorial Pacific warming delayed by aerosols and thermostat response to CO2 increase. Nat Clim Change 11(8):696-703

Heede UK, Fedorov AV, Burls NJ (2021) A stronger versus weaker Walker: understanding model differences in fast and slow tropical Pacific responses to global warming. Clim Dyn 57(9-10):2505-2522

Hu S, Fedorov AV (2019) The extreme El Niño of 2015-2016: the role of westerly and easterly wind bursts, and preconditioning by the failed 2014 event. Clim Dyn 52(12):7339-7357

Huang P, Xie S-P (2015) Mechanisms of change in ENSO-induced tropical pacific rainfall variability in a warming climate. Nat Geo-sci 8(12):922-926. https://doi.org/10.1038/ngeo2571

Jin F-F (1997) An equatorial ocean recharge paradigm for ENSO. Part I: conceptual model. J Atmos Sci 54(7):811-829. https://doi.org/10.1175/1520-0469(1997)054%3c0811:AEORPF%3e2.0.CO;2

Jin F-F, Kim ST, Bejarano L (2006) A coupled-stability index for ENSO. Geophys Res Lett. https://doi.org/10.1029/2006GL027221

Kim ST, Cai W, Jin F-F, Jin-Yi Yu (2014) ENSO stability in coupled climate models and its association with mean state. Clim Dyn 42(11):3313-3321

Kohyama T, Hartmann DL, Battisti DS (2018) Weakening of nonlinear ENSO under global warming. Geophys Res Lett 45(16):8557-8567. https://doi.org/10.1029/2018GL079085

Kug J-S, Jin F-F, Sooraj KP, Kang I-S (2008) State-dependent atmospheric noise associated with ENSO. Geophys Res Lett. https://doi.org/10.1029/2007GL032017

Liang Y and Fedorov AV (2022) Stronger westerly wind bursts under global warming: the roles of the warming pattern, madden-julian oscillation and tropical cyclones. 2022 (December): A15B-05

Lopez H, Lee S-K, Kim D, Wittenberg AT, Yeh S-W (2022) Projections of faster onset and slower decay of El Niño in the 21st century. Nat Commun 13(1):1915. https://doi.org/10.1038/s41467-022-29519-7

Maher N, Jnglin RC, Wills PD, Klavans J, Milinski S, Sanchez SC, Stevenson S, Stuecker MF, Xian Wu (2023) The future of the El Niño-southern oscillation: using large ensembles to illuminate time-varying responses and inter-model differences. Earth Syst Dyn 14(2):413-431. https://doi.org/10.5194/esd-14-413-2023

Mamalakis A, Ebert-Uphoff I and Barnes EA (2022) Explainable artificial intelligence in meteorology and climate science: model fine-tuning, calibrating trust and learning new science. In:XxAI - Beyond Explainable AI: International Workshop, Held in Conjunction with ICML 2020, July 18, 2020, Vienna, Austria, Revised and Extended Papers, edited by Andreas Holzinger, Randy Goebel, Ruth Fong, Taesup Moon, Klaus-Robert Müller, and Wojciech Samek, 315-39. Lecture Notes in Computer Science. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-031-04083-2_16

Manucharyan GE, Fedorov AV (2014) Robust ENSO across a wide range of climates. J Clim 27(15):5836-5850. https://doi.org/10.1175/JCLI-D-13-00759.1

McPhaden MJ, Santoso A, Cai W (2020) El Niño southern oscillation in a changing climate, vol 253. John Wiley & Sons

O'Neill BC, Tebaldi C, van Vuuren DP, Eyring V, Friedlingstein P, Hurtt G, Knutti R et al (2016) The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6. Geosci Model Dev 9(9):3461-3482. https://doi.org/10.5194/gmd-9-3461-2016

Philip S, Jan van Oldenborgh G (2009) Significant atmospheric non-linearities in the ENSO cycle. J Clim 22(14):4014-4028. https://doi.org/10.1175/2009JCLI2716.1

Puy M, Vialard J, Lengaigne M, Guilyardi E, DiNezio PN, Voldoire A, Balmaseda M, Madec G, Menkes C, Mcphaden MJ (2019) Influence of westerly wind events stochasticity on El Niño amplitude: the case of 2014 vs. 2015. Clim Dyn 52(12):7435-7454. https://doi.org/10.1007/s00382-017-3938-9

Seager R, Cane M, Henderson N, Lee D-E, Abernathey R, Zhang H (2019) Strengthening tropical pacific zonal sea surface temperature gradient consistent with rising greenhouse gases. Nat Clim Chang 9(7):517-522

Sevellec F, Fedorov AV (2013) The leading, interdecadal eigenmode of the atlantic meridional overturning circulation in a realistic ocean model. J Clim 26(7):2160-2183. https://doi.org/10.1175/JCLI-D-11-00023.1

Sevellec F, Fedorov AV (2015) Unstable AMOC during glacial intervals and millennial variability: the role of mean sea ice extent. Earth Planet Sci Lett 429:60-68

Trenberth KE (2012) Framing the way to relate climate extremes to climate change. Clim Change 115(2):283-290. https://doi.org/10.1007/s10584-012-0441-5

Vecchi GA, Soden BJ, Wittenberg AT, Held IM, Leetmaa A, Harrison MJ (2006) Weakening of tropical pacific atmospheric circulation due to anthropogenic forcing. Nature 441(7089):73-76

Wengel C, Dommenget D, Latif M, Bayr T, Vijayeta A (2018) What controls ENSO-amplitude diversity in climate models? Geophys Res Lett 45(4):1989-1996. https://doi.org/10.1002/2017GL076849

Wengel C, Lee S-S, Stuecker MF, Timmermann A, Chu J-E, Schloesser F (2021) Future high-resolution El Niño/southern oscillation
dynamics. Nat Clim Chang 11(9):758-765. https://doi.org/10.1038/s41558-021-01132-4

Xie S-P, Deser C, Vecchi GA, Ma J, Teng H, Wittenberg AT (2010) Global warming pattern formation: sea surface temperature and rainfall. J Clim 23(4):966-986

Yeh S-W, Cai W, Min S-K, McPhaden MJ, Dommenget D, Dewitte B, Collins M et al (2018) ENSO atmospheric teleconnections and their response to greenhouse gas forcing. Rev Geophys 56(1):185-206. https://doi.org/10.1002/2017RG000568

Yu S, Fedorov AV (2020) The role of westerly wind bursts during different seasons versus ocean heat recharge in the development of extreme El Niño in climate models. Geophys Res Lett 47(16):e2020GL088381. https://doi.org/10.1029/2020GL088381

Yun K-S, Lee J-Y, Timmermann A, Stein K, Stuecker MF, Fyfe JC, Chung E-S (2021) Increasing ENSO-rainfall variability due to changes in future tropical temperature-rainfall relationship. Commun Earth Environ 2(1):1-7. https://doi.org/10.1038/s43247-021-00108-8

Zebiak SE, Cane MA (1987) A model El Niño-southern oscillation. Mon Weather Rev 115(10):2262-2278. https://doi.org/10.1175/1520-0493(1987)115%3c2262:AMENO%3e2.0.CO;2

Zhao B, Fedorov A (2020) The Effects of Background Zonal and Meridional Winds on ENSO in a Coupled GCM. J Clim 33(6):2075-2091. https://doi.org/10.1175/JCLI-D-18-0822.1

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.