Similar North Pacific variability despite suppressed El Niño variability in the warm mid-Pliocene climate

The mid-Pliocene is the most recent geological period with similar atmospheric CO₂ concentration to the present day and similar surface temperatures to those projected at the end of this century for a moderate warming scenario. While not a perfect analogue, the mid-Pliocene can be used to study the functioning of the Earth system under similar forcings to a near future, especially regarding features in the climate system for which uncertainties exist in future projections. According to the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2), the variability in the El Niño–Southern Oscillation (ENSO) was suppressed. In this study, we investigate how teleconnections of ENSO, specifically variability in the North Pacific atmosphere, respond to a suppressed ENSO according to PlioMIP2. The multi-model mean (MMM) shows a similar sea-level pressure (SLP) variability in the Aleutian Low (AL) in the mid-Pliocene and pre-industrial, but a per-model view reveals that the change in AL variability is related to the change in ENSO variability. Even though ENSO is suppressed, the teleconnection between ENSO sea-surface temperature (SST) anomalies, tropical precipitation, and North Pacific SLP anomalies is quite robust in the mid-Pliocene. We split AL variability in a part that is ENSO-related, and a residual variability which is related to internal stochastic variability, and find that the change in ENSO-related AL variability is strongly related to the change in ENSO variability itself, while the change in residual AL variability is unrelated to ENSO change. Since the internal atmospheric variability, which is the dominant forcing of the AL variability, is largely unchanged, we are able to understand that the AL variability is largely similar even though ENSO variability is suppressed. We find that the specific change in ENSO and AL variability depends on both the model equilibrium climate sensitivity and Earth system sensitivity. Finally, we present a perspective of (extra-)tropical Pacific variability in PlioMIP2, combining our results with literature findings on changes in the tropical mean climate and in the Pacific Decadal Oscillation (PDO).

1 Introduction

The mid-Pliocene (∼3 Ma) is the most recent geological period that had a similar atmospheric CO₂ concentration to the present day (∼400 ppm; de la Vega et al., 2020). Its global surface temperatures were similar to those projected for the end of this century following the SSP2-4.5 future climate scenario (Tierney et al., 2020). Apart from elevated atmospheric CO₂, the mid-Pliocene had a similar geography – compared to earlier geological warm periods – to the present day (Haywood et al., 2011). Notable differences include a reduction in the Greenland and West Antarctic ice sheets; closure of the Bering Strait and Canadian Arctic Archipelago; and changes to topography, bathymetry, and vegetation (Dowsett et al., 2016). Because of these similarities, the mid-Pliocene climate has been called the “best analogue” to near-future climate in comparison to other past warm periods (Burke et al., 2018). Regardless of its potential as a future climate analogue, the mid-Pliocene can serve as a benchmark for climate models. The analysis of the climate system under mid-Pliocene boundary conditions can help us investigate the strength of feedbacks and the response of large-scale climate features to specific forcings.

The first coordinated effort to simulate the mid-Pliocene climate was the first phase of the Pliocene Model Intercomparison Project (PlioMIP; Haywood et al., 2010). It was followed by PlioMIP2 (Haywood et al., 2016), and the most recent phase, PlioMIP3, has just been announced (Haywood et al., 2024). PlioMIP2 was specifically designed to reduce model–proxy uncertainties. All models are supplied with consistent boundary conditions, which are focused on a specific time slice (KM5c interglacial, 3.205 Ma) where the orbital configuration was close to present-day conditions (Haywood et al., 2016; Dowsett et al., 2016).

In the past years, many studies of the mid-Pliocene climate based on PlioMIP2 have been published. In the global mean, the mid-Pliocene climate was 3.3 °C warmer than the pre-industrial climate (ranging between +1.7 and +5.2 °C) (Haywood et al., 2020). Increased CO₂ and closed gateways influenced the Arctic such that Arctic temperatures were relatively high due to Arctic amplification (De Nooijer et al., 2020). Most models simulate summer Arctic sea-ice-free conditions, and the Atlantic Meridional Overturning Circulation (AMOC) was intensified (De Nooijer et al., 2020; Weiffenbach et al., 2023). Atmospheric moisture content was increased, leading to more precipitation on average (Haywood et al., 2020) and specifically wetter conditions over the deep tropics, such as the Pacific Intertropical Convergence Zone (ITCZ, Han et al., 2021). The subtropics become drier over the ocean, but precipitation over land is generally enhanced related to enhanced monsoonal activity (Berntell et al., 2021; Feng et al., 2022). A hemispheric energy asymmetry shifts the Hadley circulation and ITCZ northwards (Han et al., 2021; Zhang et al., 2024). The Pacific Walker circulation is shifted westward (Han et al., 2021; Zhang et al., 2024), and this is associated with warmer and wetter conditions over the Indian Ocean and the Maritime Continent (Ren et al., 2023). Not all of these features are analogous to (near-)future climate projections; e.g. AMOC is projected to decrease, while the mid-Pliocene AMOC is simulated to be strengthened (Eyring et al., 2021; Weiffenbach et al., 2023).

A more puzzling feature of the mid-Pliocene climate is the behaviour of the El Niño–Southern Oscillation (ENSO). In the present day, ENSO is the dominant mode of climate variability on interannual timescales, with teleconnections to many regions of the world (Philander, 1990). Earlier work using proxy reconstructions showed a reduced zonal sea-surface temperature (SST) gradient in the tropical Pacific mean climate (Wara et al., 2005; Ravelo et al., 2006), on the basis of which it was suggested that the mid-Pliocene ENSO was settled in a “permanent El Niño” state (Fedorov et al., 2006). Proxies for variability in the mid-Pliocene do suggest, however, that there existed ENSO variability with an amplitude varying between reduced and similar to the present day (Scroxton et al., 2011; Watanabe et al., 2011; White and Ravelo, 2020). In addition, more recent reconstructions and modelling efforts suggest that the zonal SST gradient is not as reduced as previously thought and that it could be in line with model estimates (Zhang et al., 2014; Tierney et al., 2019). Modelling efforts in PlioMIP1 (Brierley, 2015) and PlioMIP2 (Oldeman et al., 2021) agree with that finding, showing that ENSO variability was reduced in the majority of the models but with considerable spread in the model ensemble. The suppression of ENSO in PlioMIP2 is explained by a series of off-equatorial processes triggered by the northward displacement of the Pacific ITCZ (Pontes et al., 2022).

