Rapid communication: Nonlinear sensitivity of the El Niño–Southern Oscillation across climate states

The El Niño–Southern Oscillation (ENSO) is the dominant mode of tropical climate variability. Understanding its sensitivity to climate states is of societal and ecosystem importance given the unabated global warming. Palaeoclimate archives and climate models suggest that ENSO activity depends on mean state conditions. However, a common mechanism that can predict ENSO variability under a range of background conditions remains elusive. Here we combine climate model simulations of past climates and future warming, performed under the Paleoclimate Modelling Intercomparison Project (PMIP) and the Coupled Model Intercomparison Project (CMIP), to evaluate ENSO activity throughout a wide range of climate states. We find that the sensitivity of ENSO to the background climate is nonlinear and tied to the climatological position of the tropical Pacific convection centres, namely the Intertropical Convergence Zone (ITCZ) and the South Pacific Convergence Zone (SPCZ). Simulations with atmospheric CO₂ lower than today display a poleward shift in the convergence zones and weakened ENSO. Moderate equatorward shifts in the convergence zones occur under moderate CO₂-induced warming increasing ENSO activity, while strong equatorward shifts reduce ENSO variability in extreme CO₂ warming scenarios, resulting in a permanent El Niño-like mean state. Our results provide a comprehensive mechanism of how tropical Pacific mean state modulates ENSO activity.

1 Introduction

The El Niño–Southern Oscillation (ENSO) is a coupled ocean–atmosphere zonal oscillation sourced in the equatorial Pacific that results in either higher (El Niño) or lower (La Niña) sea surface temperatures (SSTs) in the central eastern Pacific (EP), with strong impacts worldwide (Ropelewski and Halpert, 1987). For this reason, the scientific community has concentrated efforts in identifying changes in ENSO variability and their drivers in past climates and future projections (Cai et al., 2021; Ford et al., 2015; Pontes et al., 2022b; Thirumalai et al., 2024; White and Ravelo, 2020).

The past ∼ 3 200 000 years (∼ 3200 kyr)) encompasses key periods that have been extensively studied, as they can improve our understanding of climate dynamics and help constrain future projections (Burke et al., 2018; Thirumalai et al., 2024). Firstly, the mid-Pliocene warm period (referred to as the mid-Pliocene, 3.205 ka) featured a global mean temperature 2–4 °C higher than pre-industrial and elevated atmospheric CO₂ concentrations (∼ 400 ppmv) (Lunt et al., 2012). Over this period, proxy data do not provide a clear picture on ENSO changes and suggest the important role of the sensitivity of the equatorial thermocline in maintaining ENSO activity (Watanabe et al., 2011; White and Ravelo, 2020). Climate models, however, systematically indicate reduced variability driven by a northward shift in the Intertropical Convergence Zone (ITCZ) (Pontes et al., 2022b). Other periods of interest are the Last Interglacial (LIG; ∼ 129–116 kyr), which was the warmest interglacial of the past 800 kyr (PAGES working Group, 2016), and the Mid-Holocene Thermal Maximum (MHTM; ∼ 8–4 kyr). During both the LIG and the MHTM, atmospheric CO₂ concentrations were slightly lower than pre-industrial (PI) levels, but boreal summer insolation in the Northern Hemisphere was higher (up to 70 and 25 W m⁻², respectively). For both periods, most proxies and models agree on reduced ENSO activity (Brown et al., 2020; Emile-Geay et al., 2016), possibly related to weaker thermocline feedback (Chen et al., 2019), amplified by a stronger Walker circulation (Pausata et al., 2017). Lastly, during the Last Glacial Maximum (LGM; ∼ 21 ka), when CO₂ concentrations were ∼ 190 ppmv (Marcott et al., 2014), a deeper equatorial mixed layer is suggested to have driven reduced ENSO variability (Ford et al., 2015; Thirumalai et al., 2024).

