Global climate mode resonance due to rapidly intensifying El Niño-Southern Oscillation

Domain: earth_science
Published: 2025-10-16
Editorial outcome: major_revision

Abstract
The El Niño-Southern Oscillation (ENSO) influences climate variability globally, encompassing various other modes of variability, and thus represents a key predictable climate signal on seasonal timescales. Yet, its response to greenhouse warming remains uncertain, with models projecting a range of outcomes. Here, we demonstrate that in response to warming, a state-of-the-art high-resolution climate model simulates a rapid transition from a moderate-amplitude irregular regime, as observed in the current climate, to a highly regular oscillation with intensifying amplitude. This behaviour can be attributed to increasing air-sea feedbacks, which approach criticality in the second half of this century, and growing atmospheric noise. As ENSO intensifies in this model, it synchronizes with other prominent climate modes, such as the North Atlantic Oscillation and the Indian Ocean Dipole, thereby imprinting its regular, predictable variability on them. If realized, this global climate mode resonance would have wide-ranging whiplash impacts on regional hydroclimates.

Introduction
Despite the profound influence of the El Niño-Southern Oscillation (ENSO) on the global climate system1,2, its response to greenhouse warming remains uncertain. Climate models exhibit a wide range of possible future ENSO behaviors1,3–7, hindering confidence in regional climate projections. The dynamics of ENSO are governed by a delicate balance of positive and negative feedbacks8,9 that determine both ENSO’s instability10 and periodicity11. The relative strengths of the individual feedbacks are, in turn, determined by both model parametrizations and the structure of the climate mean state12–14. Previous research, using both simple low-order models15 and intermediate complexity models16, has demonstrated how changes in the climate mean state can affect the strength of these feedbacks (such as the zonal advective and thermocline feedbacks) and thereby ENSO characteristics (such as its growth rate, periodicity, and spatial pattern). In addition, recent studies also indicated that the interactions of ENSO with the seasonal cycle17–19 as well as with other more damped empirical modes in the climate system, such as the Indian Ocean Dipole (IOD)20, the Tropical North Atlantic (TNA) mode21, or the North Atlantic Oscillation (NAO)22 can shape the dynamics of both ENSO and these other modes23–25.

Over 20 years ago, an ENSO-resolving coupled general circulation model5,26 exhibited a very intriguing ENSO behavior. It showed a gradual increase in the linear ENSO growth rate in response to greenhouse warming and the crossing of a Hopf bifurcation in the mid-twenty-first century26. This ENSO regime shift towards supercriticality resulted in a rapid intensification in ENSO’s amplitude. Such a drastic transition in qualitative ENSO behavior has not been reported in other complex climate models.

Here we revisit the issue of anthropogenically forced rapid emergence of ENSO supercriticality and its potential repercussions on global climate using an ensemble of state-of-the-art high resolution climate model simulations (AWI-CM3, TCo319 horizonal resolution with ~31 km and 137 vertical layers in the atmosphere and ~4–25 km with 80 vertical layers in the ocean)27, subject to SSP5-8.5 greenhouse gas forcing (see “Methods”). We demonstrate that this model simulates a rapid intensification of ENSO by mid-twenty-first century, a transition to a regular, strongly seasonally-locked oscillation (with periodicities of 2, 3, 4, and 5 years), and an unprecedented resonance with other important modes of climate variability. Using a hierarchy of simplified dynamical ENSO models, with parameters estimated from the complex AWI-CM3 simulations, we study the underlying mechanisms for the qualitative change in ENSO behavior, focusing on coupled air-sea feedbacks and atmospheric noise. While models from the Coupled Model Intercomparison Project Phase 6 (CMIP6) show a wide range of possible future ENSO regularity and amplitude projections, a few models show qualitatively similar behavior to AWI-CM3. If this peculiar ENSO dynamical scenario were to materialize in the future, it could lead to both increased ENSO predictability, due to greater regularity, and, at the same time, to warming-amplified “whiplash impacts”28 on regional climates. These impacts would arise from the compounding effects of (1) more regular ENSO transitions, (2) increased ENSO sea surface temperature (SST) variance, (3) larger ENSO impacts on precipitation and the atmospheric circulation for the same SST anomaly29, and (4) synchronized fluctuations of the other climate modes.

Results
Response of SST and SLP variability to greenhouse warming
As the Earth warms in response to increasing greenhouse gas concentrations, the AWI-CM3 TCo319 model27 simulates a future increase in variability of both SST (Fig. 1e) and sea level pressure (SLP) (Fig. 1f) between the current (Period P1; 2015–2035) and end-of-the-century (Period P2; 2080–2100) climate in many regions of the globe. Further looking at the key regional aspects of these SST and SLP variability changes, we see a striking projected increase in ENSO SST variance as well as in ENSO regularity, i.e., a tendency from an intermittent to a more cyclic ENSO behavior, reminiscent of a Hopf bifurcation30 over this century (Figs. 1a, b, 2a and 3a). Interestingly, the model simulates qualitatively similar changes in the NAO index (Fig. 1c, d), the dominant empirical mode of atmospheric variability over the North Atlantic and Europe, especially in boreal winter (December–January–February: DJF). Both the dominant interannual ENSO timescale, as well as its near-annual combination tones17,18, are evident in the NAO wavelet power spectrum above a white noise background, and their variance increases with time and rising greenhouse gas concentrations (Fig. 1d). Increasing ENSO influences on the NAO are likely to have important implications for atmospheric impacts and climate predictability over Europe27.

Assessment of model fidelity
Next, using sample entropy (SampEn; see “Methods”) of Niño 3.4 SST anomalies as a metric for ENSO regularity together with the Niño 3.4 SST anomaly standard deviation (SD), we show that a 1950 control simulation of AWI-CM3 TCo319 exhibits both ENSO regularity (Fig. 2a) and ENSO amplitude (Fig. 2b) similar to the observations and that the increase of both metrics in the AWI-CM3 TCo319 future projections cannot be explained by internal variability. Furthermore, analysis of the CMIP6 model archive demonstrates that while there is a wide range of possible changes in future ENSO regularity (Fig. 2c) and amplitude (Fig. 2d), 55% of CMIP6 models show an increase in ENSO regularity (Fig. 2c) and 82% an increase in ENSO SST anomaly amplitude (Fig. 2d) by the second half of this century under the SSP5-8.5 scenario. In addition, a few models (i.e., E3SM-1-1, EC-Earth3, and EC-Earth3-Veg) show similar increased ENSO regularity and amplitude as AWI-CM3 TCo319 (Fig. 2c, d and Supplementary Fig. 1). We emphasize that we do not expect strong inter-model agreement of ENSO projections across the CMIP6 archive given that considerable biases of climate mean states and ENSO dynamics result in different ENSO regimes across models8,24.

Physical reasons for the ENSO regime shift
To understand the simulated regime shift in ENSO toward higher variance and increased regularity in AWI-CM3 TCo319, we derive two conceptual ENSO models directly from the simulation output. The first is a version of the recharge oscillator (RO) model13 that includes seasonal cycle modulations9 and encapsulates the fundamental ENSO dynamics in two coupled ordinary differential equations (ODEs) for SST anomalies in the Niño 3.4 region and the zonal-mean equatorial upper-ocean warm water volume. The second is a modified version of the extended recharge oscillator (XRO)24, which also takes into account the seasonally-modulated interactions between ENSO (as described by the RO) and the North Pacific Meridional mode (NPMM), TNA, IOD, and the Indian Ocean Basin (IOB) mode (yielding six coupled ODEs; see “Methods”). We estimate the time-evolving parameters of both conceptual models using multi-linear regression in a moving 21-year window. This allows us to compute the time-dependent linear growth rates and frequencies of ENSO via Floquet eigenanalysis31 and assess changes in the damping rates and coupling strengths of the interacting modes in response to greenhouse warming.

Both models show that ENSO’s linear growth rate increases considerably later in the twenty-first century, particularly during boreal summer (Fig. 3b). This is consistent with the increase in ENSO variance in boreal winter (Fig. 3a). We emphasize that while ENSO, in an annual mean sense, is still in the stable (subcritical) regime, the presence of atmospheric stochastic noise, can lead to a noise-induced Hopf bifurcation that occurs before criticality (i.e., the point at which its growth rate switches from negative to positive) is reached32,33. Moreover, the noise amplitude, calculated as the residual from the XRO model projection, increases substantially in boreal spring to summer (Fig. 3c), together with the growth rate, pushing ENSO into the high-variance and highly-regular regime. This is consistent with increased tropical intraseasonal atmospheric variability in the AWI-CM3 TCo319 simulation (Supplementary Fig. 2, ref. 27), which is well known to energize the ENSO mode2,9,34 and also seen in other climate model simulations35. Comparing the wavelet spectrum of the Niño 3.4 index (Fig. 1b) with the linear ENSO frequency for both the RO and XRO (Supplementary Fig. 3), we see close agreement at the dominant interannual frequencies of ~2, 3, 4, and 5 yr−1.

