# Mid-Pliocene El Niño/Southern Oscillation suppressed by Pacific intertropical convergence zone shift

Gabriel Pontes, Andréa Taschetto, Alex Sen Gupta, Agus Santoso, Ilana

Wainer, Alan Haywood, Wing-Le Chan, Ayako Abe-Ouchi, Christian

Stepanek, Gerrit Lohmann, et al.

► To cite this version:

Gabriel Pontes, Andréa Taschetto, Alex Sen Gupta, Agus Santoso, Ilana Wainer, et al. Mid-Pliocene El Niño/Southern Oscillation suppressed by Pacific intertropical convergence zone shift. Nature Geoscience, 2022, 10.1038/s41561-022-00999-y. hal-03763272

HAL Id: hal-03763272

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

Submitted on 7 Oct 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>

# Mid-Pliocene El Niño/Southern Oscillation suppressed by Pacific intertropical convergence zone shift

Gabriel M. Pontes<fnref n="1" />✉, Andréa S. Taschetto²,³, Alex Sen Gupta²,³, Agus Santoso²,³,<fnref n="4" />, Ilana Wainer<fnref n="1" />, Alan M. Haywood⁵, Wing-Le Chan⁶, Ayako Abe-Ouchi⁶, Christian Stepanek<fnref n="1" />, Gerrit Lohmann⁷, Stephen J. Hunter⁵, Julia C. Tindall⁵, Mark A. Chandler⁸, Linda E. Sohl⁸, W. Richard Peltier⁹, Deepak Chandan⁹, Youichi Kamae¹⁰, Kerim H. Nisancioglu¹¹, Zhongshi Zhang¹¹, Camille Contoux¹², Ning Tan¹²,¹³, Qiong Zhang¹⁴, Bette L. Otto-Bliesner¹⁵, Esther C. Brady¹⁵, Ran Feng<fnref n="16" />, Anna S. von der Heydt¹⁷, Michiel L. J. Baatsen¹⁷ and Arthur M. Oldeman¹⁷

The El Niño/Southern Oscillation (ENSO), the dominant driver of year-to-year climate variability in the equatorial Pacific Ocean, impacts climate pattern across the globe. However, the response of the ENSO system to past and potential future temperature increases is not fully understood. Here we investigate ENSO variability in the warmer climate of the mid-Pliocene (~3.0–3.3 Ma), when surface temperatures were ~2–3 °C above modern values, in a large ensemble of climate models—the Pliocene Model Intercomparison Project. We show that the ensemble consistently suggests a weakening of ENSO variability, with a mean reduction of 25% (±16%). We further show that shifts in the equatorial Pacific mean state cannot fully explain these changes. Instead, ENSO was suppressed by a series of off-equatorial processes triggered by a northward displacement of the Pacific intertropical convergence zone: weakened convective feedback and intensified Southern Hemisphere circulation, which inhibit various processes that initiate ENSO. The connection between the climatological intertropical convergence zone position and ENSO we find in the past is expected to operate in our warming world with important ramifications for ENSO variability.

El Niño/Southern Oscillation (ENSO) warm (El Niño) and cold (La Niña) events cause substantial changes in weather patterns and ocean circulation, impacting agriculture, fisheries, coral bleaching and cyclogenesis, among a host of other impacts<fnref n="1" />. Given its pronounced socioeconomic impacts and potential predictability of a few seasons in advance, ENSO has been under intense investigation². Whether and how ENSO changes in response to greenhouse gases and other external forcing may be studied by investigating past, present and future climates with palaeoreconstructions, instrumental records, theory and numerical simulations. There is a lack of consensus among climate models in general as to how ENSO variability will respond to future warming³,<fnref n="4" />, although models that better capture ENSO nonlinearity tend to simulate enhanced variability in the eastern equatorial Pacific⁵ and increased frequency of extreme events⁶,⁷. These changes in ENSO properties are linked to changes in the Pacific mean state marked by a weakened Walker

circulation, increased upper-ocean stratification, reduced zonal sea surface temperature (SST) gradient and equatorially enhanced warming that causes the intertropical convergence zone (ITCZ) to be displaced equatorward⁷⁻⁹.

Studies based on palaeoreconstructions have also suggested that ENSO activity is sensitive to the mean climate. A synthesis of mid-Holocene (~6 thousand years ago (ka)) records indicates a 33% reduction in ENSO amplitude in the eastern Pacific during this period¹⁰. ENSO activity over the past millennium was shown to be weaker when compared with the past half century¹¹, potentially suggesting global warming-related changes. Furthermore, there is evidence of significantly reduced ENSO variability during the Last Glacial Maximum¹² (~21 ka). Proxy data for the Pliocene (~5 to ~3 million years ago (Ma)) are controversial with regard to tropical Pacific changes¹³⁻<fnref n="16" />. A Pliocene El Niño-like mean state has been hypothesized to reduce ENSO variability¹⁷, although there is

<footnote n="1">¹Oceanographic Institute, University of São Paulo, São Paulo, Brazil. ²Climate Change Research Centre, The University of New South Wales, Sydney, New South Wales, Australia. ³ARC Centre of Excellence for Climate Extremes, The University of New South Wales, Sydney, New South Wales, Australia.</footnote>

<footnote n="4">⁴Centre for Southern Hemisphere Oceans Research CSHOR, CSIRO Oceans and Atmosphere, Hobart, Tasmania, Australia. ⁵School of Earth and Environment, University of Leeds, Leeds, UK. ⁶Centre for Earth System Dynamics, Atmosphere and Ocean Research Institute, University of Tokyo, Tokyo, Japan. ⁷Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany. ⁸NASA/GISS, Columbia CCSR, Columbia University, New York, NY, USA. ⁹Department of Physics, University of Toronto, Toronto, Ontario, Canada. ¹⁰Faculty of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Japan. ¹¹NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, Bergen, Norway. ¹²Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ Université Paris-Saclay, Gif-sur-Yvette, France. ¹³Key Laboratory of Cenozoic Geology and Environment, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China. ¹⁴Department of Physical Geography and Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden. ¹⁵National Center for Atmospheric Research, Boulder, CO, USA.</footnote>

<footnote n="16">¹⁶Department of Geosciences, College of Liberal Arts and Sciences, University of Connecticut, Storrs, CT, USA. ¹⁷Institute for Marine and Atmospheric research Utrecht (IMAU), Department of Physics, Utrecht University, Utrecht, the Netherlands. ¹⁸e-mail: gabrielpontes@usp.br</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Simulated mid-Pliocene tropical variability changes. **a**, Multimodel mean change in the amplitude (standard deviation) of SST anomaly (SSTa) variability in the PlioMIP2 models (see Extended Data Fig. 1 for PlioMIP1 models). SST anomaly is obtained through removing the mean seasonal cycle. Stippling indicates locations where there is significant model agreement (at least 70%) in the sign of the change. **b**, Change in the amplitude (standard deviation) of the Niño3 time series in each PlioMIP model. Red box in the eastern Pacific in **a** indicates the Niño3 region. Map created using the Basemap library for Python.</figcaption></figure>

evidence of significant interannual variability during this period¹⁸,¹⁹, whose magnitude could be comparable to the late Holocene²⁰. As such, tropical Pacific mean-state changes during the Pliocene and how they have impacted ENSO activity remain uncertain.

Palaeoclimate studies have suggested that the mid-Pliocene Warm Period (mPWP; ~3.3 Ma) might be a useful analogue to the end-of-century climate because of the warming magnitude²¹⁻²³. The mPWP was marked by warmer SSTs of up to 9°C and 4°C in the Northern Hemisphere and Southern Hemisphere, respectively, compared with pre-industrial times²² (~1850), with orbital forcing and elevated atmospheric CO₂ concentrations similar to present day (~400 ppm) while polar ice was reduced²³. Partly motivated by the similarities between the mPWP and scenarios of future projected warming, the Pliocene Model Intercomparison Project (PlioMIP)²⁴,²⁵ initiative was developed. Here we examine the broad PlioMIP ensemble, including phases 1 and 2 with a total of 25 models (Extended Data Tables 1 and 2), to better understand how ENSO activity might change in warmer climates.

