# Increased frequency of multi-year El Niño–Southern Oscillation events across the Holocene

Zhengyao Lu, Anna Schultze, Matthieu Carré, Chris Brierley, Peter Hopcroft, Debo Zhao, Minjie Zheng, Pascale Braconnot, Qiuzhen Yin, Johann Jungclaus, et al.

► To cite this version:

Zhengyao Lu, Anna Schultze, Matthieu Carré, Chris Brierley, Peter Hopcroft, et al. Increased frequency of multi-year El Niño–Southern Oscillation events across the Holocene. Nature Geoscience, 2025, 18 (4), pp.337-343. 10.1038/s41561-025-01670-y. hal-05075867

HAL Id: hal-05075867

https://hal.science/hal-05075867v1

Submitted on 21 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Increased frequency of multi-year El Niño–Southern Oscillation events across the Holocene

Received: 21 July 2024

Accepted: 20 February 2025

Published online: 31 March 2025

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

Check for updates

Zhengyao Lu Ⓔ¹, Anna Schultze Ⓔ¹, Matthieu Carré Ⓔ²,³, Chris Brierley Ⓔ⁴, Peter O. Hopcroft Ⓔ⁵, Debo Zhao Ⓔ⁶, Minjie Zheng⁷,⁸, Pascale Braconnot Ⓔ⁹, Qiuzhen Yin Ⓔ¹⁰, Johann H. Jungclaus Ⓔ¹¹, Xiaoxu Shi Ⓔ¹², Haijun Yang Ⓔ¹³ & Qiong Zhang Ⓔ¹⁴

El Niño–Southern Oscillation (ENSO) events, whether in warm or cold phases, that persist for two or more consecutive years (multi-year), are relatively rare. Compared with single-year events, they create cumulative impacts and are linked to extended periods of extreme weather worldwide. Here we combine central Pacific fossil coral oxygen isotope reconstructions with a multimodel ensemble of transient Holocene global climate simulations to investigate the multi-year ENSO evolution during the Holocene (beginning ~11,700 years ago), when the global climate was relatively stable and driven mainly by seasonal insolation. We find that, over the past ~7,000 years, in proxies the ratio of multi-year to single-year ENSO events increased by a factor of 5, associated with a longer ENSO period (from 3.5 to 4.1 years). This change is verified qualitatively by a subset of model simulations with a more realistic representation of ENSO periodicity. More frequent multi-year ENSO events and prolonged ENSO periods are being caused by a shallower thermocline and stronger upper-ocean stratification in the Tropical Eastern Pacific in the present day. The sensitivity of the ENSO duration to orbital forcing signals the urgency of minimizing other anthropogenic influence that may accelerate this long-term trend towards more persistent ENSO damages.

The El Niño–Southern Oscillation (ENSO) is the most prominent year-to-year variability in the climate system¹. It is characterized by transitions between two phases: positive (El Niño) and negative (La Niña) sea surface temperature (SST) anomalies in the tropical Pacific. By disrupting the global weather patterns², the hazards caused by ENSO range from droughts and forest fires to extreme rainfall and flooding in various parts of the world, damaging agriculture, ecosystems and human societies¹,³,⁴

ENSO events follow a distinct cycle with their mature phase in boreal winter and their decline in spring due to a persistence barrier (PB)⁵,⁶. However, ENSO is known to be highly variable in its evolution,

and some ENSO events tend to persist for consecutive years⁷⁻⁹; notably, the recent triple-dip La Niña persisted for three succeeding years¹⁰. It is projected that multi-year La Niña and El Niño will probably increase in the future under anthropogenic forcings¹¹,¹². This trend has already been seen in the past decades where five out of six La Niña events turned into multi-year La Niña¹³. Since 1950 CE, five multi-year El Niño events have occurred, two of which developed in the past decade¹². These multi-year ENSO variations come along with stronger, cumulative impacts compared with single-year ENSO events¹³,¹⁴. For instance, persistent wet conditions over Australia, Indonesia, tropical South America and South America, as well as dry conditions over the southern
<figure><img src="image_4.png" /></figure>

**Fig. 1 | Single-year and multi-year ENSO in proxies and reanalysis.** a–d, The temporal evolution of single-year (a and b) and multi-year (c and d) ENSO cases (red lines, El Niño; blue lines, La Niña), showing composites of selected proxy reconstruction of the coral oxygen isotope anomaly (a and c) and the Niño 3.4 SST anomaly (b and d). The solid lines and the shading show the mean and 1

s.d. of all cases in the corresponding composites. e, The locations of the fossil coral synthesis used in this study (magenta circles). The ERSST monthly SSTA variability (s.d.) during 1854–2023 CE is represented by the shaded colours. The black box is for the Niño 3.4 region.

USA, equatorial Africa, India and southeast China are climate responses linked to multi-year La Niña events¹³. Multi-year El Niño events followed by La Niña events could cause a shift from one extreme to the other in some regions, increasing their severe damages due to their prolonged states¹⁵.

The possibility of more frequent multi-year ENSO has important societal implications. Climate models still have limitations in how they represent various physical processes vital to ENSO periodicity¹⁶. This brings uncertainties to the climate projections of future occurrence of multi-year ENSO. For instance, when investigating such projections, earlier studies usually selected a subset of models based on some aspects of ENSO dynamics, such as ENSO skewness¹¹ and ENSO spatial pattern¹², but those cannot guarantee a realistic simulation of ENSO duration and timing. The chaotic nature of ENSO on annual to centennial timescales¹⁷,¹⁸, combined with the relative scarcity of multi-year events and the modeling challenges, mean that the robustness of shifts in ENSO periodicity over the instrumental record¹²,¹³,¹⁹ is still under debate.

Reconstructions of ENSO in the geological past (palaeo-ENSO archives²⁰,²¹) can greatly increase the sample size of multi-year events as well as explore their response to any external forcings (for example, seasonal insolation). Oxygen isotope records from central equatorial Pacific corals, which highly correlate with the local SST evolution, indicate the existence of multi-year ENSO in the past millennium¹⁵. Over the past decade, these monthly resolved central Pacific coral

δ¹⁸O reconstructions²²⁻²⁶ have revealed a highly variable yet overall strengthening ENSO magnitude during the Holocene, but few studies have specifically examined variations in multi-year ENSO occurrence or ENSO duration. In this study, we use a synthesis of these reconstructions and the most recent climate model ensemble of transient Holocene simulations to investigate the changes in occurrence of multi-year ENSO over the much longer timespan of the past 7,000 years (7 ka; the Holocene) and relate these to the forcing mechanisms.

## Multi-year ENSO in reconstructions and simulations