ENSO variability exerts a global influence through its oceanic and atmospheric teleconnections, which include a circumglobal connection along the tropics and links with the stratosphere (Yeh et al., 2018; Cai et al., 2019; Domeisen et al., 2019). The so-called “atmospheric bridge” explains the deterministic link between ENSO and the variability in the North Pacific atmosphere and ocean, in which SST anomalies originating from ENSO events cause extratropical atmospheric variability via tropical convection, the Hadley circulation, and atmospheric Rossby waves (Hoskins and Karoly, 1981; Mo and Livezey, 1986; Alexander et al., 2002). The variability in the North Pacific atmosphere is often referred to as Aleutian Low (AL) variability, which consists of the dominant Pacific–North American (PNA) pattern and the second-leading North Pacific Oscillation (NPO) (Barnston and Livezey, 1987; Linkin and Nigam, 2008). AL variability is in part forced by ENSO through the deterministic atmospheric bridge but also by internal stochastic variability (Di Lorenzo et al., 2013; Newman et al., 2016). The dominant mode of ocean variability in the North Pacific is the Pacific Decadal Oscillation (PDO), which is forced by AL variability through wind forcing and by ENSO through ocean waves (Newman et al., 2003, 2016; Di Lorenzo and Mantua, 2016).

What may happen to ENSO and its teleconnections to the North Pacific in the near future under global warming is unclear. It is likely that ENSO–precipitation variability will increase (Cai et al., 2021; Yun et al., 2021) and that the variability in ENSO and atmospheric teleconnections, including AL variability, will increase in the near future (Chen et al., 2018; Fredriksen et al., 2020; Cai et al., 2021). However, uncertainties are very large, in part due to internal variability, and conclusions become even less robust towards the end of this century (Fredriksen et al., 2020; Beobide-Arsuaga et al., 2021). Additionally, ENSO teleconnections can also change because mean atmospheric circulation will change, regardless of ENSO change (Yeh et al., 2018). In the long term, idealized future warming simulations under equilibrated high CO₂ forcing, however, suggest a weakening of ENSO variability (Callahan et al., 2021; Zheng et al., 2022). This is similar to what is found in PlioMIP2 (Oldeman et al., 2021; Pontes et al., 2022), implying that the mid-Pliocene ENSO response is similar to what could be expected in an equilibrated high-CO₂ future but is not similar to the near-future ENSO response. This makes the mid-Pliocene climate a valuable test case to investigate the response of North Pacific variability to a suppressed ENSO.

Although ENSO variability is reduced compared to the pre-industrial in the majority of PlioMIP2 models (Oldeman et al., 2021), the amplitude of the PDO variability does not change, at least in a subset of the PlioMIP2 ensemble, making the PDO the dominant mode of ocean variability in the Pacific sector (Katya Canal-Solis, personal communications, 2024). However, there is a considerable model spread regarding both ENSO and PDO change in PlioMIP2. In one PlioMIP2 model, it is the mid-Pliocene boundary conditions that cause an ENSO suppression and substantial suppression of the AL variability (specifically the PNA pattern), while elevated CO₂ causes increased AL variability without any change in ENSO (Oldeman et al., 2023). It could indicate a tug of war in ENSO teleconnection responses to different conditions of the mid-Pliocene. However, the model used in Oldeman et al. (2023) showed the largest ENSO reduction of PlioMIP2, so it might not be representative of the rest of the ensemble.

In this study, we aim to answer the following question: how does variability in the North Pacific atmosphere respond to a suppressed ENSO in the warm mid-Pliocene climate, according to the PlioMIP2 ensemble? Specifically, we want to know the following. (1) Is any change in AL variability related to the change in ENSO? (2) Does the ENSO–North Pacific atmosphere teleconnection strength change? (3) Are there any changes in the North Pacific variability not related to ENSO change? And (4) is the ENSO teleconnection response related to model–climate sensitivity? To answer these questions, we will investigate simulation results of the PlioMIP2 ensemble. More details on the simulations and the analysis methods are presented in Sect. 2. Section 3 presents the results and provides answers to the research questions. In Sect. 4, we discuss the results and interpret our findings. We conclude with a summary in Sect. 5.

2 Methods

2.1 Models, simulations, and data

2.1.1 The PlioMIP2 ensemble

The PlioMIP2 ensemble consists of 17 climate models (Table 1). All models performed simulations following the PlioMIP2 protocol, which includes a pre-industrial control reference simulation (E280) and a mid-Pliocene simulation (Eoi400) (Haywood et al., 2016). Pre-industrial simulations are forced with 280 ppm atmospheric CO₂, while the mid-Pliocene simulations are forced with 400 ppm CO₂. Altered mid-Pliocene boundary conditions include closed Arctic gateways (Bering Strait and Canadian Arctic Archipelago); reduced land-ice cover (Greenland ice sheet and West Antarctic Ice Sheet); shoaling of the Sunda and Sahul shelves leading to changes to straits in the Maritime Continent; and changes to vegetation, lakes, and soils (Dowsett et al., 2016; Haywood et al., 2016). Only HadGEM3 used a modern land–sea mask in the Eoi400 simulation (Williams et al., 2021). All model simulations were run for 1000 or more model years (following PlioMIP2 protocol) and can be regarded as being in climatological equilibrium. We use the last 100 years from each simulation.

Feng et al. (2020)Chandan and Peltier (2017)Baatsen et al. (2022)Feng et al. (2020)Feng et al. (2020)Stepanek et al. (2020)Zheng et al. (2019)Hunter et al. (2019)Williams et al. (2021)Tan et al. (2020)Tan et al. (2020)Lurton et al. (2020)Chan and Abe-Ouchi (2020)Kamae et al. (2016)Li et al. (2020)Li et al. (2020)

Table 1 Details of the models contributing to the PlioMIP2 ensemble, with equilibrium climate sensitivity (ECS) and Earth system sensitivity (ESS) from Haywood et al. (2020) and Williams et al. (2021). In bold are the models primarily used in this study. The asterisks (^*) show when the model has sensitivity simulations available (either E400 or Eoi280).

Additionally, we use two different sensitivity simulations that are available for a subset of PlioMIP2 models. We consider simulations with mid-Pliocene boundary conditions (BCs) but at pre-industrial CO₂ levels (Eoi280), ran by CCSM4-Utr, COSMOS and HadCM3. We also consider simulations with mid-Pliocene CO₂ levels but pre-industrial conditions otherwise (E400), ran by CCSM4-UoT, COSMOS and HadCM3.