While palaeo-modelling studies indicate reduced ENSO variability, future projections suggest enhanced ENSO variability under anthropogenic warming for the 21st century, which is attributed to enhanced upper-ocean stratification (Cai et al., 2018). Nevertheless, extended simulations of the high-emissions scenario (SSP5-8.5) until 2300 and equilibrated long-term (over 1000 years) CO₂-driven simulations (e.g. abrupt 4 and 8 times CO₂) under the Long Run Model Intercomparison Project (LongRunMIP) show a robust decrease in ENSO (Callahan et al., 2021; Peng et al., 2024; Zheng et al., 2022). ENSO decline under these scenarios has been attributed to collapsed eastern Pacific upwelling and weakened zonal SST gradient.

Studies of specific time slices and scenarios led to significant advancements in the understanding of ENSO's dynamical processes. These studies usually suggest that changes in the equatorial Pacific mean state (e.g. upwelling, zonal SST gradient, and ITCZ) are responsible for changing ENSO activity, thus contributing to the hypothesis that ENSO is modulated by the background state. Under this hypothesis, a strong mean circulation (i.e. La Niña-like mean state) inhibits the development of wind and SST anomalies, effectively suppressing ENSO variability, while a weak mean circulation (i.e. El Niño-like mean state) favours the development of anomalies that can lead to ENSO events. However, the evaluation of a single scenario provides a limited view of the broader ENSO response across climate states. To comprehensively evaluate whether and how ENSO activity can be modulated by the mean state, a consistent framework that covers a wide range of background states is needed. In this context, we combine climate model output from the Coupled Intercomparison Project phase 6 (CMIP6) and Paleoclimate Modelling Intercomparison Project phase 4 (PMIP4), covering seven different scenarios, to investigate the role of the Pacific Ocean mean state in modulating ENSO variability.

2 Data and methods

