Following Fedorov et al. ( 2015 ), we reconstruct SST anomalies from key regions (eastern Pacific, western Pacific, tropics, and midlatitudes) as well as zonal and meridional gradients. We investigate the meridional SST gradient (tropics—midlatitudes) because one mechanism for the permanent El Niño involves strong subtropical warming, which would feed warm waters into the tropical Pacific thermocline and reduce the zonal gradient (Burls & Fedorov, 2014 ; Fedorov et al., 2015 ). We focus on the mid‐Piacenzian warm period (MPWP, 3.264–3.025 Ma; Dowsett et al., 2010 ), because the Pliocene Model Intercomparison Project 1 (PlioMIP1; Haywood et al., 2013 ) targeted this specific time. We use model output from PlioMIP1 and a new simulation conducted with the Community Earth System Model (CESM1‐CAM5) to compare with proxy‐derived estimates. We aim to assess (1) whether there was a permanent El Niño, and (2) more generally, whether simulations forced with standard MPWP boundary conditions can reproduce proxy SST anomalies.

Pliocene SST anomalies (mid‐Piacenzian–preindustrial, calculated from) at the core sites used in this study. Sites are labeled with their Deep Sea Drilling Project, Ocean Drilling Project, or Integrated Ocean Drilling Project number. Layered circles show the 2.5% (outside ring), 50% (median, middle ring), and 95% confidence levels (inner circle) for each site. Shaded bands and darker boxes mark the tropical and midlatitude regions, and western and eastern Pacific regions, used for model box averages, respectively. SST = sea surface temperature.