## Reduced ENSO amplitude

The PlioMIP ensemble simulates significant reduction in the variability of SST anomalies across most of the global tropics in the mPWP compared with pre-industrial (Fig. 1a; see Extended Data Fig. 1 for PlioMIP1). Although there are notable changes in the Indian and Atlantic oceans²⁶, the most pronounced weakening occurs in the equatorial Pacific, where reduced SST variability in the eastern basin (Niño3 region) is simulated by 21 out of 23 PlioMIP models (including PlioMIP1 and 2). Considering PlioMIP2 models only, there is a multimodel mean amplitude reduction of 25% (±16% standard deviation; Fig. 1b).

Separating the Niño3 variability change into interannual (<10 yr) and longer timescale components shows that all but one model

simulation reduced amplitude in the interannual band (Extended Data Fig. 2), a timescale that is dominated by ENSO. In addition, 75% (17 out of 23) of the models suggest a shift towards lower frequencies as indicated by either an increased amplitude at low frequency (>10 yr) or a more pronounced weakening at interannual than on longer timescales. In this Article, due to data availability, our analysis is performed on the last 100 years of each model's simulation, making the decadal analysis more uncertain.

## Role of equatorial Pacific Ocean changes

ENSO dynamics are dominated by equatorial processes, which are influenced by the background state²⁷. Although the PlioMIP models simulate an amplified eastern Pacific warming (Fig. 2a), there are large intermodel differences in this pattern, as indicated by inconsistent changes in the zonal SST difference²⁸ (Extended Data Fig. 3). Of particular importance for ENSO dynamics are changes in equatorial thermal gradients in the mixed layer⁵,⁶,²⁹.

First, we evaluate changes in the thermocline slope, which plays an important role in ENSO dynamics. Stronger (weaker) westward equatorial currents are associated with increased (decreased) east-west thermocline slope⁶. The thermocline slope better represents the resultant effect of changes in zonal equatorial ocean dynamics than does the zonal SST gradient in the PlioMIP models, as reflected in a higher intermodel correlation with ENSO amplitude change (Fig. 2b and Extended Data Fig. 3). Models with a steeper mean thermocline in the mPWP are typically associated with larger ENSO amplitude reductions, while a flatter mean thermocline is associated with either a slight increase or a weak decrease in ENSO variability ($r_s = -0.43$; Fig. 2b). This indicates that an equatorial Pacific mean state with a steeper thermocline, which is associated with intensified trades and westward surface currents, is less favourable for strong ENSO variability. Under such a 'La Niña-like' mean state, stronger
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Equatorial Pacific Ocean changes. **a**, PlioMIP2 multimodel mean change in surface tropical and subsurface equatorial Pacific temperatures. The vertical profile is averaged between 5° S and 5° N. Stippling indicates significant change at the 95% level (in the SST panel, the entire basin-wide warming is significant at the 95% level). **b**, Intermodel relationship between the change in the thermocline slope between the eastern and western Pacific (Methods) and the change in the Niño3 amplitude. **c**, Intermodel relationship between the change in ocean stratification and the change in Niño3 amplitude. Ocean stratification was measured as the difference between the average temperature in the top 75 m (green box in **a**) and that at 100 m (blue line in **a**) between 150° E and 140° E.</figcaption></figure>

initial anomalies are required to weaken the climatological conditions sufficiently for El Niño development⁶.

Another mean-state factor that can affect ENSO development is the equatorial upper-ocean stratification⁵. In particular, western-central equatorial Pacific stratification influences the variability of strong ENSO events through changes in the dynamical ocean–atmosphere coupling. As such, we evaluate ocean stratification in the central-western Pacific, a region where wind anomalies trigger oceanic Kelvin and Rossby waves, which influence ENSO genesis³⁰. Models with decreased ocean stratification are typically associated with larger ENSO reductions, while weaker reductions occur in models where ocean stratification increases (Fig. 2c). Given that over half of the models show increased stratification, this relationship cannot explain the consistent decrease in ENSO across the ensemble. Similarly, the fact that many models show a decrease in thermocline slope indicates that this is not the underlying cause for ENSO amplitude reduction. In summary, while changes in the thermocline and stratification help explain intermodel differences in ENSO amplitude changes, it appears that other processes that apply across models are required for the overall weakening of ENSO variability.

## Off-equatorial Pacific changes

While ENSO development is closely related to zonal equatorial dynamics²⁷, ENSO events are also affected by a variety of other large-scale processes beyond the equatorial Pacific⁷,³¹,³². For example, changes to the mean meridional SST gradient or processes in the extratropics can play an important role in triggering ENSO events. In particular, all PlioMIP models simulate a weaker Equator-to-pole temperature gradient associated with polar amplified warming³³.

We first evaluate the role of meridional SST gradients through possible displacements of the ITCZ in the mPWP. Southward (northward) ITCZ displacements, due to reduced (increased) off-equatorial SST gradients, affect ENSO activity through increased (reduced) probability of occurrences of deep convection in the central-eastern equatorial Pacific⁸,³⁴. Here we show that a mean northward ITCZ shift during austral spring–summer, developing and mature ENSO phases, is strongly related to the ENSO weakening across models ($r_s = -0.64$; Fig. 3a). This scenario increases convergence throughout the tropical North Pacific that suppresses anomalous convergence feedback at the Equator (Extended Data Fig. 4). To further illustrate this effect, we evaluate models' performance in simulating the nonlinear relationship between SST
<figure><img src="image_5.png" /><figcaption>**Fig. 3** | ITCZ-ENSO relationship. **a**, PlioMIP2 intermodel relationship between the change in the Niño3 amplitude and mean ITCZ shift from October to February. Green star indicates values obtained from observations by comparing periods 1979-1999 and 2000-2020. The correlation coefficient was evaluated considering PlioMIP models only. **b**, Relationship between December-January-February (DJF) Niño3 SST anomalies and DJF Niño3 rainfall for the period pre-2000 (red) and post-2000 (green). **c-h**, As in **b** but for selected PlioMIP models that correctly simulate nonlinear ENSO characteristics (Methods): GISS-E2-1-G (**c**), CESM2 (**d**), MRI-CGCM2.3 (**e**), NorESM1-F (**f**), CCSM4-1deg (**g**) and CCSM4-2deg (**h**). Pre-industrial simulation in blue and Pliocene in yellow.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

**Fig. 4** | Changes to potential ENSO triggers. a, Intermodel relationship between the change in the intensity of the western Pacific trade winds (10° S–10° N; 160° E–150° W) and the amplitude (standard deviation) of its monthly variability. To ideally examine changes in the western wind bursts, we would need high-frequency output, which is not available for the PlioMIP models. b, Change in the amplitude (standard deviation) of the SPMM time series, defined as the mean SST anomaly between 15° S and 25° S and between 110° W and 120° E. c, Change in the amplitude (standard deviation) of the meridional wind variability over the SHB region (10° S–30° S; 140° E–170° E). PlioMIP2 models in b and c: a, CCSM4-UofT; b, CCSM4-2deg; c, CESM2; d, COSMOS; e, EC-EARTH3.3; f, GISS-E2-1-G; g, HadCM3; h, IPSL-CM6A-LR; i, IPSL-CM5A; j, IPSL-CM5A2; k, MIROC4m; l, MRI-CGCM2.3; m, NorESM-L; n, NorESM1-F.

anomalies and precipitation events in the eastern Pacific (Methods and Extended Data Fig. 5). Six models correctly simulate this characteristic and indicate that the further north the mean ITCZ migrates, the less probable are occurrences of deep convection events in the eastern Pacific associated with ENSO SST anomalies (Fig. 3). The ITCZ shift can fully explain ENSO weakening across these wix models ($r_s = 0.94$; Extended Data Fig. 6).