Multi-year ENSO events are evident in scattered coral δ¹⁸O anomalies from the equatorial central Pacific (Methods), with the development and decay of each event showing strong similarities with the Niño 3.4 (5°N to 5°S, 170°W to 120°W) SST reanalysis data of 1854–2023 CE (Fig. 1). The anomalies of multi-year ENSO events start to decline after the first-year’s peak, reaching a trough after -6 months, but strengthen again and peak for a second year after -12 months. Compared with single-year ENSO, multi-year events are preceded by stronger opposing anomalies; for instance, a stronger El Niño is followed by a multi-year La Niña (Fig. 1a,b, blue line versus Fig. 1c,d, blue line). Over the 1,265-years-worth of monthly records spanning the Holocene, we count 87 and 108 multi-year El Niño and La Niña events and 307 and 260 single-year El Niño and La Niña events, respectively.
<figure><img src="image_5.png" /></figure>

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

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

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

**Fig. 2 | Multi-year ENSO in proxies and simulations during the past 7 ka.**
a–h, The evolution of the ratio of multi-year to single-year El Niño + La Niña $r_{\text{ENSO}}$ (a and e), El Niño $r_{\text{EN}}$ (b and f), La Niña $r_{\text{LN}}$ (c and g) and ENSO-dominant period (d and h) from proxy reconstructions (a–d) and model simulations (e–h). In a–d, each black circle represents a fossil coral slice, and the purple circles show

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

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

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

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

the values for ERSST reanalysis of 1854–2023 CE. The blue circles in a provide reference diameters of the circles (proportional to data length). The dashed blue lines are weighted (by data length) linear fits. In e–h, the black curves show the selected model ensemble mean and the grey shading shows the model ensemble spread of 1 s.d.

Individual fossil coral slices reveal that multi-year ENSO occurrence is common throughout the Holocene, with at least one multi-year event happening in almost every coral slice longer than 20 years (Fig. 2a and Extended Data Fig. 1a). Most records show a ratio of multi-year to single-year events below 1, suggesting that single-year events have dominated throughout the Holocene (Fig. 2a–c). Also, the multi-year La Niña events are generally more frequent than the multi-year El Niño events (Fig. 2b versus Fig. 2c). These features are consistent with the instrumental record (Fig. 2a–c, purple circles).

The multimodel ensemble of the Holocene transient simulations (Methods) reproduces considerable multi-year ENSO events during the Holocene (Fig. 2e–g). The multimodel ensemble usually provides a more reliable representation for long-term ENSO evolution, potentially cancelling out the intrinsic decadal or centennial variability of ENSO in individual simulations¹⁷,¹⁸. We set a criterion to select a subset of four simulations to exclude the models biased by too-regular ENSO cycles²⁷ or those that show limited performance on critical ENSO

metrics²⁸ considered unsuitable for this study (Methods). However, the multimodel simulations still tend to underestimate the occurrence of consecutive ENSO, and the simulated ratio of multi-year to single-year ENSO is about half the level suggested by the proxy data and observation (see more details in the Methods).

## The Holocene evolution of multi-year ENSO occurrence

Despite potential uncertainties due to data length or chronology (Supplementary Discussion), a robust linear increase in the ratio of multi-year to single-year ENSO events over the past 7,000 years emerges from the proxy reconstructions (Fig. 2a–c, dashed blue line), driven by an increasing (decreasing) occurrence of multi-year (single-year) ENSO (Extended Data Fig. 1). Modern ratios obtained from instrumental data (Fig. 2a–c, purple circles) are very close to those obtained in recent (long) fossil corals, which supports the reliability of the palaeoclimate reconstruction. However, in short coral records, the large internal
<figure><img src="image_13.png" /><figcaption>b Corr ($r_{\text{ENSO, ENSO amplitude}}$)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>c Correlation</figcaption></figure>

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

**Fig. 3 | Multi-year ENSO, ENSO period and amplitude relationship in proxies and selected simulations.** a,b, The relationship between the ratio of multi-year to single-year ENSO $r_{\text{ENSO}}$ on the x axis and the ENSO-dominant period (a) or ENSO amplitude (b) on the y axis. The dashed blue lines show the linear fit. The blue circles in a provide reference diameters of the circles (proportional to data length). The purple circle in a is for ERSST reanalysis of 1854–2023 CE, but it is not shown in b because it is not appropriate to compare the variability of SST and $\delta^{18}\text{O}$ directly. In a and b, the *P* values in the Pearson correlation (corr) are shown in the figure legend. c, A comparison of corr($r_{\text{ENSO, ENSO period}}$) (blue) and corr($r_{\text{ENSO, ENSO amplitude}}$) (red) in proxies and four selected simulations.

variability in ENSO introduces noise (outliers in Fig. 2b,c and Extended Data Fig. 1b,c), which can be minimized by considering El Niño and La Niña together as this increases the total number of events (Fig. 2a). Nonetheless, we estimate an increasing trend in the ratio of combined multi-year events to single-year events ($r_{\text{ENSO}}$), and those of multi-year to single-year El Niño ($r_{\text{EN}}$) and La Niña ($r_{\text{LN}}$) are all statistically significant (*P* < 0.05; samples weighted by data length; Methods). The ratio $r_{\text{ENSO}}$ increases from 0.09 during the earliest 1,000 years of the proxy compilation to 0.44 during the last millennium (compared with 0.42 in the instrumental period) by a factor of 5.

The four-model ensemble mean also shows an increasing proportion of multi-year events during the Holocene (Fig. 2e–g; *P* = 0.064, *P* = 0.827 and *P* < 0.05), which coincides qualitatively with the proxy reconstructions. The increase in the ratios is more modest, but the ratios themselves are also smaller. The uncertainties due to the sliding window size when analysing the simulations can be almost neglected when window size is larger than 50 years (Supplementary Discussion).

Notably, the proxies and model simulations consistently show that this increasing occurrence of consecutive ENSO events is contributed by more frequent multi-year La Niña as well as multi-year El Niño, with a greater contribution from the former (as seen in its larger slope of linear fit). The averaged $r_{\text{EN}}$ and $r_{\text{LN}}$ of all proxies are 0.28 and 0.42, respectively. These features are in agreement with the observation that multi-year La Niña cases are more common than multi-year El Niño$^{7,19}$ (Fig. 2b,c, purple circle).

## Prolonged ENSO cycles and their cause

The Holocene annual mean climate is relatively stable compared with the glacial period$^{29}$. Nonetheless, changes in the distribution of incoming solar radiation across the year (caused by orbital precession) still drove important regional climate variations$^{30,31}$. This seasonal insolation change caused a marked reduction in ENSO amplitude in the eastern and central Pacific$^{23,24,32,33}$ that has been linked to systematic changes in the upper-ocean structure and climate background in the tropical Pacific$^{21,34}$. They could further affect the ENSO periodicity.

The fact that multi-year El Niño and La Niña consistently became more frequent is associated with a prolonged ENSO period, as the anomalies remain in the warm or cold phase for longer. The ENSO-dominant period (Methods) is found to significantly (*P* < 0.05) increase in the proxy data by roughly 6 months (Fig. 2d). This is consistent with an increase in ENSO period in the model ensemble (although this is not statistically significant, *P* = 0.058). The correlation between the $r_{\text{ENSO}}$ and the ENSO period is 0.58 as estimated from the reconstructions (Fig. 3a) and verifies the relationship between consecutive ENSO occurrence and ENSO period.