Next, we use the Bjerknes stability analysis (see “Methods”) to determine which feedback changes are responsible for ENSO’s growing instability. Comparing the difference (Fig. 3f) between period 2 (P2: 2080–2100; Fig. 3e) and period 1 (P1: 2015–2035; Fig. 3d) of the individual ENSO feedbacks that make up ENSO’s net growth rate, we find that the increased growth rate primarily results from a reduced damping from thermocline adjustment (ε) and a modest enhancement of the Bjerknes feedback (R). The changes in R can be explained by enhanced Ekman feedback (EK) and reduced dynamic damping (DD), which are slightly offset by a reduced thermocline feedback (TH) and increased thermodynamic damping (TD). The enhanced Ekman feedback is primarily driven by enhanced stratification and an amplified anomalous vertical current response to SST anomalies (Supplementary Fig. 4d, i). The reduced dynamic damping is linked to the weakening of the climatological trade winds and surface currents. For the zonal advective feedback changes, we see a cancellation effect between the mean state and feedback changes, with reduced dT/dx but strongly enhanced anomalous zonal current response to SST anomalies (Supplementary Fig. 4a, d, h). The reduced thermocline feedback is explained by a weakening of the climatological trade winds and climatological upwelling (Supplementary Fig. 4b, f).

We attribute the reduced damping from thermocline adjustment (ε) to the widening of ENSO wind stress anomalies (Fig. 3g, h and Supplementary Fig. 5a), which leads to increased wind stress curl off the equator, exciting longer Rossby waves. These are less effectively reflected at the western boundary and thus lead to reduced damping from thermocline adjustment (see “ENSO wind stress structure” in “Methods”). Importantly, the projected changes in the ENSO-associated zonal wind stress in the AWI-CM3 TCo319 simulations (Supplementary Fig. 5a) are consistent with the CMIP6 multi-model mean of the response36. The changes of R and ε together lead to the largest growth rate increase in the boreal spring season (Fig. 3f).

The AWI-CM3 model explicitly resolves Tropical Instability Waves (TIWs), which also play a key role in ENSO thermodynamics4,37, mainly as a negative feedback38. With ocean background currents changing in response to the simulated greenhouse warming, the statistics of TIWs will also change, and their contribution to the ENSO heat budget4 will also change. This effect has not been explicitly accounted for in our calculations of the ENSO instability index, but it is captured, among other effects, in the residual term (Fig. 3d–f).

Global climate mode resonance
For climate modes, which are characterized by the variations in SST, such as the TNA mode (Fig. 4b), the IOB (Fig. 4c), the IOD (Fig. 4d), and the NPMM (Fig. 4e), we find a considerable intensification of their amplitude, ranging from 40–75%. All these modes are known to interact with ENSO18,23,24. To further elucidate the apparent emergence of synchronization between ENSO and the other climate modes (i.e., NPMM, TNA, IOD, and IOB) in response to greenhouse warming (Fig. 4), we calculate the power spectrum for each mode for P1 and P2 respectively (Supplementary Fig. 6). We see a clear indication of statistically significant and intensifying spectral peaks for the NPMM, TNA, IOD, and IOB indices at the dominant ENSO frequencies (fE) as well as at the near-annual ENSO combination tone frequencies (1±fE) above a Lorentz/Hasselmann background spectrum39, consistent with theoretical expectations17,18. The emergence of interannual and near-annual ENSO timescales in the spectra of other climate modes provides the first clear evidence that their coupling with ENSO strengthens in response to greenhouse warming in this model. In addition, we see a clear intensification of both the ENSO and ENSO combination mode17 wind stress variability (Supplementary Fig. 7).

To better quantify the dynamical linkage between different climate modes, we calculate the phase difference (see “Methods”) between the Niño 3.4 index (representing ENSO) and the other indices (Fig. 4f–j). A bounded phase difference indicates stronger phase-coupling between the climate modes under consideration (Supplementary Fig. 8). Here, we focus on 1:1 phase synchronization, meaning that the oscillators maintain a (near-) constant phase difference.

Starting around 2060, a strong preference for bounded phase differences emerges between ENSO and each of the other modes (i.e., blue curves in Fig. 4f–j have high values), thus providing further evidence of an intensifying connection between ENSO and the other modes as time progresses. The time-evolving probability distribution of phase differences (red shading in Fig. 4f–j) also documents that the range of phase differences narrows considerably for all modes, supporting the evidence for strengthened phase synchronization with ENSO.

Under present-day conditions, the influence of ENSO on the NAO is detectable but generally weak and indirect. Consequently, current seasonal forecasts for the European climate do not treat ENSO as a primary driver of NAO-related variability or impacts. According to our simulations, however, this could change in the future. The amplitude of the NAO increases by about 20% during the simulation (Fig. 4a), but the overall seasonal amplitude modulation with a peak in February remains relatively stable. The ±π preferred phase difference to ENSO for the NAO is consistent with the well-known out-of-phase relationship between ENSO and the NAO. As greenhouse gas concentrations increase, the preferred phase difference between ENSO and the NAO becomes more robust (Fig. 4f), with dominant peaks between −π and −π/2, consistent with an increasingly negative correlation in boreal winter. Overall, we expect these dynamics to potentially influence the seasonal predictability of the NAO, as illustrated in Supplementary Fig. 9, which shows the projected extended persistence of both ENSO and NAO. It is also illustrated in Fig. 1, which shows that after 2060, 2-year-long La Niña events typically create two subsequent strong NAO events.

Whether the intensification of the NAO linkage to ENSO emerges, simply because the ENSO SST amplitude increases during the simulation, or whether large-scale reorganizations of the extratropical atmospheric circulation facilitate the coupling between the Pacific-North American (PNA) pattern and the NAO in wintertime (DJF) is quantified by a regression analysis for both periods (Fig. 5). For P1, the regression of precipitation and 500 hPa geopotential height (Z500) anomalies onto Niño 3.4 SST anomalies and Niño 4 precipitation shows no substantial impact in Europe. During P2, for the regressions with both the Niño 3.4 SST index in degrees Kelvin (Fig. 5e) and with the normalized Niño 3.4 SST index (Fig. 5d), we see increases in the regression coefficient amplitudes and a clear negative NAO pattern emerging. The strongest signal can be found for the normalized Niño 3.4 index, which indicates that the amplification of the SST signal itself plays a vital role in strengthening the teleconnection. However, even a 1 K SST change in the Niño 3.4 region, the NAO pattern emerges much more strongly in P2 than in P1, indicating a higher sensitivity of the atmospheric response to the same amplitude of warming in the eastern equatorial Pacific. As a result of the increased synchronization between ENSO and the NAO pattern in boreal winter, we also find an enhanced precipitation response over western Europe27 with considerably wetter conditions occurring over the Iberian Peninsula during El Niño events (Fig. 5). Recent studies have suggested that the ENSO–NAO teleconnection may strengthen in response to greenhouse warming, either due to amplified ENSO forcing or increased extratropical sensitivity. Reference 40 showed that CMIP5 models simulate both enhanced tropical Pacific precipitation variability and a more robust ENSO–NAO correlation in the future, with evidence for stronger stratospheric and tropospheric pathways. Reference 41, using a pacemaker approach with fixed ENSO SST anomalies, attributed the strengthened NAO response to increased extratropical sensitivity linked to a more zonally extended Pacific jet. Our study confirms and extends these results by showing that in a high-resolution coupled model, both mechanisms occur simultaneously: ENSO amplitude increases due to stronger air–sea feedbacks and enhanced atmospheric noise, while normalized regressions reveal that the extratropical circulation becomes more sensitive to ENSO forcing. In addition, we identify a regime shift toward a highly regular, seasonally phase-locked ENSO that entrains other modes like the NAO, suggesting a possible global synchronization of interannual variability under climate change.

The increase in ENSO variance and regularity directly leads to increased variance and regularity of the other SST climate modes. In addition, we find that both a reduction of their individual damping rates (Supplementary Fig. 10a–d) and increased ENSO teleconnection strength (Supplementary Fig. 10e–l) also contribute to the global climate mode resonance. The enhanced variance of the NPMM is primarily driven by reduced damping from February to July (Supplementary Fig. 10a), likely associated with SST warming in the subtropical northeastern Pacific and enhanced wind–evaporation–SST (WES) feedback42. The reduction in damping rate and the strengthening of ENSO teleconnection are essential for increasing IOD and TNA variance (Supplementary Fig. 10c, d, g, h). Interestingly, the IOB exhibits an enhanced early variance peak in February–March and a secondary peak in May–June by the end of the twenty-first century (Fig. 4c). This change is related to reduced damping in addition to enhanced ENSO teleconnection from November to April (Supplementary Fig. 10b, j).