2.1.2 Observational data

We compare the pre-industrial results to observational products of the present day and the historical period. For SSTs, we use the NOAA Extended Reconstructed SST (ERSST) v5 dataset (Huang et al., 2017). For atmospheric variables, we use the NOAA 20th Century Reanalysis (20CR) project v3 (Slivinski et al., 2021). We will refer to both products as “NOAA observations”. For consistency, we use 100 years, namely 1916–2015. These NOAA products have data available from the 19th century, which might be more like the pre-industrial, but the caveat is that the data become less reliable, so we chose to use recent data with higher reliability and consistency instead.

2.2 Analysis methods

We use 100 years of monthly SST, sea-level pressure (SLP), and total precipitation fields. For anomalies, we remove the climatology and a linear trend. For NOAA observations, we instead perform a LOWESS filtering (50-year running mean) to remove the trend from anthropogenic climate change. For model data, performing LOWESS filtering instead of removing a linear trend did not make any difference to the results. When computing multi-model means (MMMs) in space, we interpolate the model data on a rectilinear ∼1° grid.

We study variability through defined climate indices. We study ENSO variability through SST anomalies in the Niño 3.4 region (5° S–5° N, 150–90° W), which was also shown to capture ENSO variability well in the mid-Pliocene (Oldeman et al., 2021). To study Aleutian Low (AL) variability, we use SLP anomalies in the AL region (30–65° N, 160° E–140° W), also known as the North Pacific Index (Trenberth and Hurrell, 1994; Chen et al., 2020). This region also captures AL variability in the mid-Pliocene simulations, as shown in the Results section. We study precipitation variability in the western equatorial Pacific (WEP), which is the region that has the strongest ENSO-related precipitation anomalies (Deser and Wallace, 1990; Williams et al., 2024), defined here as 6° S–6° N, 120–180° E following Oldeman et al. (2023). In all cases we take area averages using grid weights based on the cosine of the latitude. Amplitudes of variability are generally defined as the standard deviation (SD) of the associated time series.

In order to connect ENSO teleconnections to other modes of variability, we compute the linear regression (linear slope) and correlation between the climate indices, where we quantify the teleconnection strength by means of the linear regression (e.g. Williams et al., 2024). Correlations are assumed statistically significant only if the corresponding p-value is below 0.05. Since the ENSO signal is strongest in boreal winter, we focus on the DJF mean Niño 3.4 index. The present-day ENSO leads AL variability, and we focus on the AL variability in DJFM. These months are chosen because the atmospheric response can take ∼1–2 weeks (e.g. Trenberth and Hurrell, 1994; Alexander et al., 2002; Newman et al., 2016). We test whether this assumption is valid by computing the lead–lag correlations between ENSO and the AL time series, and we find that ENSO leads AL with a 0–1 month lag with strongest correlations around January Niño 3.4 for NOAA observations and the majority of the models (Fig. S1 in the Supplement).

The deterministic teleconnection linking ENSO and AL variability, the atmospheric bridge, consists of several steps. ENSO SST anomalies leading to tropical convection anomalies (i.e. precipitation anomalies) are an important precursor to extratropical SLP anomalies in the North Pacific. Examining this step is relevant considering the substantial changes in the Indo-Pacific mean hydrological cycle in PlioMIP2 (e.g. Han et al., 2021; Ren et al., 2023; Zhang et al., 2024). Therefore, we also consider the regression between the Niño 3.4 index and tropical precipitation and the regression between WEP precipitation and North Pacific SLP anomalies.

Since we are interested in the ENSO–AL teleconnection change in the mid-Pliocene, we check whether the PlioMIP2 models are able to simulate that connection well in the pre-industrial. We find that IPSL-CM6A, MIROC4m, and MRI-CGCM2.3 do not show statistically significant correlations between the Niño 3.4 and AL indices in the E280 simulation in any relevant combination of months or in any subsection of the AL region. Hence, we do not use the results of IPSL-CM6A, MIROC4m, and MRI-CGCM2.3 in this study. A more detailed justification of this omission is included in the Supplement (text and Figs. S1–S3).

To separate the AL variability in a part that is related to ENSO and in a part that is not related to ENSO, we split the total AL variability (A_tot) by following a linear regression model (LRM; see e.g. Chiang and Vimont, 2004; Deser et al., 2017):

$$\begin{array}{}\text{(1)}& A_{\text{tot}}\left(t\right)=A_{\mathrm{N}}\left(t\right)+A_{\text{res}}\left(t\right),\end{array}$$

in hPa, where A_N represents the part of the AL variability that linearly regresses with the Niño 3.4 index and A_res represents any residual variability:

$$\begin{array}{}\text{(2)}& A_{\mathrm{N}}\left(t\right)={\mathit{\beta }}_{\mathrm{N},A_{\text{tot}}}\cdot N\left(t\right),\end{array}$$

where N is the Niño 3.4 index in °C and ${\mathit{\beta }}_{\mathrm{N},A_{\text{tot}}}$ is the linear regression between the Niño 3.4 index and the AL index in hPa °C⁻¹. We compute A_tot and A_N, and A_res then follows from the LRM. Ordinary least squares ensures that the LRM is constructed such that the time variance of the total AL variability is the sum of the variance of the part of the AL variability that regresses with the Niño 3.4 index (Niño-regr. part) and the variance of the residual:

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{A_{\text{tot}}}^{\mathrm{2}}={\mathit{\sigma }}_{A_{\mathrm{N}}}^{\mathrm{2}}+{\mathit{\sigma }}_{A_{\text{res}}}^{\mathrm{2}},\end{array}$$

in hPa². The Niño-regr. part represents the part of the AL variability that covaries linearly with ENSO and can be seen as the part of the variability that is explained or caused by ENSO variability. By definition, the Niño-regr. part of the AL variability and the residual AL variability are uncorrelated. The residual thus represents that part of the AL variability that is either nonlinearly related to ENSO or does not covary with ENSO at all. This last part could be any internal stochastic variability (e.g. variability related to the jet streams, sea-ice cover, or Arctic Oscillation).