We propose that extended periods of ENSO provide a broad explanation for more frequent multi-year El Niño as well as multi-year La Niña cases, with its root cause in the changes of upper-ocean structure (detailed below). This is based on the classic linear oscillatory theories of ENSO, rather than mechanisms tied to surface and/or nonlinear processes such as the cross-equatorial wind$^{35}$, North Pacific variability$^{12,36}$, trade wind anomalies$^{11}$ or teleconnections from other ocean basins$^{37}$. These theories are used to explain recent trends in multi-year La Niña or El Niño alone and tend to focus on anthropogenic influence within the past few decades$^{12,13}$ or the near future$^{11}$. We next demonstrate how the linear framework that encompasses slowly varying insolation forcing links the mean climate changes and prolonged ENSO periods.

We adopt the conceptual recharge oscillator (RO) model$^{38}$ (Methods) to investigate the relationship between the climate background, sensitivity in key ENSO feedbacks, and ENSO periods. By perturbing one parameter for a set of sensitivity experiments with RO, we elucidate the dependence of ENSO period on the thermocline depth (Methods and Extended Data Fig. 2). The RO sensitivity experiments show a series of steady-state ENSO conditions: when the zonal mean thermocline depth begins to decrease artificially in RO until 40% of it is removed (parameter $\lambda$ changes from 0 to 0.4), the period of ENSO increases from 3.82 years to 4.63 years, or prolonged by 21% (Extended Data Table 1).

The thermocline depth–ENSO period relationship revealed by RO is in line with the transient Holocene simulations. Three models have monthly subsurface ocean data available, and they all show a decreased vertical temperature gradient in the central and eastern equatorial Pacific from January to August and, thus, on an annual mean basis (Fig. 4c), implying an overall deeper and/or warmer thermocline during the mid-Holocene, consistent with a shorter ENSO period and fewer multi-year ENSO events. Although cooler surface temperature has been recorded in the easternmost Pacific$^{39,40}$, subsurface temperature estimates from the western equatorial Pacific$^{41}$, and eastern Pacific$^{42,43}$ suggest a warmer and/or deeper thermocline in the mid-Holocene, in agreement with the mid-Holocene (6 ka) multimodel time-slice ensemble$^{44}$ and with our transient experiments.

The deepened tropical Pacific thermocline weakens the thermocline feedback, because the wind stress anomalies are less effective
<figure><img src="image_16.png" /><figcaption>**Fig. 4 | Influence of orbital forcing on the upper ocean in the Eastern Pacific.**

a–c, The seasonal cycle of the difference between 6 to 5.5 ka and 1 to 0.5 ka of insolation (W m⁻²) (a), zonal mean SST of the Eastern Pacific (180–100° W) (b) and zonal mean subsurface ocean (at -150 m depth) temperature of the Eastern Pacific (180–100° W) (c). b shows the ensemble mean of four selected simulations, and c shows the ensemble mean of three selected simulations (LOVECLIM data not included due to availability).</figcaption></figure>

in reinforcing the thermocline anomalies that are linked by the mean upwelling⁵,⁴⁵. Towards the late Holocene, increased thermocline feedback can contribute to more efficient ocean–atmosphere coupling and reduced ENSO PB⁵ (Extended Data Fig. 3). The PB for modern ENSO is present typically during the Northern Hemisphere spring when SST anomaly (SSTA) sharply weakens. When PB becomes weaker, ENSO anomalies after the first-year’s peak are more likely to persist and lead to prolonged ENSO, regardless of El Niño or La Niña. Similarly in principle, a deeper thermocline and enhanced thermocline feedback are found to be associated with a pronounced onset rate of multi-year La Niña events following super El Niño events in the past century¹³.

Alternatively, stronger ENSO activities may also contribute to more multi-year ENSO. Multiple earlier studies suggest reduced ENSO amplitude throughout most of the Holocene compared with the twentieth century²¹,⁴⁶. In particular, ENSO activity in the tropical central Pacific appears to reach a minimum around 5 ka and gradually increases after that²⁴ (Extended Data Fig. 4). In theory, a stronger ENSO case creates larger ocean heat content variation lagging the surface temperature anomalies, and subsequently it can take longer time for an opposite condition to neutralize the preceding strong anomalies through the slow recharge–discharge ocean heat content process⁴⁷. This view is supported by the weak (*P* = 0.07) correlation of 0.39 between *r*<sub>ENSO</sub> and ENSO amplitude in proxies (Fig. 3b and Methods) and selected simulations (Fig. 3c), which suggests that the precession-forced intensifying of ENSO across the Holocene also partly leads to a favourable condition for more persistent ENSO. Nevertheless, a consensus among proxies and simulations (Fig. 3c) indicates that *r*<sub>ENSO</sub> is more closely linked to ENSO period rather than amplitude and underscores the primary role of ENSO period in driving multi-year ENSO variability.

## Orbital forcing of the equatorial Pacific thermocline

Finally, we investigate how orbital forcing drove the tropical Pacific thermocline structure changes, and the resultant prolonged ENSO cycles towards the late Holocene. The forcing mechanism for deeper and warmer thermocline in the mid/early Holocene was first explored in earlier modelling studies²¹,³³,⁴⁸. By comparing the subsurface sea temperature, we can identify the warmer temperature anomalies throughout most of the year in the Tropical Eastern Pacific (Fig. 4c and Extended Data Fig. 5). This change is in contrast to the SST anomalies in the tropics which are colder during the first half year but warmer during the second half year (Fig. 4b and Extended Data Fig. 6). The SST largely follows the local insolation forcing, but with a phase-lagged response of 1–2 months (Fig. 4a,b) due to large heat capacity of sea water. Conversely, the subsurface ocean warming can be explained by the subduction process. Increased solar radiation in the Southern

Hemisphere late winter/early spring warms the surface water in the subtropical South Pacific (Fig. 4b). This warm signal then propagates to the deeper Tropical Eastern Pacific thermocline through subduction when the local thermocline is the deepest and can be ventilated into the tropics⁴⁹. Eventually, this warmer water reaches the equatorial thermocline, heating the subsurface and deepening the thermocline throughout most of the year (Fig. 4c and Extended Data Fig. 5). With the subsurface warming, the surface cooling⁵,⁴⁰ further weakens the vertical gradient. In short, due to the ocean subduction, the precessional forcing in the subtropical South Pacific led to weaker ocean stratification and deeper thermocline depth in the tropical Pacific during the early/mid-Holocene epoch.