This PlioMIP result is analogous to the reduced ENSO activity in the past two decades, which corresponded with a more northward position of the ITCZ³⁵ (Fig. 3a). The multi-decadal period pre-2000 was marked by enhanced ENSO variability while post-2000, it has been reduced by 21%²⁶,³⁷, resulting in weaker rainfall events in the eastern Pacific (Fig. 3b). This reduction has been attributed to a negative phase of the Interdecadal Pacific Oscillation³⁸ with enhanced trade winds and surface ocean currents, which resembles a La Niña-like mean state with the Pacific ITCZ displaced northwards. Consistently, the PlioMIP models indicate a larger reduction in ENSO activity when shifted towards a La Niña-like mean state (Fig. 2b).

We also evaluate possible changes to other processes that are favourable for initiating ENSO events, such as the reversal of the easterly trade winds in the western Pacific³⁹. In the PlioMIP models, the annual mean intensification of the western Pacific trade winds corresponds with weaker wind variability over this region (Fig. 4a). Climatologically stronger easterly trades tend to inhibit: (1) the stochastic forcing of westerly wind bursts in the western Pacific⁴⁰ that triggers the positive thermocline feedback; (2) southward shifts of the ITCZ through positive wind–evaporation–SST feedback³⁴, which cools the equatorial Pacific Ocean, increasing the meridional SST gradient; and (3) eastward displacements of the Walker circulation.

Further, we evaluate patterns of variability that promote wind anomalies in the western Pacific and contribute to the development of El Niño events. First, the South Pacific meridional mode

(SPMM), analogue to the North Pacific meridional mode (NPMM), is initiated by the weakening of off-equatorial southeast trade winds in the eastern Pacific. This alters the latent heat flux, triggering a wind–evaporation–SST feedback that propagates wind anomalies into the tropics³². We find that all but two PlioMIP2 models simulate decreased SPMM variability in the mPWP (Fig. 4b). Equivalent changes in the NPMM are not consistent across models and do not help explain ENSO changes (Extended Data Fig. 7).

Finally, extreme El Niño events are amplified by an anomalous zonal pressure dipole in the Southern Hemisphere³¹. In such condition, an anomalous high pressure over Australia facilitates cold surges through the Coral Sea (the Southern Hemisphere booster (SHB))³¹, which promote westerly wind bursts in the western Pacific conducive for El Niño development. This meridional wind variability in the SHB region also decreases in 10 out of 12 PlioMIP2 models (Fig. 4c). All these changes are associated with reduced probability of El Niño initiation, which results in weaker ENSO activity. It is important to note that a northward ITCZ shift probably had a major effect on ENSO triggers from the Southern Hemisphere due to changes in the large-scale atmospheric circulation, as we evaluate next.

## Large-scale forcing

The Pacific ITCZ–ENSO relationship demonstrated in the previous section can be a result of either a large-scale global ITCZ shift modulating ENSO or a local response of the Pacific ITCZ to changes in ENSO activity. The PlioMIP models indicate that the northward ITCZ shift during the mPWP occurs in all basins, as indicated by anomalous meridional dipoles in rainfall across the global tropics (Fig. 5a; see Extended Data Fig. 8 for PlioMIP1). In addition, the PlioMIP models systematically simulate asymmetric polar amplified warming in both hemispheres (Fig. 5b), which can give rise to large-scale changes in the meridional temperature gradient and affect the ITCZ position through changes in atmospheric heat fluxes⁴¹.
a

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

b

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

c

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

d

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

e

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

f

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

**Fig. 5** | Energetic constraints for the ITCZ position. **a,** DJF precipitation change in the PlioMIP2 models (mPWP minus pre-industrial). Stippling indicates where the change is significant at the 95% level. **b,** Multimodel mean change zonally averaged SST for PlioMIP1 (magenta) and PlioMIP2 (red). Banding indicates standard deviation range. **c,** Cam4: changes in DJF atmospheric energy flux, computed as the residual between the total top-of-the-atmosphere and surface energy fluxes, in the AGCM experiments forced with PlioMIP1 and PlioMIP2 climatological SST and sea ice (Methods). Banding indicates standard deviation range of a five-member ensemble. Negative anomalies in the Northern Hemisphere indicate weakening of the northward heat transport, while negative anomalies in the Southern Hemisphere indicate intensification of the southward heat transport. **d,** CAM4: changes in the meridional stream function in the AGCM experiment forced with PlioMIP2 SST and sea ice. Contours indicate pre-industrial stream function (zero contour in bold). Colours indicate change. **e,** Intermodel relationship between changes in the intensity of the zonal western Pacific trades and ITCZ shift during austral summer. **f,** Changes in low-level (850 hPa) winds and stream function in the PlioMIP2 models. Wind changes are shown only where changes are significant at the 95% level. Maps created using the Basemap library for Python.

It is important to note that increased rainfall south of the Equator in the eastern Pacific may be a result of double-ITCZ bias in the PlioMIP models⁹,⁴². A more consistent northward ITCZ shift across the tropical North Pacific is evident through increased low-level wind convergence (Extended Data Fig. 4), which indicates that increased precipitation in the eastern Pacific is probably a result of the thermodynamic effect over the double-ITCZ region¹³.

The ITCZ northward shift is not consistent with the equatorial warming (Fig. 2a), which would otherwise tend to shift the ITCZ southwards. To assess the role of the large-scale SST warming

patterns in the ITCZ shift, we performed sensitivity experiments using an atmospheric general circulation model (AGCM; the National Center for Atmospheric Research Community Atmospheric Model version 4 (CAM4)). Here the AGCM is forced with PlioMIP climatological SSTs, which allow us to isolate changes in atmospheric circulation from changes in ocean-atmosphere variability, such as ENSO. Note that the mPWP climatological-mean warming pattern, used to force the atmospheric model, may still contain some nonlinear influence of ENSO changes, but this effect is negligible (Methods).
<figure><img src="image_13.png" /><figcaption>**Fig. 6** | Schematic of the drivers of suppressed ENSO activity in the mPWP. A northward ITCZ shift reduces the probability of occurrence of deep convection in the central-eastern Pacific. Energetic constraints for the ITCZ position indicate that higher rates of warming in the Northern Hemisphere drive a northward ITCZ shift and enhanced Southern Hemisphere Hadley circulation. These changes are also associated with intensified surface subtropical high and western Pacific trades. Enhanced trade winds suppress wind variability in the western Pacific, which is important for El Niño initiation. An intensified subtropical high is thought to impede zonal pressure anomalies across the tropical South Pacific and, thus, suppress the activity of the SPMM and SHB, which are important for the development of strong El Niño events. Map created using MATLAB®.</figcaption></figure>

In the present climate, during austral summer, increased insolation in the Southern Hemisphere results in intensification of the Northern Hemisphere Hadley circulation, northward energy flux across the Equator (energy-flux-Equator), and southward ITCZ shift¹¹. In the mPWP, the AGCM simulates decreased northward energy-flux-Equator during the austral summer (Fig. 5c). Due to the mutual relationship between changes in the energy-flux-Equator and ITCZ position, a decreased northward energy-flux-Equator is accompanied by a northward ITCZ shift in agreement with a recent PlioMIP2 study⁴⁴. Higher rates of warming in the Northern Hemisphere drive an intensification and northward expansion of the Southern Hemisphere Hadley cell and weaker circulation in the Northern counterpart (Fig. 5d; see Extended Data Fig. 8 for PlioMIP1), which reduces the atmospheric energy input from the Southern Hemisphere to the Northern Hemisphere during the austral summer.

The AGCM experiments suggest that the meridional displacement of the ITCZ is a global feature of the PlioMIP simulations and occurs due to the extratropical warming. The experiments indicate an overall decrease in the northward atmospheric heat transport in the Northern Hemisphere and a slight increase in the southward heat transport in the Southern Hemisphere (Fig. 5c), which initially points to changes in pole-to-pole temperature gradient. One of the most robust features of the mPWP simulations is the asymmetric polar amplified warming (Fig. 5b), which increases the inter-hemispheric temperature gradient and was caused by reduced sea-ice volume⁴⁵. However, whether the mPWP ITCZ shift was driven by sea-ice changes needs to be further investigated.