Regarding model–climate sensitivities, we use both equilibrium climate sensitivity (ECS) and Earth system sensitivity (ESS). ECS is defined as the global mean surface temperature response to a doubling of CO₂ with pre-industrial boundary conditions once the energy balance has reached equilibrium. ESS is defined as the temperature response to a CO₂ doubling and to other forcing changes – in other words, including responses to feedbacks with long timescales such as those involving ice sheets. ESS is relevant in the context of (past) climates where there are more changes in forcings than elevated CO₂. We obtain the values for ECS and ESS from Haywood et al. (2020) and Williams et al. (2021), which are included in Table 1.

Using the sensitivity simulations, we define a fraction of the total ENSO–AL response (FoR) in the mid-Pliocene, which is due to elevated CO₂, as follows:

$$\begin{array}{}\text{(4)}& {\text{FoR}}_{{\mathrm{CO}}_{\mathrm{2}}}=d_{{\mathrm{CO}}_{\mathrm{2}}}/(d_{{\mathrm{CO}}_{\mathrm{2}}}+d_{\text{BCs}}),\end{array}$$

where $d_{{\mathrm{CO}}_{\mathrm{2}}}$ and d_BCs are the distances in terms of relative ENSO change and relative AL change due to elevated CO₂ and due to the mid-Pliocene BCs, respectively. Using sensitivity simulation E400 (elevated CO₂ with pre-industrial BCs), we compute the distances as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{d}_{{\mathrm{CO}}_{\mathrm{2}}}=\\ {\displaystyle }& {\displaystyle }\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\text{N, E400}}-{\mathit{\sigma }}_{\text{N, E280}}}{{\mathit{\sigma }}_{\text{N, E280}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\text{A, E400}}-{\mathit{\sigma }}_{\text{A, E280}}}{{\mathit{\sigma }}_{\text{A, E280}}}}\right)}^{\mathrm{2}}}\end{array}\text{(6)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{d}_{\text{BCs}}=\\ {\displaystyle }& {\displaystyle }\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\text{N, Eoi400}}-{\mathit{\sigma }}_{\text{N, E400}}}{{\mathit{\sigma }}_{\text{N, E400}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\text{A, Eoi400}}-{\mathit{\sigma }}_{\text{A, E400}}}{{\mathit{\sigma }}_{\text{A, E400}}}}\right)}^{\mathrm{2}}},\end{array}\end{array}$$

where σ_N,x and σ_A,x are the amplitudes (standard deviation) of the Niño 3.4 index and AL index, respectively, in simulation “x”. In the case of sensitivity simulation Eoi280, the $d_{{\mathrm{CO}}_{\mathrm{2}}}$ and d_BCs are determined using the differences between Eoi400−Eoi280 and Eoi280−E280, respectively. If the change in AL would be zero throughout all simulations, the distances simply reduce to the relative ENSO change, and the FoR becomes a fraction of the ENSO change between the set of simulations.

3 Results

3.1 AL variability and ENSO change

Figure 1 shows the DJFM SLP SD for NOAA observations (a), the E280 multi-model mean (MMM; b), the Eoi400 MMM (c), and the difference between Eoi400 and E280 (d). The pre-industrial MMM reproduces the spatial pattern of the SLP SD of the NOAA observations well but overestimates the amplitude substantially, which is a known model bias (Chen et al., 2018; Eyring et al., 2021). The mid-Pliocene MMM is not substantially different from the pre-industrial MMM, implying that in all cases almost all of the North Pacific atmospheric variability is captured in the AL region. The MMM of mid-Pliocene minus pre-industrial differences is small; furthermore, more than half of the models in the ensemble do not agree on the sign of change in the AL region. The largest and most consistent change across the ensemble is a reduction in SLP variability along the North American west coast. Even there, the maximum change corresponds to approximately −10 %.

Figure 1 Standard deviation (SD) of the DJFM mean sea-level pressure (SLP) for (a) NOAA observations, (b) the E280 multi-model mean (MMM), (c) the Eoi400 MMM, and (d) the difference between the Eoi400 MMM and the E280 MMM. The cyan rectangle indicates the Aleutian Low region. The arcing in panel (d) shows when more than 7 out of 14 models agree with the sign of change.

While the MMM suggests no substantial change in AL variability, a per-model look reveals a variety of responses which are related to the ENSO amplitude change. Figure 2 shows scatter plots of the ENSO amplitude (defined as Niño 3.4 SD) and AL amplitude in the pre-industrial E280 (a), the mid-Pliocene Eoi400 (b), and as the difference between the simulations (c). The ENSO amplitude in the NOAA observations is reproduced well by the E280 ensemble, while the AL amplitude is overestimated by the majority of the models. The Eoi400 ensemble shows a large spread in both ENSO and AL amplitudes compared to the E280. In the E280 and Eoi400, the AL amplitude does not correlate well with the ENSO amplitude (i.e. no statistically significant ensemble correlation, p-value>0.05). However, the change in AL amplitude in the ensemble is related to the change in ENSO amplitude with a statistically significant correlation coefficient of 0.76. Generally, models with slight ENSO change show a similar or increased AL variability in the mid-Pliocene, while the models with substantial ENSO reduction show a similar or reduced AL variability. There is a considerable model spread regarding this relation, though, implying that the change in AL variability is not just related to ENSO change.

Figure 2 Scatter plots of DJF Niño SD versus DJFM AL SD in the E280 (a) and Eoi400 (b) and a scatter plot of Eoi400−E280 differences in Niño SD and AL SD (c). The red triangles show the multi-model mean (MMM), and the blue triangle shows the NOAA observations. The dashed red line in panel (c) is a linear fit through the points, with a significant correlation coefficient of 0.76.

3.2 ENSO teleconnection with the North Pacific atmosphere

Figure 3 shows the linear regression, representing the ENSO–AL teleconnection strength, for the E280 and Eoi400 (a) and the relative change between the two (b). Results per model are included in the Supplement (Fig. S3). The majority of the models overestimate the regression in the pre-industrial when compared to NOAA observations, which is expected, since most models overestimate AL variability, while ENSO variability is similar to observations and correlations are similar too (Fig. S3). The large spread in modelled ENSO–AL teleconnection strength compared to observations is also reported in CMIP5 generation models (e.g. Deser et al., 2017). The MMM shows a similar regression in the pre-industrial and mid-Pliocene. The Eoi400 regression for CCSM4-Utr and EC-Earth3.3 is set to zero because the correlation between ENSO and AL becomes statistically insignificant in the mid-Pliocene simulation (i.e. p-value>0.05). Not doing this would not change the MMM substantially because the correlations are weak nonetheless. Figure 2b shows that the MMM (dashed red line) and median (red line in boxplot) values are slightly positive. However, the confidence intervals around the median (boxplot notches) encompass 0 %, meaning that we cannot statistically speak of an increase in regression. An ensemble end-member is COSMOS, showing the largest increase in regression strength, whereas CCSM4-Utr and EC-Earth3.3 show the largest decrease. The change in regression is weakly but significantly correlated to the change in ENSO and the change in AL across the ensemble (Fig. S4).