Here, we analyse four past climate scenarios part of PMIP4: the mid-Pliocene (KM5c interglacial time slice; Haywood et al., 2020), the Last Interglacial (LIG, 127 ka; Otto-Bliesner et al., 2021), the Mid-Holocene Thermal Maximum (MHTM; 6 ka; Brierley et al., 2020), and the Last Glacial Maximum (LGM; 21 ka; Kageyama et al., 2021). In addition, we incorporate in our analysis extended future projections of the high- (SSP5-8.5) and low-emissions (SSP1-2.6) scenarios (O'Neill et al., 2016). In SSP5-8.5, anthropogenic radiative forcing continues to increase beyond the year 2100, reaching a level of ∼ 12 W m⁻² in 2250 and stabilising afterwards. In SSP1-2.6, radiative forcing reaches 2.6 W m⁻² in 2100 and then slowly declines, stabilising at around 2.0 W m⁻² onwards from 2200. Finally, we also include the hypothetical scenario that abruptly increases atmospheric CO₂ to 4 times pre-industrial levels (abrupt-4xCO2). These data total 83 simulations covering seven different scenarios (Table S1 in the Supplement). We use each model's pre-industrial (PI) control simulation (CO₂ = 284 ppmv) as reference for all scenarios. For all palaeoclimate experiments (including abrupt4xCO2 and piControl), we use the last 100 years of each model's simulation. For the future scenarios, SSP5-8.5 and SSP1-2.6, we use the time slice from years 2251 to 2300 (referred to as SSP585-23 and SSP126-23, respectively) as representative of a near-equilibrium state.

Given that ENSO is a complex phenomenon of ocean–atmosphere interaction, not all climate models accurately simulate key processes of ENSO dynamics, resulting in poor ENSO simulations. To avoid being misled by results from models that do not capture ENSO dynamics properly, we apply a selection criterion based on the nonlinear convective feedback, which contributes to the observed positive skewness of SST anomalies in the eastern Pacific (Appendix A) (Dommenget et al., 2012). Here, we use data from the Global Precipitation Climatology Project (Adler et al., 2003) and the Hadley Centre Sea Surface Temperature (Rayner et al., 2003) datasets to calculate the strength of this feedback in observations. This results in 61 selected simulations of past climates and future scenarios (Table S1).

We use the standard deviation of SST anomalies in the Niño3 region (5° S–5° N, 150–90° W) as a measure of ENSO amplitude (Taschetto et al., 2014). SST time series were linearly detrended to remove possible model drift in all simulations. SST anomalies are computed by removing the mean annual cycle at a monthly resolution.

Lastly, we evaluate the nonlinear ENSO sensitivity to the background state through the nonlinear coefficient (a) of the quadratic regression model ($y\left(x\right)=a{x}^{\mathrm{2}}+bx+g$), while the precision of the nonlinear relationship is assessed through the coefficient of determination (R²).

3 Results and discussion

3.1 Mean state changes and ENSO activity

All three past warm climates (mid-Holocene, Last Interglacial, and mid-Pliocene) share similar conditions in the tropical Pacific Ocean, though for different reasons (Brierley et al., 2020; Otto-Bliesner et al., 2021; Pontes et al., 2020). Firstly, the ITCZ shifts northward in all warm palaeoclimate scenarios considered due to increased differences in inter-hemispheric warming. In the mid-Pliocene simulations, this is attributed to enhanced albedo feedback in the Northern Hemisphere due to large sea ice and land ice losses and enhanced vegetation area, while increased insolation in the Northern Hemisphere is the main driver in the MHTM and LIG (Brierley et al., 2020; Otto-Bliesner et al., 2021; Pontes et al., 2020). An ITCZ displaced northward allows intensified equatorial trade winds in the central western Pacific (Fig. 1a–c), which cause SSTs to be lower in the equatorial Pacific than in the surrounding area due to the wind–evaporation–SST feedback (Vimont et al., 2001), thus strengthening the equatorial upwelling (Fig. 1a–c). Such climatic anomalies resemble a La Niña-like mean state, associated with reduced ENSO variability (Fig. 1g; Table S1).

Figure 1 Climate anomalies and ENSO variability change. Multi-model annual mean changes in SST and wind field for the mid-Holocene (a), Last Interglacial (b), mid-Pliocene (c), SSP126-23 (d), abrupt-4xCO2 (e), and SSP585-23 (f). The colour bar varies for each panel to highlight regions where SST changes are higher (warm colours) and lower (cold colours) than the mean tropical Pacific warming in each ensemble, making changes in SST gradients easily identifiable. Arrows represent changes in wind stress in all panels but panel (c), where they represent changes in surface winds due to data availability. Arrows are plotted where there is a significant change in either zonal or meridional component (70 % model agreement in the sign of the change). Purple lines indicate the multi-model mean 8 mm d⁻¹ climatological rainfall contours for each time slice as an estimate of the ITCZ and SPCZ positions. (g) Change in Niño3 standard deviation for each simulation used in this study. The pre-industrial control simulation is used as reference when quantifying changes. Analyses shown in this figure include results from all models, with no selection criteria being applied. Values for each simulation can be found in Table S1. Maps in this figure were plotted using the Cartopy Python library.

In contrast, CO2-driven scenarios (SSP126-23, SSP585-23, and abrupt4xCO2) resemble an El Niño-like mean state (Fig. 1d). It is interesting to note that, as the SST warming intensifies, both the ITCZ and the South Pacific Convergence Zone (SPCZ) tend to move equatorward (Fig. 1d–f), progressively weakening the equatorial trades and eastern upwelling. In moderate CO₂-driven warming scenarios (e.g. SSP126-23), these changes allow enhanced ENSO variability (Fig. 1g). Interestingly, however, ENSO variability tends to decrease under strong warming scenarios (e.g. SSP585-23), when the ITCZ and SPCZ are nearly collapsed at the Equator (Fig. 1e–g), thus resulting in a nonlinear ENSO response to the background state. As the position of the convergence zones can be directly modulated by changes in the atmospheric meridional heat transport (Pontes et al., 2020; Schneider et al., 2014), which provides a linkage between equatorial Pacific and large-scale changes (e.g. orbital forcing and ice sheets), we use them as a proxy for changes in the tropical Pacific mean state.

3.2 Role of the Pacific convergence zones

To evaluate the effect of the position of the convergence zones on ENSO activity, we develop an index that captures their combined displacement during the developing and mature ENSO phases (austral spring and summer) (Pontes et al., 2022a). The ITCZ and SPCZ positions are computed as the precipitation-weighted average over latitudes in which precipitation is greater than 70 % of the maximum zonally averaged precipitation in each hemisphere (0–20° N and 20–0° S). This methodology captures migrations of the ITCZ and the SPCZ independently from one another. As our objective is to quantify their overall displacement relative to the Equator (e.g. equatorward shift in the ITCZ and SPCZ represent an El Niño-like mean state), we consider their absolute values. The index D is obtained through quantifying their combined displacement relative to the model's pre-industrial position:

$$\begin{array}{}\text{(1)}& D=\left[\left|{\mathrm{ITCZ}}_{\mathrm{S}}\right|-\left|{\mathrm{ITCZ}}_{\mathrm{PI}}\right|\right]+\left[\left|{\mathrm{SPCZ}}_{\mathrm{S}}\right|-\left|{\mathrm{SPCZ}}_{\mathrm{PI}}\right|\right],\end{array}$$

where the subscript “S” indicates the position of the ITCZ and SPCZ in the perturbation scenarios and the subscript “PI” denotes their position in their respective pre-industrial simulation. Hereafter, the index D is referred to as the convection centre index and is, by definition, positive for poleward displacements of the convergence zones. It is important noting that, in all scenarios, the Pacific ITCZ lies in the Northern Hemisphere; therefore a poleward movement reflects a northward shift, whereas, for the SPCZ, a poleward displacement reflects a southward shift.

It is important to note that double-ITCZ biases may affect the calculated SPCZ position. The double-ITCZ bias is an artificial feature produced by most climate models that overestimates the tropical precipitation south of the Equator in the central eastern Pacific. This biased precipitation artificially induces the index D to capture an apparent SPCZ displaced northwards. Nevertheless, although being an artificial feature, the double ITCZ dynamically influences the tropical Pacific climate through simulating weaker trade winds and warmer SSTs than in observations, thereby affecting the model's ENSO response. We thus do not differentiate between the southern branch of the double ITCZ and the SPCZ in our evaluation.

We find that the ENSO–convection centre relationship is nonlinear (Fig. 2 and Table S1). This relationship indicates that there are two mean states that tend to inhibit ENSO development and that there is an optimal distance at which the convection centres must be from the Equator to maximise ENSO variability. Overall, poleward (D > 0°) and strong equatorward (D < −8°​​​​​​​) shifts in the convection centres, relative to their position in the PI climate, are associated with reduced ENSO variability, whereas moderate equatorward shifts (−8° < D < 0°) enhance ENSO activity. An analysis of how this relationship affects central Pacific (CP) and eastern Pacific ENSO types and the 21st-century transient ENSO response is included in Appendix B.

Figure 2 Nonlinear ENSO sensitivity across climate states. Relationship between the convection centre index and the change in Niño3 amplitude, measured by its standard deviation (SD). The solid black line indicates the quadratic fit based on the least-squares method. Banding indicates the 95 % confidence interval based on a 1000-sample bootstrap. The mean displacement of the convection centres during boreal spring and summer is considered (i.e. encompassing developing and mature ENSO phases). R² indicates the coefficient of determination, and a indicates the nonlinear coefficient of the quadratic regression model. Error estimates for R² and a were calculated as 1 standard deviation based on 1000 bootstrap realisations. The convection centre index is, by definition, positive for poleward movements of the convergence zones referenced at each model's pre-industrial position. Values for each model are presented in Table S1.

In exploring the causes of the nonlinear sensitivity of ENSO to background state, we find that the position of the convergence zones modulates both atmospheric and oceanic processes important for ENSO development. Firstly, we investigate the rainfall response to the convection centre migration through evaluating the frequency of extreme rainfall events (> 5 mm d⁻¹) in the eastern Pacific (Niño3 region) associated with the displacement of the convergence zones (Fig. 3a). We find that poleward shifts in the convergence zones reduce the frequency of extreme events from one event per 15.3 ± 7.5 years in the PI to one event per 26.9 ± 18.3 years, indicating suppressed convective feedback (Fig. 3a). The convective feedback is further suppressed as the convergence zones shift poleward until a full La Niña-like mean state is reached and no extreme El Niño events occur, which is illustrated by the CCSM4-Utrecht mid-Pliocene simulation (Fig. 3a; Table S1). Moderate equatorward shifts in the convection centres allow more intense convective feedback and increased frequency of extreme events (one event per 2.9 ± 0.4 years; Fig. 3a). This is consistent with studies which suggested that, as ENSO's frequency approaches the annual cycle frequency (represent here by one event per year), ENSO variability tends to increase (Liu, 2002; Timmermann et al., 2007), a process known as “ENSO frequency entrainment”. Strong equatorward displacements of the convection centres (−12° < D < −8°) result in more intense convective feedback, associated with reduced SST variability (Fig. 3a). In this case, ENSO frequency is nearly fully entrained in the annual cycle at 1.2 ± 0.06 events per year, thus making El Niño events more predictable. Finally, full equatorward shifts in the convergence zones (D < −12°) indicate a “permanent El Niño” situation, and ENSO's frequency is fully entrained into the annual cycle (one event per year), thus indicating these events are not related to intensified anomalous convective feedback but are due to the convergence zones being within the equatorial band.

Figure 3 Equatorial Pacific climate and ENSO feedbacks. (a) Nonlinear convective feedback: relationship between DJF Niño3 SST anomalies and DJF Niño3 rainfall. The first 100 years of 1000 bootstrap realisation is shown to illustrate the relationship in each group: pre-industrial control and moderate equatorward (−8° < D < 0°), strong equatorward (−12° < D < −8°), full equatorward (D < −12°), poleward (D > 0°), and strong poleward (represented by CCSM4-Utrecht mid-Pliocene simulation) displacements. The frequency of extreme events is indicated by the bootstrap mean frequency of events that exceeded 5 mm d⁻¹ (dashed orange line) across all realisations. Error estimates are indicated by the standard deviation of the realisations. (b) Inter-model relationship between the convection centre index and the change in easterly wind variability in the western equatorial Pacific (5° S–5° N, 160–210° E). (c) Inter-model relationship between the convection centre index and the change in easterly wind variability. (d) Inter-model relationship between the convection centre index and the wind–thermocline coupling coefficient. R² indicates the coefficient of determination of the quadratic regression model. Error estimates are given by the standard deviation of 1000 bootstrap realisations. Results shown in this figure were obtained through analyses of the subset models that best simulate the equatorial nonlinear convective feedback (Table S1).

Another important process in ENSO dynamics is the easterly wind variability in the western Pacific (Niño4 region: 160° E–150° W). The intensity of the winds in this region has been shown to determine the amplitude of their variability (Pontes et al., 2022b), which is related to the intensity and frequency of westerly wind bursts, important for ENSO initiation (Chen et al., 2015). In this context, we find that the position of the convection centres also modulates the wind variability through a quadratic relationship (R² = 0.51; Fig. 3c). The two scenarios that show reduced wind variability are consistent with reduced ENSO activity. Poleward movement of the convection centres increases the horizontal scale of the wind flow that reaches the western Pacific, generating wind anomalies in the South Pacific Subtropical High region, indicating that this new regime is more geostrophic and therefore more linearly balanced and deterministic (Fig. 1a–c). On the other hand, strong equatorward shifts in the convection centres require weak horizontal wind intensity within the equatorial band, thus suppressing the development of wind anomalies (Fig. 3b). Despite the consistency between changes in easterly wind variability and the position of the convergence zones, changes in wind variability cannot solely explain the nonlinear sensitivity of ENSO across climate states (Fig. 3b). The relationship between changes in easterly wind variability and ENSO is linear and does not distinguish between mean states of reduced (e.g. palaeo) and amplified (e.g. future) equatorial warming.

Finally, to investigate the possible modulation of the position of the convergence zones on the ocean–atmosphere dynamical coupling, we analyse the response of the wind–thermocline coefficient (Appendix C) (Jin et al., 2006). This parameter measures the sensitivity of thermocline slope anomalies to wind stress anomalies, which, during El Niño events, result in eastward heat advection by downwelling equatorial Kelvin waves (Timmermann et al., 2018). Our results indicate that the displacement of the convection centres also modulates the wind–thermocline coefficient through a quadratic relationship (R² = 0.45; Fig. 3d). Through affecting the horizontal wind field, the position of the convergence zones underpins the momentum transfer from the atmosphere, influencing oceanic stratification and the thermocline slope in the equatorial Pacific. As strong equatorward shifts in the convergence zones drive weak trades and collapsed eastern upwelling, they reduce the effectiveness of wind anomalies in generating swings of the thermocline. Furthermore, the propagation of Kelvin waves in a flat thermocline does not effectively promote oscillations between El Niño and La Niña anomalies. Poleward shifts in the convection centres increase the climatological thermocline slope, requiring stronger wind anomalies to promote thermocline oscillation and thus reducing the effectiveness of dynamical air–sea coupling.

4 Conclusion

Our findings suggest that meridional shifts in the convergence zones modulate key processes for ENSO development, such as equatorial convection, trade wind variability, and zonal thermocline oscillations. The interplay of these processes results in a nonlinear ENSO sensitivity to mean states, encompassing three distinct background conditions deviating from the pre-industrial climate (Fig. 4): (1) poleward migrations of the convection centres intensify the equatorial trades, weaken the convective feedback, increase momentum transfer to the upper ocean, reduce ocean stratification, and ultimately result in a weaker dynamical coupling, thus dampening ENSO variability and resembling a La Niña-like mean state (Fig. 4c); (2) moderate equatorward shifts in the convection centres (−8° < D < 0°) reduce the intensity of the equatorial trades and increase upper-ocean stratification, consequently amplifying the dynamical air–sea coupling and permitting equatorial trades and the thermocline slope to be rapidly reversed, thus enhancing ENSO activity and extreme rainfall events (Fig. 4b); (3) strong equatorward shifts in the convection centres create a permanent El Niño-like mean state in the eastern Pacific (Fig. 4a). The fact that the convection centres lie at the Equator does not allow momentum transfer to the ocean, resulting in a highly stratified ocean and dampened dynamical coupling. Nonetheless, this scenario is associated with intense warming in the eastern Pacific that enhances the climatological thermodynamical coupling, where climatological high SSTs allow continuous intense rainfall.

Figure 4 Mechanisms for changes in ENSO activity across climate states. Schematic of mean state features associated with ENSO variability. (a) A strong equatorward shift in the convection centres results in a flat thermocline, a highly stratified upper ocean, and weak dynamical coupling, thus reducing ENSO activity. (b) A moderate equatorward shift in the convection centres weakens the equatorial coupled circulation and increases the dynamical coupling, increasing thermocline swings and ENSO variability. (c) A poleward shift in convection centres enhances equatorial coupled circulation and reduces the dynamical coupling and ENSO activity. All comparison statements are relative to the pre-industrial climate. Maps in this figure were plotted using the Cartopy Python library.

The proposed mechanism can reconcile divergences in ENSO responses in proxy data and modelling studies. Firstly, modelling results from the idealised abrupt4xCO2 simulations indicate inconsistent inter-model ENSO response. We showed that this is related to the sensitivity of ITCZ and SPCZ migrations in each model. This argument can be expanded for meltwater-induced AMOC weakening modelling studies and proxy data covering Heinrich stadials and the Younger Dryas, which often do not agree on the ENSO response (Glaubke et al., 2024; Liu, 2002; Timmermann et al., 2007). This could be related to the different responses of the Pacific convergence zones to the intensity and duration of meltwater pulses. Finally, our results provide novel insights into the relationship between ENSO frequency and the annual cycle (Liu, 2002), suggesting that a moderate equatorward shift in the convergence zones relaxes both mean horizontal and mean vertical circulations, allowing nonlinear interactions between the annual cycle and ENSO's frequency, thus increasing ENSO's frequency towards the annual cycle. Strong equatorward displacements of the convergence zones result in a strong annual cycle associated with an El Niño-like mean state, thus fully entraining ENSO's frequency into the annual cycle and reducing ENSO's interannual variance.

Our modelling assessment suggests that the projected increase in ENSO variability has likely not occurred in the past ∼ 3.3 Myr. While there are significant uncertainties in both palaeoclimate simulations and palaeo-proxy records, both tend to suggest reduced ENSO variability during the LGM, the LIG, and the mid- to late Holocene (Brown et al., 2020; Emile-Geay et al., 2016; Sachs et al., 2009). While these results agree with our understanding of past ITCZ changes (Sachs et al., 2018; Schneider et al., 2014), past dynamics of the SPCZ need to be better understood and constrained to enhance our ability to predict future ENSO behaviour.

Appendix A

The criteria used to select models that correctly capture the strength of the convective feedback are based on an essential definition of extreme ENSO-related rainfall events, which are defined as precipitation events greater than 5 mm d⁻¹ in the Niño3 region. In observations, the 5 mm d⁻¹ rainfall rate is achieved at an SST anomaly of 2 °C, which gives a convective feedback of 2.5 mm d⁻¹ °C⁻¹ (Fig. A1 and Table S1). To ensure that models capture the observed strength, they are required to simulate convective feedback greater than 2 mm d⁻¹ °C⁻¹ in their pre-industrial runs. Models that do not correctly simulate this feedback tend to simulate SSTs well below or above the convective threshold of 26–28 °C (Johnson and Xie, 2010), resulting in unrealistic convective feedback strength. The model's ability to properly simulate the convective skewness filters out models that systematically simulate overly wet (i.e. double ITCZs) or dry (i.e. overly strong cold tongue) conditions in the eastern equatorial Pacific.

Figure A1 Model selection criteria: convective feedback. Relationship between DJF Niño3 SST anomalies and DJF Niño3 rainfall in observations, selected models, and non-selected models. The nonlinear convective criterion is used to select models for the analyses shown in Figs. 2 and 3 in the main article.

Appendix B

To evaluate the transient response of ENSO activity to the position of the convergence zones during the 21st century, we use the SSP3-7.0 scenario. As we compare ENSO variability with the climatological position of the convergence zones, this scenario shows a less intense transient forcing than the SSP5-8.5, thus potentially making this comparison less inaccurate. To reduce this uncertainty, we use the time slice from 2051 to 2100. Results show that the climatological position of the convergence zones also modulates ENSO activity during a transient warming (Fig. B1). Future anthropogenic warming under the SSP3-7.0 scenario tends to displace the convection centres equatorward by up to 4.2°, lying in the group of “enhanced ENSO activity”. Consistently, all models but one simulate increased ENSO activity in the second half of the 21st century under the SSP3-7.0 scenario (Fig. B1).

Figure B1 Relationship between changes in ENSO variability and the convection centre index. The solid black line indicates the quadratic fit based on the least-squares method. Banding indicates the 95 % confidence interval based on a 1000-sample bootstrap. R² indicates the coefficient of determination, and a indicates the nonlinear coefficient of the quadratic regression model. Error estimates for R² and a are calculated as 1 standard deviation based on 1000 bootstrap realisations.

Figure B2 Model selection criteria evaluation: Bjerknes feedback. Relationship between the first and second principal components of SST anomalies in the tropical Pacific in observation (a), selected (b), and non-selected (c) models. The nonlinear coefficient of the quadratic fitted model is indicated by a.

To account for ENSO complexity, we evaluate the central Pacific (CP) and eastern Pacific (EP) ENSO types (Takahashi et al., 2011), which have distinct anomaly centres in the equatorial Pacific and distinct impacts on remote areas. CP ENSO variability is characterised by low-intensity warming in the central western Pacific, while EP ENSO events are recognised by strong warming in the eastern Pacific. CP and EP types are distinguished by the degree of nonlinearity of SST advection from the western to eastern Pacific, where enhanced nonlinearity leads to the development of intense SST anomalies in the eastern Pacific, associated with EP events. Given this complexity, not all models are able to correctly simulate the two ENSO types due to weak nonlinear SST advection feedback (Takahashi et al., 2011). For this reason, in assessing the response of the CP and EP to the position of the convergence zones, in addition to applying the nonlinear convective feedback criterion for model selection, we also require models to be able to simulate the nonlinear SST advection feedback, known as “Bjerknes feedback”. The model's ability to simulate the Bjerknes feedback is assessed through the nonlinear relationship between the first two principal components of monthly SST anomalies in the tropical Pacific. Models are required to simulate the parameter a, given by the nonlinear coefficient of the fitted quadratic model, at greater than half of the observed value (a_obs = 0.32, Fig. B2).

Figure B3 Relationship between C (a) and E (b) indices and the convection centre index. The solid black line indicates the quadratic fit based on the least-squares method. Banding indicates the 95 % confidence interval based on 1000 bootstrap realisations. R² indicates the coefficient of determination, and a is the nonlinear coefficient of the quadratic regression model. Error estimates for R² and a were calculated as 1 standard deviation of 1000 bootstrap realisations.

As we evaluate ENSO dynamics in climate models, the simulated anomaly centres of CP and EP ENSO may vary across models. As such, indices for CP ENSO (C index) and EP ENSO (E index) were computed through combining the first two empirical orthogonal functions of SST anomalies in tropical Pacific (15° S–15° N, 140° E–80° W). The first two principal component time series are combined to obtain the E ($(\mathrm{PC}\mathrm{1}-\mathrm{PC}\mathrm{2})/\sqrt{\mathrm{2}}$) and C ($(\mathrm{PC}\mathrm{1}+\mathrm{PC}\mathrm{2})/\sqrt{\mathrm{2}}$) indices. The standard deviations of the principal component time series of both C and E indices are used as a proxy of their amplitudes.

The relationship between ENSO amplitude change and the displacement of the convection centres reveals important differences between CP and EP ENSO types. Firstly, CP ENSO activity is likely less sensitive to changes in the mean state, as indicated by a wider shape of the quadratic fit (Fig. B3a; a = −0.79) and a weaker relationship (R² = 0.33) compared to EP ENSO (Fig. B3b; a = −2.31; R² = 0.78). EP ENSO is strongly modulated by the position of the convection centres (R² = 0.78) and exhibits a high nonlinear sensitivity to climate states (a = −2.31​​​​​​​; Fig. B3b). According to this finding, the climate system supports a maximum increase in EP ENSO variance of approximately 55 %, which is achieved with an overall equatorward displacement in the position of the convection centres of ∼ 4.7° (Fig. B3b). It is important to note that the Niño3 index likely captures a combination of CP and EP ENSO variabilities, as its region encompasses both CP and EP anomalies (Takahashi et al., 2011).

Appendix C

The efficiency of the dynamical coupling between the ocean and the atmosphere is measured through the intensity of the wind–thermocline coupling coefficient (Jin et al., 2006). This coefficient measures the sensitivity of the tilt mode of thermocline slope anomalies to wind stress anomalies, which, during El Niño events, results in eastward temperature advection by downwelling equatorial Kelvin waves:

$$\begin{array}{}\text{(C1)}& \langle h{\rangle }_{\mathrm{E}}-\langle h{\rangle }_{\mathrm{W}}={\mathit{\beta }}_h\langle {\mathit{\tau }}_x\rangle ,\end{array}$$

where h indicates the thermocline depth, β_h indicates the wind–thermocline coupling coefficient, and τ_x indicates the zonal wind stress. Subscripts “E” and “W” denote the area average in the eastern (5° S–5° N, 150–90° W) and western (5° S–5° N, 160° E–150° W) equatorial Pacific, respectively. The thermocline depth is computed from the mean temperature profile in each of the boxes indicated above as the weighted average depth, based on depths in which the temperature gradients are greater than 70 % of its maximum (Pontes et al., 2022b). The wind–thermocline coefficient is computed from monthly anomalies, which capture the evolution of the thermocline slope within each single ENSO event.