The large-scale changes in the meridional circulation probably induce changes in horizontal circulation. In the Pacific Ocean, the PlioMIP models indicate that a northward ITCZ shift is significantly related to intensified western Pacific trades (Fig. 5d), which is analogous to synchronized shifts of the Walker and Hadley circulations during different ENSO phases⁴⁶. An analysis of the low-level circulation indicates that the anomalously stronger western trades in the mPWP are sourced at the subtropical South Pacific due to an intensified circulation of the South Pacific subtropical-high system (Fig. 5e,f; see Extended Data Fig. 7 for PlioMIP1). These changes

are not exclusive to the South Pacific but occur in all ocean basins (Fig. 5f). The synchronized changes in the meridional and zonal atmospheric circulation are probably a result of global changes in atmospheric heat fluxes during the warmer mPWP. This illustrates a possible influence of changes in global atmospheric dynamics on ENSO in a warmer climate.

## Implications for past and future climates

The results presented here suggest a link between reduced ENSO amplitude and the northward shift of the ITCZ in the mPWP, associated with stronger climatological circulation in the Southern Hemisphere (Fig. 6). The northward shift of the ITCZ reduces the probability of ENSO-related rainfall events in the eastern Pacific. Northward ITCZ shift and intensified Southern Hemisphere Hadley and subtropical circulations are a response to enhanced Northern Hemisphere warming via energetic constraints for the ITCZ position⁴¹ (Fig. 6). This intensified Southern Hemisphere circulation reduces wind variability in the western Pacific, may suppress zonal sea-level pressure anomalies imposed by the SPMM and the SHB and weakens and shifts the South Pacific convergence zone polewards⁴⁷, reducing its interaction with equatorial processes (Fig. 6). As such, the climatological stability imposed by intensified tropical Southern Hemisphere circulation acts to increase ENSO stability, as ENSO, by definition, is a deviation from the mean climate, and thus stronger climatological circulations can be viewed as unfavourable to ENSO-induced changes¹⁰.

In addition to the reduced ENSO amplitude, SST variability in other tropical basins also decreases (Fig. 1a). This may also contribute to weakened ENSO variability via pan-tropical teleconnections related to a delayed and weaker negative feedback⁴⁸, although reduced variability in other tropical basins itself might also be a consequence of reduced ENSO variability. Ref. ²⁶ reported that all PlioMIP1 models simulate reduced tropical North Atlantic variability associated with a warming of this basin and a northward Atlantic ITCZ shift. Taken together, these results suggest that a northward shift of the global ITCZ may mute tropical Pacific and Atlantic SST variability.
Our results are subject to a number of uncertainties in the
simulations tied to sparse and limited proxy data, which are used
to constrain the PlioMIP experiments, and systematic climate
model biases⁴⁹. Changes in the inter-hemispheric SST gradient, for
example, could be affected by uncertainties in the extension of the
mPWP ice sheets⁵⁰, poor representation of certain polar feedbacks⁵¹
(interactive land-ice), climate sensitivity⁵² and biases in tropical
convection and SST of the climate models, such as double ITCZs⁴²
and an overly strong cold tongue. Despite different model biases,
we show that the current generation of climate models simulate
consistent changes to ENSO related to shifts in the ITCZ position
in the mPWP.

Palaeoclimate states may have particular relevance as analogues for the future climate. However, our findings indicate that, although the mPWP warming is comparable to that projected for the end of the twenty-first century under a 'business as usual' scenario (~3°C) (ref. ²¹), the simulated mPWP ENSO response is the opposite of that projected⁸,¹¹. An important factor in this difference is that the mPWP represents an equilibrium climate, albeit with similar CO₂ levels as today. If equilibrium conditions are of particular importance, this suggests that a more Pliocene-like climate might be possible if present-day CO₂ concentrations were to be maintained constant and a steady state is reached. However, the current rate of atmospheric CO₂ rise is unprecedented in Earth's history, resulting in warming trends that differ compared with past regimes. Thus, relating past and future warmings is not straightforward. Here evaluating the empirically based mPWP warming, we find that a northward ITCZ shift suppresses ENSO activity. If this relationship can be applied to the twenty-first-century projections where a southward shift of the Pacific ITCZ is projected⁷, then an increase in ENSO variability⁵ appears to be a potential outcome.

Online content

Any methods, additional references, Nature Research report-
ing summaries, source data, extended data, supplementary in-
mation, 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-022-00999-y.

Received: 7 April 2021; Accepted: 29 June 2022;
Published online: 11 August 2022

References

1. McPhaden, M. J., Zebiak, S. E. & Glantz, M. H. ENSO as an integrating concept in Earth science. *Science* **314**, 1740–1745 (2006).

2. Santoso, A. et al. Dynamics and predictability of El Niño-Southern Oscillation: an Australian perspective on progress and challenges. *Bull. Am. Meteorol. Soc.* **100**, 403–420 (2019).

3. Collins, M. et al. The impact of global warming on the tropical Pacific Ocean and El Niño. *Nat. Geosci.* **3**, 391–397 (2010).

4. Taschetto, A. S. et al. Cold tongue and warm pool ENSO Events in CMIP5: mean state and future projections. *J. Clim.* **27**, 2861–2885 (2014).

5. Cai, W. et al. Increased variability of eastern Pacific El Niño under greenhouse warming. *Nature* **564**, 201–206 (2018).

6. Santoso, A. et al. Late-twentieth-century emergence of the El Niño propagation asymmetry and future projections. *Nature* **504**, 126–130 (2013).

7. Cai, W. et al. ENSO and greenhouse warming. *Nat. Clim. Change* **5**, 849–859 (2015).

8. Cai, W. et al. Increasing frequency of extreme El Niño events due to greenhouse warming. *Nat. Clim. Change* **4**, 111–116 (2014).

9. Mamalakis, A. et al. Zonally contrasting shifts of the tropical rain belt in response to climate change. *Nat. Clim. Change* https://doi.org/10.1038/s41558-020-00963-x (2021).

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

11. Grothe, P. R. et al. Enhanced El Niño–Southern Oscillation variability in recent decades. *Geophys. Res. Lett.* https://doi.org/10.1029/2019gl083906 (2019).

12. Ford, H. L., Ravelo, A. C. & Polissar, P. J. Reduced El Niño-Southern Oscillation during the Last Glacial Maximum. *Science* **347**, 255–258 (2015).

13. Wara, M. W., Ravelo, A. C. & Delaney, M. L. Permanent El Niño-like conditions during the Pliocene Warm Period. *Science* **309**, 758–761 (2005).

14. Tierney, J. E., Haywood, A. M., Feng, R., Bhattacharya, T. & Otto-Bliesner, B. L. Pliocene warmth consistent with greenhouse gas forcing. *Geophys. Res. Lett.* **46**, 9136–9144 (2019).

15. Wycech, J. B., Gill, E., Rajagopalan, B., Marchitto, T. M. & Molnar, P. H. Multiproxy reduced-dimension reconstruction of Pliocene equatorial Pacific sea surface temperatures. *Paleoceanogr. Paleoclimatol. 35*, e2019PA003685 (2020).

16. O'Brien, C. L. et al. High sea surface temperatures in tropical warm pools during the Pliocene. *Nat. Geosci.* **7**, 606–611 (2014).

17. Brierley, C. M. et al. Greatly expanded tropical warm pool and weakened Hadley circulation in the early Pliocene. *Science* **323**, 1714–1718 (2009).

18. Watanabe, T. et al. Permanent El Niño during the Pliocene Warm Period not supported by coral evidence. *Nature* **471**, 209–211 (2011).

19. Scroxton, N. et al. Persistent El Niño-Southern Oscillation variation during the Pliocene epoch. *Paleoceanography* **26**, PA2215 (2011).