Using multiple approaches, we demonstrate an increasing occurrence of multi-year ENSO from the early/mid to late Holocene mainly under the orbital forcing. This increasing trend is robust in a synthesis of proxy reconstructions. It is also seen, albeit weaker, across the four climate models that adequately reproduce multi-year ENSO. This trend, caused by both increasingly frequent multi-year La Niña and multi-year El Niño events, is associated with a prolonged ENSO period. We attribute the lower frequency of multi-year events during the mid/early Holocene to a weaker vertical temperature gradient in the central and eastern equatorial Pacific, due to local surface cooling and subduction of warmer SST in the subtropical region of the Southern Hemisphere.

ENSO flavour change may interact with different processes of ENSO excitation and feedbacks¹² and lead to changes in ENSO duration. Proxy records from both the central Pacific and the eastern Pacific suggest that the present-day mix of flavours has occurred since 5 ka (ref. 24). However, constraining the shifts in the most active region of ENSO between the central Pacific and the eastern Pacific is challenging because of the scarcity of reconstructions and model inconsistency (Extended Data Fig. 7). These uncertainties highlight the need for further research with high-temporal-resolution datasets covering the Pacific basin to solidify our findings. Nonetheless, the reconstructions within the Niño 3.4 region are adequate to effectively record the history of multi-year ENSO (Supplementary Discussion), even if they cannot distinguish that of its separate flavours.

Historical observations and simulations have already demonstrated that the probability of multi-year La Niña and El Niño events has increased in the past decades¹¹⁻¹³. The research here shows that the frequency of multi-year ENSO events has been trending upwards for the past 7,000 years, hinting that the recent probability may be unprecedented since the early/mid Holocene, particularly as the slowly varying upper-ocean stratification background at orbital timescales has become more favourable for prolonged ENSO. This implies that climate models with more realistic CO₂-induced ocean response⁵⁰ can be key for robust projections of ENSO periodicity. The shallower
tropical Pacific thermocline, in addition to other emerging nonlinear and/or surface processes under anthropogenic forcings<sup>11,12</sup> (although not evident at the Holocene timescale; Extended Data Fig. 8), will probably drive increased prevalence of multi-year ENSO events in the twenty-first century and will combine with the stronger amplitude of ENSO<sup>24,46</sup> to have consequences across the globe. More climate mitigation efforts are urgently needed to tackle these impacts.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-025-01670-y.

## References

1. McPhaden, M. J., Santoso, A. & Cai, W. *El Niño Southern Oscillation in a Changing Climate*, Vol. 253 (Wiley, 2020).

2. Liu, Z. & Alexander, M. Atmospheric bridge, oceanic tunnel, and global climatic teleconnections. *Rev. Geophys*. https://doi.org/10.1029/2005RG000172 (2007).

3. Cai, W. et al. Changing El Niño-Southern Oscillation in a warming climate. *Nat. Rev. Earth Environ*. https://doi.org/10.1038/s43017-021-00199-z (2021).

4. Wu, M. et al. Regional responses of vegetation productivity to the two phases of ENSO. *Geophys. Res. Lett.* **51**, e2024GL108176 (2024).

5. Jin, Y., Liu, Z., Lu, Z. & He, C. Seasonal cycle of background in the tropical Pacific as a cause of ENSO Spring Persistence Barrier. *Geophys. Res. Lett.* **46**, 13371–13378 (2019).

6. Torrence, C. & Webster, P. J. The annual cycle of persistence in the El Niño/Southern Oscillation. *Q. J. R. Meteorol. Soc.* **124**, 1985–2004 (1998).

7. Okumura, Y. M. & Deser, C. Asymmetry in the duration of El Niño and La Niña. *J. Clim.* **23**, 5826–5843 (2010).

8. Cole, J. E., Overpeck, J. T. & Cook, E. R. Multiyear La Niña events and persistent drought in the contiguous United States. *Geophys. Res. Lett.* **29**, 25-1-25-4 (2002).

9. Lee, S. K. et al. Spring persistence, transition, and resurgence of El Niño. *Geophys. Res. Lett.* **41**, 8578–8585 (2014).

10. Jeong, H., Park, H.-S., Chowdary, J. S. & Xie, S.-P. Triple-dip La Niña contributes to Pakistan flooding and Southern China drought in summer 2022. *Bull. Am. Meteorol. Soc.* **104**, E1570–E1586 (2023).

11. Geng, T. et al. Increased occurrences of consecutive La Niña events under global warming. *Nature* **619**, 774–781 (2023).

12. Ding, R. et al. Multi-year El Niño events tied to the North Pacific Oscillation. *Nat. Commun.* **13**, 3871 (2022).

13. Wang, B. et al. Understanding the recent increase in multiyear La Niñas. *Nat. Clim. Change* **13**, 1075–1081 (2023).

14. Okumura, Y. M., DiNezio, P. & Deser, C. Evolving impacts of multiyear La Niña events on atmospheric circulation and U.S. drought. *Geophys. Res. Lett.* **44**, 11,614–11,623 (2017).

15. Sanchez, S. C. & Karnauskas, K. B. Diversity in the persistence of El Niño events over the last millennium. *Geophys. Res. Lett.* **48**, e2021GL093698 (2021).

16. Guilyardi, E., Capotondi, A., Lengaigne, M., Thual, S. & Wittenberg, A. T. in *El Niño Southern Oscillation in a Changing Climate* (eds McPhaden, M.J. et al.) 199–226 (Wiley, 2020).

17. Wittenberg, A. T. Are historical records sufficient to constrain ENSO simulations? *Geophys. Res. Lett.* https://doi.org/10.1029/2009GL038710 (2009).

18. Braconnot, P. et al. Impact of multiscale variability on last 6,000 years Indian and West African monsoon rain. *Geophys. Res. Lett.* **46**, 14021–14029 (2019).

19. Timmermann, A. et al. El Niño-Southern Oscillation complexity. *Nature* **559**, 535–545 (2018).

20. Lu, Z., Liu, Z., Zhu, J. & Cobb, K. M. A Review of Paleo El Niño-Southern Oscillation. *Atmosphere* **9**, 130 (2018).

21. Liu, Z. et al. Evolution and forcing mechanisms of El Niño over the past 21,000 years. *Nature* **515**, 550–553 (2014).

22. Cobb, K. M. et al. Highly variable El Niño-Southern Oscillation throughout the Holocene. *Science* **339**, 67–70 (2013).

23. Emile-Geay, J. et al. Links between tropical Pacific seasonal, interannual and orbital variability during the Holocene. *Nat. Geosci.* **9**, 168–173 (2016).

24. Carré, M. et al. High-resolution marine data and transient simulations support orbital forcing of ENSO amplitude since the mid-Holocene. *Quat. Sci. Rev.* **268**, 107125 (2021).

25. Lawman, A. E. et al. Unraveling forced responses of extreme El Niño variability over the Holocene. *Sci. Adv.* **8**, eabm4313 (2022).

26. Grothe, P. R. et al. Enhanced El Niño-Southern oscillation variability in recent decades. *Geophys. Res. Lett.* **47**, e2019GL083906 (2020).

27. Deser, C., Capotondi, A., Saravanan, R. & Phillips, A. S. Tropical Pacific and Atlantic climate variability in CCSM3. *J. Clim.* **19**, 2451–2481 (2006).