Figure 3 (a) Regression (linear slope) between DJF Niño and DJFM AL in the E280 and Eoi400 per model. The red triangle is the multi-model mean (MMM), and the dash-dotted blue line shows the NOAA observations. The Eoi400 value of CCSM4-Utr and EC-Earth3-LR is set to zero because there is no significant correlation between the two variables (i.e. p-value>0.05). (b) Relative change in this regression per model, with the dashed red line indicating the multi-model mean. A boxplot is included (the red line is the median).

Figure 4 shows the E280 MMM regression between the Niño 3.4 index and tropical Pacific precipitation (b) and the WEP precipitation index and SLP in the North Pacific (a). Regressions are only shown when the majority of the models have a statistically significant correlation (i.e. p-value<0.05) in that grid cell. Pre-industrial ENSO variability leads to a strong precipitation signal in the WEP (b; box drawn), and precipitation anomalies in that region lead to a strong SLP signal in the AL region (a; box drawn). The shape and amplitude of the regression patterns look similar to NOAA observations (Fig. S5) and present-day simulations (Williams et al., 2024). Figure 4c shows the change in ENSO–WEP precipitation regression and the change in WEP precipitation–AL regression. The MMM indicates no change in both regressions. It agrees with the previous finding that the ENSO–AL regression does not change, and it shows that this is not because of concealed counteracting changes in the teleconnection processes. Figure 4c also shows that the strong increase in ENSO–AL regression in COSMOS is due to both the ENSO–WEP precipitation and the WEP precipitation–AL regression strengthening. For CCSM4-Utr, the weaker regression is mainly due to the WEP precipitation–AL link weakening, while the weaker regression in EC-Earth3.3 is mainly related to a weakening of the ENSO–WEP precipitation link.

Figure 4 (a) E280 multi-model mean (MMM) regression between DJF precipitation in the western equatorial Pacific (WEP) and DJFM SLP in the North Pacific. The cyan rectangle indicates the Aleutian Low region. (b) E280 MMM regression between DJF Niño 3.4 and DJF precipitation. The cyan rectangle indicates the WEP region. For both panels (a) and (b), values are only shown if more than 7 out of 14 models have a statistically significant correlation (p-value<0.05). (c) Scatter plot of Eoi400−E280 change in regression between Niño 3.4 and WEP precipitation versus change in regression between WEP precipitation and AL. The red triangle is the MMM. Regressions in either E280 or Eoi400 are set to 0 when correlations are not statistically significant (i.e. p-value>0.05)

3.3 Separating ENSO- and non-ENSO-related AL variability

We split the AL variability in a part that regresses with ENSO (Niño-regr. part) and a residual following the LRM as explained in Sect. 2. The variances of both parts add up to the total AL variance. Figure 5 shows the Niño-regressing AL variance, the residual AL variance, and the total AL variance for all models and simulations. In the MMM, the Niño-regr. AL variance decreases in the mid-Pliocene, which can be understood, since the ENSO amplitude reduces (Fig. 2c), while the ENSO–AL teleconnection strength does not change (Fig. 3). The AL residual variance does not change in the MMM. The pre-industrial MMM overestimates both parts of the AL variance compared to NOAA observations. It also overestimates the AL variance fraction related to ENSO (0.30 over 0.18), which can be explained, since the MMM regression is overestimated (Fig. 3a), whereas the ENSO amplitude is similar (Fig. 2a). CMIP6 generation models have been shown to generally underestimate the tropical influence on variability in the North Pacific (Zhao et al., 2021), which in this case is true for 4 out 14 models. In the MMM, the fraction of the AL variance related to ENSO decreases from 0.30 to 0.22 in the mid-Pliocene, and 10 out of 14 models agree with that sign of change in variance fraction. For 9 models this is because the ENSO–AL variance decreases. For NorESM1-F the variance fraction decreases because the residual variance increases, while for NorESM-L the variance fraction increases despite a slight reduction in the ENSO-related variance, caused by a stronger reduction in the residual variance. The ENSO-related AL variance is zero for CCSM4-Utr and EC-Earth3.3 because the ENSO–AL regression is set to zero.

Figure 5 Aleutian Low (AL) variance, split (following a linear regression model; LRM) into the AL variance that regresses with Niño (Niño-regr. part; hatched) and the residual AL variance. Results per model for the E280 (blue) and Eoi400 (red), including the multi-model mean (MMM) and the result from NOAA observations. The full bar length represents the total AL variance. Values below are the variance fraction of the Niño-regr. part for E280 (blue, top) and Eoi400 (red, bottom) and are in bold when the variance fraction is lower in the Eoi400 compared to in the E280.

Figure 6 shows the change in the AL amplitude (in terms of SD), split into the part that regresses with ENSO (a) and the residual variability (b) as a function of the change in ENSO amplitude (similar to Fig. 2c). By using the LRM to split the AL variability, we are able to separate the change in AL variability into a part regressing with ENSO where the change is strongly dependent on the change in ENSO variability (corr. coef. of 0.89 over 0.76 for the total AL change; Fig. 2c) and a residual where the change is not related to ENSO change (corr. becomes insignificant). We could expect the ensemble correlation in Fig. 6a to be higher than the ensemble correlation in Fig. 2c if the linear regression between ENSO and the AL were the same between the pre-industrial and the mid-Pliocene. While the MMM regression is largely unchanged (Fig. 3b), the regression change per model can be substantial, implying that the correlation in Fig. 6a is not necessarily higher merely by construct. In the MMM, there is a slight reduction in the ENSO-related AL variability, where 10 of the 14 models agree with that sign of change. The MMM residual AL variability shows no change, with 7 models showing a reduction and 7 showing an increase. The results indicate that the change in total AL variability is primarily driven by a change in ENSO. The residual AL variability does change slightly per model (mean absolute error (MAE) of 0.47 hPa) but on average less than the ENSO-related part (MAE of 0.79 hPa). Since the residual variability is similar and dominates the total AL variability (see Fig. 5), the total AL variability does not seem to change much in the mid-Pliocene along the ensemble, even though ENSO is suppressed.