20. White, S. M. & Ravelo, A. C. Dampened El Niño in the early Pliocene Warm Period. *Geophys. Res. Lett.* https://doi.org/10.1029/2019GL085504 (2020).

21. Burke, K. D. et al. Pliocene and Eocene provide best analogs for near-future climates. *Proc. Natl Acad. Sci. USA* **115**, 13288–13293 (2018).

22. McClymont, E. L. et al. Lessons from a high-CO₂ world: an ocean view from ~3 million years ago. *Clim. Past* **16**, 1599–1615 (2020).

23. Haywood, A. M., Dowsett, H. J. & Dolan, A. M. Integrating geological archives and climate models for the mid-Pliocene Warm Period. *Nat. Commun.* **6**, 10646 (2016).

24. Haywood, A. M. et al. Pliocene Model Intercomparison Project (PlioMIP): experimental design and boundary conditions (experiment 2). *Geosci. Model Dev.* **4**, 571–577 (2011).

25. Haywood, A. M. et al. The Pliocene Model Intercomparison Project (PlioMIP) phase 2: scientific objectives and experimental design. *Clim. Past* **12**, 663–675 (2016).

26. Pontes, G. M., Wainer, I., Prado, L. & Brierley, C. Reduced Atlantic variability in the mid-Pliocene. *Climatic Change* **160**, 445–461 (2020).

27. Jin, F. F., Kim, S. T. & Bejarano, L. A coupled-stability index for ENSO. *Geophys. Res. Lett.* **33**, L23708 (2006).

28. Oldeman, A. M. et al. Reduced El-Niño variability in the mid-Pliocene according to the PlioMIP2 ensemble. *Clim. Past* **17**, 2427–2450 (2021).

29. Wang, B. et al. Historical change of El Niño properties sheds light on future changes of extreme El Niño. *Proc. Natl Acad. Sci. USA* **116**, 22512–22517 (2019).

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

31. Hong, L.-C., LinHo & Jin, F.-F. A Southern Hemisphere booster of super El Niño. *Geophys. Res. Lett.* **41**, 2142–2149 (2014).

32. Zhang, H., Clement, A. & Di Nezio, P. The South Pacific meridional mode: a mechanism for ENSO-like variability. *J. Clim.* **27**, 769–783 (2014).

33. Haywood, A. M. et al. The Pliocene Model Intercomparison Project phase 2: large-scale climate features and climate sensitivity. *Clim. Past* **16**, 2095–2123 (2020).

34. Okajima, H., Xie, S.-P. & Numaguti, A. Interhemispheric coherence of tropical climate variability: effect of the climatological ITCZ. *J. Meteorol. Soc. Jpn. Ser.* **2** **81**, 1371–1386 (2003).

35. Hu, S. & Fedorov, A. V. Cross-equatorial winds control El Niño diversity and change. *Nat. Clim. Change* **8**, 798–802 (2018).

36. Lübbecke, J. F. & McPhaden, M. J. Assessing the twenty-first-century shift in ENSO variability in terms of the Bjerknes stability index. *J. Clim.* **27**, 2577–2587 (2014).

37. Zhao, B. & Fedorov, A. The effects of background zonal and meridional winds on ENSO in a coupled GCM. *J. Clim.* https://doi.org/10.1175/jcli-d-18-0822.1 (2019).

38. England, M. H. et al. Recent intensification of wind-driven circulation in the Pacific and the ongoing warming hiatus. *Nat. Clim. Change* **4**, 222–227 (2014).

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

40. Chen, D. et al. Strong influence of westerly wind bursts on El Niño diversity. *Nat. Geosci.* **8**, 339–345 (2015).

41. Schneider, T., Bischoff, T. & Haug, G. H. Migrations and dynamics of the intertropical convergence zone. *Nature* **513**, 45–53 (2014).

42. Tian, B. & Dong, X. The double-ITCZ bias in CMIP3, CMIP5 and CMIP6 models based on annual mean precipitation. *Geophys. Res. Lett.* https://doi.org/10.1029/2020GL087232 (2020).

43. Held, I. M. & Soden, B. J. Robust responses of the hydrological cycle to global warming. *J. Clim.* **19**, 5686–5699 (2006).
44. Han, Z. et al. Evaluating the large-scale hydrological cycle response within the Pliocene Model Intercomparison Project phase 2 (PlioMIP2) ensemble. *Clim. Past* **17**, 2537–2558 (2021).

45. Hill, D. J. et al. Evaluating the dominant components of warming in Pliocene climate simulations. *Clim. Past* **10**, 79–90 (2014).

46. Yun, K. S., Timmermann, A. & Stuecker, M. Synchronized spatial shifts of Hadley and Walker circulations. *Earth Syst. Dyn.* **12**, 121–132 (2021).

47. Pontes, G. M. et al. Drier tropical and subtropical Southern Hemisphere in the mid-Pliocene Warm Period. *Sci. Rep.* **10**, 13458 (2020).

48. Cai, W. et al. Pantropical climate interactions. *Science* **363**, eaav4236 (2019).

49. Fedorov, A. V. et al. Patterns and mechanisms of early Pliocene warmth. *Nature* **496**, 43–49 (2013).

50. Koenig, S. J. et al. Ice sheet model dependency of the simulated Greenland ice sheet in the mid-Pliocene. *Clim. Past* **11**, 369–381 (2015).

51. Fischer, H. et al. Palaeoclimate constraints on the impact of 2°C anthropogenic warming and beyond. *Nat. Geosci.* **11**, 474–485 (2018).

52. Zhu, J., Poulsen, C. J. & Otto-Bliesner, B. L. High climate sensitivity in CMIP6 model not supported by paleoclimate. *Nat. Clim. Change* **10**, 378–379 (2020).

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

Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
## Methods

**Models and data.** A total of 9 PlioMIP1 and 16 PlioMIP2 models were analysed. See Extended Data Table 1 for a list of the models included in our analysis. The number of models used in each analysis varies according to data availability in the PlioMIP1 and PlioMIP2 databases. The last 100 years of each model's simulation is used. We additionally use observational SST and precipitation from the Extended Reconstructed Sea Surface Temperature version 5 and Global Precipitation Climatology Project datasets, respectively.

**PlioMIP1 and 2 protocols.** PlioMIP phases 1 and 2 apply rather similar boundary conditions (Extended Data Table 2). Nonetheless, there were significant differences at some regional locations, which can potentially affect large-scale climate⁵⁴,⁵⁵. Both phases applied an mPWP land-sea mask, but PlioMIP2 land-sea mask accounts for glacial isostatic adjustments and changes dynamic topography. This resulted in a subaerial Canadian Archipelago, Bering Strait and emergence of Sunda and Sahul shelves in the Indonesia and Australia region. Phase 2 models also applied soils and lakes reconstructions and a newer reconstruction of the Greenland ice sheet that now accounts for a 70% retreat, instead of the 50% retreat applied in phase 1. These reconstructions were derived from the US Geological Survey PRISM dataset, specifically, the most recent and fourth iteration of the reconstructions (PRISM)⁵³. In PlioMIP2, modelling groups could use either the reconstructed vegetation⁵⁴, same as PlioMIP1, or the dynamic model vegetation option. COSMOS was the only model among PlioMIP2 participants to use dynamic vegetation. CO₂ concentrations were set to 405 and 400 ppm in phases 1 and 2, respectively. For a detailed description of each model's implementation, see refs. ⁵⁵⁻⁷¹ (Extended Data Table 1). Beyond differences in boundary conditions, there are also fundamental differences in the conception of PlioMIP2 versus PlioMIP1. Although both phases have not applied changes in orbital parameters, phase 1 was designed to simulate a time-averaged global SST reconstruction between 3.0 and 3.3 Ma, while phase 2 focuses on a narrower time slice (Marine Isotope Stage KM5c at 3.205 Ma) with almost identical orbital parameters to modern.