28. Planton, Y. Y. et al. Evaluating climate models with the CLIVAR 2020 ENSO metrics package. *Bull. Am. Meteorol. Soc.* **102**, E193–E217 (2021).

29. Kaufman, D. S. & Broadman, E. Revisiting the Holocene global temperature conundrum. *Nature* **614**, 425–435 (2023).

30. Claussen, M., Dallmeyer, A. & Bader, J. Theory and modeling of the African humid period and the green Sahara. *Oxford Research Encyclopedia of Climate Science* https://doi.org/10.1093/acrefore/9780190228620.013.532 (2017)

31. Chen, J., Zhang, Q., Kjellström, E., Lu, Z. & Chen, F. The contribution of vegetation-climate feedback and resultant sea ice loss to amplified Arctic warming during the mid-Holocene. *Geophys. Res. Lett.* **49**, e2022GL098816 (2022).

32. Carré, M. et al. Holocene history of ENSO variance and asymmetry in the eastern tropical Pacific. *Science* **345**, 1045–1048 (2014).

33. Luan, Y., Braconnot, P., Yu, Y., Zheng, W. & Marti, O. Early and mid-Holocene climate in the tropical Pacific: seasonal cycle and interannual variability induced by insolation changes. *Climate* **8**, 1093–1108 (2012).

34. Lu, Z., Liu, Z., Chen, G. & Guan, J. Prominent precession band variance in ENSO intensity over the last 300,000 years. *Geophys. Res. Lett.* **46**, 9786–9795 (2019).

35. Lee, C. W., Tseng, Y. H., Sui, C. H., Zheng, F. & Wu, E. T. Characteristics of the prolonged El Niño events during 1960–2020. *Geophys. Res. Lett.* **47**, e2020GL088345 (2020).

36. Kim, J.-W., Yu, J.-Y. & Tian, B. Overemphasized role of preceding strong El Niño in generating multi-year La Niña events. *Nat. Commun.* **14**, 6790 (2023).

37. Wu, X., Okumura, Y. M. & DiNezio, P. N. What controls the duration of El Niño and La Niña events? *J. Clim.* **32**, 5941–5965 (2019).

38. Jin, F.-F. An equatorial ocean recharge paradigm for ENSO. Part I: conceptual model. *J. Atmos. Sci.* **54**, 811–829 (1997).

39. Koutavas, A., Lynch-Stieglitz, J., Marchitto, T. M. Jr & Sachs, J. P. El Niño-like pattern in ice age tropical Pacific sea surface temperature. *Science* **297**, 226–230 (2002).

40. Carré, M. et al. Mid-Holocene mean climate in the south eastern Pacific and its influence on South America. *Quat. Int.* **253**, 55–66 (2012).

41. Dang, H. et al. Pacific warm pool subsurface heat sequestration modulated Walker circulation and ENSO activity during the Holocene. *Sci. Adv.* **6**, eabc0402 (2020).
42. White, S. M., Ravelo, A. C. & Polissar, P. J. Dampened El Niño in the early and mid-Holocene due to insolation-forced warming/ deepening of the thermocline. *Geophys. Res. Lett.* **45**, 316–326 (2018).

43. Sadekov, A. Y. et al. Palaeoclimate reconstructions reveal a strong link between El Niño-Southern Oscillation and Tropical Pacific mean state. *Nat. Commun.* **4**, 2692 (2013).

44. Chen, L., Zheng, W. & Braconnot, P. Towards understanding the suppressed ENSO activity during mid-Holocene in PMIP2 and PMIP3 simulations. *Clim. Dyn.* **53**, 1095–1110 (2019).

45. Iwakiri, T. & Watanabe, M. Mechanisms linking multi-year La Niña with preceding strong El Niño. *Sci. Rep.* **11**, 17465 (2021).

46. Lu, Z. & Liu, Z. Examining El Niño in the Holocene: implications and challenges. *Natl Sci. Rev.* **5**, 807–809 (2018).

47. Wyrtki, K. Water displacements in the Pacific and the genesis of El Niño cycles. *J. Geophys. Res. Oceans* **90**, 7129–7132 (1985).

48. Saint-Lu, M., Braconnot, P., Leloup, J. & Marti, O. The role of El Niño in the global energy redistribution: a case study in the mid-Holocene. *Clim. Dyn.* **52**, 7135–7152 (2019).

49. Stommel, H. Determination of water mass properties of water pumped down from the Ekman layer to the geostrophic flow below. *Proc. Natl Acad. Sci. USA* **76**, 3051–3055 (1979).

50. Lee, S. et al. On the future zonal contrasts of equatorial Pacific climate: perspectives from observations, simulations, and theories. *npj Clim. Atmos. Sci.* **5**, 82 (2022).

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

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

© The Author(s) 2025

<footnote n="1">¹Department of Physical Geography and Ecosystem Science, Lund University, Lund, Sweden. ²IPSL/Laboratoire d'Océanographie et du Climat: Expérimentations et Approches Numériques (CNRS-IRD-MNHN-Sorbonne Universités), Paris, France. ³CIDIS, Facultad de Ciencias y Filosofía, Universidad Peruana Cayetano Heredia, Lima, Peru. ⁴Department of Geography, University College London, London, UK. ⁵School of Geography, Earth and Environmental Sciences, University of Birmingham, Birmingham, UK. ⁶Key Laboratory of Marine Geology and Environment, Institute of Oceanology, Chinese Academy of Sciences, Qingdao, China. ⁷Department of Geology, Lund University, Lund, Sweden. ⁸Institute for Atmospheric and Climate Science, ETH Zürich, Zurich, Switzerland. ⁹Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. ¹⁰Earth and Climate Research Center, Earth and Life Institute, Université catholique de Louvain, Louvain-la-Neuve, Belgium. ¹¹Max Planck Institut für Meteorologie, Hamburg, Germany. ¹²Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai, China.</footnote>

<footnote n="13">¹³Department of Atmospheric and Oceanic Sciences, Fudan University, Shanghai, China. ¹⁴Department of Physical Geography and Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden. ✉e-mail: zhengyao.lu@nateko.lu.se</footnote>
# Methods

## Proxy reconstructions

We synthesize the monthly resolved fossil coral oxygen isotope data from Kiritimati and Fanning Atolls in the equatorial central Pacific from several previous studies<sup>24–26</sup>. The Atolls are located within or near the Niño 3.4 region (5° S to 5° N, 170° W to 120° W, an ENSO-active region) (Fig. 1e) and collectively have over 1,200 years of coral proxy data spanning the last -7,000 years of the Holocene. We use only the coral slice records longer than 20 years (22 out 38 records selected; Supplementary Table 1), while the mean (median) length of them is 57.5 (38.5) years. The oxygen isotopic composition of a coral’s carbonate skeleton is a proxy that jointly reflects SST and the oxygen isotopic composition of seawater, which is linearly related to sea-surface salinity. On interannual timescales, El Niño (La Niña) events bring warmer and wetter (cooler and drier) conditions to Kiritimati and Fanning, which collectively yield negative (positive) coral δ<sup>18</sup>O anomalies. Previous work (for example, Cobb et al.<sup>22</sup>) demonstrates the high correlation between modern coral δ<sup>18</sup>O at these sites and the Niño 3.4 SST index (a practical ENSO index). Bivalve δ<sup>18</sup>O records from the equatorial eastern Pacific<sup>32</sup>, although monthly resolved, are not included in our synthesis owing to the short length of data.