Figure 6 Scatter plots of the change in the Eoi400−E280 DJF Niño 3.4 SD versus (a) the change in DJFM AL variability that regresses with Niño (Niño-regr. part) and (b) the change in the residual DJFM AL variability. The red triangle is the multi-model mean (MMM). The dashed red line in panel (a) is a linear fit through the points, with a significant correlation coefficient of 0.89. Multi-model correlation in panel (b) is not statistically significant (i.e. p-value>0.05).

3.4 ENSO and AL variability response in relation to climate sensitivities

Oldeman et al. (2023) use sensitivity simulations to show that CCSM4-Utr's Eoi400 North Pacific variability response is largely dictated by the response to the mid-Pliocene BCs (e.g. closed Arctic gateways, reduced ice sheets) and not by the response to elevated CO₂. They hypothesize that this is related to the relatively high sensitivity of the model to the mid-Pliocene BCs compared to its sensitivity to elevated CO₂ in terms of global mean temperature response and compared to the PlioMIP2 ensemble. In this section we will investigate this issue in more detail by looking at the ENSO change and the AL change in response to mid-Pliocene boundary conditions and elevated CO₂ using sensitivity simulations from a subset of the PlioMIP2 models.

Figure 7a shows the relative change in ENSO variability and AL variability in response to elevated CO₂ (in red), mid-Pliocene BCs (in blue), and the “total” mid-Pliocene (i.e. the Eoi400 simulation, in black). There is a clear separation of responses to both forcings. In response to elevated CO₂, five out of six simulations show increased AL variability, of which four also show slightly increased ENSO variability. In response to the mid-Pliocene BCs, five out of six simulations show reduced ENSO variability, of which four also show reduced AL variability. It indicates a “tug of war” between the opposite ENSO–AL responses to either mid-Pliocene BCs or elevated CO₂, where the total mid-Pliocene (Eoi400) response differs per model, explaining the ensemble spread in Fig. 2c. We argue that whichever forcing “wins” this tug of war is related to the model sensitivity. Figure 7b shows the fraction of the ENSO–AL response which is due to elevated CO₂ as a function of the ratio between the model sensitivity to elevated CO₂ (i.e. ECS) and the total mid-Pliocene climate sensitivity (ESS). It indicates that a positive relation exists between the $\text{ECS}/\text{ESS}$ ratio and the fraction of the ENSO–AL response due to elevated CO₂. The larger the relative model sensitivity to CO₂, the more of the ENSO–AL response is related to elevated CO₂, and likewise the larger the relative model sensitivity to the mid-Pliocene BCs, the more of the ENSO–AL response is related to the mid-Pliocene BCs. This relationship is not necessarily intuitive; both ECS and ESS are a measure of the annual global mean surface temperature change in response to a certain (combination of) forcing(s), while the fraction of ENSO–AL response is related to the change in ENSO and AL variability in boreal winter due to different forcings. It demonstrates that the (relative) sensitivity of a climate model to a specific forcing is connected with more responses of the climate system than just global mean temperature change.

Figure 7 (a) Scatter plot of relative change in DJF Niño 3.4 SD vs. relative change in DJFM AL variability for a subset of the PlioMIP2 ensemble and using sensitivity studies. The total mid-Pliocene response is in black (Eoi400−E280), the response to elevated CO₂ is in red (white markers for E400−E280 and filled markers for Eoi400−Eoi280), and the response to mid-Pliocene BCs is in blue (white markers for Eoi400−E400 and filled markers for Eoi280−E280). The Eoi400−E280 responses of remaining PlioMIP2 models are denoted by small crosses. Note the different markers used per model compared to previous figures. (b) Fraction of ENSO–AL response (FoR) due to elevated CO₂ as a function of the ratio between ECS and ESS. FoR is defined using the distances (in terms of ENSO change and AL change) between the different simulations (see Methods section for details). The ECS and ESS are listed in Table 1.

4 Discussion

4.1 Performance of PlioMIP2 with regard to observations and ensemble spread

In this section we discuss the performance of the PlioMIP2 ensemble in assessing changes to the ENSO–AL teleconnection in the mid-Pliocene in terms of comparing pre-industrial results with NOAA observations and the within-ensemble performance. The PlioMIP2 models are generally better at reproducing the amplitude of ENSO variability compared to NOAA observations (Fig. 2a; for more detail see Oldeman et al., 2021). The model that best reproduces the NOAA observations in the pre-industrial is NorESM-L, which also shows changes to ENSO and AL variability in the mid-Pliocene that are very close to the MMM (Fig. 2c). Most models overestimate the ENSO–AL regression strength (Fig. 3a), which is not surprising considering that the amplitude of ENSO variability is similar, that the amplitude of AL variability is overestimated, and that the correlation coefficient between ENSO and AL indices is similar for most E280 simulations and NOAA observations (Fig. S3). The shape and amplitude of the pre-industrial MMM ENSO–precipitation relation (Fig. 4b) is similar to NOAA observations and is similar to the response in present-day simulations with higher-resolution models (Williams et al., 2024). The fraction of AL variance related to ENSO is overestimated in the E280 MMM compared to NOAA observations (Fig. 5). An important feature is that the AL variance fraction of the residual is larger than the fraction of AL variance related to ENSO, which is captured by all E280 simulations except for CCSM4-Utr. Finally, it should be noted that any discrepancies between NOAA observations and the pre-industrial simulations can also arise from the fact that we are comparing equilibrated pre-industrial simulations with historical observations. Even though the observations have been detrended, there will be an anthropogenic signal present which is absent from the simulations considered.