Methods
Earth system model experiments
This study employs the Alfred Wegener Institute Climate Model (AWI-CM3), which couples the OpenIFS (Open Integrated Forecasting System) atmosphere model with the FESOM2 (Finite-volume Sea Ice-Ocean Model) ocean model. The atmospheric component operates at a horizontal resolution of approximately 31 km (TCo319) with 137 vertical pressure levels. It contains the WAM (Wave Model) surface gravity wave model and the hydrology model H-TESSEL (Hydrology in the Tiled ECMWF Scheme for Surface Exchange over Land). The ocean model features a variable horizontal resolution ranging from 5 to 27 km, depending among others on latitude, and includes 80 vertical layers. The reader is referred to a more detailed description of AWI-CM3 and the simulations (except the additional ensemble members)27.

A full transient simulation covering the period from 1950 to 2100, branched off from a 100-year-long spin-up simulation, was conducted using historical forcing from 1950 to 2014 CE, followed by the high-emissions shared socio-economic pathway (SSP) 5-8.5 scenario thereafter. To assess the robustness of our results, three additional ensemble simulations spanning the period from 2055 to 2100 were conducted. In line with the micro-perturbation initialization strategy (e.g., ref. 35), a small perturbation was added to the wave model restart files.

A 150-year TCo319 control simulation with fixed 1950 forcing (CTL1950) was also analyzed to characterize internal climate variability and to serve as a baseline comparison. The CTL1950 simulation reproduces the observed SST and precipitation mean state and SST variability reasonably well (Supplementary Fig. 11). The CTL1950 simulation also reproduces the observed characteristics of ENSO, including its seasonal synchronization, spectral characteristics, and teleconnections (Supplementary Fig. 12).

Observational and CMIP6 data
We use four observational SST reconstructions/reanalysis: the Hadley Centre Sea Ice and Sea Surface Temperature dataset v. 1.1 (HadISST45), the Extended Reconstructed Sea Surface Temperature v. 5 (ERSSTv546), the Centennial in situ Observation-Based Estimates of Sea Surface Temperature v. 2 (COBE247) for 1871–2024, and the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis 5 (ERA548) for 1940–2024. In addition, ERA5 monthly precipitation and sea level pressure (SLP) fields are used for 1940–2024.

We also analyze simulations from 49 CMIP6 models, including historical runs and SSP5-8.5 projections, providing monthly SST and SLP data (Supplementary Table 1). The simulations are forced by historical anthropogenic and natural forcings up to 2014, followed by future greenhouse-gas forcing under the SSP5-8.5 scenario through 2100, covering the period 1850–2100. All model output was re-gridded to a common 1°×1° horizontal resolution using bilinear interpolation.

Definitions of climate variability modes
All monthly fields are first quadratically detrended, and then anomalies are computed by subtracting the 21-year running monthly climatology. We then calculate the climate mode indices with the detrended anomalies. The ENSO SST index is represented as SST anomalies averaged over the cold tongue region (180°–90°W, 6°S–6°N). The Niño 3.4 SST index is defined as SST anomalies averaged over the Niño 3.4 region (170°–120°W, 5°S–5°N). The Niño 4 precipitation index is defined as precipitation anomalies averaged over the Niño 4 region (160°E–150°W, 5°S–5°N). The North Pacific Meridional Mode (NPMM) index is defined as SST anomalies averaged over 160°–120°W, 10°–25°N49. The NPMM SST anomaly index reproduces the original maximum covariance analysis (MCA)-based NPMM index with a correlation of ~0.9549. The Tropical North Atlantic (TNA) index is defined as SST anomalies averaged over 55°–15°W, 5°–25°N50. The Indian Ocean Basin (IOB) mode index is defined as SST anomalies averaged over 40°–100°E, 20°S–20°N51. The Indian Ocean Dipole (IOD) mode index is defined as SST anomalies averaged over 50°–70°E, 10°S–10°N minus those averaged over 90°–110°E, 10°S–0°N20. The North Atlantic Oscillation (NAO) index is defined as the principal-component (PC) time series of the leading EOF of sea level pressure (SLP) anomalies over the North Atlantic sector (80°W–40°E, 20°–80°N)52.

Spectral analysis: wavelets
Wavelet analysis was used to explore the different climate modes’ time-frequency characteristics and identify dominant periodicities53. A continuous wavelet transform was performed using the Morlet wavelet, which offers an optimal balance between frequency and time localization, making it ideal for capturing oscillatory behavior. The wavelet power spectrum was calculated for each ensemble member, followed by calculating the ensemble-mean of the spectra (Fig. 1b, d) to highlight changes in dominant frequencies over time. Statistical significance was assessed using Monte Carlo simulations: Niño 3.4 SST was tested against an AR(1) process with a lag(-1 month) autocorrelation coefficient of 0.9, while the NAO index was tested against a white-noise process.

Spectral analysis: power spectral density
The power spectral density (PSD) of the various climate mode indices (the four ensemble members for P2 were concatenated, leading to 84 years of monthly data for P2 but only 21 years of monthly data for P1) was calculated using the Multi-Taper Method (MTM) with 3 tapers and nfft=256 for P1 and nfft=1024 for P254. A Lorentzian background spectrum null hypothesis was chosen for the SST modes. The significance level of spectral peaks was therefore determined from the respective PSD percentile at each frequency of PSDs calculated with the same MTM method for 10,000 discrete AR(1) processes with (i) the same lag(-1 month) autocorrelation, (ii) the same variance, and (iii) the same data length as the respective index that is tested against. A white noise null hypothesis was chosen for the NAO atmospheric mode. Here, the significance level of spectral peaks was determined from the respective PSD percentile at each frequency of PSDs calculated with the same MTM method, generated from 10,000 realizations of white noise time series with (i) the same variance and (ii) the same data length as the index that is tested against.

ENSO regularity calculation using sample entropy
We assessed ENSO regularity by calculating the sample entropy (SampEn) of the Niño 3.4 SST anomaly index, a nonlinear metric that measures how predictable a time series is55. SampEn is an alternative to approximate entropy (ApEn56) for quantifying the randomness of a time series; unlike ApEn, it excludes self-matches, thereby reducing bias and improving consistency across different data segments57. Compared to an ENSO regularity metric that is defined by the sharpness of the ENSO spectral peak58, SampEn is less sensitive to record length, spectral windowing, and other frequency-domain constraints. Lower SampEn values indicate more regular (periodic) behavior, whereas higher values indicate more irregular behavior. Given a time series {x1, x2, ..., xN}, SampEn estimates the negative natural logarithm of the conditional probability that two sequences of length m that match within tolerance r will also match for m+1 points, while excluding self-matches to reduce bias and improve consistency. We used an embedding dimension m=2 (the dimensionality of the ENSO RO model) and a tolerance r=0.2σ, where σ is the standard deviation of the time series. The SampEn metric is defined as:

SampEn(m,r,N) = -ln(Am(r)/Bm(r)) (1)

where Am(r) the probability that two sequences of length m+1 are similar (matches), Bm(r) is the probability that two sequences of length m are similar (candidates). Explicitly,

Am(r) = 1/(N-m-1) * sum_{i=1}^{N-m} (number of times that d(u^{m+1}(j), u^{m+1}(i)) < r) (2)

Bm(r) = 1/(N-m) * sum_{i=1}^{N-m} (number of times that d(u^m(j), u^m(i)) < r) (3)

where u^m(i) is the subsequence {xi, xi+1, ..., xi+m-1}, d(u^m(j), u^m(i)) = max_{k=1,...,m} |xi+k-1 - xj+k-1| is the maximum norm distance. Since the number of matches is always less than or equal to the number of possible vectors, the ratio Am(r)/Bm(r) is a conditional probability less than unity, ensuring SampEn is positive. This method provides a single diagnostic value characterizing the regularity of the ENSO time series, facilitating consistent comparison across observations and simulations under different climate scenarios. The observed SampEn effectively captures the decadal variations in ENSO forecast skill, including the modest decrease in predictability after the 2000s (black curve in Fig. 2a, refs. 59–61).

Recharge Oscillator (RO) model formulation and fit
The Recharge-Oscillator (RO) model is a widely recognized framework for diagnosing and understanding ENSO dynamics in observations and climate models9,13. The RO model captures the oscillatory behavior of El Niño and La Niña events through two coupled equations describing the evolution of ENSO SST anomalies (T) and equatorial Pacific zonal-mean thermocline depth anomalies (h). In its linear form, the model is expressed as:

d/dt [T_ENSO, h_ENSO] = L_ENSO * [T_ENSO, h_ENSO] + [ξ_ENSO, ξ_h] (4)

In Eq. (4), L_ENSO = [R, F1; -F2, -ε], where R represents the SST anomaly growth rate, collectively describing the Bjerknes feedback, ε denotes ocean damping rate related to the energy leakage at the western boundary and mixing, F1 represents the effectiveness of the heat content discharge–recharge in controlling SST anomalies, and F2 represents the discharge–recharge efficiency; ξ_ENSO and ξ_h are the residual terms encompassing non-resolved nonlinear processes and stochastic noise. The complex eigenvalue of the RO model operator L_ENSO is the Bjerknes-Wyrtki-Jin (BWJ) index for the ENSO growth rate and periodicity and can be written as:

BWJ = (R-ε)/2 + i * sqrt(F1*F2 - (R+ε)^2/4) (5)

where the real part of the index is referred to as the ENSO Bjerknes stability (BJ) index in ref. 10, and the imaginary part as the reciprocal of the Wyrtki periodicity (WF) index in ref. 11.

Due to the strong seasonal dependence of ENSO, we explicitly incorporate seasonality by estimating the parameters as:

L = L0 + Lc*cos(ωt) + Ls*sin(ωt) (6)

where ω = 2π/(12 months), and the subscripts 0 and 1 indicate the mean and annual cycle components, respectively. The parameters of the operators are estimated by multivariate linear regression. We applied the RO model fitting to each ensemble member of the AWI-CM3 simulation using a 21-year moving window.

Extended Recharge Oscillator (XRO) model formulation and fit
To understand the climate mode interactions, we use a modified version of a recently developed extended Recharge Oscillator (XRO), which allows for two-way interactions between ENSO and the other modes24, which consists of a recharge oscillator model for ENSO coupled to seasonally-modulated stochastic-deterministic models for the other climate modes25:

d/dt [T_ENSO, T_NPMM, T_IOB, T_IOD, T_TNA, h] = L * [T_ENSO, T_NPMM, T_IOB, T_IOD, T_TNA, h] + [ξ_ENSO, ξ_NPMM, ξ_IOB, ξ_IOD, ξ_TNA, ξ_h] (7)

The linear dynamical operator L contains four submatrices, organized as follows:

L = [L_ENSO, C1; C2, L_M] (8)

where the linear operator submatrix L_ENSO describes the ENSO internal recharge-discharge dynamics, L_M represents the internal processes and interactions among the other climate modes; C are coupling submatrices, with C2 describing the impact of ENSO on other climate modes and C1 describing the feedback of other modes on ENSO. Similar to the RO model fit, we explicitly incorporate seasonality by estimating the XRO parameters with Eq. (7). The noise parameters are determined from the residuals of the XRO fit. There are a total of 12 noise parameters, i.e., a noise amplitude and decorrelation timescale for each of the 6 state variables in the system. The noise amplitudes σ_ξ are estimated from the standard deviations of the residuals of the XRO fit. The decorrelation timescales τ_ξ are estimated as τ_ξ = Δt / ln(1 - δ_ξ), where δ_ξ is the autocorrelation of the residuals of the XRO fit. The range of observed noise time scales τ_ξ are between 0.25~0.70 months.

Floquet theory is used to determine the stability of the coupled system with an annual cycle basic state. We determine the ENSO growth rate and periodicity using the RO model operator L_ENSO or XRO model operator L via Floquet exponent analysis31. The Floquet exponents are determined numerically by integrating dQ/dt = L*Q over 1 year with a time step of 3.65 days from an initial condition Q(0)=I to form a monodromy matrix M = Q(T), T = 1 year. The Floquet exponents (σ_j) are calculated from the least damped complex eigenvalues α_j of M by using σ_j = ln(α_j)/T.

Atmospheric noise definitions
To quantify changes in atmospheric noise amplitude, we define atmospheric noise as surface zonal wind variability that is independent of SST-related variability, following ref. 62. Specifically, an EOF analysis was applied to monthly tropical SST anomalies to derive the leading modes of variability and their associated PCs. Multivariate linear regression was then performed, regressing the monthly surface zonal wind stress anomalies onto the first 15 leading PCs, which together account for approximately 75% of the variance in tropical SST anomalies. The time series of atmospheric noise was then defined as the residual wind stress, obtained from removing the SST-forced signals. Finally, the amplitude of noise is defined as the standard deviation of the residual wind stress.

Heat budget and Bjerknes stability analysis
To assess possible physical processes that control the changes of ENSO growth rate, we calculate the ocean mixed layer heat budget using a partial flux form63,64:

∂T'/∂t = Q'_0 - ∂(u*T')/∂x - ∂(v*T')/∂y - ∂(w*T')/∂z - u'*∂T_bar/∂x - v'*∂T_bar/∂y - w'*∂T_bar/∂z - u'*∂T'/∂x - v'*∂T'/∂y - w'*∂T'/∂z + Q_Res (9)

where the overbars denote the climatological monthly-mean seasonal cycle, and the primes denote anomalies; T is the ocean temperature, u, v, and w the ocean zonal, meridional, and vertical current velocities, Q_Res the unresolved surface heat flux and the effects of the unresolved surface heat flux. The net surface heat flux Q'_0 is calculated as:

Q'_0 = Q'_net / (ρ*Cp*H) (10)

where H is the ocean mixed layer depth (50 m), ρ is the seawater density (1025 kg m−3), Cp is the specific heat capacity of seawater (3994 J kg−1 K−1), and Q'_net represents the net surface heat flux, consisting of four components:

Q'_net = Q'_SW + Q'_LW + Q'_LH + Q'_SH (11)

representing shortwave radiation, longwave radiation, latent heat flux, and sensible heat flux anomalies, respectively. The shortwave radiation transmitted through the bottom of the mixed layer (Q'_bot) is parameterized following ref. 65 as:

Q'_bot = Q'_SW * (0.58*e^(-H/0.35) + 0.42*e^(-H/23)) (12)

Taking the volume-average by integrating both hand sides of Eq. (9) from the ocean mixed layer depth to the ocean surface and spatially over the eastern equatorial Pacific box where the SST variability of ENSO is dominant, the heat budget equation can be written as follows:

∂<T'>/∂t = <Q'_0> - <∂(u*T')/∂x> - <∂(v*T')/∂y> - <∂(w*T')/∂z> - <u'*∂T_bar/∂x> - <v'*∂T_bar/∂y> - <w'*∂T_bar/∂z> - <u'*∂T'/∂x> - <v'*∂T'/∂y> - <w'*∂T'/∂z> + <Q_Res> (13)

where the angle brackets denote volume/area-average variables over the eastern Pacific box (180°–90°W, 6°S–6°N), T'_bot and w_bar the anomalous temperature and climatological vertical velocity at the bottom interface of the ocean mixed layer, respectively. We group the various feedbacks on the right-hand side of Eq. (9) as thermal damping by the net surface heat flux (TD), dynamic damping by mean horizontal currents (DD), thermocline feedback (TH), zonal advective feedback (ZA), meridional advective feedback and vertical advective feedback (combined as Ekman feedback, EK), nonlinear dynamic heating (NDH), and sub-grid scale contributions (SG), respectively. Our definitions of DD and TH represent the advection of anomalous temperature by mean currents via lateral and bottom interfaces, respectively.

Next, we apply the following linear relationship between various oceanic processes with T_e and h9 and replacing <T'> with T_e for simplicity:

TH = R_TH * T_e + F_1TH * h (14a)
ZA = R_ZA * T_e + F_1ZA * h (14b)
EK = R_EK * T_e + F_1EK * h (14c)
DD = R_DD * T_e + F_1DD * h (14d)
TD = R_TD * T_e (14e)
NDH = R_NDH * T_e + F_1NDH * h (14f)
Res = R_Res * T_e + F_1Res * h (14g)

Substituting these into the ENSO mixed-layer temperature equation, the equation reduces to:

∂T_e/∂t = (R_TH + R_ZA + R_EK + R_DD + R_TD + R_NDH + R_Res) * T_e + (F_1TH + F_1ZA + F_1EK + F_1DD + F_1NDH + F_1Res) * h (15)

From this, the total SST growth rate R and phase-transition rate F1 are expressed as:

R = R_TH + R_ZA + R_EK + R_DD + R_TD + R_NDH + R_Res
F1 = F_1TH + F_1ZA + F_1EK + F_1DD + F_1NDH + F_1Res

To estimate these parameters, we account for seasonality explicitly outlined in Eq. (6) using multivariate linear regression. Ensemble mean parameters are calculated by performing regression on concatenated results from all ensemble members. These derived parameters capture the contributions of each process and provide a comprehensive understanding of ENSO feedbacks and their contribution to the ENSO growth rate. The R_NDH and R_Res terms are combined as R_AllRes, which consists of tropical instability waves and other unresolved processes.

We have also tested our analysis using different ocean mixed layer depths (30 m, 40 m, and 50 m) and alternative definitions of the eastern Pacific box (e.g., the Niño 3.4 region). The results remain quantitatively similar and do not affect our conclusions.

ENSO wind stress structure
The ENSO wind stress structure affects both the growth rate and periodicity of ENSO9,36. Here, we determine the ENSO wind stress pattern using linear regression of monthly zonal wind stress anomalies (τ'_x) onto the monthly ENSO SST anomaly index. Following ref. 36, the meridional width of τ'_x is defined as the latitudinal distance between the meridional boundaries, the nearest local minima of τ'_x on either side of the equator, based on the meridional profile averaged over 160°E–150°W (Fig. 2g). The zonal length is defined as the zonal distance between the zonal boundaries of τ'_x, determined from the zonal profile averaged over 5°S to 5°N across the tropical Pacific. The eastern boundary is identified as the zero-crossing point in the zonal profile; the western boundary is defined as the local minimum within the region from 140°E to 190°E; and zonal centroid location as the midpoint between zonal boundaries (Supplementary Fig. 5b).

Under the assumption that the zonal wind stress anomaly τ'_x has a Gaussian meridional structure, the RO h equation can be reformulated to explicitly incorporate the wind stress structure36:

ε = ε_m + α^2 / x_c (16)
F2 / α / x_c (17)

where ε_m is a linear mechanical damping coefficient (2.5 yr)−1, the parameter α characterizes the meridional structure of τ'_x, such that the meridional width of τ'_x is proportional to 1/α, and x_c denotes the nondimensional location of the center of wind stress anomalies along the equator (Supplementary Fig. 5b). From the meridional widening of τ'_x (with x_c ~ 0.35), a reduction in α, thereby weakening the damping from thermocline adjustment (ε, Fig. 2h) and recharge/discharge efficiency (F2, Fig. 2i). These changes enhance the growth rate and increase the periodicity of ENSO following the RO eigenvalue Eq. (5).

Phase synchronization
Phase synchronization characterizes how interactions between two different processes can lead to synergistic dynamical phase behavior (e.g., the temporal alignment of respective local maxima or minima). A special case of phase-synchronization is frequency-locking31, where two systems with normally independent oscillatory behavior develop a joint frequency and locked dynamics. In general, for phase-synchronized systems, the difference between the phases of the interacting processes, or multiples thereof, is statistically distinguishable from a uniform white noise process. In particular, phase synchronization can emerge in nonlinear dynamical systems, which include product terms between two signals.

Practically, phase relationships between two processes, which are represented by their time series x(t) and y(t) can be calculated using their complex analytical signals66,67 z1(t) = A1(t)e^(iφ1(t)) = x(t) + ix_H(t), and z2(t) = A2(t)e^(iφ2(t)) = y(t) + iy_H(t), where x_H(t), y_H(t), represent the Hilbert transforms of x(t), y(t) and φ1(t), φ2(t) are the respective phases of the signal. The only requirement for the two signals x(t), y(t) is that they satisfy Bedrosian’s Product Theorem68. This can be achieved by appropriate band-pass filtering. To quantify m:n phase synchronization (m, n are integers), we calculate δφ_{n,m} = nφ1(t) - mφ2(t). For two un-synchronized (decoherent) processes, we would expect (1) the time evolution of the unwrapped phases δφ_{n,m} to essentially increase, without any long periods of phase “friction”, which can be characterized as consistent plateaus of near constant δφ_{n,m}; (2) the probability distribution of δφ_{n,m} (without unwrapping procedure) to be indistinguishable from a uniform distribution between 0–360 degrees66,69. Respectively, for m:n synchronized signals, (1) we would expect to see long periods of little growth in unwrapped phase differences; (2) the frequency histogram of δφ_{n,m} (without unwrapping procedure) across all angles (0–360°) should deviate statistically (e.g., Kolmogorov-Smirnov test) from a uniform distribution.

We use this approach to calculate the 1:1 difference of unwrapped phases between time series characterizing different climate modes simulated by the AWI-CM3 model. All time series were band pass filtered between 2.5 months and 10 years using a Butterworth filter prior to the analysis to focus on annual to interannual variability and to satisfy Bedrosian’s theorem. We use the standard deviation of the phase difference PDF, σ[PDF(δΦ_{1,1})], as a measure of the synchronization strength. If there is no synchronization, the phase differences are uniformly distributed, yielding a flat PDF and σ[PDF(δΦ_{1,1})] ≈ 0 (Supplementary Fig. 8a, b). In contrast, if there is strong synchronization, the phase differences concentrate around a preferred value (−π/2 for the examples shown in Supplementary Fig. 8c, d), producing a sharply peaked PDF, resulting in a larger σ[PDF(δΦ_{1,1})] (Supplementary Fig. 8c, d).

Statistical significance test
The Fisher z-transformation was used to test statistical significance of the correlation differences as follows:

Z = 0.5 * ln((1+r1)/(1-r1)) - ln((1+r2)/(1-r2)) / sqrt(1/(n1-3) + 1/(n2-3)) (18)

where r1 and r2 are the correlation coefficients, n1 and n2 are the sample sizes of the first and second group samples. The absolute value |Z| is then compared against a critical value from the t-distribution for a two-tailed test. We rejected the null hypothesis that the two correlations are not significantly different at a 90% confidence level if |Z| exceeds the critical value. The Student’s t test was used to test statistical significance of the regression coefficient.

Given small sample sizes, we apply a bootstrap method to assess whether changes in ENSO regularity are statistically significant. Specifically, the 21-year moving-window values of ENSO regularity are randomly resampled with replacement to generate 10,000 synthetic realizations. From these realizations, we compute the standard deviation of the resulting distribution of mean values.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of “Global climate mode resonance due to rapidly intensifying El Nino-Southern Oscillation” submitted to Nature Communications. This paper examines how greenhouse-gas–driven warming leads to highly regular, large-amplitude ENSO cycles by mid-21st century, in response to strengthened air–sea feedbacks and weakened thermocline damping. The paper also shows that this intensification of ENSO phase-locks and amplifies variability in the North Atlantic, Indian, and Pacific climate modes of variability (e.g. NAO, TNA, IOD), suggesting both enhanced predictability and heightened global hydroclimate risks under alternating El Niño/La Niña events.\nThis is a very good paper, with great implications for global climate risk forecasting on sub-seasonal to seasonal scales. I also really like the methodology, combining high-resolution coupled climate models (AWI-CM3) with reduced-order dynamical diagnostics. I am only concerned about the use of a single coupled climate models, and the absence of a control simulations with change in greenhouse gas emissions. Can we find the same patterns in lower resolution CMIP6 ScenarioMIP and/or HighResMIP models? Is the intensification of ENSO variability exceeding the maximum range of variability that could be simulated in piControl experiments? I was also wondering how the method used to analyse the amplitude and phase synchronization would compare with the methodology developed by Jajcay et al. (2018; https://doi.org/10.1038/s41612-018-0043-7).\nI would therefore recommend major revisions.", "reviewer_2": "The authors present a well-written, timely study on future ENSO intensification and global synchronization under greenhouse warming, using a high-resolution climate model. The abstract clearly communicates key findings, especially the shift to a highly regular, high-amplitude ENSO regime and its impacts on other climate modes. This addresses an important unresolved question, and the proposed mechanism linking air-sea feedbacks, atmospheric noise, and teleconnection synchronization is compelling. I recommend the manuscript for publication, provided the authors address some important points.\n• A major limitation is the reliance on a single climate model (AWI-CM3) to infer future changes in ENSO amplitude, teleconnection strength, and global hydroclimate impacts. Given well-known inter-model diferences in ENSO characteristics, the lack of multi-model comparison raises concerns about the robustness and generalizability of results.\n• Additionally, the study lacks validation of modeled trends against observational or reanalysis data. Without benchmarking AWI-CM3’s current ENSO behavior and atmospheric responses, confidence in future projections is reduced. Incorporating multi-model analyses and observational validation would greatly enhance credibility and applicability.\n• The suggestion that future ENSO will become more regular, intense, and predictable—synchronizing other major climate modes like the North Atlantic Oscillation and Indian Ocean Dipole—is speculative. ENSO predictability as well as teleconnections exhibit substantial variability influenced by internal dynamics and background state changes. Assuming a more regular ENSO leads to improved predictability or consistent control over other modes oversimplifies the complex climate system. Without supporting observational or multi-model evidence, this conclusion should be framed more cautiously, clearly distinguishing model results from real-world implications. The authors may consider referencing relevant literature on ENSO predictability for discussion or introduction.\nhttps://www.nature.com/articles/s41612-024-00845-5", "reviewer_3": "Stuecker et al. investigates climate change experiments in a high resolution coupled GCM under a high-CO2-emissions scenario and finds a strong increase in the amplitude and periodicity in ENSO. Consequently, other climate modes phase lock onto the strong and periodic ENSO signal and increase in amplitude. They use a hierarchy of idealized models to diagnose the changes in ENSO amplitude and periodicity. The societal impacts of such a transition are manifest, but it remains unclear whether this feature is a peculiarity of the chosen model or whether it reflects reality.\nThe paper is clearly written, the results are well-presented, and the importance of the findings are clear. I consider the paper suitable for publication after minor revisions. However, the paper would be improved if the authors provide some additional commentary and analysis on why this model exhibits this unique ENSO behavior, but not any other models. I offer some small suggestions to improve the manuscript.\nMajor Comments\nThe results of this paper are very interesting, but one important point is not addressed: why does the AWI-CM3 model simulate a rapid switch to a highly regular oscillation with high amplitude, while other climate models do not? Since this model is an outlier, the impact of the paper will be increased if an argument can be made that this model’s ENSO response is in some respects more realistic than other models.\nA few ideas are hinted at throughout the paper, such as intensifying tropical intraseasonal variability (Ln. 111), or the relaxation of the mean SST gradient and trade winds (Ln.130). But these changes are seen in other climate models as well (especially the weakening of the east-west SST gradient), so what makes AWI-CM3 an outlier? Is its high horizontal and vertical resolution important? If so, do we see this increase in other HiResMIP models?\nIf the authors can discuss possible reasons conventional CMIP6 models do not simulate the response, it would be great. Furthermore, if they can speculate on why this particular model might better reflect reality compared to other CMIP models, that would further improve the impact. But overall it is still a very nice result that is worth being published.\nMinor Comments\nLn. 52-53. It could be interesting to show Fig 1B for some selected CMIP6 models for comparison, to demonstrate that this behavior of this model is unique.\nLn. 123 It’s not obvious to me that ZA is very significant, it looks like it averages about zero over the annual cycle.\nLn. 174. I’m confused why, if the range of phase differences narrows, the standard deviation of the PDFs trends upwards (blue curve)? Intuitively, I would expect the standard deviation to become smaller.\nFig.4: The features on the figure are very small. Perhaps three rows and two columns will allow each sub-panel to be larger."}, "author_response": "Manuscript NCOMMS-25-44516\nPoint-by-point response to the reviewers\nWe thank the three reviewers for their valuable comments and suggestions, which have\nimproved the manuscript significantly. We have included a point-by-point response to the\ncomments below (Reviewer comments are in black font, our responses in blue font, changes\nin the revised text in Times New Roman italics).\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nReview of “Global climate mode resonance due to rapidly intensifying El Nino-Southern\nOscillation” submitted to Nature Communications. This paper examines how greenhouse-\ngas–driven warming leads to highly regular, large-amplitude ENSO cycles by mid-21st\ncentury, in response to strengthened air–sea feedbacks and weakened thermocline\ndamping. The paper also shows that this intensification of ENSO phase-locks and amplifies\nvariability in the North Atlantic, Indian, and Pacific climate modes of variability (e.g. NAO,\nTNA, IOD), suggesting both enhanced predictability and heightened global hydroclimate\nrisks under alternating El Niño/La Niña events.\nThis is a very good paper, with great implications for global climate risk forecasting on sub-\nseasonal to seasonal scales. I also really like the methodology, combining high-resolution\ncoupled climate models (AWI-CM3) with reduced-order dynamical diagnostics.\nThank you for the assessment and valuable feedback. We address your specific comments\nbelow.\nI am only concerned about the use of a single coupled climate models, and the absence of\na control simulations with change in greenhouse gas emissions. Can we find the same\npatterns in lower resolution CMIP6 ScenarioMIP and/or HighResMIP models? Is the\nintensification of ENSO variability exceeding the maximum range of variability that could be\nsimulated in piControl experiments?\nFollowing your suggestion, we have added comparisons of the AWI-CM3 projections with a\n150-yr control simulation with the same model as well as with observational and reanalysis\ndata (new Fig. 2a,b in the revised manuscript). The AWI-CM3 model simulates ENSO\namplitude and regularity within the observational range. Both ENSO amplitude and\nregularity show a significant increase when comparing the AWI-CM3 SSP5-8.5 projections\nwith the internal variability of these metrics in the AWI-CM3 control simulation.\nWe also analysed ENSO regularity and amplitude in all 49 available CMIP6 models (Fig.\n2c,d). As expected, the projected changes in ENSO regularity and amplitude exhibit a wide\nspread; however, AWI-CM3 is not an outlier. Several CMIP6 models produce similar results\n(Fig. 2c,d; Supplementary Fig. 1). The ENSO changes across all these models are shown in\nFig. R1 below for your reference.\nFig. 2| ENSO regularity and amplitude changes in observations, AWI-CM3, and CMIP6\nmodels. a,b Moving 21-year changes of sample entropy (SampEn) (a) and standard deviation (SD)\n(b) of the Niño3.4 SST anomaly index in observations (black), the 150-year AWI-CM3 1950\ncontrol simulation (CTL1950, blue), and the AWI-CM3 SSP5-8.5 scenario simulation (red). Solid\nlines and shading indicate the multi-product/multi-member mean and the 1 SD spread of 10,000\ninter-realizations using a bootstrap method, respectively. The error bars indicate the mean and 1\nSD of values across all overlapping 21-year windows for the observations (1950-2024),\nCTL1950 (150-yr), and SSP5-8.5 (2050-2100). c,d Violin plots showing the probability density distribution\nof ENSO regularity change (c) and amplitude change (d) across 49 CMIP6 models. The change is\ncalculated as the difference between 2050–2100 and 1900–2000. Dots represent individual model\nresults, with E3SM-1-1, EC-Earth3, and EC-Earth3-Veg highlighted in distinct colours; the red\nstar marks the AWI-CM3 result (2050-2100 minus CTL1950). In a and c, the y-axis is inverted to\nemphasize decreasing SampEn indicating increasing regularity.\nSupplementary Fig. 1| Increased ENSO regularity in selected CMIP6 models. a Time\nevolution of the ENSO Niño3.4 SSTA index [K] for EC-Earth3 (3-month running-mean applied).\nb Wavelets of the normalized Niño3.4 SSTA index for EC-Earth3. Black contours enclose regions\nabove 90% confidence level significance tested against a Lorentzian spectrum. c,d and e,f Same\nas a,b but for EC-Earth3-Veg and E3SM-1-1, respectively.\nFig. R1| ENSO time series in CMIP6 models. Time evolution of ENSO Niño3.4 SST anomaly\nindex [K] (3-month running-mean applied) for ERSSTv5 (1900-2024), the 150-year AWI-CM3\nCTL1950 simulation, as well as the historical+SSP5-8.5 simulations from 49 CMIP6 models\n(1900-2100).\nCorresponding changes to the revised manuscript text (L92–106, a new section titled\n“Assessment of model fidelity.”):\nAssessment of model fidelity\nNext, using sample entropy (SampEn; see Methods) of Niño3.4 SST anomalies as a metric for\nENSO regularity together with the Niño3.4 SST anomaly standard deviation (SD), we show that a\n1950 control simulation of AWI-CM3 TCo319 exhibits both ENSO regularity (Fig. 2a) and ENSO\namplitude (Fig. 2b) similar to the observations and that the increase of both metrics in the AWI-\nCM3 TCo319 future projections cannot be explained by internal variability. Furthermore, analysis\nof the CMIP6 model archive demonstrates that while there is a wide range of possible changes in\nfuture ENSO regularity (Fig. 2c) and amplitude (Fig. 2d), 55% of CMIP6 models show an increase\nin ENSO regularity (Fig. 2c) and 82% an increase in ENSO SST anomaly amplitude (Fig. 2d) by\nthe second half of this century under the SSP5-8.5 scenario. In addition, a few models (i.e., E3SM-1-1,\nEC-Earth3, and EC-Earth3-Veg) show similar increased ENSO regularity and amplitude as\nAWI-CM3 TCo319 (Fig. 2c,d, Supplementary Fig. 1). We emphasize that we do not expect strong\ninter-model agreement of ENSO projections across the CMIP6 archive given that considerable\nbiases of climate mean states and ENSO dynamics result in different ENSO regimes across models\n8,24.\nI was also wondering how the method used to analyse the amplitude and phase\nsynchronization would compare with the methodology developed by Jajcay et al.\n(2018; https://doi.org/10.1038/s41612-018-0043-7).\nI would therefore recommend major revisions.\nThank you for raising this point. The methodology of Jajcay et al. (2018) primarily examines\nintra-series phase interactions by analysing how diferent time scales within the same time\nseries (e.g., Niño3.4 at diferent wavelet periods) interact using a phase-dynamics\nframework combined with conditional mutual information and statistical randomization.\nThis approach is well-suited for identifying scale–scale interactions or cross-frequency\ncoupling within a single process.\nIn contrast, our analysis targets inter-series phase synchronization — that is, how two\ndistinct processes (e.g., Niño3.4 and other climate mode indices) exhibit joint phase\nlocking. We employ Hilbert-transform–based analytic signals and directly test for m:n\nphase synchronization through the evolution and statistical distribution of phase\ndiferences, rather than through scale-resolved mutual information. Our method\nemphasizes dynamical phase locking between separate time series, whereas the Jajcay et\nal. approach emphasizes interactions among diferent oscillatory components of the same\nsignal. In short, Jajcay et al. investigate cross-scale interactions within a single process,\nwhile our method examines phase synchronization between diferent processes. Both\nmethods are complementary but designed to address diferent scientific questions.\nReviewer #2 (Remarks to the Author):\nThe authors present a well-written, timely study on future ENSO intensification and global\nsynchronization under greenhouse warming, using a high-resolution climate model. The\nabstract clearly communicates key findings, especially the shift to a highly regular, high-\namplitude ENSO regime and its impacts on other climate modes. This addresses an\nimportant unresolved question, and the proposed mechanism linking air-sea feedbacks,\natmospheric noise, and teleconnection synchronization is compelling. I recommend the\nmanuscript for publication, provided the authors address some important points.\nThank you for the assessment and valuable feedback. We address your specific comments\nbelow.\n• A major limitation is the reliance on a single climate model (AWI-CM3) to infer future\nchanges in ENSO amplitude, teleconnection strength, and global hydroclimate impacts.\nGiven well-known inter-model diferences in ENSO characteristics, the lack of multi-model\ncomparison raises concerns about the robustness and generalizability of results.\n• Additionally, the study lacks validation of modeled trends against observational or\nreanalysis data. Without benchmarking AWI-CM3’s current ENSO behavior and\natmospheric responses, confidence in future projections is reduced. Incorporating\nmulti-model analyses and observational validation would greatly enhance credibility and\napplicability.\nThank you for raising these important points. To address these issues, we included in the\nrevised version:\n(i) A detailed model evaluation, showing that a 150-yr AWI-CM3 1950 control\nsimulation reproduces both the observed climate mean state and SST variability\nreasonably well (new Supplementary Fig. 11) as well as reasonably reproduces\nkey ENSO metrics, such as its seasonal synchronization, spectral\ncharacteristics, and its teleconnections (new Supplementary Fig. 12). These\nresults are presented in the revised manuscript (L294–299):\nA 150-year TCo319 control simulation with fixed 1950 forcing (CTL1950) was also\nanalysed to characterize internal climate variability and to serve as a baseline\ncomparison. The CTL1950 simulation reproduces the observed SST and precipitation\nmean state and SST variability reasonably well (Supplementary Fig. 11). The CTL1950\nsimulation also reproduces the observed characteristics of ENSO including its\nseasonal synchronization, spectral characteristics, and teleconnections\n(Supplementary Fig. 12).\n(ii) A comparison of future changes of both ENSO regularity and amplitude in the\nAWI-CM3 model with the CMIP6 model archive (new Fig. 2). While the CMIP6\nmodels show a large range of projected changes (new Fig. 2c,d), a few models\nexhibit similar behaviour as AWI-CM3 (new Supplementary Fig. 1), providing\nsupport that the AWI-CM3 results are not spurious. Corresponding changes to\nthe revised manuscript text (L92–106, a new section titled “Assessment of model\nfidelity.”):\nAssessment of model fidelity\nNext, using sample entropy (SampEn; see Methods) of Niño3.4 SST anomalies as a\nmetric for ENSO regularity together with the Niño3.4 SST anomaly standard deviation\n(SD), we show that a 1950 control simulation of AWI-CM3 TCo319 exhibits both ENSO\nregularity (Fig. 2a) and ENSO amplitude (Fig. 2b) similar to the observations and that\nthe increase of both metrics in the AWI-CM3 TCo319 future projections cannot be\nexplained by internal variability. Furthermore, analysis of the CMIP6 model archive\ndemonstrates that while there is a wide range of possible changes in future ENSO\nregularity (Fig. 2c) and amplitude (Fig. 2d), 55% of CMIP6 models show an increase\nin ENSO regularity (Fig. 2c) and 82% an increase in ENSO SST anomaly amplitude\n(Fig. 2d) by the second half of this century under the SSP5-8.5 scenario. In addition, a\nfew models (i.e., E3SM-1-1, EC-Earth3, and EC-Earth3-Veg) show similar increased\nENSO regularity and amplitude as AWI-CM3 TCo319 (Fig. 2c,d, Supplementary Fig.\n1). We emphasize that we do not expect strong inter-model agreement of ENSO\nprojections across the CMIP6 archive given that considerable biases of climate mean\nstates and ENSO dynamics result in different ENSO regimes across models 8,24.\nIn addition, we revised the text in the abstract, introduction, and discussion making it clear\nthat the results are one possible scenario within a wide range of possible ENSO future\nbehaviours:\nAbstract (with tracked changes):\nChanges to the text in the in introduction and discussion sections in italics:\nL65-70:\nWhile models from the Coupled Model Intercomparison Project Phase 6 (CMIP6) show a wide\nrange of possible future ENSO regularity and amplitude projections, a few models show\nqualitatively similar behaviour to AWI-CM3. If this peculiar ENSO dynamical scenario were to\nmaterialize in the future, it could lead to both increased ENSO predictability, due to greater\nregularity, and, at the same time, to warming-amplified “whiplash impacts” 28 on regional\nclimates.\nL272-274:\nWhile similar increases in ENSO regularity and amplitude can be found in some CMIP6 model\nprojections (Fig. 2), future work needs to assess the detailed dynamics across these models and\nassess the likelihood of the ENSO regime changes seen in AWI-CM3.\nSupplementary Fig. 11| Fidelity of AWI-CM3 in simulating the mean state and SST\nvariability. a,b Annual mean SST climatology in ERA5 (1940-2024) and the AWI-CM3 CTL1950\nsimulation; c,d Annual mean precipitation climatology in ERA5 (1940-2024) and the AWI-CM3\nCTL1950 simulation; e,f Monthly SST anomaly standard deviation in ERA5 (1940-2024) and the\nAWI-CM3 CTL1950 simulation.\nSupplementary Fig. 12| Fidelity of AWI-CM3 in simulating ENSO and its teleconnections. a\nSeasonal standard deviation of the Niño3.4 SST anomaly index in ERA5 (gray bars) and the AWI-\nCM3 CTL1950 simulation (red line); b Power spectral density (PSD) estimated using the Multi-\nTaper Method, of the Niño3.4 SSTA index for ERA5 (1940-2024, black curve) and the AWI-CM3\nCTL1950 simulation (red curve). The confidence levels (CL) for the spectral peaks were derived\nfrom the percentile distributions at each frequency across 10,000 synthetic power spectra generated\nfrom discrete AR(1) processes for the Niño3.4 index with the same lag(-1 month) autocorrelation,\nvariance, and data length to the corresponding time series (See Methods). The red shading in a and\nb is one standard deviation spreads across five non-overlapping 30-yr segments of the 150-yr\nCTL1950 simulation. c,d Regressed December-January-February (DJF) anomalies onto the\nnormalized Niño3.4 index for ERA5 (1940-2024) and the AWI-CM3 CTL1950 simulation for SST\n(°C), land precipitation (mm day-1), SLP (contour interval of 0.3 hPa) anomalies.\n• The suggestion that future ENSO will become more regular, intense, and predictable—\nsynchronizing other major climate modes like the North Atlantic Oscillation and Indian\nOcean Dipole—is speculative. ENSO predictability as well as teleconnections exhibit\nsubstantial variability influenced by internal dynamics and background state changes.\nAssuming a more regular ENSO leads to improved predictability or consistent control over\nother modes oversimplifies the complex climate system. Without supporting observational\nor multi-model evidence, this conclusion should be framed more cautiously, clearly\ndistinguishing model results from real-world implications. The authors may consider\nreferencing relevant literature on ENSO predictability for discussion or introduction.\nhttps://www.nature.com/articles/s41612-024-00845-5\nThank you for these comments. We revised the text (see reply to your comments above) to\nemphasize that the results discussed here are a possible scenario but that we do not\nassess its likelihood. We demonstrate in the revised manuscript using an analysis of ENSO\nregularity and amplitude changes in the CMIP6 archive that a few other models exhibit an\nincrease in both ENSO amplitude and regularity. In addition, we now also include a\nreference to a previous study (Berner et al. 2020), which found an increase in ENSO\nregularity and amplitude in the CMIP5 model CESM1.\nThank you for pointing out the Ehsan et al. (2024) study. We now cite this paper in addition\nto two other papers on decadal variations of ENSO predictability when we discuss how\nsample entropy captures observed decadal variations in ENSO forecast skill:\nL381-383:\nThe observed SampEn effectively captures the decadal variations in ENSO forecast skill,\nincluding the modest decrease in predictability after the 2000s (black curve in Fig. 2a, refs 59-61).\nReviewer #3 (Remarks to the Author):\nStuecker et al. investigates climate change experiments in a high resolution coupled GCM\nunder a high-CO2-emissions scenario and finds a strong increase in the amplitude and\nperiodicity in ENSO. Consequently, other climate modes phase lock onto the strong and\nperiodic ENSO signal and increase in amplitude. They use a hierarchy of idealized models\nto diagnose the changes in ENSO amplitude and periodicity. The societal impacts of such a\ntransition are manifest, but it remains unclear whether this feature is a peculiarity of the\nchosen model or whether it reflects reality.\nThe paper is clearly written, the results are well-presented, and the importance of the\nfindings are clear. I consider the paper suitable for publication after minor revisions.\nHowever, the paper would be improved if the authors provide some additional commentary\nand analysis on why this model exhibits this unique ENSO behavior, but not any other\nmodels. I ofer some small suggestions to improve the manuscript.\nThank you for the assessment and valuable feedback. We address your specific comments\nbelow.\nMajor Comments\nThe results of this paper are very interesting, but one important point is not addressed: why\ndoes the AWI-CM3 model simulate a rapid switch to a highly regular oscillation with high\namplitude, while other climate models do not? Since this model is an outlier, the impact of\nthe paper will be increased if an argument can be made that this model’s ENSO response is\nin some respects more realistic than other models.\nA few ideas are hinted at throughout the paper, such as intensifying tropical intraseasonal\nvariability (Ln. 111), or the relaxation of the mean SST gradient and trade winds (Ln.130).\nBut these changes are seen in other climate models as well (especially the weakening of\nthe east-west SST gradient), so what makes AWI-CM3 an outlier? Is its high horizontal and\nvertical resolution important? If so, do we see this increase in other HiResMIP models?\nIf the authors can discuss possible reasons conventional CMIP6 models do not simulate\nthe response, it would be great. Furthermore, if they can speculate on why this particular\nmodel might better reflect reality compared to other CMIP models, that would further\nimprove the impact. But overall it is still a very nice result that is worth being published.\nWe added an analysis of the CMIP6 archive to address this important point (see new Fig. 2\nin the revised manuscript and also our response to reviewer #1). We find while models\nshow a large range of possible projections, a few models show qualitatively similar\nbehaviour of future projected ENSO regularity and amplitude changes to AWI-CM3.\nWe conjecture that multiple criteria likely need to be fulfilled for a model to exhibit a similar\nENSO regime shift:\n(i) relatively weak biases in the climate mean state and its variability (new\nSupplementary Figs. 11 and 12)\n(ii) relatively weak biases in ENSO’s interaction with other climate modes. We\nshowed in the supplementary material of Zhao et al. (2024) using the XRO\nframework that CMIP models tend to have large biases in these interaction\nterms.\nSee the following new additions to the results (including a new “assessment of model\nfidelity” section) summarizing these points (changes marked in italics):\nL65-70:\nWhile models from the Coupled Model Intercomparison Project Phase 6 (CMIP6) show a wide\nrange of possible future ENSO regularity and amplitude projections, a few models show\nqualitatively similar behaviour to AWI-CM3. If this peculiar ENSO dynamical scenario were to\nmaterialize in the future, it could lead to both increased ENSO predictability, due to greater\nregularity, and, at the same time, to warming-amplified “whiplash impacts” 28 on regional\nclimates.\nL92-106:\nAssessment of model fidelity\n“Next, using sample entropy (SampEn; see Methods) of Niño3.4 SST anomalies as a metric for\nENSO regularity together with the Niño3.4 SST anomaly standard deviation (SD), we show that\na 1950 control simulation of AWI-CM3 TCo319 exhibits both ENSO regularity (Fig. 2a) and\nENSO amplitude (Fig. 2b) similar to the observations and that the increase of both metrics in the\nAWI-CM3 TCo319 future projections cannot be explained by internal variability. Furthermore,\nanalysis of the CMIP6 model archive demonstrates that while there is a wide range of possible\nchanges in future ENSO regularity (Fig. 2c) and amplitude (Fig. 2d), 55% of CMIP6 models\nshow an increase in ENSO regularity (Fig. 2c) and 82% an increase in ENSO SST anomaly\namplitude (Fig. 2d) by the second half of this century under the SSP5-8.5 scenario. In addition, a\nfew models (i.e., E3SM-1-1, EC-Earth3, and EC-Earth3-Veg) show similar increased ENSO\nregularity and amplitude as AWI-CM3 TCo319 (Fig. 2c,d, Supplementary Fig. 1). We emphasize\nthat we do not expect strong inter-model agreement of ENSO projections across the CMIP6\narchive given that considerable biases of climate mean states and ENSO dynamics result in\ndifferent ENSO regimes across models 8,24.”\nWe emphasize that more work is needed to investigate inter-model diferences that are\nbeyond the scope of this paper. We added the following statement in the discussion\nsection:\nL272-274:\nWhile similar increases in ENSO regularity and amplitude can be found in some CMIP6 model\nprojections (Fig. 2), future work needs to assess the detailed dynamics across these models and\nassess the likelihood of the ENSO regime changes seen in AWI-CM3.\nMinor Comments\nLn. 52-53. It could be interesting to show Fig 1B for some selected CMIP6 models for\ncomparison, to demonstrate that this behavior of this model is unique.\nFollowing your suggestion, we have selected the few models that show a qualitatively\nsimilar behaviour (new Supplementary Fig. 1). Most of the models do not exhibit this\nbehaviour (see Fig. R1 in response to reviewer #1).\nSupplementary Fig. 1| Increased ENSO regularity in selected CMIP6 models. a Time\nevolution of the ENSO Nino3.4 SSTA index [K] for EC-Earth3 (3-month running-mean applied).\nb Wavelets of the normalized Nino3.4 SSTA index for EC-Earth3. Black contours enclose regions\nabove 90% confidence level significance tested against a Lorentzian spectrum. c,d and e,f Same\nas a,b but for EC-Earth3-Veg and E3SM-1-1, respectively.\nLn. 123 It’s not obvious to me that ZA is very significant, it looks like it averages about zero\nover the annual cycle.\nThank you for catching this. We revised the text accordingly.\nLn. 174. I’m confused why, if the range of phase diferences narrows, the standard\ndeviation of the PDFs trends upwards (blue curve)? Intuitively, I would expect the standard\ndeviation to become smaller.\nWe apologize for the unclear description. We revised the methods and main text to clarify\nthis point (L552-561; note that the previous version of the manuscript incorrectly stated 2:1\ninstead of 1:1 phase synchronization, which has been corrected in the revised version). In\naddition, we added a new supplementary figure (Supplementary Fig. 8) to illustrate the\nphase synchronization metric.\nL184-186:\nA bounded phase difference indicates stronger phase-coupling between the climate modes under\nconsideration (Supplementary Fig. 8). Here we focus on 1:1 phase synchronization, meaning\nthat the oscillators maintain a (near-)constant phase difference.\nL552-561:\nWe use this approach to calculate the 1:1 difference of unwrapped phases between timeseries\ncharacterizing different climate modes simulated by the AWI-CM3 model. All timeseries were\nbandpass filtered between 2.5-months and 10-years using a Butterworth filter prior to the\nanalysis to focus on annual to interannual variability and to satisfy Bedrosian’s theorem. We use\nthe standard deviation of the phase difference PDF, 𝜎[𝑃𝐷𝐹(𝛿𝷿 )], as a measure of the\n!,!\nsynchronization strength. If there is no synchronization, the phase differences are uniformly\ndistributed, yielding a flat PDF and 𝜎[𝑃𝐷𝐹(𝛿𝷿 )]) ≈ 0 (Supplementary Fig. 8a,b). In\n!,!\ncontrast, if there is strong synchronization, the phase differences concentrate around a preferred\nvalue (−π/2 for the example shown in Supplementary Fig. 8c,d), producing a sharply peaked\nPDF, resulting in a larger 𝜎[𝑃𝐷𝐹(𝛿𝷿 )] (Supplementary Fig. 8c,d).\n!,!\nSupplementary Fig. 8| Illustration of phase synchronization analysis. a,b No phase\nsynchronization case: Two independent white-noise time series (a) and the probability density\nfunction (PDF) of their phase difference (b). c,d Strong phase synchronization case: two idealized\n4-year cosine time series with small white-noise perturbation (c) and the PDF of their phase\ndifference (d). In b and d, the synchronization strength defined by the standard deviation of the\nphase difference PDF (𝜎[PDF(𝛿Φ )]).\n!,!\nFig.4: The features on the figure are very small. Perhaps three rows and two columns will\nallow each sub-panel to be larger.\nRevised as suggested. Thanks for the helpful suggestion."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I believe the authors have now addressed all my comments. The paper is significantly stronger, and deserves to be accepted for publication in Nature Communications.", "reviewer_3": "I am satisfied with the authors’ careful reply to my comments and consider the paper ready for publication. Congratulations on an excellent paper."}, "author_response": "Manuscript NCOMMS-25-44516A\nPoint-by-point response to the reviewers\nWe have included a point-by-point response to the comments below (Reviewer comments\nare in black font, our responses in blue font).\nReviewer #1 (Remarks to the Author):\nI believe the authors have now addressed all my comments. The paper is significantly\nstronger, and deserves to be accepted for publication in Nature Communications.\nWe sincerely thank you for the positive assessment and supportive recommendation.\nReviewer #3 (Remarks to the Author):\nI am satisfied with the authors’ careful reply to my comments and consider the paper ready\nfor publication. Congratulations on an excellent paper.\nThank you again for the constructive comments throughout the review process.\nReviewer #3 (Remarks on code availability):\nCurrently, the Github links to the code do not work.\nThe code is now publicly available, and the links along with permanent DOIs have been\nupdated in the revised manuscript."}]

Decision letter