## Models and simulations

We analyse the Holocene ENSO evolution in the Holocene transient simulation intercomparison contributed by eight research groups worldwide (Supplementary Table 2), which provides an opportunity to examine the continuous climate variations to orbitally induced changes in multiple model realizations. These simulations (except TRACE; more information on TRACE is given in Supplementary Table 2) were spun up from an early/mid Holocene initial condition (range from 10 ka to 6 ka, equilibrium or transient) and continued to run till the present day, driven primarily by the forcing of seasonal and latitudinal distribution of insolation (for example, Fig. 4a). Basic monthly climate output from the simulations, such as surface air temperature, SST are regridded to 1° by 1° for analysis.

The simulated multi-year versus single-year ENSO ratio ($r_{\text{ENSO}}$) using Niño 3.4 region SSTA (or surface air temperature anomalies used as SSTA owing to data availability) was calculated for each simulation to evaluate the model performance of multi-year ENSO events (Supplementary Fig. 1). We then use a criterion of $r_{\text{ENSO}}$ averaged during the last 1 ka, given $r_{\text{ENSO}} = 0.44$ in the fossil coral data of the last 1 ka and 0.42 in Extended Reconstructed Sea Surface Temperature (ERSST) (1854–2023 CE)<sup>51</sup>. Four out of eight models were selected with the best match of $r_{\text{ENSO}}$ for further analysis and calculating ensemble mean, namely EC-Earth3-LR<sup>52</sup>, IPSL<sup>53</sup>, MPI-ESM<sup>54,55</sup> and LOVECLIM1.3 (ref. 56). Other models simulating too regular and short ENSO cycles and too few multi-year ENSO events were excluded from the main results. For instance, TRACE (CCSM3 model) simulation is characterized by notable quasi-biannual ENSO, almost unable to reproduce any consecutive ENSO (Supplementary Fig. 1a, TRACE). These selected models can realistically represent observed SSTA evolution of single-year and multi-year El Niño and La Niña events (Supplementary Fig. 2). They also demonstrate the highest capability in realistically capturing some key ENSO assessment metrics (ENSO period, ENSO skewness and SST zonal structure)<sup>28</sup> (Supplementary Fig. 1).

Nevertheless, we notice an overall underestimation of multi-year ENSO occurrence, but the sources of model biases need further investigation. The common ENSO periodicity biases in climate models can be model dependent and come from biases in the background climatology<sup>16,57</sup>, from improper combined ENSO feedback processes<sup>16</sup> or, partly, from coarse model resolution<sup>27</sup>.

## Definition of multi-year ENSO

Due to different amplitudes of ENSO simulated by models, we do not use a fixed threshold of specific Niño 3.4 SSTA but 0.5 standard

deviation (s.d.)<sup>11</sup> which is calculated every 100 years (comparable to the length of a long coral data, and for more stable results; Supplementary Discussion). For every 100 years of monthly Niño 3.4 SSTA, after applying for a running mean filter to retrieve interannual variability of 1.5–7 years, and a 3-month smoothing<sup>11</sup>, we first identify the peak ENSO seasons. For instance, if the ENSO peak is in December (ENSO phase locking<sup>58</sup>), then the ENSO season is from -2 to +2 months (October (0) to February (1)). Then, if the SST is 0.5 s.d. above/below average for a whole ENSO season, it is categorized as an El Niño/La Niña event. Multi-year events are identified when the El Niño (or La Niña) criterion is met for 2 years in a row during the whole ENSO season, in contrast to neutral conditions or an opposite phase in the second year. Note that in this way 3-year events are simplified as two 2-year events. The calculation of multi-year ENSO occurrence for coral δ<sup>18</sup>O anomalies is similar but using the s.d. of each coral slice depending on its length (not a fixed 100-year window) (Fig. 1a,c). It should be pointed out that the ENSO peak month in coral data could be highly variable and did not always occur in December throughout the Holocene (Supplementary Fig. 3). We also tested different thresholds of ENSO (for example, 0.4× s.d. or 0.6× s.d.) for the coral data, and the results are not affected (figure not shown).

## ENSO-dominant period (spectrum analysis) and amplitude

ENSO periods in proxy reconstructions and model simulations (Niño 3.4 SSTA) are estimated by spectrum analysis. We first apply a running mean filter to exclude the variability longer than 7 years or shorter than 1.5 years in the time series of δ<sup>18</sup>O anomalies and SSTA. The rest ENSO-band time series are then processed by power spectrum analysis and a calculation of its s.d. for ENSO amplitude. Due to the irregularity of ENSO, the power spectra show high values of variability at a broad interannual band with multiple peaks. Rather than the period associated with maximum power, we instead use the weighted-averaged (by power) period within the band of 1.5–7 years to determine the ENSO-dominant period (as a centroid).

## Significance test for the linear trend

We apply the commonly used Mann–Kendall (MK) statistical test<sup>59,60</sup> for the significance test of linear trends in time series of the occurrence of multi-year ENSO and ENSO periods (Fig. 2, dashed blue lines). The MK is rank-based and non-parametric and has the advantage that the input variable does not have to be normally distributed. The null hypothesis (H0) of the two-tailed MK suggests that there is no change of trend, while the alternative hypothesis (H1) implies a positive or negative trend over time. The H0 was rejected when the P value fell below 0.05. The weight of sample data (proportional to length) is considered when calculating the linear trend of coral slice data and its significance.

## The conceptual RO

We use the classic conceptual RO model<sup>38</sup> to verify the relationship between the depth of the mean thermocline in the tropical Pacific and the period of ENSO. It is based on coupling SST evolution to ocean adjustment processes to provide a simplified framework with which to understand the oceanic and atmospheric roles in driving variations in ENSO cycles. RO does not incorporate multi-year events, but it allows for changes in the period of ENSO. We modify the model to give us the ability to introduce changes in the zonal mean thermocline depth in the equatorial Pacific. Compared with complex GCMs, this approach allowed us to explicitly isolate and examine the specific impact of thermocline depth on ENSO periodicity.

The RO model can be expressed as