In terms of ENSO and AL change (e.g. Figs. 2c, 6, and 7a), the clear end-members of the ensemble are CCSM4-Utr on the one hand, showing the largest reduction in ENSO and AL variability, and COSMOS on the other hand, showing the largest increase in ENSO and AL variability. Our results show that this spread in responses is related to the relative model sensitivities of these two models to either the Pliocene BCs or the elevated CO₂ (Fig. 7b). COSMOS is one of the coarsest models in terms of ocean and atmosphere resolution (Table 1). Its E280 simulation performs well in terms of AL variability compared to NOAA observations, but the ENSO variability is greatly overestimated. Williams et al. (2024) show improved simulation of the ENSO–AL teleconnection with increased model resolution. CCSM4-Utr, on the other hand, performs well in terms of ENSO variability but substantially overestimates AL variability compared to NOAA observations. However, as shown in Oldeman et al. (2023), the patterns of North Pacific atmospheric variability (i.e. the spatial patterns, amplitudes, and variance fractions of the PNA pattern and NPO) are very well reproduced, even though the total amplitude is overestimated.

4.2 Changes in the residual Aleutian Low variability

We find that, in the mid-Pliocene, a change in AL variability is related to a change in ENSO variability (Fig. 2c). This relation originates from a change in AL variability which linearly covaries with the ENSO signal (Fig. 6a). The residual AL variability also changes, but the change is model-dependent and not related to a change in ENSO variability (Fig. 6b). The change in residual AL variability is also not related to a change in residual WEP precipitation (Fig. S6). In this section, we explore the residual Aleutian Low variability in more detail and hypothesize what its change might be related to.

4.2.1 Nonlinear ENSO interactions

We separated the AL variability in an ENSO-related and residual part following an LRM. This implies that the residual AL variability will also contain nonlinear ENSO contributions. Any nonlinear atmospheric response to ENSO SST variability can be regarded as a sum of nonlinear responses to a linear ENSO and of linear responses to a nonlinear ENSO (Frauen et al., 2014). ENSO variability is known to be nonlinear, originating from diversity in ENSO events and from asymmetry in pattern and duration between El Niño and La Niña events (e.g. Ashok et al., 2007; Yeh et al., 2009; Okumura, 2019; Cai et al., 2021). In the PlioMIP2 MMM, both the skewness of the Niño 3.4 index and the ratio between central Pacific and eastern Pacific ENSO events are unchanged compared to the pre-industrial (Oldeman et al., 2021). However, there is a considerable model spread, both in terms of change in the mid-Pliocene and of differences in observations. The atmospheric response to ENSO variability is complex and known to encompass nonlinearity originating from a variety of factors (e.g. Yeh et al., 2018; Domeisen et al., 2019; Jiménez-Esteve and Domeisen, 2019), the main factor being the nonlinear relationship between ENSO variability and tropical precipitation (Deser and Wallace, 1990). This nonlinear ENSO–precipitation relationship changes in a warming climate (Yun et al., 2021). Only about half of the PlioMIP2 models capture the nonlinear nature of the ENSO–precipitation relationship (Pontes et al., 2022). In conclusion, since there is considerable model spread in changes in both ENSO skewness and kurtosis (Oldeman et al., 2021) and in the ENSO–precipitation relation (Pontes et al., 2022, and Fig. 4c of this study), nonlinearity in the atmospheric response to ENSO could explain some of the residual AL variability, but the exact contribution is likely model-dependent. Considering ENSO diversity and its nonlinear interactions, removing ENSO influence from North Pacific variability using linear regression might not always be the best option (Zhao et al., 2021). Improving the separation of (nonlinear) ENSO influence from the North Pacific variability response in PlioMIP2 is outside the scope of this work.

4.2.2 Stochastic internal variability

Apart from the deterministic link with the tropics via the atmospheric bridge, AL variability is forced by internal stochastic variability originating from the extra-tropical atmosphere (Di Lorenzo et al., 2013; Newman et al., 2016). Indeed, the existence of modes of winter variability in the Northern Hemisphere atmosphere can largely be explained by internal or stochastic variability (Branstator, 2002; Branstator and Selten, 2009). Idealized models with no ENSO variability or ocean dynamics show that the dominant SLP variability in the North Pacific is located in the Aleutian Low region, resulting from internal atmospheric dynamics (Pierce, 2001; Alexander, 2010). Modes of variability in the Northern Hemisphere extratropics such as the NAO and Arctic Oscillation, but also Pacific modes such as the PNA pattern, are related to variability in the jet streams (e.g. Ambaum et al., 2001; Linkin and Nigam, 2008) and are related to each other through a circumglobal teleconnection (e.g. Barnston and Livezey, 1987; Branstator, 2002). It suggests that changes in atmospheric circulation originating from the mid-Pliocene boundary conditions (e.g. the reduced Greenland ice sheet) can ultimately lead to changes in atmospheric variability in the North Pacific. Research on the Last Glacial Maximum atmosphere indeed shows a distorted pattern of variability in the North Pacific related to jet stream changes originating from changes in ice sheet extent and not related to ENSO (Hu et al., 2020). Apart from the internal atmospheric dynamics, North Pacific winter atmospheric variability has been linked to variations in sea-ice extent and North Pacific SSTs (e.g. Linkin and Nigam, 2008; Hurwitz et al., 2012; Simon et al., 2022). Across the PlioMIP2 ensemble, changes in mid-Pliocene sea-ice extent (via De Nooijer et al., 2020) or North Pacific SSTs do not relate to changes in the residual AL variability (Fig. S7). Garfinkel et al. (2020) show that patterns in the North Pacific wintertime atmosphere can be explained by a sum of forcings related to the land–sea contrast, heat fluxes in the ocean, and topography. In the mid-Pliocene simulations, the land–sea contrast is different due to changes in the land–sea mask and in vegetation and lakes, and the topography is different, albeit not substantially. Differences in ocean heat fluxes could originate from the closure of the Arctic gateways, which has been shown to cause a strengthening of the mid-Pliocene AMOC (Weiffenbach et al., 2023). Changes in the AMOC and changes in the Arctic surface temperatures do not relate to changes in the residual AL variability (Fig. S7). Considering that every climate model will have different sensitivities of its internal atmospheric variability response to the different mid-Pliocene boundary conditions, it might be hard, if not impossible, to find one common driver that is able to explain changes to the residual AL variability in the mid-Pliocene.

4.3 Synthesis of tropical–North Pacific variability in PlioMIP2

In the last part of the Discussion, we present a synthesis of the results presented in this paper on the mid-Pliocene changes in North Pacific atmosphere variability combined with other published results using PlioMIP2 data, specifically regarding tropical–North Pacific variability (Oldeman et al., 2021; Pontes et al., 2022) and Indo-Pacific tropical mean climate (Han et al., 2021; Ren et al., 2023; Zhang et al., 2024). Figure 8 presents a summary view of the interactions of mid-Pliocene tropical and extratropical variability and mean state changes according to the PlioMIP2 MMM.