Here, we recalculate Pliocene patterns of warmth and SST gradients using a single proxy and a probabilistic approach. Given controversies associated with the Mg/Ca and TEX 86 paleothermometers (O'Brien et al., 2014 ; Ravelo et al., 2014 ) we base our reconstruction solely on alkenone ( ) paleothermometric data, an established proxy for SSTs. Eighty‐three percent of the SST records from the Pliocene have measurements, providing global coverage (Figure 1 ). The main limitation of is that it has an upper bound of ∼30 °C, challenging inference in warm waters. However, a new Bayesian calibration, BAYSPLINE, improves prediction of SSTs near the saturation point (Tierney & Tingley, 2018 ) and allows us to generate probabilistic representations of Pliocene SSTs that fully incorporate uncertainties in the proxy.

As an equilibrated climate state under CO 2 levels (300–450 ppm; Foster et al., 2017 ) similar to the current value (410 ppm), the Pliocene Epoch (5.3–2.6 Ma) offers useful insights into the climate dynamics of the near future. One intriguing feature of the Pliocene is evidence of a “permanent El Niño” in the tropical Pacific (Fedorov et al., 2006 ; Wara et al., 2005 ), that is, a severe reduction to near elimination of the average east‐west temperature gradient along the equator. However, the magnitude of this reduction is uncertain, with some studies suggesting only a moderate change (O'Brien et al., 2014 ; Zhang et al., 2014 ), while others confirm the original estimate (Brierley et al., 2015 ; Ravelo et al., 2014 ). Resolving this debate is important, because CO 2 forcing alone cannot explain the Pliocene permanent El Niño. Changes not simulated by current state‐of‐the‐art models, such as enhanced tropical cyclone activity (Fedorov et al., 2010 ), and extreme changes in cloud albedo (Burls & Fedorov, 2014 ) are needed to produce a strong reduction in the Pacific sea surface temperature (SST) gradient. If such shifts took place in the Pliocene, they could become active in the warmer future, with implications for climate prediction.

To create a field reconstruction of Pliocene SST anomalies, we follow the reduced space methodology of Gill et al. ( 2016 ) with modifications to propagate proxy uncertainties and accommodate Bayesian inference. The reduced space method uses the leading EOFs of historical SST variability—including the global warming signature and patterns associated with El Niño (Figure S1)—in order to reconstruct the Pliocene Pacific SST field. We posit that this is a reasonable assumption, as Pliocene continental configuration is similar to present day, and the leading EOFs in the CESM1 preindustrial and Pliocene SST simulations are similar (Figure S2). Briefly, the method proceeds as follows: Pacific SST anomalies from the ERSSTv5 product (Huang et al., 2017 ) are subject to singular value decomposition, alongside a limited field that consists only of the locations where there are Pliocene proxy data. The leading principal components (PCs) of the limited field are used to predict the PCs of the full field using Bayesian linear regression. The Pliocene proxy SST anomalies are then converted to PC space, and the regressions are used to predict full field PCs. The full field for the Pliocene is computed from eigenvalue expansion. Mathematical details are provided in Text S3.

We analyzed simulated SSTs from nine models participating in PlioMIP1 Experiment 2, run under 405 ppm CO 2 (Haywood et al., 2011 ). The other boundary conditions for the experiment are based on the PRISM3D reconstruction (Dowsett et al., 2010 ) and described in Haywood et al. ( 2010 , 2011 ); see also Text S2. We also use results from a new Pliocene simulation run with CESM1 (Hurrell et al., 2013 ), under the same boundary conditions as PlioMIP1 Experiment 2 (see Text S2), and the extended Representative Concentration Pathway (RCP) 2.6 simulation conducted with CESM1 as part of the CMIP5 effort (Meinshausen et al., 2011 ; Taylor et al., 2012 ). RCP2.6 is an aggressive mitigation scenario, under which CO 2 concentrations peak at ∼450 ppm by 2050 and then decline to 360 ppm by 2300. These CO 2 levels are similar to Pliocene values (300–450 ppm; Foster et al., 2017 ). For all experiments, we analyze the last 50 years of the simulation. To compare model output with the regional alkenone reconstructions, we take box averages across the regions of interest, including the western Pacific warm pool (10°S–10°N, 130–170°E), the tropics (15°S–15°N), the eastern equatorial Pacific cold tongue (2°S–2°N, 140–90°W) and the midlatitudes (50–35°S, 35–50°N) (Figure 1 ). Median anomalies and 2 σ standard errors (based on bootstrap estimates of annual data) are used to compare the model output with the proxy data. Box averages are preferred over using output from the individual site locations because they are less sensitive to model‐specific biases in the Pacific cold tongue, warm pool, and midlatitude fronts; however, our results are not sensitive to this choice (see Text S2).

Reconstruction of SST anomalies and gradients for the past 4 Ma fromdata and comparison to model simulations. Each panel shows anomalies relative to preindustrial. A Gaussian smoothing of 400 ka was applied to the reconstructions; note that smoothed data do not equal 0 at 0 Ma because they average glacial‐interglacial cycles. Shaded areas represent the 95% confidence interval. Gray bar encloses the mid‐Piacenzian warm period. Blue dots show PlioMIP1 model simulations; orange and red dots show CESM Pliocene and RCP2.6 simulations, respectively. Model estimates are offset from one another for clarity and are plotted with 2error bars that account for interannual variability. Gray dots show previous estimates from Wara et al. () and Fedorov et al. (). The “Permanent El Niño” threshold of −3 °C in anomaly space is equivalent to an absolute west‐east gradient of 1.2 °C (see main text). SST = sea surface temperature; CESM = Community Earth System Model; RCP = Representative Concentration Pathway.

Our SST anomaly reconstructions are based on published alkenone data from 28 sites (Figure 1 and supporting information Table S1). We use these data to reconstruct regional SST anomalies and the Pacific SST field (see section 2.3 ). The regions targeted include the midlatitudes, following the definition in Fedorov et al. ( 2015 ; 30–50°N or °S, excludes subtropical upwelling zones, includes Site 982 at 57°N), the eastern equatorial Pacific, the western Pacific warm pool, and the tropics (see Table S1 for a list of sites assigned to each region). Our analysis is limited to the last 4 Ma because prior to this time, only one site is available in the western Pacific (ODP 806) and this biases the regional average. We calibrated data to SSTs using BAYSPLINE (Tierney & Tingley, 2018 ) with a prior standard deviation of 5 °C, a loose prior that only minimally impacts the posterior. For each site, anomalies were calculated relative to a preindustrial estimate of SST, which, wherever possible, was based on a coretop or latest Holocene value (see Text S1). In order to equally weight the contribution of each site, data were interpolated to a common time step of 5 ka between 0 and 4 Ma prior to taking a regional average. To emphasize long‐term (supra‐orbital) variability, a 400‐ka (full‐width, 6 σ ) Gaussian smoothing is applied to these estimates in Figure 2 . To assess whether the proxy data indicate a permanent El Niño, we formally define this as a zonal gradient of less than 1.2 °C, the average gradient during the strongest historically observed El Niño events (1877, 1982, 1997, and 2015). This limit is calculated from the ERSSTv5 product (Huang et al., 2017 ), using the same regional averages as the model simulations (see section 2.2 ).

The reduced space reconstruction shows enhanced warming in the eastern equatorial Pacific cold tongue and the Peru and California margin upwelling zones and a horseshoe pattern of warming that extends up the Alaskan margin and southwestward across the northern subtropical Pacific (Figure 3 a). Compared to the PRISM3 reconstruction, our SST anomalies are 1–2 °C warmer in the tropics and 1–3 °C cooler in the subtropics (Figure S3). For model comparison, we focus on CESM1, which simulates a zonal gradient reduction within range of the proxy estimates (Figure 2 c). The same spatial features—enhanced warming in the northeast Pacific and eastern equatorial region—are clearly seen in the simulations of both the MPWP and RCP 2.6 (Figures 3 b and 3 c). The magnitudes of the simulated anomalies are smaller than the median field, but well within uncertainty of the reconstruction (Figure S4a). The model simulations show different responses in the northwest Pacific and southern Ocean (Figure 3 ); however, our reconstruction has limited skill (and data) in these regions (Figure S4b). Across the tropics and the subtropics (35°S to 35°N), where most proxy data lie and the reconstruction has the strongest skill, the correlation between the reconstructed anomalies and the CESM Pliocene anomalies is 0.71.

Our reconstructed regional SST anomalies indicate that the western Pacific warm pool was 1.3 °C (95% confidence interval (CI) = 0.6 to 2.0 °C, Figure 2 a) warmer than preindustrial during the MPWP, similar to the tropics as a whole (1.4 °C, 95% CI = 1.0 to 2.0 °C, Figure 2 b). The magnitude of warming in the eastern Pacific and midlatitudes is larger, such that there is a moderate but significant decrease in both the tropical Pacific zonal gradient (−0.9 °C, 95% CI = −1.9 to −0.1 °C) and the meridional gradient (−1.2 °C, 95% CI = −1.7 to −0.6 °C; Figures 2 c and 2 d). Our estimate for the change in the zonal gradient is well above the threshold for a permanent El Niño (Figure 2 c). The proxy‐inferred anomalies compare favorably with simulations of MPWP climate from PlioMIP and CESM1 (Figure 2 ). The majority of the models predict warm pool, tropical, and zonal gradient anomalies within the 95% CI of the proxy estimates (8/10). Half of the models produce a meridional gradient reduction within the 95% range (5/10). CESM predictions for years 2250–2300 under RCP 2.6 show a warming in the western Pacific and tropics and a decrease in the zonal gradient that are similar to the Pliocene data and simulations (1.1 °C, 1.3 °C, and −0.5 °C, respectively; Figure 2 ).

4 Discussion

Overall, our new reconstructions of Pliocene SSTs show good agreement with model simulations during the MPWP. Most models can reproduce the magnitude of warm pool and tropical warmth and moderate zonal gradient reduction (Figure 2). Some models do underestimate the decrease in the zonal or meridional gradient, but the mismatch is not as large as previously thought (Figures 2c and 2d). Notably, the RCP2.6 simulation with CESM1 results in Pliocene‐like warming and zonal gradient reduction (Figure 2). Although Pliocene changes in orography, ice, and vegetation contribute to a rise in global mean air temperature (Lunt et al., 2012), the impact of these effects on SSTs is smaller, especially in the tropics and subtropics. Our results with CESM1 suggest that for the regions of interest here, Pliocene CO 2 forcing is dominant.

Our results support the inferences based on the TEX 86 paleothermometer that the Pliocene Pacific zonal gradient was moderately lower than preindustrial (∼−1 °C; O'Brien et al., 2014; Zhang et al., 2014) and disagree with studies that have argued for a −3 to −4 °C change (Fedorov et al., 2015; Wara et al., 2005) consistent with a permanent El Niño (Figure 2c). The difference between our result and the latter work can be attributed to several factors. The first is our use of BAYSPLINE, which ameliorates the tendency of previous calibrations (e.g., Müller et al., 1998) to underestimate tropical SSTs (Tierney & Tingley, 2018), produces warmer estimates for tropical sites during the Pliocene (Figure S5), and therefore yields smaller zonal and meridional gradient anomalies. We acknowledge that, at Sites 1143 and 806, is close to, and occasionally reaches, its limit (a value of 1). In the BAYSPLINE calibration, the limit does not correspond to a single value of SST but rather is probabilistically estimated to be 30.4 °C ± 5.5 °C (2σ). However, comparison between ‐inferred SSTs and two other paleothermometers (TEX 86 and Mg/Ca with a seawater adjustment) at these sites shows good agreement in the MPWP (Figure S5). This confirms that our choice to focus on a single proxy does not skew our results.

In the western Pacific warm pool, previous work relied on analysis of faunal assemblages (Dowsett et al., 2009, 2011) or the Mg/Ca paleothermometer (Fedorov et al., 2015; Wara et al., 2005) and found either no evidence for warmer temperatures than preindustrial or a cooling relative to preindustrial (e.g., Figures 2a and S2). However, faunal assemblage methods cannot by design predict temperatures above their calibration range (Dowsett, 2007), and Mg/Ca inferences did not account for the increase in Mg/Ca of seawater (Mg/Ca sw ) over the last several million years (Evans et al., 2016). Both factors cause systematic underestimation of western Pacific warming during the Pliocene and a greater inferred reduction in the zonal gradient. Since Mg/Ca sensitivity to temperature is nonlinear, the change in Mg/Ca sw has a bigger impact on the warmer western Pacific sites than the eastern Pacific sites. This results in artificially small gradient calculations, even if only Mg/Ca data are used. The sensitivity of Mg/Ca to other environmental factors, including salinity, pH, and dissolution (Gray et al., 2018; Regenberg et al., 2014), may further confound SST inference (O'Brien et al., 2014).

Another factor contributing to the more moderate results in our study is that we consider anomalies relative to preindustrial conditions for the MPWP specifically so that we can make a proportionate comparison to the model simulations. Previous work estimated gradient changes using a running mean that smoothed over Late Pleistocene glacial cycles (Brierley et al., 2015; Fedorov et al., 2015). This includes the impact of glaciation and therefore yields larger anomalies. Previous work also compared proxy data spanning all of the Pliocene, 5.3–2.6 Ma, to the model simulations, even though CO 2 levels may have changed over this period. This said, if we extend our proxy average to 4–2.6 Ma, we obtain similar results: the zonal gradient changes by −1.1 °C (95% CI = −2.2 to −0.28 °C) and the meridional gradient by −1.3 °C (95% CI = −1.8 to −0.60 °C).

A final difference is our choice to include all available sites in the western and eastern Pacific regional averages. In our view, this is the most methodologically defensible approach, as it ensures that random noise associated with a single proxy site, and/or genuine but localized variations in SST, have less impact on the averages. In particular, there has been debate over whether Site 1143, located in the South China Sea (Figure 1), is representative of the western Pacific (Brierley et al., 2015; O'Brien et al., 2014; Zhang et al., 2014). Based on analysis of the ERSSTv5 product, Site 1143 has a mean annual preindustrial SST of 28.1 °C, and SST variability is strongly correlated with the open‐water warm pool (10°S–10°N, 130–170°E; r=0.66, p≪0.0001), and Site 806 (r=0.60, p≪0.0001). This site also sits underneath the ascending branch of the Walker cell. For these reasons, we retain Site 1143 in our regional average; however, our results are not sensitive to this choice (Figure S6).

Previous work suggested that the reduction in the Pacific zonal gradient during the Pliocene is closely tied to changes in the meridional gradient; that is, entrainment of warm subtropical waters into the Pacific thermocline drives a relaxation of the zonal gradient (Burls & Fedorov, 2014; Fedorov et al., 2015). However, we do not find a significant relationship between the meridional and zonal SST gradients across models (Figure 4a). Likewise, while the proxy data indicate that both the zonal and meridional gradients decrease back in time, the changes in the zonal gradient are larger (by about a factor of 2; Figures 2c and 2d). In the model simulations, the change in the sea level pressure gradient across the Pacific—a proxy for the Walker Circulation—is a stronger predictor of zonal gradient change than the meridional gradient (Figure 4b). This underscores the importance of tropical rather than extratropical ocean‐atmosphere dynamics in driving a reduced zonal gradient during the Pliocene.

Figure 4 Open in figure viewer PowerPoint The Pacific zonal SST gradient change versus the meridional SST gradient change (a) and sea level pressure gradient (a proxy for Walker circulation) change across the equatorial Pacific (5°S to 5°N, 80–140°W—5°S to 5°N, 155–180°E; b) across model simulations. Error bars represent the 2σ bootstrap uncertainties. Pink box in panel (a) encloses the proxy‐inferred range of values for the mid‐Piacenzian warm period. The MIROC model is not plotted in (b) because sea level pressure data were not available. Spearman ρ correlations are shown in the lower right corner; the value in panel (a) is not significant (p=0.13). SST = sea surface temperature.

There is a notable diversity in the zonal gradient response across the models (Figure 4b), with two models predicting an increase during the Pliocene. This could reflect cold tongue bias, which impacts both the direction and magnitude of zonal gradient change in response to elevated CO 2 (Li et al., 2016). Models with an exaggerated cold tongue tend to simulate smaller reductions in the gradient or even a strengthened gradient, because they underestimate precipitation in the western Pacific and therefore overestimate warming in that area (Li et al., 2016). This may explain the behavior of the COSMOS model, which has a cold tongue that extends to Indonesia. It predicts overly strong warming in the western Pacific with a high standard deviation (2.7 ± 0.27 °C; Figure 2a) and a strengthened gradient. However, across the remaining PlioMIP models there is no relationship between western Pacific warming and the zonal gradient change (ρ=−0.08). There is also no significant relationship between bias in cold tongue temperature (assessed relative to ERSSTv5, 1854–1900), nor bias in the longitudinal extent of the cold tongue (assessed as the location of the 28 °C isotherm at the equator), and the simulated zonal gradient change, respectively (ρ=0.01, ρ=−0.21, p>0.5 in both cases). This suggests to us that cold tongue bias is not responsible for the model diversity.

Rather, to a first order the magnitude of zonal gradient change in each model is directly related to its Walker circulation change (Figure 4b). We suspect that the diversity of Walker response across models reflects different parameterizations that impact convection, cloud feedbacks, evaporative damping, and ocean‐atmosphere coupling in the tropics. Comparison between CCSM4 and CESM1 is illustrative in this respect, as these models differ only in the generation of atmospheric model used. CESM1, which incorporates updated cloud parameterizations and microphysics (Gettelman, 2015), simulates a much stronger weakening of the Walker circulation (Figure S7).