$$ 
\begin{aligned}
\frac{dh_W}{dt} &= -rh_W - \alpha b T_E \\
\frac{dT_E}{dt} &= RT_E + \gamma h_W - e_n(h_W + b T_E)^3.
\end{aligned}
\qquad (1) $$
These two equations describe the subsurface ocean adjustment dynamics and the SST dynamics, respectively. In equation (1), $h_W$ denotes the thermocline depth anomaly in the tropical Western Pacific, and $T_E$ denotes the SST anomaly in the Tropical Central to Eastern Pacific; $r$ describes the collective damping rate of the upper ocean through mixing and the equatorial energy loss to the west and east boundary layer currents; $\alpha$ and $b$ are coupling coefficients of wind stress and Sverdrup transport, and wind stress and SST anomalies, respectively; $R$ is a collective Bjerknes feedback strength, and $\gamma$ is the thermocline feedback coefficient; the cubic term in the SST equation introduces nonlinear effects to the coupled system with $e_n$ to quantify its strength.

To study the sensitivity of the oscillation system to zonal mean equatorial thermocline depth (denoted $\dot{h}$ and $\dot{h} = (h_W + h_E)/2$), we update the RO model with variable $\dot{h}$ (anomaly)$^{35}$. Due to the quasi-balance between the thermocline depth tilt (and, thus, the pressure gradient force) along the Equator and the wind stress over the equatorial band ($\tau$), there is a relationship $h_E - h_W = \tau = \mu T_E$ (refs. 61,62), in which $\mu$ is the wind stress coupling coefficient. Also, $b = b_0\mu$, where $b$ is the high-end estimation of $b_0$. If we introduce a parameter $\lambda$ for the change in the zonal mean thermocline anomaly and replace $h_W$ by $h_W - \lambda\dot{h}$ in the SST equation (the second) in equation (1), so the original nonlinear RO system is updated as

$$ 
\begin{aligned}
\frac{dh_W}{dt} &= -rh_W - \alpha b T_E \\
\frac{dT_E}{dt} &= (R - \gamma b_0 \mu \lambda / 2) T_E + \gamma(1 - \lambda)h_W - e_n(h_W + bT_E)^3.
\end{aligned}
\quad (2) $$

Other coefficients and parameters in this RO with variable thermocline depth are given values according to the original study$^{38}$ as $r = 0.25$, $\alpha = 0.125$, $b = 2.22$, $b_0 = 2.5$, $R = 0.667$, $\gamma = 0.75$, $\mu = 8/9$ and $e_n = 3$. We then conduct a series of sensitivity experiments focusing on different levels of $\lambda$ (Extended Data Fig. 2 and Supplementary Table 1). Note that, when $\lambda$ equals 1, the $h_W$ term is removed from the SST equation, which becomes decoupled from the ocean dynamics equation; when $\lambda$ equals 0, the system is not changed. The relationship between the thermocline depth and ENSO period is still robust, although slightly weakened, if we introduce a seasonal cycle to $\lambda$ (Supplementary Table 1).

## Data availability

The Niño 3.4 index data (NetCDF files) from all eight simulations are available via Zenodo at https://doi.org/10.5281/zenodo.14007727 (ref. 63). All fossil coral data and metadata are compiled and archived by Grothe et al. at the NCDC at https://www.ncdc.noaa.gov/paleo/study/22415. The ERSST (anomaly) reanalysis is available at https://iridl.Ideo.columbia.edu/SOURCES/.NOAA./NCDC.ERSST./version5./anom/. The TRACE simulation is available at https://www.earthsystem-grid.org/project/trace.html.

## Code availability

The MATLAB code for generating Figs. 1–4 is available via Zenodo at https://doi.org/10.5281/zenodo.14007727 (ref. 63).

## References

51. Huang, B. et al. Extended reconstructed sea surface temperature, version 5 (ERSSTv5): upgrades, validations, and intercomparisons. *J. Clim.* **30**, 8179–8205.

52. Zhang, Q. et al. Simulating the mid-Holocene, last interglacial and mid-Pliocene climate with EC-Earth3-LR. *Geosci. Model Dev.* **14**, 1147–1169 (2021).

53. Braconnot, P., Zhu, D., Marti, O. & Servonnat, J. Strengths and challenges for transient mid- to late Holocene simulations with dynamical vegetation. *Climate* **15**, 997–1024 (2019).

54. Bader, J. et al. Global temperature modes shed light on the Holocene temperature conundrum. *Nat. Commun.* **11**, 4726 (2020).

55. Dallmeyer, A. et al. Holocene vegetation transitions and their climatic drivers in MPI-ESM1.2. *Climate* **17**, 2481–2513 (2021).

56. Yin, Q. Z., Wu, Z. P., Berger, A., Goosse, H. & Hodell, D. Insolation triggered abrupt weakening of Atlantic circulation at the end of interglacials. *Science* **373**, 1035–1040 (2021).

57. Guilyardi, E. El Niño–mean state–seasonal cycle interactions in a multi-model ensemble. *Clim. Dyn.* **26**, 329–348 (2006).

58. Neelin, J. D., Jin, F.-F. & Syu, H.-H. Variations in ENSO phase locking. *J. Clim.* **13**, 2570–2590 (2000).

59. Mann, H. B. Nonparametric tests against trend. *Econometrica* **13**, 245–259 (1945).

60. Kendall, M. G. *Rank Correlation Methods* (Griffin, 1948).

61. Schneider, E. K., Huang, B. & Shukla, J. Ocean wave dynamics and El Niño. *J. Clim.* **8**, 2415–2439 (1995).

62. Deser, C. & Wallace, J. M. Large-scale atmospheric circulation features of warm and cold episodes in the tropical Pacific. *J. Clim.* **3**, 1254–1281 (1990).

63. Lu, Z. Data and analysis code for “increased frequency of multi-year El Niño-Southern Oscillation events across the Holocene”. *Zenodo* https://doi.org/10.5281/zenodo.14007727 (2024).

64. Jin, Y., Lu, Z. & Liu, Z. Controls of spring persistence barrier strength in different ENSO regimes and implications for 21st century changes. *Geophys. Res. Lett.* **47**, e2020GL088010 (2020).

## Acknowledgements

