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

Article https://doi.org/10.1038/s41467-025-64619-0 Global climate mode resonance due to rapidly intensifying El Niño-Southern Oscillation Malte F. Stuecker 1,2,9 , Sen Zhao 3,9 , Axel Timmermann 4,5 , Rohit Ghosh 6 , Tido Semmler 7 ,Sun-SeonLee 4,5 , Ja-Yeon Moon 4,5 , Fei-Fei Jin 2,3 & Thomas Jung 6,8 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 moderateamplitude 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, whichapproachcriticalityinthesecond 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. 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,hinderingconfidence in regional climate projections. The dynamics of ENSO are governed by a delicate balance of positive and negative feedbacks8,9that determine both ENSO’sinstability 10 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 Received: 10 June 2025 Accepted: 23 September 2025 Check for updates 1 Department of Oceanography, University of HawaiʻiatMānoa, Honolulu, USA. 2 International Pacific Research Center, University of HawaiʻiatMānoa, Honolulu, USA. 3 Department of Atmospheric Sciences, University of HawaiʻiatMānoa, Honolulu, USA. 4 Center for Climate Physics, Institute for Basic Science, Busan, Republic of Korea. 5 Pusan National University, Busan, Republic of Korea. 6 Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany. 7 Met Éireann, Dublin, Ireland. 8 Department of Physics and Electrical Engineering, University of Bremen, Bremen, Germany. 9 These authors contributed equally: Malte F. Stuecker, Sen Zhao. e-mail: st[email protected];[email protected] Nature Communications | (2025) 16:9013 1 1234567890():,; 1234567890():,; 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 2020 2030 2040 2050 2060 2070 2080 2090 2100 −3 −2 −1 0 1 2 3 SST anomaly (K) aENSO SST index (member=0) 2020 2030 2040 2050 2060 2070 2080 2090 2100 0.5 1 2 4 8 16 Period (years) bEnsemble-mean wavelet of ENSO SSTA index 2020 2030 2040 2050 2060 2070 2080 2090 2100 −2 −1 0 1 2 Normalized cNAO index (member=0) 2020 2030 2040 2050 2060 2070 2080 2090 2100 0.5 1 2 4 8 16 Period (years) dEnsemble-mean wavelet of NAO index 0.062 0.125 0.25 0.5 1 2 4 8 16 32 64 128 0.062 0.125 0.25 0.5 1 2 4 8 16 e SSTA SD change (P2 - P1) f SLPA DJF SD change (P2 - P1) −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 K −135 −90 −45 045 90 135 Pa Fig. 1 | Oceanic and atmospheric variability in response to greenhouse warming. a,cTime evolution of the El Niño-Southern Oscillation (ENSO) sea surface temperature (SST) index [K] and the normalized North Atlantic Oscillation (NAO) index [n.u.] for ensemble member 0 (3-month running-mean applied). b,dEnsemble-mean (four members post year 2055) wavelets of the ENSO SST and NAO indices. Black contours enclose regions above 90% confidence level significance tested against a Lorentzian (b) and a white noise (d)spectrum,respectively. e,fSpatial maps of the ensemble-mean standard deviation changes between period 2 (P2, 2080–2100) and period 1 (P1, 2015–2035) of the monthly SST anomaly (SSTA) (e) and December–January–February (DJF) sea level pressure (SLP) anomaly (SLPA) (f). Dots indicate regions where all four ensemble members show disagreement inthe sign ofchanges. Boxes in (e) indicate the cold tongue ENSO region as well as the key regions associated with SST variability of the North PacificMeridional Mode (NPMM), Indian Ocean Dipole (IOD), Indian Ocean Basin (IOB) mode, and the Tropical North Atlantic (TNA) mode. Boxes in (f) indicate the SLP centers of the NAO. Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 2 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, 2aand3a). 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 nearannual 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ño3.4 SST anomalies as a metric for ENSO regularity together with the Niño3.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 thatwhile there is a wide range of possible changes in future ENSO regularity (Fig. 2c) and amplitude (Fig. 2d),55%ofCMIP6modelsshowanincreaseinENSOregularity (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 modulations9and encapsulates the fundamental ENSO dynamics in two coupled ordinary differential equations (ODEs) for SST anomalies in the Niño3.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 PacificMeridionalmode(NPMM),TNA,IOD,andtheIndian 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.Thisallowsustocomputethetime-dependentlinear growth rates and frequencies of ENSO via Floquet eigen analysis31 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 1900 1925 1950 1975 2000 2025 2050 2075 Center year of 21-yr moving window 0.6 0.8 1.0 1.2 1.4 Sample entropy (SampEn) aENSO regularity in Observation & AWI-CM3 Observation AWI-CM3 CTL1950 AWI-CM3 SSP585 1900 1925 1950 1975 2000 2025 2050 2075 Center year of 21-yr moving window 0.6 0.8 1.0 1.2 1.4 1.6 Standard deviation (°C) bENSO amplitude in Observation & AWI-CM3 Observation AWI-CM3 CTL1950 AWI-CM3 SSP585 Observation (1950-2024) CTL1950 SSP585 (2050-2100) Observation (1950-2024) CTL1950 SSP585 (2050-2100) −0.3 0.0 0.3 2050-2100 minus 1900-2000 (SampEn) 55% Consensus cENSO regularity changes in CMIP6 −0.5 0.0 0.5 1.0 2050-2100 minus 1900-2000 (°C) 82% Consensus dENSO amplitude changes in CMIP6 AWI-CM3 E3SM-1-1 EC-Earth3 EC-Earth3-Veg Fig. 2 | El Niño-Southern Oscillation (ENSO) regularity and amplitude changes in observations, Alfred Wegener Institute Climate Model (AWI-CM3), and Coupled Model Intercomparison Project Phase 6 (CMIP6) models. Moving 21year changes of sample entropy (SampEn) (a) and standarddeviation (SD) (b)ofthe Niño3.4 SST anomaly index in observations (black), the 150-year AWI-CM3 1950 control simulation (CTL1950, blue), and the AWI-CM3 shared socio-economic pathway (SSP)5-8.5 scenario simulation (red). Solid lines and shading indicate the multi-product/multi-member mean and the 1 SD spread of 10,000 inter-realizations using a bootstrap method, respectively. The error bars indicate the mean and 1 SD of values across all overlapping 21-year windows for the observations (1950–2024), CTL1950 (150-yr), and SSP5-8.5 (2050–2100). Violin plots showing the probability density distribution of ENSO regularity change (c) and amplitude change (d)across 49 CMIP6 models. The change is calculated as the difference between 2050–2100 and 1900–2000. Dots represent individual model results, with E3SM-1-1, EC-Earth3, and EC-Earth3-Veg highlighted in distinct colors; the red star marks the AWI-CM3 result (2050–2100 minus CTL1950). In aand c,they-axis is inverted to emphasize decreasing SampEn, indicating increasing regularity. Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 3 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 highvariance 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ño3.4 index (Fig. 1b) with thelinearENSO 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’sgrowing 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’snetgrowthrate,wefind 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 Rcan 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 Fig. 3 | El Niño-Southern Oscillation (ENSO) amplitude, Bjerknes stability, and wind stress structure changes. a Ensemble-mean moving 21-year seasonally varying (shading) and monthly (blue curve) standard deviation (SD) changes of the ENSO sea surface temperature anomaly (SSTA) index relative to the monthly SD during 2015–2035 (ratio). bGrowth rate of the leading eigenmode (ENSO) obtained via Floquet analysis for the Recharge Oscillator (RO) (blue curve) and eXtended Recharge Oscillator (XRO) (magenta curve) models in a moving 21-year window. The seasonal RO growth rate is indicated in shading. cNoise amplitude defined as the seasonally varying SD (shading) and monthly SD (blue curve) of the XRO fit residual. Different ENSO feedbacks (as a function of calendar month) obtained via Bjerknes stability analysis for period 1 (P1) (d, 2015–2035; one ensemble member) andperiod2(P2)(e, 2080–2100; ensemble-mean of four members), as well as their difference (f). The horizontal axes in d–fshow the Bjerknes (BJ) index, SSTA growth rate ðRÞ, ocean damping rate (ε), thermocline feedback (TH), zonal advective feedback (ZA), Ekman feedback (EK), thermal damping by the net surface heat flux (TD), dynamic damping by mean horizontal currents (DD), all residuals from nonlinear dynamic heating and unresolved processes (AllRes), as well as the dynamic damping components by mean zonal (DDu) and meridional (DDv) currents and the thermal damping components by shortwave radiation (Q SW ), longwave radiation (Q LW ), latent heat flux (Q LH ), and sensible heat flux (Q SH ). gLatitudinal profiles (averaged of 160°E–150°W) of regressed zonal wind stress anomalies τ0 xonto the monthly ENSO SST index in a moving 21-year window. Dashed curves and shadings indicate the ensemble mean and one SD spread of meridional boundaries defined by the local minima. h,iBasin thermocline adjustment rate (ε) and Recharge/discharge efficiency (F2)inamoving21-year window. In hand i, red curves denote the meridional width of ENSO zonal wind stress anomalies τ0 x(see “Methods”). Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 4 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 Rand ε 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 budget4will 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 2030 2040 2050 2060 2070 2080 2090 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec aAmplitude of NAO 2030 2040 2050 2060 2070 2080 2090 −π −π 2 0 π 2 π δΦ1, 1 fPhase-synchronization of NAO with Niño3.4 2030 2040 2050 2060 2070 2080 2090 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec bAmplitude of TNA 2030 2040 2050 2060 2070 2080 2090 −π −π 2 0 π 2 π δΦ1, 1 gPhase-synchronization of TNA with Niño3.4 2030 2040 2050 2060 2070 2080 2090 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec cAmplitude of IOB 2030 2040 2050 2060 2070 2080 2090 −π −π 2 0 π 2 π δΦ1, 1 hPhase-synchronization of IOB with Niño3.4 2030 2040 2050 2060 2070 2080 2090 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec dAmplitude of IOD 2030 2040 2050 2060 2070 2080 2090 −π −π 2 0 π 2 π δΦ1, 1 iPhase-synchronization of IOD with Niño3.4 2030 2040 2050 2060 2070 2080 2090 Center year of 21-yr moving window Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec eAmplitude of NPMM 2030 2040 2050 2060 2070 2080 2090 Center year of 21-yr moving window −π −π 2 0 π 2 π δΦ1, 1 jPhase-synchronization of NPMM with Niño3.4 0 10 20 σ change ratio (%) 0 20 40 60 80 σ change ratio (%) 0 25 50 75 σ change ratio (%) 0 20 40 60 80 σ change ratio (%) 0 20 40 σ change ratio (%) −160 −120 −80 −40 0 40 80 120 160 σ change ratio (%) relative to 2015-2035 0.06 0.08 0.10 0.12 σ[PDF(δΦ1, 1)] 0.10 0.15 0.20 σ[PDF(δΦ1, 1)] 0.15 0.20 0.25 σ[PDF(δΦ1, 1)] 0.06 0.08 0.10 σ[PDF(δΦ1, 1)] 0.06 0.08 0.10 0.12 σ[PDF(δΦ1, 1)] 0.15 0.25 0.35 0.45 0.55 0.65 PDF(δΦ1, 1) Fig. 4 | Amplitude and phase synchronization changes of climate modes. a–eEnsemble-mean moving 21-year seasonally varying (shading) and monthly (blue curve) standard deviation (SD) changes of each climate mode relative to its monthly SD during 2015–2035 (ratio) for the North Atlantic Oscillation (NAO), Tropical North Atlantic (TNA) mode, Indian Ocean Basin (IOB) mode, Indian Ocean Dipole (IOD), and North Pacific Meridional Mode (NPMM) indices, respectively. f–jEnsemble-mean moving 21-year phase synchronization quantified by the histogram of phase differences (shading) between each climate modeand the Niño3.4 index for NAO, TNA, IOB, IOD, and NPMM, respectively. The blue curves in f–jindicate synchronization strength defined by the SD of the histogram density (PDF) of phase differences. The blue shading in a–jindicates the one SD spread among the four ensemble members. Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 5 amplitude, ranging from 40–75%. All these modes are known to interact with ENSO18,23,24. To further elucidate the apparent emergent 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 (f E )as well as at the near-annual ENSO combination tone frequencies (1±f E ) above a Lorentz/Hasselmann background spectrum39, consistent with theoretical expectations17,18. The emergence of interannual and nearannual ENSO timescales in the spectraof 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ño3.4 index (representing ENSO) and the other indices (Fig. 4f–j). A bounded phase difference indicates stronger phasecoupling 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, 2year-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 PacificNorth 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ño3.4 SST anomalies and Niño4 precipitation shows no substantial impact in Europe. During P2, for the regressions with a Reg PrA/Z500A onto normalized Niño3.4 SSTA over P1 (2015-2035) d Reg PrA/Z500A onto normalized Niño3.4 SSTA over P2 (2080-2100) b Reg PrA/Z500A onto Niño3.4 SSTA over P1 (2015-2035) e Reg PrA/Z500A onto Niño3.4 SSTA over P2 (2080-2100) c Reg PrA/Z500A onto Niño4 PrA over P1 (2015-2035) f Reg PrA/Z500A onto Niño4 PrA over P2 (2080-2100) −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 mm day−1 −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 mm day−1 K−1 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 mm day−1 (mm day−1)−1 Fig. 5 | Changes of the El Niño-Southern Oscillation (ENSO) Northern Hemisphere teleconnection during boreal winter. a Regression of normalized Niño3.4 sea surface temperature anomalies (SSTA) with precipitation anomalies (PrA) (shading, [mm day−1]) and 500 hPa geopotential height anomalies (Z500A) (contours, ±5, ±15, …[m]) during December–January–February (DJF) 2015–2035 (P1); bRegression of Niño3.4 SST anomalies with precipitation anomalies (shading, [mm day−1K−1]) and Z500 anomalies (contours, ±5, ±15, …[m K−1]) during DJF 2015–2035; cRegression of Niño4 precipitation anomalies with precipitation anomalies (shading, [mm day−1(mm day−1)−1]) and Z500 anomalies (contours, ±5, ±15, …[m (mm day−1)−1]) during DJF 2015–2035. d–fSame as a–cbut for ensemblemean results during DJF 2080–2100 (P2), respectively. Stippling and dark contours indicate statistical significance at the 95% level using a one-sided Student’sttest. Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 6 both the Niño3.4 SST index in degrees Kelvin (Fig. 5e) and with the normalized Niño3.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ño3.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ño3.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. Reference40 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. Reference41,usingapacemaker 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). Discussion Previous work has shown that the atmosphere is “ringing”with a range of combination tones in response to ENSO forcing, causing an “ENSO frequency cascade”43. Here we show that ENSO not only causes these deterministic signals in the atmosphere, but its intensification can also energize both atmospheric and air-sea coupled climate modes across the globe in response to greenhouse warming. This global climate resonance towards the strengthening ENSO signal is apparent in the raw timeseries (Fig. 1), their amplitudes, spectra, and phase synchronization characteristics (Fig. 4and Supplementary Fig. 6). ENSO variability increases in the AWI-CM3 simulation because of an intensification of atmospheric noise in boreal summer (Fig. 3candSupplementary Fig. 2c), reduced damping from thermocline adjustment, and the modest increase of the Bjerknes feedback (Fig. 3f). Its global synchronization towards the extratropics (e.g., with the NAO) is boosted by the growing ENSO amplitude, but also by an overall reorganization of atmospheric teleconnection patterns (Fig. 5), due to a different future atmospheric mean flow. The changes in the mean flow —particularly a stronger, more zonally extended Pacificjetand increased Atlantic baroclinicity—can enhance the propagation, amplification, and momentum deposition of ENSO-forced stationary waves into the North Atlantic sector, thereby strengthening the teleconnection to the NAO. Increasing ENSO amplitude and teleconnection patterns imply that remote extratropical precipitation responses, such as in Southern California and the Iberian Peninsula (Fig. 5), will also become stronger between alternating El Niño and La Niña events. Even though the recurrence of these events may become more predictable due to the increased regularity/periodicity and amplitude (Fig. 1)44,futureENSO teleconnections might generate enhanced whiplash effects on hydroclimate, which requires additional planning and management strategies to minimize the costs of climate damage. While similar increases in ENSO regularity and amplitude can be found in some CMIP6 model projections (Fig. 2), future work needs to assess the detailed dynamics across these models and assess the likelihood of the ENSO regime changes seen in AWI-CM3. Methods Earth system model experiments This study employs the Alfred Wegener Institute Climate Model (AWICM3), which couples the OpenIFS (Open Integrated Forecasting System) atmosphere model with the FESOM2 (Finite-volume Sea IceOcean Model) ocean model. The atmospheric component operates at a horizontal resolution of approximately 31km (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 AWICM3 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 highemission 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)for1871–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 Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 7 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 SSP58.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ño3.4 SST index is defined as SST anomalies averaged over the Niño3.4 region (170°–120°W, 5°S–5°N). The Niño4 precipitation index is defined as precipitation anomalies averaged over the Niño4 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 definedasSSTanomaliesaveraged 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 associated with the leading Empirical Orthogonal Function (EOF) of area-weighted SLP anomalies over the North Atlantic sector 90°W–40°E, 20°–80°N52. 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 analysis was conducted on normalized Niño3.4 SST and NAO indices, with the resulting wavelet power spectra computed 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ño3.4 SST was tested against an AR(1) process with a lag(−1month) 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 =256forP1andnfft =1024 for P2 54. 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(−1month) 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ño3.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,SampEnislesssensitivetorecord 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 fx1,x2,...,xNg, SampEn estimates the negative natural logarithm of the conditional probability that two sequences of length m thatmatchwithintolerancerwill 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 ofthetimeseries.TheSampEnmetricisdefined as: SampEn m,r,NðÞ=ln AmrðÞ BmrðÞ  ð1Þ where AmrðÞthe probability that two sequences of length m+1 are similar (matches), BmrðÞis the probability that two sequences of length mare similar (candidates). Explicitly, AmrðÞ=1 Nm+1 1 NmX Nm i=1 X Nm j=1, j≠i number of times thatdu m+1 jðÞ,um+1 iðÞ  <r  ð2Þ BmrðÞ=1 Nm+1 1 NmX Nm i=1 X Nm j=1, j≠i number of times that du mjðÞ,umiðÞ  <r  ð3Þ where umiðÞ is the subsequence fxi,xi+1,... ,xi+m1g, du mjðÞ,umiðÞ  =max k=1,...,mjxi+k1xj+k1jis the maximum norm distance. Since the number of matches is always less than or equal to the number of possible vectors, the ratio AmrðÞ=BmrðÞ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 (TENSO) and equatorial Pacific zonal-mean thermocline depth anomalies (h). In its linear form, the model is expressed as: d dt TENSO h  =LENSO TENSO h  +ξENSO ξh  ð4Þ In Eq. (4), LENSO =RF 1 F2ε  ,whereRrepresents the SST anomaly growth rate, collectively describing the Bjerknes feedback, ε Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 8 denotes ocean damping rate related to the energy leakage at the western boundary and mixing, F1represents the effectiveness of the heat content discharge–recharge in controlling SST anomalies, and F2 represents the discharge–recharge efficiency; ξENSO and ξhare the residual terms encompassing non-resolved nonlinear processes and stochastic noise. The complex eigenvalue of the RO model operator LENSO is the Bjerknes-Wyrtki-Jin (BWJ) index for the ENSO growth rate and periodicity and can be written as: BWJ = Rε 2+iffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F1F2R+εðÞ 2 4 sð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 1cos ωt+Ls 1sin ω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 TENSO h TNPMM TIOB TIOD TTNA 0 B B B B B B B B @ 1 C C C C C C C C A =L TENSO h TNPMM TIOB TIOD TTNA 0 B B B B B B B B @ 1 C C C C C C C C A + ξENSO ξh ξNPMM ξIOB ξIOD ξTNA 0 B B B B B B B B @ 1 C C C C C C C C A ð7Þ The linear dynamical operator Lcontains four submatrices, organized as follows: L=LENSO C1 C2LM  ð8Þ where the linear operator submatrix LENSO describes the ENSO internal recharge-discharge dynamics, LMrepresents the internal processes and interactions among the other climate modes; Care coupling submatrices, with C2describing the impact of ENSO on other climate modes and C1describing 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 time scale 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 time scales are estimated as rξ=lnða1Þ=δt, where a1are the lag(−1 month) autocorrelations of the residual of the XRO fit. The range of observed noise time scales r1 ξ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 Floquet ENSO growth rate and periodicity using the RO model operator (LENSO) or XRO model operator Lvia 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 Q0ðÞ=Ito form a monodromy matrix M=QTðÞ,T= 1 year. The Floquet exponents (σj) are calculated from the least damped complex eigenvalues αjof Mby 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: ∂T0 ∂t=Q0∂uT0 ðÞ ∂x∂vT0 ðÞ ∂y∂wT0 ðÞ ∂zu0∂T ∂xv0∂T ∂yw0∂T ∂z u0∂T0 ∂x+v0∂T0 ∂y+w0∂T0 ∂z  +QRes,ð9Þ where the overbars denote the climatological monthly-mean seasonal cycle, and the primes denote anomalies; Tis the ocean temperature, u, v,andwthe ocean zonal, meridional, and vertical current velocities, Q the net surface heat flux effect on the ocean mixed layer, and QRes the unresolved residual. The anomalous heat flux term Q0is calculated as: Q0=Q0 net Q0 bot ρCpHð10Þ where His the ocean mixed layer depth (50 m), ρis the seawater density (1025 kg m−3), Cpis the specificheatcapacityofseawater (3994 J kg−1K−1), and Q0 net represents the net surface heat flux, consisting of four components: Q0 net =Q0 SW +Q0 LW +Q0 LH +Q0 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 (Q0 bot)isparameterized following ref. 65 as: Q0 bot =Q0 SW 0:58eH 0:35 +0:42eH 23  ð12Þ Taking the volume-average by integrating both hand sides of Eq. (9) from the oceanmixed layer depth to the ocean surface and spatially over the eastern equatorial Pacific box where the SST variability of Article https://doi.org/10.1038/s41467-025-64619-0 Nature Communications | (2025) 16:9013 9