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

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 1pctCO2 and abrupt 4xCO2 experiments. Several models, however, do predict no change or ENSO weakening, especially in the idealized experiments. Critically, the strongest forcing (abrupt-4xCO2) 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.


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 1 3 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. 2009Guilyardi et al. , 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. 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;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 . 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  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 2 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 2 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 2 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.

CMIP6 archive
We analyze two types of experiments from the Climate Model Intercomparison Project Phase 6 (CMIP6) archive ). The first type consists of two hypothetical or idealized CO 2 -only experiments: abrupt-4xCO 2 rise and 1pctCO 2 gradual CO 2 rise (1% per year) where CO 2 -increases are relative to the pre-industrial level of 280 ppm. The second type includes more realistic, fullforcing 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 2 and 1pctCO 2 ), where ensemble members are 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.

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-scipy 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 , 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 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

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 −1 ) 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: Here, CD is damping by the mean currents given by: Angle brackets indicate a spatial average for the centraleastern Pacific (180-280 °E) and bars indicate a time-mean over the chosen period of the experiment. u, v and 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 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: Here, 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: To have units of yr −1 , is further multiplied by 1 where 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: where a and * a are the coefficients of linear regression of wind stress onto temperature anomalies. u , w , and h are the coefficients of linear regression of zonal current, vertical current and thermocline anomalies, respectively, onto wind stress anomalies. T x is the mean zonal temperature gradient and T z is the mean vertical stratification. w is the mean vertical velocity. All regression coefficients are for the centraleast equatorial Pacific (180-280 °E), except * 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.

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: where (x, y, t) denotes stochastic forcing by random wind stress variations, 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 over the equatorial Pacific region (130-280 °E, 5 °S-5 °N) and compute the standard deviation over time.

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 2 scenario is also well correlated with the SSP5-8.5 response, but it is more muted. On average, the abrupt-4xCO 2 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 2 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 2 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 2 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 2 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 2 experiment (R = − 0.56, p = 0.05). Other correlations are much smaller and not significant. Using the Nino3.4 index, E-index or C-index did not improve these correlations (Supplementary Figs. 1 and 2). Correlations between changes in ENSO rainfall variability and the mean zonal SST gradient in the gradual experiments (1pctCO 2 and SSP5-8.5) are generally higher (Fig. 5a-c) than those for ENSO SST variability, except for the abrupt-4xCO 2 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 of extreme events versus changes in ENSO SST amplitude in the 1pctCO 2 , abrupt4xCO 2 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 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 2 , 1pctCO 2 , SSP5-8.5 and SSP1-2.6 experiments, respectively. Thin errorbars for SSP experiments represent the maximum and minimum values across three ensemble members 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.

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 2 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 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 2 , abrupt-4xCO 2 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-ent relaxes in the 1pctCO 2 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 2 experiment (R = − 0.56)

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 2 , abrupt-4xCO 2 and SSP5-8.5 scenarios, respec-tively. 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 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 2 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 mag-nitude of the Bjerknes Index. All feedback terms and the BJ index are in units of year −1 . Noise is in N m −2 and ENSO SST amplitude is in o 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 nonlinearities that could counteract increased linear stability, we calculate changes to wind stress-SST coupling ( a ) and thermodynamic damping ( ) for three different ranges of temperature anomalies (corresponding to large positive, moderate, and large negative values, respectively). We find that the thermodynamic damping , 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 , 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 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, 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 2 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.

Long-term warming experiments
Lastly, we investigate three abrupt-4xCO 2 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.

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 2 increase exhibit a stronger ENSO response to warming on average than the abrupt-4xCO 2 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.
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 equatorial 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 disappears in other types of experiments considered. Likewise, we find some correlation between ENSO SST amplitude and ocean stratification in the abrupt 4xCO 2 experiment (R = − 0.56, p = 0.05), but it decreases in the SSP5-8.5 scenario (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 consequences for adaptation and mitigation because changes in drought or flooding conditions may cause more damage than changes in SST and because ENSO remote teleconnections 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 precipitation to SST anomalies (Yun et al. 2021). However, we find that changes in the background state alone cannot fully 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 2 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 2 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 andPhilander 2000, 2001;An and Jin 2000), and ocean GCMs (Sévellec andFedorov 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 2 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 2 scenarios, such as the 1pctCO 2 or abrupt-4xCO 2 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 2 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 2 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.
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 FINESST 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:// creat iveco mmons. org/ licen ses/ by/4. 0/.