Z.L. and A.S. acknowledge the Swedish Research Council Vetenskapsrådet (grant no. 2022-03617), the Swedish Research Council Formas (grant no. 2020-02267), and the Strategic Research Area MERGE (Modeling the Regional and Global Earth System). Q.Z. is funded by the Swedish Research Council Vetenskapsrådet (grant nos. 2017-04232 and 2022-03129). X.S. is supported by the Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai) (grant no. SML2023SP204), the National Natural Science Foundation of China (NSFC) (grant no. 42206256) and the Ocean Negative Carbon Emissions (ONCE) Program. The EC-Earth simulations were performed on resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at Linköping University and the ECMWF’s computing and archive facilities. The LOVECLIM1.3 simulations were performed with the support of the Fonds de la Recherche Scientifique-FNRS research grant no. T.0246.23. The IPSL simulation was performed as part of the GENCI computing project gen2212 and gen12006 on the TGCC computing centre. They were supported by JPI-Belmont project PACMEDY (ANR-15-JCLI-0003-01). The HadCM3 simulations were run on the University of Birmingham’s BlueBEAR HPC service (http://www.birmingham.ac.uk/bear) and supported through a Birmingham Fellowship to P.O.H. Z.L. thanks W. Cai, M. McPhaden and Z. Liu for their helpful comments.

## Author contributions

Z.L. conceived this study. Z.L. analysed the data with input from A.S. and Q.Z. Holocene transient simulations were performed by P.B., P.O.H., H.Y., J.H.J., X.S., Q.Y. and Q.Z. All authors interpreted and discussed the results. Z.L. and A.S. wrote the paper with contributions from all authors. All authors reviewed the manuscript.

## Funding

Open access funding provided by Lund University.
**Competing interests**

The authors declare no competing interests.

**Additional information**

**Extended data** is available for this paper at
https://doi.org/10.1038/s41561-025-01670-y.

**Supplementary information** The online version
contains supplementary material available at
https://doi.org/10.1038/s41561-025-01670-y.

**Correspondence and requests for materials** should be addressed to
Zhengyao Lu.

**Peer review information** *Nature Geoscience* thanks Giovanni Liguori, Weiyi Sun and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: James Super, in collaboration with the *Nature Geoscience* team.

**Reprints and permissions information** is available at
www.nature.com/reprints.
Extended Data Table 1 | Simulated ENSO dominant periods in the RO

<table><thead><tr><th>r=0.25, α=0.125, b=2.22, b_0=2.5, R=0.667, γ=0.75, μ=8/9, and e_n=3</th><th>λ</th><th>Period (years)</th><th>λ *(1+sin(2πt/T))</th><th>Period (years)</th></tr></thead><tbody><tr><td></td><td>0</td><td>3.82</td><td>0</td><td>3.82</td></tr><tr><td></td><td>0.1</td><td>3.93</td><td>0.1</td><td>4.01</td></tr><tr><td></td><td>0.2</td><td>4.08</td><td>0.2</td><td>4.01</td></tr><tr><td></td><td>0.3</td><td>4.29</td><td>0.3</td><td>4.21</td></tr><tr><td></td><td>0.4</td><td>4.63</td><td>0.4</td><td>4.48</td></tr><tr><td></td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

The simulated ENSO dominant periods in the RO with variable thermocline depth (modulated by λ). A larger λ means a shallower thermocline depth. The seasonally-varying λ can approximate deeper thermocline particularly during the first half of year in the early/mid Holocene due to orbital forcing (Extended Data Fig. 5). The results show that the deeper thermocline-longer ENSO period relationship generally exists irrespective of whether λ has a seasonal cycle or not.
<figure><img src="image_17.png" /><figcaption>**Extended Data Fig. 1** | Multi-year ENSO cases in proxies. Same as Fig. 2, but for the occurrence (number per 100 years) of (a)–(c) multi-year and (d)–(f) single year El Nino + La Nina, El Nino and La Nina.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 2 | ENSO cycles in two idealised RO simulations.** The simulated $T_E$ and $H_W$ in the RO model (the last 30-year of 500-year simulations are shown).

(a) and (b) are for different $\lambda$ values (without seasonal cycle), and the thermocline depth is shallower in (b). Simulated ENSO periods in (b) are extended.
<figure><img src="image_19.png" /></figure>

**Extended Data Fig. 3** | ENSO persistence barrier in proxies and selected simulations. (a) The evolution of ENSO persistence barrier (PB) strength from proxy reconstructions. The relationship between *r*<sub>ENSO</sub> on x-axis, and PB strength on y-axis. The blue line shows the linear fit. Blue circles in (b) provides reference diameters of the circles (proportional to data length), and the p value in Pearson correlation is shown in the legend. (c) The correlation between *r*<sub>ENSO</sub> and PB strength in proxies and 4 selected simulations over the Holocene. PB is derived

from the autocorrelation function (ACF) of Niño 3.4 SST anomalies which is a function of initial months *t* and lag months *τ*. Following Jin et al. (2020)⁶⁴, for a calendar month *t*, we identify *τ*<sub>B</sub> (*t*) as the specific lag of maximum ACF decline, which is calculated as the lag gradient in the time step of 1 month as PB(*t*) = {[*t*<sub>*t*,τ*<sub>B</sub>(*t*)</sub>−1]−*r*<sub>*t*,τ*<sub>B</sub>(*t*)+1</sub>} / 2 = max<sub>*t*</sub> {[*t*<sub>*t*,*τ*<sub>B</sub>(*t*)</sub>−1]−*r*<sub>*t*,*τ*<sub>B</sub>(*t*)+1</sub>} where PB(*t*) is the maximum gradient for each month. The intensity of the PB is then estimated using the sum of monthly PB(*t*) as PB = Σ<sub>*t*=1</sub><sup>12</sup> PB(*t*).
<figure><img src="image_20.png" /><figcaption>**Extended Data Fig. 4** | **The reconstructed Holocene evolution of ENSO amplitude.** The Holocene evolution of ENSO amplitude from the reconstruction synthesis, as indicated by the standard deviation of δ¹⁸O anomalies (subtracting seasonal cycle) from each Central Pacific coral slice longer than 20-years. The blue dashed line is weighted (by data length) linear fit.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Extended Data Fig. 5** | Influence of orbital forcing on the Eastern Pacific upper ocean. (a)–(d) Same as Fig. 4(c) but for all three model simulations. (e)–(g) the seasonal cycle of vertical temperature gradient (surface-subsurface) differences between 6ka and 1ka for 5°N–5°S, 180°W–100°W. The blue dashed lines are for the annual mean.</figcaption></figure>
<figure><img src="image_22.png" /></figure>

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

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

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

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

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

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

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

Extended Data Fig. 6 | Influence of orbital forcing on the Eastern Pacific SST. Same as Fig. 4(b) but for all 8 model simulations.
<figure><img src="image_30.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_37.png" /><figcaption>**Extended Data Fig. 7** | Simulated ENSO spatial pattern in the Holocene. Hovmöller diagrams of the equatorial Pacific (5° N-5° S) SSTA interannual std in all 8 simulations. The Nino 3.4 region (5° N-5° S, 170° W-120° W) is marked by the black line near the bottom of each panel.</figcaption></figure>
<figure><img src="image_38.png" /></figure>

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

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

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

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

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

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

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

**Extended Data Fig. 8** | Simulated trend of annual mean temperature change during the last 7 ka. The trend of annual mean temperature change during the last 7 ka in all Holocene transient simulations. These long-term changes can be compared to the warming spatial pattern over the past century or in the future projections, for example, stronger warming in the western tropical Pacific (Wang et al.¹³), or in the subtropical northeast Pacific (Geng et al.¹¹; see their Fig. 4a). The

general warming patterns in Holocene simulations show the maximum warming region varies widely across individual simulations. Particularly, the stronger warming in the subtropical northeast Pacific or warm pool region is not evident in selected simulations. These results suggest at the Holocene timescale, the key mechanisms for multi-year ENSO enhancement may differ from it under the recent/future anthropogenic warming conditions.