**Statistical significance of the changes.** This is measured through model agreement on the sign of the change. This method is based on a binomial distribution of equal probability (P = q = 0.5). Here we consider that all models have an equal probability of simulating positive and negative changes in the mPWP simulation. As such, the cumulative probability distribution function of a binomial distribution of N=9 (PlioMIP1) and N=16 (PlioMIP2) models shows that the 95% probability level is reached when there is a model agreement on the sign of the change of 7 and 11 models, respectively. In addition, we use the non-parametric Spearman rank-correlation test (r_s) to determine whether there is a monotonic relation between two variables. Note that the assumption of sample independence may not be completely satisfied, given that climate models share common components and physical equations. In addition, the CESM family of models may be overrepresented in the PlioMIP ensemble; however, the differing results obtained among their simulations may allow us to consider these models independent. To illustrate that, we performed a sensitivity analysis where each model from the CESM family was considered at a time when computing the Spearman rank correlation for the relationship shown in Fig. 3a. The coefficients ranged from -0.55 (P=0.01) to -0.63 (P=3×10⁻³).

**ENSO amplitude.** The standard deviation of Niño3 index is used to represent ENSO amplitude. The Niño3 index is calculated from SST anomalies averaged over the eastern Pacific region between 5°N and 5°S latitude and between 150°W and 90°W longitude. SST anomalies were computed by removing the mean annual cycle.

**Frequency separation.** The amplitude of low-frequency variability (>10 yr) is evaluated through the variance of the 11 yr running mean Niño3 time series in each model. The amplitude of the interannual period is estimated as the variance of the residual time series, that is, original Niño3 time series subtracted from the Niño3 decadal time series.

**Thermocline slope.** Thermocline slope is the difference between mean eastern (5°S–5°N; 150°W–90°W) and western (5°S–5°N; 160°E–150°W) Pacific thermocline depths. The thermocline depths are computed from the mean temperature profile in each of the boxes indicated in the preceding. This is the weighted average depth, based on depths in which the temperature gradients are greater than 50% of its maximum.

**Equatorial Pacific Ocean stratification.** This is the difference between the mean temperature in the top 75 m and the temperature at 100 m from 150°E to 140°W (as indicated in Fig. 2a).

**Pacific ITCZ position.** The ITCZ position⁴⁵ is taken as the average latitudes over which precipitation in the tropical Pacific Ocean (20°S–20°N) is greater than 50% of the maximum zonally averaged precipitation over 120°E–90°W. This method may consider the double-ITCZ bias⁴² if the double-ITCZ associated precipitation is greater than 50% of the maximum. The double-ITCZ bias is an artificial feature

produced by most climate models that overestimates the tropical precipitation south of the Equator. Here we define the ITCZ bias as the difference between simulated pre-industrial Pacific ITCZ position and the observed position averaged from 1979 to 2020. Although the PlioMIP models suffer from double-ITCZ bias (mean bias: -4.1°±2.1° s.d.), we do not find a statistically significant relationship with ENSO amplitude changes (r_s = -0.16; P = 0.45). Note that the ITCZ bias is evaluated by comparing the pre-industrial model simulations with modern climate, which is already under the influence of global warming.

**Criteria for model selection.** Models were selected according to their ability to simulate ENSO nonlinear characteristics⁴¹. Models were required to be able to simulate December–January–February Niño3 precipitation greater than 5 mm d⁻¹ and Niño3 precipitation skewness greater than 1 in the pre-industrial control run. These criteria underscore the essential definition of an extreme El Niño¹, which is fundamental to the ENSO system in observations⁷⁵. Out of 14 PlioMIP2 models, 6 models met these criteria (Extended Data Fig. 5). The skewness criterion filters out models that systematically simulate overly wet and dry conditions in the eastern equatorial Pacific. Such biases tend to reduce rainfall skewness in the models as they simulate SSTs well below or above the convective threshold of 26–28°C (ref. ⁷⁶), affecting Niño3 precipitation variability.

**Atmospheric subtropical high systems.** Quantifying the intensity of the subtropical highs is not a simple task when dealing with different climate backgrounds (+2–3 K) as the global pressure weakens in a warmer atmosphere. To overcome this pressure issue, we compute the stream function at 850 hPa to identify the position and intensity of the subtropical high systems. The stream function can be derived from the geostrophic balance:

$$f \times \mathbf{v} = -\frac{1}{\rho} \nabla_H p$$

where $\nabla_H$, $\mathbf{v} = (u_g, v_g)$, p, f and ρ denotes the horizontal gradient, velocity vector, pressure, Coriolis function and density, respectively. Knowing that for a fluid of horizontally uniform density, the geostrophic flow in an f-plane is non-divergent, that is

$$\frac{\partial u_g}{\partial x} + \frac{\partial v_g}{\partial y} = 0 \text{ for } \rho = \rho_0 \text{ (g) and } f = f_0 = 2\Omega \sin\theta,$$

we can define a stream function that yields to

$$u_g = -\frac{\partial\psi}{\partial y} = -\frac{1}{\rho_0 f_0} \frac{\partial p}{\partial y} \text{ and } v_g = \frac{\partial\psi}{\partial x} = \frac{1}{\rho_0 f_0} \frac{\partial p}{\partial x}$$

Note that in the Southern Hemisphere, increased pressure gradients over geostrophic flows result in intensified anticyclonic circulation (negative stream function).

**SPMM amplitude.** SPMM amplitude is computed as the amplitude (standard deviation) of mean SST anomalies from 15°S to 25°S and from 110°W to 120°W.

**NPMM amplitude.** NPMM amplitude is computed as the amplitude (standard deviation) of mean SST anomalies from 20°S to 25°S and from 142°W to 138°W.

**SHB amplitude.** SHB amplitude is computed as the amplitude (standard deviation) of meridional wind anomalies from 10°S to 30°S and from 140°E to 170°E.

**CAM4 experiments.** We undertook four experiments, with multiple ensemble members, using the National Center for Atmospheric Research CAM4: (1) mean mid-Pliocene SST and sea-ice forcing from PlioMIP1. PlioMIP1 SST and sea-ice were time and ensemble averaged to force the CAM4 model; (2) mean pre-industrial SST and sea ice as simulated by PlioMIP1 models for comparison; experiments 3 and 4 consisted of repeating experiments 1 and 2 but with PlioMIP2 SST and sea ice. For each experiment, five ensemble members were integrated with slightly different initial conditions: each ensemble member was initialized from a different day of the year. The CO₂ forcing was kept as pre-industrial at 280 ppm, and no changes over continental areas were made in all experiments. Each experiment was run for 31 years. The first year of each simulation was discarded due to the atmospheric spin-up. To check whether nonlinearities in ENSO affected the mean SST change, we compared the multimodel mean mPWP warming during all years and during non-ENSO years only. Differences in the tropical Pacific were approximately two orders of magnitude (<0.05 K) lower than the mean tropical Pacific warming (~2 K).

## Data availability