Figure 8 Summary view of interactions and mechanisms of mid-Pliocene tropical and extratropical variability in the North Pacific according to the PlioMIP2 multi-model mean. Blue boxes are ocean features, and yellow boxes are atmosphere features. Black arrows are relationships with variability of which we know the direction of influence (through lead–lag), black connector lines are relationships with the mean state (through ensemble correlations) of which we do not know the direction, and ★ indicates the dominant driver out of the two main drivers. The grey box indicates features primarily assessed in the present study. PDO results are obtained from Katya Canal-Solis and Julia Tindall through personal communications (2024) (PC).

In PlioMIP2, ENSO variability is suppressed (Oldeman et al., 2021), which is caused by a series of off-equatorial processes triggered by a northward displacement of the Pacific ITCZ (Pontes et al., 2022). In the tropical Indo-Pacific mean climate, the rising branch of the Pacific Walker circulation (PWC) is shifted westwards, both in the annual mean (Han et al., 2021) and in boreal winter (Zhang et al., 2024). Furthermore, the Indian Ocean, the Maritime Continent, and the WEP warm, but the Indian Ocean warms more than the WEP (Ren et al., 2023). We hypothesize that the westward shift of the PWC, the Indo-Pacific warming asymmetry, and the ENSO suppression are related, and we find that there is indeed a significant ensemble correlation between the change in these three variables (Fig. S8). Explaining the direction and causality of these links is out of the scope of this work.

In the present-day extratropics, ENSO forces AL variability and PDO variability, AL variability is also forced by internal stochastic variability which is dominant, and PDO variability is also forced by the AL variability which is dominant over ENSO forcing (e.g. Trenberth and Hurrell, 1994; Di Lorenzo et al., 2013; Newman et al., 2016). In PlioMIP2, the amplitude of PDO variability is similar to the pre-industrial, and ENSO leads PDO variability with a similar correlation (Fig. S9). In this study, we find that the amplitude of AL variability is also similar and that the leading regression with ENSO is largely unchanged. Furthermore, the internal stochastic atmospheric variability forcing AL variability is similar. So, despite ENSO being suppressed, AL variability is similar because the dominant driver (internal variability) is similar as in the pre-industrial. Likewise, PDO variability is similar despite ENSO being suppressed, because the dominant driver (AL variability) is similar to that in the pre-industrial. We also confirm that, in PlioMIP2, AL variability leads the PDO, with a higher correlation than with ENSO (Fig. S1).

The benefit of the view presented in Fig. 8 is that we are able to combine previously published results on changes to (extra-)tropical variability and mean climate in PlioMIP2 to the results of this study. It provides a view that shows the PlioMIP2 MMM changes tendency, with which a majority of the models (generally) agree. A weakness, however, is that information on spread in the ensemble is not included. For example, while a clear majority of the PlioMIP2 ensemble shows a westward shift in the PWC (11 out of 13; Han et al., 2021) and reduced ENSO variability (15 out of 17; Oldeman et al., 2021), the change in AL variability is similar according to the MMM, but most models show a substantial change which is either positive or negative. In this specific example, we have described these opposing responses as a “tug of war” between different responses to different boundary conditions (Fig. 7). In conclusion, the MMM view can provide useful information, but it can conceal a range of changes in opposing directions that can relate to the amplitude of a more obvious change in one direction.

5 Conclusions

The mid-Pliocene is the most recent geological period with similar atmospheric CO₂ concentrations to the present day and similar surface temperatures to those projected at the end of this century for a moderate warming scenario. According to the mid-Pliocene modelling ensemble PlioMIP2, ENSO variability was suppressed (Oldeman et al., 2021; Pontes et al., 2022). In this study, we investigate how variability in the North Pacific atmosphere responds to a suppressed ENSO in the warm mid-Pliocene climate according to the PlioMIP2 ensemble.

We find limited changes to the variability in the Aleutian Low (AL) in the MMM compared to the pre-industrial. Models with a similar ENSO variability show an increase in AL variability in the mid-Pliocene, while models with a suppressed ENSO variability also show suppressed AL variability. The majority of the ensemble shows relatively little change in ENSO–AL teleconnection strength. We separate the AL variability into a part that regresses with ENSO and into a residual, which mainly represents internal variability. We find, agreeing with observations, that most of the AL variance is explained by the residual and that a smaller part is explained by the variability that regresses with ENSO (in both climates). The change in ENSO-related AL variability is strongly related to the change in ENSO itself, while the change in the residual AL variability is not related to ENSO change. A brief investigation does not reveal one change that is able to explain changes in the residual AL variability for the whole ensemble. We find that the specific change in ENSO and AL variability depends on the relative model sensitivity either to elevated CO₂ or to the other mid-Pliocene boundary conditions, which include closed Arctic gateways and reduced ice sheet extent. Specifically, models that are relatively sensitive to elevated CO₂ generally show ENSO variability that is similar or increased compared to the pre-industrial and AL variability that is increased, while models that are relatively sensitive to the mid-Pliocene boundary conditions generally show reduced ENSO and similar or reduced AL variability.

We present a summary perspective of tropical–North Pacific variability in PlioMIP2, combining our results with the literature. Changes in the tropical Pacific mean climate result in a suppression of the mid-Pliocene ENSO variability, but since the dominant internal variability in the North Pacific extratropics is largely unchanged, both AL variability in the atmosphere and PDO variability in the ocean are similar to the pre-industrial and are not necessarily suppressed. Evaluating the effect of past climate boundary conditions, including changes to ice sheets and Arctic gateways, on the internal variability of the extratropical atmosphere is a topic for future research.

While the mid-Pliocene is not a perfect analogue for near-future climate, investigating atmosphere–ocean interactions in the mid-Pliocene provides a useful view of the functioning of the Earth system under different forcings. Our results show that teleconnections of a suppressed ENSO in a warmer past climate are quite robust. Furthermore, we are able to explain how suppressed ENSO variability does not have to lead to a suppression of its connected modes of variability in the North Pacific. Lastly, we demonstrate that, in addition to equilibrium climate sensitivity, we need Earth system sensitivity in order to explain the spread in simulated climate variability responses in the mid-Pliocene.