PlioMIP2 data (with exception of IPSL-CM6A and GISS2.1 G) are available upon request to Alan M. Haywood (a.m.haywood@leeds.ac.uk). PlioMIP2 data from CESM2, EC-Earth3.3, NorESM1-F, IPSL-CM6A and GISS2.1 G can be obtained
directly through the Earth System Grid Federation repository (ESGF; https://esgf-node.llnl.gov/search/cmpip6/). Observational SST and precipitation data can be found in the NOAA-USA (NOAA Extended Reconstructed SST version 5) and NCAR-USA (Global precipitation climatology project) online repositories, respectively.

## Code availability

Computer codes are available at https://github.com/gmpontes/Nature_Geoscience_ENSO_ITCZ_PlioMIP.git or upon request to the corresponding author.

## References

53. Dowsett, H. et al. The PRISM4 (mid-Piacenzian) paleoenvironmental reconstruction. *Clim. Past* **12**, 1519–1538 (2016).

54. Salzmann, U., Haywood, A. M., Lunt, D. J., Valdes, P. J. & Hill, D. J. A new global biome reconstruction and data-model comparison for the Middle Pliocene. *Glob. Ecol. Biogeogr.* **17**, 432–447 (2008).

55. Rosenbloom, N. A., Otto-Bliesner, B. L., Brady, E. C. & Lawrence, P. J. Simulating the mid-Pliocene Warm Period with the CCSM4 model. *Geosci. Model Dev.* **6**, 549–561 (2013).

56. Stepanek, C. & Lohmann, G. Modelling mid-Pliocene climate with COSMOS. *Geosci. Model Dev.* **5**, 1221–1243 (2012).

57. Zheng, W., Zhang, Z., Chen, L. & Yu, Y. The mid-Pliocene climate simulated by FGOALS-g2. *Geosci. Model Dev.* **6**, 1127–1135 (2013).

58. Chandler, M. A., Sohl, L. E., Jonas, J. A., Dowsett, H. J. & Kelley, M. Simulations of the mid-Pliocene Warm Period using two versions of the NASA/GISS ModelE2-R Coupled Model. *Geosci. Model Dev.* **6**, 517–531 (2013).

59. Bragg, F. J., Lunt, D. J. & Haywood, A. M. Mid-Pliocene climate modelled using the UK Hadley Centre Model: PlioMIP experiments 1 and 2. *Geosci. Model Dev.* **5**, 1109–1125 (2012).

60. Contoux, C., Ramstein, G. & Jost, A. Modelling the mid-Pliocene Warm Period climate with the IPSL coupled model and its atmospheric component LMDZ5A. *Geosci. Model Dev.* **5**, 903–917 (2012).

61. Chan, W.-L., Abe-Ouchi, A. & Ohgaito, R. Simulating the mid-Pliocene climate with the MIROC general circulation model: experimental design and initial results. *Geosci. Model Dev.* **4**, 1035–1049 (2011).

62. Kamae, Y. & Ueda, H. Mid-Pliocene global climate simulation with MRI-CGCM2.3: set-up and initial results of PlioMIP experiments 1 and 2. *Geosci. Model Dev.* **5**, 793–808 (2012).

63. Zhang, Z. et al. Pre-industrial and mid-Pliocene simulations with NorESM-L. *Geosci. Model Dev.* **5**, 523–533 (2012).

64. Peltier, W. R. & Vettoretti, G. Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: a “kicked” salt oscillator in the Atlantic. *Geophys. Res. Lett.* **41**, 7306–7313 (2014).

65. Chandan, D. & Peltier, W. R. Regional and global climate for the mid-Pliocene using the University of Toronto version of CCSM4 and PlioMIP2 boundary conditions. *Clim. Past* **13**, 919–942 (2017).

66. Feng, R., Otto-Bliesner, B. L., Brady, E. C. & Rosenbloom, N. Increased climate response and Earth system sensitivity from CCSM4 to CESM2 in mid-Pliocene simulations. *J. Adv. Model. Earth Syst.* **12**, e2019MS002033 (2020).

67. Stepanek, C., Samakinwa, E., Knorr, G. & Lohmann, G. Contribution of the coupled atmosphere-ocean-sea ice-vegetation model COSMOS to the PlioMIP2. *Clim. Past* **16**, 2275–2323 (2020).

68. Zheng, J., Zhang, Q., Li, Q., Zhang, Q. & Cai, M. Contribution of sea ice albedo and insulation effects to Arctic amplification in the EC-Earth Pliocene simulation. *Clim. Past* **15**, 291–305 (2019).

69. Hunter, S. J., Haywood, A. M., Dolan, A. M. & Tindall, J. C. The HadCM3 contribution to PlioMIP phase 2. *Clim. Past* **15**, 1691–1713 (2019).

70. Lurton, T. et al. Implementation of the CMIP6 forcing data in the IPSL-CM6A-LR model. *J. Adv. Model. Earth Syst.* **12**, e2019MS001940 (2020).

71. Tan, N. et al. Modeling a modern-like $p_{CO2}$ warm period (Marine Isotope Stage KM5c) with two versions of an Institut Pierre Simon Laplace atmosphere-ocean coupled general circulation model. *Clim. Past* **16**, 1–16 (2020).

72. Chan, W.-L. & Abe-Ouchi, A. Pliocene Model Intercomparison Project (PlioMIP2) simulations using the Model for Interdisciplinary Research on Climate (MIROC4m). *Clim. Past* **16**, 1523–1545 (2020).

73. Kamae, Y., Yoshida, K. & Ueda, H. Sensitivity of Pliocene climate simulations in MRI-CGCM2.3 to respective boundary conditions. *Clim. Past* **12**, 1619–1634 (2016).

74. Li, X., Guo, C., Zhang, Z., Otterå, O. H. & Zhang, R. PlioMIP2 simulations with NorESM-L and NorESM1-F. *Clim. Past* **16**, 183–197 (2020).

75. Santoso, A., Mcphaden, M. J. & Cai, W. The defining characteristics of ENSO extremes and the strong 2015/2016 El Niño. *Rev. Geophys.* **55**, 1079–1129 (2017).

76. Johnson, N. C. & Xie, S. P. Changes in the sea surface temperature threshold for tropical convection. *Nat. Geosci.* **3**, 842–845 (2010).

## Acknowledgements

This work was supported by the São Paulo Research Foundation (FAPESP-Brazil grants no. 2016/23670-0, no. 2019/0882-1 and no. 2021/11035-6), the Australian Research Council (ARC FT160100495) including the ARC Centre of Excellence for Climate Extremes (CE110001028) and the NCI National Facility, Canberra. A.S. is supported by CSHOR, a joint research centre between QNLM and CSIRO, and the Australian Government's National Environmental Science Program. PlioMIP2 experiments were supported by FP7 Ideas Programme: European Research Council, Past Earth Network, CEMAC—University of Leeds, JSPS, Earth Simulator at JAMSTEC, Helmholtz Climate Initiative REKLIM, Alfred Wegener Institute's research programme Marine, Coastal and Polar Systems, Swedish Research Council, Swedish National Infrastructure for Computing, Canadian Innovation Foundation, UNINETT Sigma2—the National Infrastructure for High Performance Computing and Data Storage in Norway, Très Grand Centre de calcul du CEA—GENCI, National Science Foundation (NSF—USA), SURFsara Dutch National Computing and Netherlands Organisation for Scientific Research, Exact Sciences. This research is sponsored by National Science Foundation Grants 2103055 to R.F. and 1418411 to B.O.-B. The CESM project is supported primarily by the National Science Foundation. This material is based on work supported by NCAR, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement 1852977. Computing and data storage resources, including the Cheyenne supercomputer (https://doi.org/10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR.

## Author contributions

G.M.P., A.S.T., A.S.G. and A.S. designed the study. G.M.P., A.S.T., A.S.G., A.S. and I.W. contributed to the interpretation of the data and discussions. G.M.P. conducted the analysis, prepared the figures and wrote the original manuscript. A.S.G. produced the schematic in Fig. 6. A.M.H., W.-L.C., A.A.-O., C.S., G.L., S.J.H., J.C.T., M.A.C., L.E.S., W.R.P., D.C., Y.K., K.H.N., Z.Z., C.C., N.T., Q.Z., B.L.O.-B., E.C.B., R.F., A.S.v.d.H. and M.L.J.B. performed the PlioMIP2 simulations. G.M.P. and A.S.T. performed the CAM4 experiments. All authors commented and reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Extended data is available for this paper at https://doi.org/10.1038/s41561-022-00999-y.

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

**Correspondence and requests for materials** should be addressed to Gabriel M. Pontes.

**Peer review information** *Nature Geoscience* thanks Matthew Collins, Soon-II An 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.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 1 | Change in amplitude of Sea Surface Temperature anomalies in the PlioMIP1 models.** SST anomalies were computed by removing the mean annual cycle. Amplitude is defined as the standard deviation of the SST anomaly timeseries. Stippling indicates where the change is significant at the 95% level. Map created using the Basemap library for Python.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 2 | Mid-Pliocene changes in eastern Pacific SST variability by frequency bands.** Percentage change in Niño3 (5°N-5°S; 150°-90°W) standard deviation separated by interannual (dark blue) and low-frequency (>10 yr; light blue) variability. The amplitude of low-frequency oscillations is evaluated through the variance of the 11-year running mean Niño3 time series in each model. The amplitude of the interannual period is estimated as the variance of the residual time series, that is, original Niño3 timeseries subtracted from the Niño3 decadal timeseries.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 3** | Relationship between ENSO amplitude change and zonal SST difference. Zonal SST change evaluated as the difference between the cold tongue (5°S-5°N; 120°W-100°W) and warm pool regions (5°S-5°N; 150°E-170°E) in the equatorial Pacific. PlioMIP1 models are represented by magenta squares while PlioMIP2 models are represented by red circles.
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 4** | Changes in the low-level wind (850hPa) divergence. a PlioMIP1 and b PlioMIP2 from preindustrial. Stippling indicates where changes are significant at the 95% level. A consistent increased convergence in the tropical North Pacific indicates a northward ITCZ shift across the Pacific Ocean. Maps created using the Basemap library for Python.
<figure><img src="image_18.png" /><figcaption>CCSM4-2deg*</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>UofT-CCSM4</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>CCSM4-1deg*</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>GISS-E2-1-G*</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>IPSL-CM6A-LR</figcaption></figure>

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

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

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

<figure><img src="image_26.png" /><figcaption>MRI-CGCM2.3*</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>NorESM-L</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>NorESM1-F*</figcaption></figure>

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

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

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

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

Extended Data Fig. 5 | See next page for caption.
**Extended Data Fig. 5** | ENSO non-linear characteristics for each PlioMIP2 model. Relationship between Niño3 SST anomalies (x-axis, °C) and Niño3 rainfall (y-axis, mm/day) for the last 100 years of the pre-industrial control simulation of each model. Observed relationship was computed from GPCP and ERSSSTv5 datasets from 1979 to 2020. Models that simulate rainfall skew greater than 1 and rainfall anomalies greater than 5 mm/day are marked with a 'red star'.
<figure><img src="image_33.png" /></figure>

**Extended Data Fig. 6** | ITCZ-ENSO relationship for selected models. Relationship between the change in the Niño3 amplitude and the mean October-to-February ITCZ shift for the models that correctly captured ENSO non-linear characteristics.
NPMM amplitude change [%]

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

Extended Data Fig. 7 | North Pacific Meridional Mode amplitude change. Simulated change in the amplitude of the North Pacific Meridional Mode.
a

DJF precipitation change

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

b

CAM4 forced with PlioMIP1 SST

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

c

DJF stream function and wind changes at 850 hPa

1 m/s

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

**Extended Data Fig. 8** | Atmospheric circulation changes for PlioMIP1 models. a multi-model mean DJF precipitation change (mPWP minus pre-industrial). Stippling indicates where changes are significant at the 95% level. b Changes in the meridional streamfunction in the AGCM experiment forced with climatological PlioMIP1 SST and sea-ice (see Methods). c multi-model mean change in low-level (850 hPa) winds and streamfunction. Wind changes are only plotted where there is a significant change at the 95% level. Map created using the Basemap library for Python.
<table><caption>Extended Data Table 1 | Models used in the present study</caption><thead><tr><th></th><th>Institute</th><th>Ref.</th></tr></thead><tbody><tr><td colspan="3"><strong>PlioMIP1</strong></td></tr><tr><td>CCSM4</td><td>National Center for Atmospheric Research - USA</td><td>52</td></tr><tr><td>COSMOS</td><td>Alfred Wegener Institute - Germany</td><td>53</td></tr><tr><td>FGOALS-g2</td><td>LASG - China</td><td>54</td></tr><tr><td>GISS-E2-R</td><td>Goddard Institute for Space Studies - USA</td><td>55</td></tr><tr><td>HadCM3</td><td>University of Leeds - United Kingdom</td><td>56</td></tr><tr><td>IPSL-CM5A</td><td>Laboratoire des Sciences du Climat et de l'Environnement - France</td><td>57</td></tr><tr><td>MIROC4m</td><td>University of Tokyo, NIES, FRCGC - Japan</td><td>58</td></tr><tr><td>MRI-CGCM2.3</td><td>University of Tsukuba - Japan</td><td>59</td></tr><tr><td>NorESM-L</td><td>Norwegian Research Centre - Norway</td><td>60</td></tr><tr><td colspan="3"><strong>PlioMIP2</strong></td></tr><tr><td>CCSM4-UofT</td><td>University of Toronto - Canada</td><td>61,62</td></tr><tr><td>CCSM4-2deg</td><td>Utrecht University - The Netherlands</td><td>-</td></tr><tr><td>CCSM4-1deg</td><td>National Center for Atmospheric Research - USA</td><td>63</td></tr><tr><td>CESM1.2</td><td>National Center for Atmospheric Research - USA</td><td>63</td></tr><tr><td>CESM2.1</td><td>National Center for Atmospheric Research - USA</td><td>63</td></tr><tr><td>COSMOS</td><td>Alfred Wegener Institute - Germany</td><td>64</td></tr><tr><td>EC-EARTH3.3</td><td>Stockholm University - Sweden</td><td>65</td></tr><tr><td>GISS-E2-1-G</td><td>Goddard institute for Space Studies - USA</td><td>-</td></tr><tr><td>HadCM3</td><td>University of Leeds - United Kingdom</td><td>66</td></tr><tr><td>IPSL-CM6A-LR</td><td>Laboratoire des Sciences du Climat et de l'Environnement - France</td><td>67</td></tr><tr><td>IPSL-CM5A</td><td>Laboratoire des Sciences du Climat et de l'Environnement - France</td><td>68</td></tr><tr><td>IPSL-CM5A2</td><td>Laboratoire des Sciences du Climat et de l'Environnement - France</td><td>68</td></tr><tr><td>MIROC4m</td><td>University of Tokyo, NIES, FRCGC - Japan</td><td>69</td></tr><tr><td>MRI-CGCM2.3</td><td>University of Tsukuba - Japan</td><td>70</td></tr><tr><td>NorESM-L</td><td>Norwegian Research Centre - Norway</td><td>71</td></tr><tr><td>NorESM1-F</td><td>Norwegian Research Centre - Norway</td><td>71</td></tr></tbody></table>

For more information on how each modelling group has applied the PlioMIP boundary conditions and performed the numerical simulation, please see their respective referenced study.
<table><caption>Extended Data Table 2 | PlioMIP phases 1 and 2 boundary conditions</caption><thead><tr><th></th><th>LSM</th><th>Topo.</th><th>soils</th><th>lakes</th><th>Ice<br>Greenl.</th><th>Ice<br>Antarct.</th><th>vegetation</th><th>CO<sub>2</sub></th><th>Orbital<br>Parameters</th></tr></thead><tbody><tr><td>Pre-<br>industrial</td><td>Modern</td><td>Modern</td><td>Modern</td><td>Modern</td><td>Modern</td><td>Modern</td><td>Modern</td><td>280</td><td>Modern</td></tr><tr><td>PlioMIP1<br>(experiment<br>2)</td><td>Plio</td><td>Plio</td><td>Modern</td><td>Modern</td><td>Plio<br>(~50%<br>retreat)</td><td>Plio</td><td>Plio</td><td>405</td><td>Modern</td></tr><tr><td>PlioMIP2<br>(Eoi400)</td><td>Plio –<br>enhanced*</td><td>Plio –<br>enhanced*</td><td>Plio</td><td>Plio</td><td>Plio<br>(~70%<br>retreat)</td><td>Plio</td><td>Plio</td><td>400</td><td>Modern</td></tr></tbody></table>

LSM: land-sea mask. Topo: topography. Ice: land-ice. *MRI-CGCM2.3 was the only model that could not apply the enhanced reconstruction from PRISM4 and, thus, applied the previous, standard, reconstruction.