Articles | Volume 18, issue 4
Research article
06 Apr 2022
Research article |  | 06 Apr 2022

Warm mid-Pliocene conditions without high climate sensitivity: the CCSM4-Utrecht (CESM 1.0.5) contribution to the PlioMIP2

Michiel L. J. Baatsen, Anna S. von der Heydt, Michael A. Kliphuis, Arthur M. Oldeman, and Julia E. Weiffenbach

We present the Utrecht contribution to the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2), using the Community Earth System Model version 1.0.5 (CCSM4-Utr). Using a standard pre-industrial configuration and the enhanced PlioMIP2 set of boundary conditions, we perform a set of simulations at various levels of atmospheric pCO2 (280, 400, and 560 ppm). This allows us to make an assessment of the mid-Pliocene reference (Eoi400) climate versus available proxy records and a pre-industrial control (E280), as well as determine the sensitivity to different external forcing mechanisms.

We find that our simulated Pliocene climate is considerably warmer than the pre-industrial reference, even under the same levels of atmospheric pCO2. Compared to the E280 case, our simulated Eoi400 climate is on average almost 5 C warmer at the surface. Our Eoi400 case is among the warmest within the PlioMIP2 ensemble and only comparable to the results of models with a much higher climate sensitivity (i.e. CESM2, EC-Earth3.3, and HadGEM3). This is accompanied by a considerable polar amplification factor, increased globally averaged precipitation, and greatly reduced sea ice cover with respect to the pre-industrial reference. In addition to radiative feedbacks (mainly surface albedo, CO2, and water vapour), a major contribution to the enhanced Pliocene warmth in these simulations is the warm model initialisation followed by a long spin-up, as opposed to starting from pre-industrial or present-day conditions. Added warmth in the deep ocean is partly the result of using an altered vertical mixing parameterisation in the Pliocene simulations, but this has a negligible effect at the surface. We find a stronger and deeper Atlantic meridional overturning circulation (AMOC) in the Eoi400 case, but the associated meridional heat transport is mostly unaffected. In addition to the mean state, we find significant shifts in the behaviour of the dominant modes of variability at annual to decadal timescales. The Eoi400 El Niño–Southern Oscillation (ENSO) amplitude is greatly reduced (−68 %) versus the E280 one, while the AMOC becomes more variable. There is also a strong coupling between AMOC strength and North Atlantic sea surface temperature (SST) variability in the Eoi400, while North Pacific SST anomalies seem to have a reduced global influence with respect to the E280 through the weakened ENSO.

1 Introduction

The Pliocene Model Intercomparison Project Phase 2 (PlioMIP2; Haywood et al.2016) aims to simulate the mid-Piacenzian Warm Period (mPWP; 3.264 to 3.025 Ma; Haywood et al.2013a) of the mid-Pliocene using a suite of global numerical climate models. During this interval, the Earth's climate saw warm and relatively stable conditions with atmospheric CO2 levels similar to those seen today. Being relatively recent in the planet's geologic history, geographical boundary conditions during the mid-Piacenzian were similar to the present. Therefore, the simulated climatic conditions for this specific time interval may serve as a suitable analogue of what to expect in the next century (Burke et al.2018).

Previous modelling efforts within PlioMIP1 resulted in globally averaged near-surface air temperatures between 1.8 and 3.6 C warmer than those seen today (Haywood et al.2013b) for the mid-Pliocene. The pattern of warming showed clear polar amplification, resulting in a reduced Equator-to-pole temperature gradient. However, this reduction was still considerably smaller than suggested by various temperature proxy records, and models mostly failed to reproduce the strong warming pattern over the North Atlantic Ocean. The simulated mid-Pliocene warming response was found to be mostly a result of increased atmospheric CO2, with surface albedo feedbacks the primary driver behind polar amplification.

In addition to an extended model ensemble, PlioMIP2 offers the use of model boundary conditions from the Pliocene Research, Interpretation and Synoptic Mapping version 4 (PRISM4; Dowsett et al.2016). This includes updated estimates of vegetation cover, coastlines, and topography, as well as reduced ice sheet cover, compared to PRISM3 (Dowsett et al.2010). As a result, models now show an average warming response of 1.7 to 5.2 C to these mid-Pliocene boundary conditions. In addition, most models have a mid-Pliocene temperature pattern that better resembles that of proxy reconstructions, particularly in warmer high-latitude and North Atlantic regions, compared to PlioMIP1 (Haywood et al.2020).

Here, we present the results of a set of simulations using the Community Earth System Model (CESM) version 1.0.5 within the PlioMIP2 effort. A detailed description of our model configuration, as used in these experiments, is provided in Sect. 2. This is followed by a general discussion of the PlioMIP2 experimental design and naming conventions in Sect. 3, together with an overview of our specific set of model experiments. The results of these are presented in Sect. refs:res, in which we assess climate sensitivity (Sect. 4.1), the general mPWP conditions in our simulations (Sect. 4.2), meridional overturning and heat transports (Sect. 4.3), a comparison to proxy records (Sect. 4.4), sensitivity to the applied PlioMIP2 model boundary conditions (Sect. 4.5), and modes of internal variability in the simulated climate (Sect. 4.6).

2 Model description

2.1 The CESM 1.0.5

The Community Earth System Model (Hurrell et al.2013) is a fully coupled atmosphere–land–ice–ocean general circulation model (GCM) that was developed at the National Center for Atmospheric Research (NCAR) in Boulder, Colorado. For use in palaeoclimate modelling, version 1.0.5 of the CESM is a suitable choice motivated by a trade-off between increasing model complexity and computational cost. This version of the CESM is equivalent to the last version of its predecessor, the Community Climate System Model (CCSM4; Blackmon et al.2001; Gent et al.2011). In other PlioMIP2 studies, our model simulations are therefore referred to as CCSM4-Utr.

2.2 Atmosphere

The atmospheric component of the CESM is the Community Atmosphere Model (CAM4; Neale et al.2013), which uses a finite-volume dynamical core. The model grid has a nominal horizontal resolution of 2 (2.5× 1.9; 144 × 96 grid cells). A total of 26 vertical levels extend upward to 2 hPa using a hybrid sigma vertical coordinate. In this configuration, the model has a warming response of 3.17 C per doubling of CO2 starting from pre-industrial conditions (Baatsen et al.2020), which is very similar to the reported value of 3.14 C by Bitz et al. (2012), and higher than the  2.5 C in CCSM3 (Kiehl et al.2006).

In accordance to the PlioMIP2 protocol (Haywood et al.2016), atmospheric concentrations other than that of CO2 are kept at their pre-industrial levels: 671 ppb CH4, 270 ppb N2O, and no CFCs. The solar constant is kept at 1361.27 W m−2 in all of our model experiments. Astronomical orbital parameters are set to their present-day configurations for both pre-industrial and Pliocene simulations, again as suggested by the PlioMIP2 protocol. Atmospheric aerosols are fixed using a pre-industrial climatology for the pre-industrial cases. An adjusted climatology is used for the Pliocene cases, resulting from a bulk aerosol model simulation to incorporate the effect of the PRISM4 boundary conditions.

2.3 Land

The physical, chemical, and biological processes taking place on land are represented in the Community Land Model (CLM4; Oleson et al.2010; Lawrence et al.2011). A static rather than dynamic vegetation model is used here to avoid runaway feedback effects, which can become an issue, especially in warm greenhouse climates (e.g. dieback of vegetation at low latitudes; Loptson et al.2014; Herold et al.2014). Either the pre-industrial biomes or the PRISM4 megabiomes (see also Supplement Table S1) are translated into fractions of the corresponding CLM4 plant functional types (PFTs), from which a set of monthly forcing files is ultimately used in the model. Fresh water runoff is treated by a simple river-routing scheme, in which all runoff is transported to one of the surrounding eight model grid cells until the ocean is reached. The direction is determined by the local topography gradient and manually adjusted where the runoff scheme would otherwise form closed loops. Within the CLM4, land ice is implemented as bare soil with a given surface elevation and its specific surface properties (e.g. albedo, evaporation and run-off). As suggested by Haywood et al. (2016), the static land ice configuration used here is based on the results of previous modelling efforts in PlioMIP1. This includes a greatly reduced Greenland Ice Sheet compared to the present day, as well as an absent West Antarctic Ice Sheet.

2.4 Ocean

The CESM1 uses the Los Alamos National Laboratory (LANL) Parallel Ocean Program version 2 (POP2; Smith et al.2010) for the ocean model component. The standard configuration is applied here, with a nominal 1 (1.25× 0.9) horizontal resolution on a curvilinear grid placing the northern pole over Greenland. The POP2 model is set up with 60 vertical layers of varying thickness, i.e. between 10 m near the surface and 250 m at greater depth. Horizontal viscosity is considered anisotropic (Smith and McWilliams2003), and horizontal tracer diffusion follows the parameterisation of Gent and McWilliams (1990). The model further uses the KPP scheme to determine vertical mixing coefficients (Large et al.1994). More information and discussion on the ocean model physics and parameterisations can be found in Danabasoglu et al. (2008, 2012).

The sea ice component consists of the LANL Community Ice Code version 4 (CICE4; Hunke and Lipscomb2008). For simplicity, sea ice only forms when the sea surface cools down to −1.8C, after which its dynamical behaviour (e.g. melt and advection) is treated by the model specifically.

In all of our model experiments considered here, tidal mixing and overflow parameterisations are switched off. The pre-industrial reference uses a uniform background vertical diffusivity, which is set at 0.16 cm2 s−1. In contrast, the Pliocene simulations have a vertically varying background vertical diffusivity determined by


where vdc1=0.524 cm2 s−1, vdc2=0.313 cm2 s−1, dpth=1000 m (the reference depth), and linv=4.5×10-3 m−1 (the inverse scaling length). κw thus varies between 0.1 cm2 s−1 at the surface and 1 cm2 s−1 at the bottom of the ocean. The latter is done to be consistent with other palaeo-climate simulations using the same model set-up, e.g. Baatsen et al. (2020). To study the effect of the different mixing parameter settings, we carry out an additional pre-industrial simulation using the same oceanic configuration as for the Pliocene cases (see Sect. 3, Figs. S4, S7, and S9 in the Supplement). The enhanced vertical mixing in the deep ocean results in an overall warming of the deep ocean at the expense of a slight cooling at upper levels. Near the surface, globally averaged temperatures are near identical, but some regional changes are seen as a result of the different vertical mixing schemes. For completeness, another pre-industrial simulation is carried out using the standard configuration in which the tidal mixing and overflow parameterisations are switched on. The results show no significant temperature differences seen at any depth level and thus suggests that the combined effect of these parameterisations is negligible within the scope of this study and the specific model configuration used here.

3 PlioMIP2 experiments

3.1 PRISM4 boundary conditions

Following the methodology outlined by Haywood et al. (2016), we perform a set of pre-industrial and Pliocene simulations using the CESM1.0.5. As we carry out a number of different model experiments, we will refer to these as either pre-industrial or Pliocene (rather than mid-Piacenzian) simulations regardless of atmospheric pCO2. For Pliocene cases we incorporate the set of enhanced boundary conditions based on the PRISM4 (Dowsett et al.2016). These include altered topography and bathymetry, coastlines, land surface properties (i.e. vegetation, soil type, and ice sheet coverage) and atmospheric composition with respect to pre-industrial conditions. Our Pliocene model geography is implemented by applying the PRISM4 topography and bathymetry anomalies to the pre-industrial reference after re-gridding onto the specific model resolution used here. The PRISM4 land–sea mask is then also interpolated onto the model grid and applied to the Pliocene model geography. A manual check of the entire model grid is done, altering land and ocean cells where needed and making sure that any marine passages are either at least two cells wide when open or closed if needed. The resulting model geography for the Pliocene simulations is shown in Fig. 1 and compared to that of the pre-industrial reference. Corresponding vegetation cover and aerosol optical depth can be found in Fig. S1 of the Supplement.

Figure 1Model geography as applied in our different (a) Pliocene and (b) pre-industrial cases. Both the bathymetry and topography are shown at the native model grid, with the CAM4 rectangular grid superimposed on the finer POP2 curvilinear grid wherever the land fraction is at least 0.5 (also indicated by the thick grey contour line). The light blue contour line shows the extent of land ice, as imposed by the land use distribution. (c, e) Northern and (d, f) southern polar stereographic views of Pliocene and pre-industrial model geography.

Some of the main aspects of the Pliocene model geography with respect to the pre-industrial reference include a closure of both the Bering Strait and Northwest Passage, making the Arctic Ocean much more isolated from the other ocean basins. Prior to large-scale glaciation on the Northern Hemisphere during the Quaternary period, a considerable part of the northern Eurasian continental shelf and the Hudson Bay area were exposed. Furthermore, large parts of the Maritime Continent are also considered to have been above sea level during the mid-Pliocene. Finally, there was significantly less land ice cover compared to the present, covering only part of southeastern Greenland and East Antarctica. Note that the Strait of Gibraltar is open while the Panama Gateway is closed in all of our simulations (this is obscured in the figures by the superposition of the atmospheric model grid onto the bathymetry).

3.2 Experimental design and model spin-up

In accordance with the PlioMIP2 guidelines we performed both a pre-industrial reference and Pliocene control core simulation, which are referred to as E280 and Eoi400, respectively (following the naming conventions of Haywood et al.2016). In addition to the pre-industrial reference we added two climate sensitivity simulations, with a doubling (E560) and quadrupling (E1120) of atmospheric CO2. Two more sensitivity experiments were added to the Pliocene control, in which we applied pre-industrial (Eoi280) and doubled (Eoi560) CO2 levels. This enables us to make a complete assessment of the model's sensitivity to either Pliocene boundary conditions or radiative forcing for the full range of possible reference states. Besides altered model boundary conditions, the Pliocene simulations differ from the pre-industrial ones as a result of the oceanic vertical mixing parameters. Another mixing sensitivity experiment is thus added as a continuation from the E280 simulation using the same ocean mixing configuration for the ocean as in the Pliocene cases, referred to as E280,P. Results of an additional (500-year) pre-industrial simulation with tidal mixing and overflow parameterisations switched on (E280,S) are not used besides a technical robustness check and are therefore not considered. See also Table 1 for an overview of the different simulations.

Figure 2Time series of globally averaged temperatures for the entire length of our three different Pliocene simulations: Eoi400 (black), Eoi280 (blue), and Eoi560 (red). Shown are the (a) upper (dark), deep (medium), and full-depth (light) ocean temperature; (b) near-surface air temperature; and (c) globally averaged top-of-model (TOM) net radiative flux. Thick lines in (b) and (c) show the corresponding time series after applying a 100-year smoothing mask. The estimated equilibrium temperatures are indicated at the end in (a) and (b) using large markers and the same colour convention. The globally averaged mean temperature (Tg in b) and net radiative flux (R in c) over the last 100 years are added in the legends (bracketed values for the estimated equilibrium).


The pre-industrial reference simulation is initialised using present-day (2000 CE) temperature and salinity fields from the PHC2 dataset (Steele et al.2001). As these are slightly warmer than pre-industrial conditions, a long spin-up of >3000 years is carried out to equilibrate the full ocean. The remaining drift of volume-weighted average ocean temperature at the end of this simulation is brought down to 10-4 K yr−1 (Table 1). In contrast, the Pliocene control simulation is started with a highly idealised ocean temperature and salinity distribution. The initial temperature is horizontally homogeneous, while salinity is set to 35 psu globally. The vertical profile of temperature decreases downward between 15 C at the surface and 4 C in the deep ocean. We thus apply similar initial conditions to those used in the Eocene simulations of Baatsen et al. (2020) but with cooler deep-ocean temperatures. The latter is done to match the volume-integrated global ocean temperatures of Rosenbloom et al. (2013), who used the CCSM4 to model the mid-Pliocene within the PlioMIP1. Thus, we start the model with a total ocean heat content similar to previously found Pliocene conditions using a very similar model configuration. Still, we perform 2000 model years of spin-up to acquire a well-equilibrated oceanic state for our Pliocene control. An overview of the temperature evolution during this spin-up, along with the top-of-model radiative imbalance can be found in Fig. 2, which shows substantial temperature adjustments. After an initial warming phase, an abrupt change in temperature trends is seen at  1000 years into the Eoi400 simulation. This is the result of a shift in the oceanic circulation, associated with the formation of deep overturning cells in the North Atlantic and Southern Ocean. These act to mix heat from the upper ocean into the deep ocean but also increase high-latitude temperatures through enhanced meridional heat transport (see also Fig. 7). Similar temperature trends were seen by Chandan and Richard Peltier (2017), especially considering their POP1 type vertical mixing simulations, in which they implement an identical ocean set-up to our mid-Pliocene cases (and E280,P) but with a cooler initial state. Additional information on the spin-up behaviour of the set of model simulations presented here can be found in Figs. S2 and S3 of the Supplement.

All three pre-industrial sensitivity experiments are started from the equilibrated state at the end of the E280 control simulation. The E560 and E1120 cases are continued for another 1000 and 2000 model years, respectively, while the E280,P is continued for 2300 model years (see also Fig. S4 of the Supplement). The E560 and E1120 simulations still have substantial drifts in deep-ocean temperature, but this is much lower near the surface. Using the transient behaviour, however, we can estimate the actual equilibrium temperatures (see Sect. 3.3) and therefore the climate sensitivity of the model within this specific set-up. In the global average, the effect of altered ocean mixing parameters is quite clear: there is no temperature change at the surface, but the upper ocean cools slightly while the deep ocean warms. The more efficient downward mixing of heat leads to an average warming of the total water column (by ∼0.8C) in the extrapolated equilibrium. The Eoi280 and Eoi560 simulations are initialised from the Eoi400 control at year 2000 and continued for 1000 and 1500 years, respectively. As their initial radiative forcing is relatively small compared to a full doubling of CO2, it does not take as long for the temperature drifts to become small. The Eoi560 simulation is longer compared to the Eoi280 one as it was showing substantial internal variability after about 700 years owing to large (i.e. >10 Sv) fluctuations in the Atlantic meridional overturning circulation strength (Fig. S5 in the Supplement).

3.3 Data analysis and methods

Within the PlioMIP2 database, model fields from the last 100 model years of the E280, E560, Eoi280, Eoi400, and Eoi560 cases are publicly available. We use the same 100 years here for our analyses and to calculate time means in order to best match the results of other studies. Some more extensive data time series are considered as well but only when specifically noted. The overview of modelled atmospheric and oceanic conditions given in Sect. 4.2 is based on climatologies using the same 100 years at the end of the E280 and Eoi400 runs. From the E280,P simulation we use 200-year averages instead (i.e. 4400–4600) as it exhibits some long-term fluctuations (see Fig. S4 in the Supplement). In addition to the model data presented here, we make use of satellite-based passive microwave sea ice concentration data (Peng et al.2013). We consider the 1988–2000 mean sea ice concentration, excluding more recent years in which a rapid decline of Arctic sea ice was observed. We also use mid-Pliocene (more specifically the KM5c interglacial period) sea surface temperature (SST) proxies from McClymont et al. (2020). This record includes absolute SST estimates using the newly calibrated U37K and Mg/Ca proxies, as well as SST anomalies, which are calculated with respect to the 1870–1899 ERSST reanalysis (Huang et al.2017).

In addition to direct climatological or annual mean fields at the end of each simulation, we also present a number of estimated equilibrium values using the complete time series of the according temperatures and radiative fluxes. Based on the method of Gregory et al. (2004), we linearly extrapolate the trend of a specific temperature measure towards a net zero top-of-model (TOM) radiative flux. It should be noted that the original technique of Gregory et al. (2004) is designed to estimate equilibrium temperatures from relatively short (∼50 years) simulations following an initial perturbation. Using a much longer time series, we do include the effects of slower (e.g. ocean thermal adjustment; Baatsen et al.2020; Farnsworth et al.2019) and non-linear (e.g. clouds; Bloch-Johnson et al.2015; Knutti and Rugenstein2015) feedbacks and therefore can acquire a better estimate of the actual equilibrium temperature (but also require a much longer simulation time). The net TOM radiation is calculated by subtracting the (surface-weighted) globally averaged net outgoing longwave flux from the net incoming shortwave flux. A visualisation of the standard procedure for globally averaged near-surface air temperature in the Eoi280, Eoi400, and Eoi560 cases can be found in Fig. S2 of the Supplement. An extension towards upper (<1 km), deep (>2 km), and full-depth ocean temperatures from all of our simulations is also provided in Fig. S3. As shown by Gregory et al. (2004), the equilibration of the globally averaged near-surface air temperature from an initial shock in radiative forcing becomes linear as a function of the net TOM radiative flux. However, their simulations use a slab ocean rather than a full-depth ocean model component. Using our model configuration, the equilibration is slightly non-linear instead, especially at the start of the simulation (see also Baatsen et al.2020). Therefore, we exclude the first 100 and 250 years of data from this analysis for atmospheric and oceanic temperatures, respectively.

To get a quantitative assessment of the temperature contributions from different components in the radiative balance between our sensitivity simulations, we repeat the analysis of Hill et al. (2014) using an energy balance framework similar to Lunt et al. (2021); Heinemann et al. (2009). We adopt the radiative balance equation as follows:


where S is the incoming solar radiation, αp the planetary albedo, H the meridional heat transport, Fcloud the cloud radiative forcing, FCO2 the radiative forcing resulting from a doubling of atmospheric pCO2, ϵ the emissivity, σ the Stefan–Boltzmann constant, and τ the surface temperature. This equation can be rearranged to determine the surface temperature, using the annual mean model fields from each simulation. The temperature effect of a single component between two different simulations can then be estimated by the temperature difference by adjusting this component in the energy balance equation. Note that in the reference temperature both Fcloud and FCO2 are 0. In contrast to Lunt et al. (2021), we consider the full 2D fields and then take the zonal average of the resulting temperature difference. We also include the effects of surface albedo, shortwave and longwave cloud forcing, and CO2 as components in the energy balance.

In our analysis of the ocean circulation, we consider the barotropic stream function (BSF) and meridional overturning stream function (MSF). Both of these are calculated from the monthly averaged horizontal flow fields and integrated as a function of latitude and either depth (BSF) or longitude (MSF), starting from the southern pole located on Antarctica. The last 100 years of each simulation are again used for the calculation of the BSF, but the last 500 years are used for the MSF (to avoid a possible influence of centennial-scale variability). From the MSF, we also determine the strength of the Atlantic meridional overturning circulation (AMOC). This is defined as the maximum of the overturning stream function below 1000 m and north of 30 S using only the flow field in the Atlantic Ocean.

We look at the temporal variability of SSTs through El Niño–Southern Oscillation (ENSO), Atlantic/Pacific (Multi-)Decadal Variability (AMV/PMV) and the AMOC time series. ENSO is characterised by the Niño indices, taking the monthly average SST anomaly over 0–10 S, 90–80 W (Niño 1+2) and 5 N–5 S, 170–120 W (Niño 3.4). To capture ENSO variability on decadal timescales as well, we use 200 rather than 100 years of monthly SST data. We subtract the 12-monthly climatological mean and apply a linear detrending to each model grid point. The resulting Niño time series are then treated with a 5-month running mean.

The calculation of AMV/PMV indices requires the use of empirical orthogonal functions (EOFs), using 500 years of annual mean SST data as these processes occur on longer timescales (see, e.g. Deser et al.2012; Trenberth and Shea2006, and; last access: 17 August 2021). We apply a local detrending to annual SST time series at each grid point by subtracting a 200-year smoothing spline. We then determine the first three spatial EOFs of the resulting SST anomalies for the North Atlantic (10–70 N) and North Pacific (20–60 N) oceans. The AMV mode is then determined as the EOF with which the time series correlates best with the North Atlantic average SST (rather than the first EOF, as the variability may be dominated by, e.g. ENSO or PMV). The PMV is simply the dominant EOF in the North Pacific Ocean obtained from the annual SST anomalies. Patterns of SST variability corresponding to these modes can then be obtained by correlating the local SST anomalies to each of the specific mode time series.

We perform a spectral analysis of the different modes of variability using a multi-tapered method (MTM) as it is more suitable for climatic time series of limited length (Ghil et al.2002). This method is based on a standard Fourier analysis but uses tapers to limit cut-off effects at the edges of the respective time series. We apply this technique using three tapers and a bandwidth parameter of 2. All of the resulting power spectra are tested against a red-noise null hypothesis, using 10 000 random AR1 surrogate time series. The median of these red-noise power spectra is used to normalise those of the corresponding modes of variability, making the latter easier to interpret as well as consider their statistical significance using 90%, 95%, and 99% confidence levels. Finally, we also correlate the time series of different modes: AMV, PMV, AMOC, and ENSO (now using 500 years of annual SST anomalies for the latter as well, rather than 200 years of monthly SST anomalies). These correlations are considered significant when their corresponding p value is less than 0.05 (i.e. white-noise null hypothesis).

4 Results

4.1 Global averages and climate sensitivity

An overview of globally averaged temperatures, drifts, and radiative balance at the end of each of our model simulations is provided in Table 1. These values, including some additional metrics of the vertical thermal distribution of the ocean, are visualised in Fig. 3. Throughout the different cases, we see a consistent globally averaged equilibrium warming of ∼3C per CO2 doubling in the near-surface atmosphere. This warming response reduces slightly to ∼2.5C in the ocean. The globally averaged sea surface temperature (SST) increases by ∼2.1C per CO2 doubling. In addition to different surface feedbacks over land versus sea, the inhomogeneous distribution of land and sea (i.e. more land at high latitudes) acts to further differentiate between the average temperature over the land and sea surface. There is some non-linear response towards higher CO2 concentrations, consistent with, e.g. Lunt et al. (2021); Baatsen et al. (2020); Caballero and Huber (2013), which is closely related to the corresponding radiative forcing in the model. Starting from the pre-industrial reference, a single CO2 doubling results in an initial perturbation of 3.49 W m−2, while a quadrupling of CO2 induces 7.93 W m−2, meaning that the radiative forcing of 4 × CO2 is 2.27 (rather than 2) times that of 2 × CO2 (Baatsen et al.2020; Etminan et al.2016). Using the corresponding deviation of the globally averaged near-surface air temperature (2.93 and 6.34 C in the E560 and E1120, respectively), we find a climate sensitivity parameter of 0.80–0.84 K W−1 m2. Using the E1120 extrapolated temperature, we report an estimated equilibrium climate sensitivity (ECS) of 3.17 C per CO2 doubling, noting the non-linearity towards higher pCO2. Comparing the Eoi560 Eoi280 temperatures, we find a difference of 3.2 C, indicating a similar ECS between the modelled pre-industrial and Pliocene states. The observation that the Pliocene ECS tends towards the higher regime of pre-industrial (i.e. E1120 E560 rather than E560 E280) suggests that the non-linearity in ECS is a function of the reference temperature rather than pCO2, as the Eoi280 temperatures are comparable to the E560 ones.

Table 1Overview of globally averaged observables and equilibration measures for all of our simulated cases, including simulation length (+ denotes a continuation from the above control case), mean annual 2 m air temperature (MAT), average tropical temperature (MATT), average polar temperature (MATP), mean annual precipitation (MAP), top-of-model net radiative flux (RTOM), sea surface temperature (SST), full-depth ocean temperature (OT), and ocean temperature drift. Subscript “e” denotes estimated equilibration values.

Download Print Version | Download XLSX

Strikingly, the applied mid-Pliocene model boundary conditions result in an average warming similar to a doubling of atmospheric pCO2 regardless of background CO2 level. Looking at the upper-ocean (<1000 m) versus deep-ocean (>2000 m) average temperatures in Fig. 3, the offset between pre-industrial and mid-Pliocene conditions seems depth-dependent. While the upper ocean is about 2 C warmer in our mid-Pliocene simulations, the deep ocean is over 3.5 C warmer compared to the pre-industrial cases. Much of this warming pattern can be explained by taking into account the extrapolated equilibrium temperatures of the E280,P simulation as well (crosses in Fig. 3, see also Fig. S4). The vertical redistribution of heat can be mostly attributed to the vertical mixing parameters in the mid-Pliocene model set-up. After correcting the pre-industrial temperatures with the E280,P E280 difference, the average offsets between pre-industrial and mid-Pliocene simulations again closely resembles that of a CO2 doubling (i.e. ∼2.5C; bracketed values in Fig. 3). A more detailed look into the distribution of different temperature responses to only a change in atmospheric pCO2 or other model boundary conditions is presented in Sect. 4.5.

Compared to the E280 pre-industrial reference, the simulated climate of our Eoi400 mid-Pliocene control is on average almost 5 C warmer at the surface. This is accompanied by a 13% increase in mean annual precipitation and a considerable polar amplification factor. The latter is found by comparing the change in average polar (> 66.5 N/S) to tropical (< 23.5 N/S) temperature, i.e. MATP/MATT, between different cases. This gives a polar amplification factor of 3.11 and 3.17 for the E560 and E1120, respectively, compared to the E280. A slightly smaller value of 2.85 is found between the Eoi560 and Eoi280, but a much larger polar amplification factor of 4.75 is seen between the Eoi400 and E280, including the effect of changing the model boundary conditions. This value increases further to 5.78 (Eoi280 E280) or 5.42 (Eoi560 E560) between Pliocene and pre-industrial cases at equal atmospheric pCO2, indicating a strongly enhanced temperature response towards high latitudes to the applied Pliocene model boundary conditions.

As shown in Fig. 3 and Table 1, the full-depth average of the oceans is about 2.5 C warmer in our Pliocene simulations compared to the equivalent pre-industrial case at equilibrium. Differences in land ice cover and vegetation type but also snow and sea ice coverage provide a substantial contribution to the modelled Pliocene warmth (see also Table S2 in the Supplement). Shortwave surface fluxes alone account for a  6 W m−2 net forcing in the globally averaged radiative balance, comparing Pliocene to pre-industrial cases at equal pCO2 (Eoi280 E280: 6.6 W m−2; Eoi560 E560: 5.4 W m−2). The responsible surface albedo feedback thus plays a primary role, decreasing slightly towards a warmer reference state. Although small, shortwave cloud feedbacks (Eoi280 E280: −1.2 W m−2; Eoi560 E560: 0 W m−2) counteract the reduced effect of surface albedo towards warmer states, helping to explain why we see a similar climate sensitivity across all of our model cases. Longwave radiative fluxes play a similar role to enhance Pliocene warmth in our simulations, mainly through the lapse rate and water vapour feedback. A warmer ocean surface and atmosphere lead to an increase in the total column water vapour, which acts as a greenhouse gas. In contrast to shortwave fluxes, this effect becomes larger towards warmer conditions (Eoi280 E280: 11.4 W m−2; Eoi560 E560: 12.1 W m−2) but is again mitigated by small longwave cloud feedbacks (Eoi280 E280: 0.7 W m−2; Eoi560 E560: 0.2 W m−2).

Figure 3Overview of global average temperature in all of our different pre-industrial and Pliocene cases compared to the E280. Open markers connected by thin lines show the time mean at the end of each simulation. Filled markers connected by thick lines show the corresponding equilibrium values. Cool colours are used for pre-industrial (PI) cases, and warm colours are used for the Pliocene (Plio) ones. We consider four different temperature measures: near-surface air temperature, upper-ocean, full-ocean, and deep-ocean area/volume weighted averages. Crosses indicate extrapolated temperatures from the E280,P case. Pliocene–pre-industrial offset values are given, determined by the mean temperature difference between cases at equal CO2 (using a log interpolation between 280 and 560 ppm to estimate the missing E400). A similar offset, corrected for the difference between the E280 and E280,P, is shown in brackets.


4.2 Simulated mid-Pliocene conditions

4.2.1 Atmosphere

Looking at the annual mean near-surface air temperature in the Eoi400 versus E280 case (Fig. 4a, b), we can identify the direct influence of many of the changes made in the PlioMIP2 mid-Pliocene boundary conditions (see Fig. S7 for a side-by-side comparison between the E280, E280,P, and Eoi280 cases). The strongest temperature deviations are seen over parts of Antarctica and Greenland in the absence of land ice. Relatively cool temperatures in the Eoi400 can be found over East Antarctica, where we have implemented a thicker ice sheet compared to the pre-industrial conditions. We see similar temperatures between the E280 and Eoi400 over much of the land at low latitudes, which are mainly a result of surface properties in the Pliocene boundary conditions that counteract the overall warming (e.g. lakes in Africa, vegetation in Australia; see Supplement Fig. S1). Warmest month mean temperatures of >40C occur over many low-latitude and midlatitude continental regions, while coldest month mean temperature contours steadily migrate poleward in the Eoi400. As shown in Table 1, we see a considerable polar amplification in the overall warming pattern between the Eoi400 and E280, with mostly small temperature differences in the tropics increasing to > 10 C over both the Arctic and Southern Ocean.

Figure 4Annual mean atmospheric model fields from our Eoi400 case (a, c, e) and the Eoi400 E280 difference. (a, b) Near-surface (2 m) air temperature (shading), warmest-month mean maximum (yellow contour at 40 C), and coldest-month mean minimum (blue contours at 0, −20, and −40C). Contours in (a) are for Eoi400, and those in (b) are for E280. Panels (c, d) show precipitation (m yr−1); differences in (d) are smoothed over land. Panel (e) shows the height of the 500 hPa pressure level (shading) and mean sea level pressure (MSLP; contours every 5 hPa; black:  1000 hPa; white: < 1000 hPa; thick line every 20 hPa). Panel (f) shows the 500 hPa height difference with respect to the global average Eoi400 E280 value (i.e. 96 m); contours are used for MSLP (every 2 hPa, solid: >0: dashed: <0; thick white line at 0 hPa).

Besides being warmer, the Eoi400 is substantially wetter compared to the E280 (Fig. 4c, d), with a global average annual precipitation of 1180 mm versus 1042 mm (i.e. a 13%; see Table 1). Polar regions are overall wetter, but this is most noticeable over the Southern Ocean. Reduced ice sheet cover over much of Antarctica and warmer surface temperatures both act to increase precipitation over much of the coastal region. An overall poleward migration of storm tracks is seen, which is most pronounced over the North Atlantic Ocean. At lower latitudes, changes in the precipitation pattern are dominated by a shift towards the Eastern Hemisphere. Regions surrounding the Indian and West Pacific oceans are much wetter in the Eoi400, while those around the East Pacific and Atlantic oceans become drier. This is related to a combined westward shift and expansion of the Walker circulation in our mid-Pliocene simulations (shown in most PlioMIP2 simulations by Han et al.2021). Like many other comparable model studies, our E280 shows a double Intertropical Convergence Zone (ITCZ) over the Pacific Ocean (e.g. Bellucci et al.2010; Haywood et al.2020). The southern branch of the ITCZ is weaker over the Pacific Ocean in the Eoi400, while shifting westward and southward towards a pronounced South Pacific Convergence Zone (SPCZ). A more active Indian Monsoon can explain increased rainfall over the Middle East and South Asia, as well as some of the cooler temperatures in northern India. Higher precipitation rates in the Eoi400 are also seen over Australia and much of northern Africa, which are directly related to the altered land surface properties (see also Fig. S1 in the Supplement).

In addition to the direct effects of the PlioMIP2 boundary conditions, some more general changes in the atmospheric circulation are seen between the E280 and Eoi400 cases (Fig. 4e, f – note that the global average change in 500 hPa height is subtracted in f). A reduction in the Equator-to-pole gradient is not only present in surface temperatures but also in mid-tropospheric height surfaces. This translates to an overall reduction of baroclinic instability, despite some of the increases in precipitation seen at midlatitudes and high latitudes. Reduced surface elevation over much of West Antarctica and Greenland leads to a lower surface pressure in the Eoi400, while the opposite is seen over parts of East Antarctica. These changes in the surface topography also likely influence the circulation, as several stationary ridges are evident in the 500 hPa height difference. Preferred midlatitude ridging is evident over the North Pacific Ocean in the Eoi400, which is in turn reflected in the temperature and precipitation differences with the E280. In addition to the advective adjustment (i.e. warmer and wetter west of the high-pressure and geopotential anomaly), the radiative feedback and subsidence act to also increase temperatures towards the centre of the anomaly.

4.2.2 Ocean

Annual mean sea surface temperatures (SSTs) reflect the patterns and changes seen in the Eoi400 and E280 atmospheric temperatures (Fig. 5a, b; a similar comparison using the E280,P instead can be found in Supplement Fig. S8). In contrast to near-surface air temperature, the Arctic Ocean SSTs show the smallest difference between the Eoi400 and E280 cases despite large reductions in sea ice cover. Another area with relatively cool temperatures in the Pliocene is the eastern tropical Pacific Ocean, as it is still  2 C warmer compared to the E280. Much warmer SSTs are seen over much of the North Atlantic Ocean, as well as the Northwest Pacific Ocean. Apart from the West Antarctic coastal region, the entire Southern Ocean is  4–8 C warmer in our Pliocene control.

Figure 5Annual mean oceanic model fields from our Eoi400 case (a, c, e) and the Eoi400 E280 difference: (a, b) sea surface temperature (shading) and sea ice fraction (contours; white: 0.15; black: 0.5 and 0.9). Contours in (a) show Eoi400, contours in (b) show E280. Panels (c, d) show sea surface salinity (shading) and mixed-layer depth (thin line: 100 m; thick black line: 250 m; thick white line: 500 m). Panels (e, f) show the barotropic stream function (BSF; shading) and zonal wind stress (contours every 5×10-2 Pa; dashed: <0; thick white line at 0 hPa). Contours in (e) show Eoi400, and contours in (f) show E280.

Sea surface salinity (SSS) patterns in the Eoi400 are mostly driven by precipitation patterns and gateway changes (Fig. 5c, d) between the E280. High precipitation amounts over most of Asia and the Indian Ocean cause low SSS through surface fluxes and river runoff. The Arctic Ocean is less saline compared to the E280, partly as a result of increased runoff but also due to the closed Bering Strait and Northwest Passage. The outflow of these low-salinity Arctic surface waters can be seen to the east and south of Greenland, resulting in a rather complex interaction with much more saline waters in the northern Atlantic Ocean. The Pacific Ocean is in general more saline in the Eoi400, with regional changes from the E280 that are the combined effect of precipitation patterns, enhanced river runoff, and the Bering Strait closure.

The depth-averaged flow, represented by the barotropic stream function (BSF) in Fig. 5e, f is very similar between the Eoi400 and E280 cases. Regardless of the reduced meridional gradient in surface temperatures, the Antarctic Circumpolar Current (ACC) is slightly stronger in our Pliocene control. In contrast to what is seen at the surface, the depth-averaged density gradient is steeper in the Eoi400 and therefore enhances the density-driven component of the ACC. Southern Hemisphere subtropical gyres are slightly weaker, while those in the Northern Hemisphere show a poleward extension in the Eoi400. North Atlantic subtropical and subpolar gyres are both stronger compared to the E280.

4.2.3 Sea ice

Even with the relatively simple model set-up used here, our E280 simulation has realistic sea ice cover compared to late 20th century observations from Peng et al. (2013), as shown in Fig. 6. The mean sea ice maximum extent corresponds well with the observed modern-day maximum in both hemispheres, but the minimum is overestimated, especially in the Southern Hemisphere. The disagreement may, however, be mostly explained by the difference between pre-industrial and modern conditions.

Figure 6(a) Eoi400 mean maximum monthly sea ice fraction over northern and southern polar regions. Coloured contours show coldest-month mean SST, and the black contour indicates the late 21st century observed September (NH) or March (SH) sea ice edge at an 0.15 fraction (indicated on the colour bar; Fetterer et al.2017). Panel (b) is the same as (a) but for the E280 case using the same contour intervals and sea ice edge data. Panels (c, d) are similar to (a, b) but show the monthly sea ice minimum, warmest-month mean SST and late 21st century observed March (NH) or September (SH) sea ice edge.

Sea ice cover is drastically reduced in our Eoi400 simulation compared to the E280 one. During the late summer minimum extent, sea ice concentrations drop to nearly zero across both polar regions, with only some remaining over exposed waters in West Antarctica. Changes in the maximum sea ice cover are less dramatic over the Arctic region, while the Southern Ocean still has much less sea ice cover in the Eoi400 compared to the E280 case. This corresponds well with the  4–8 C warmer SSTs seen across the Southern Ocean in our Pliocene simulations (Fig. 5b). The isolated nature of the Arctic Ocean allows its surface waters to remain relatively cool, with some sea ice persisting until late summer. Despite substantial refreezing during wintertime, open waters in the early winter months prevent the air from cooling down further. This explains the much larger near-surface air temperature anomalies over the Arctic region between our Eoi400 and E280 simulations (Fig. 4b).

4.3 Meridional overturning and heat transport

4.3.1 Global

Our Eoi400 case is characterised by a stronger and deeper global meridional overturning circulation (MOC) compared to the E280. This is reflected in the global meridional overturning stream function (MSF) in Fig. 7a, b by the dominant northern overturning cell north of 40 S. The deep southern (i.e. negative) overturning cell linked to Antarctic bottom waters is slightly stronger in the Eoi400 as well. Note that the appearance of this overturning cell is mostly masked in the Southern Ocean by strong Ekman upwelling linked to the ACC. The Eoi560 MOC is similar to the Eoi400 one but is slightly weaker overall, while the opposite holds for the Eoi280. This could be partly due to the transient effect of the simulations but suggests a tendency for the MOC to get weaker towards higher atmospheric pCO2.

Figure 7(a) Global oceanic meridional overturning stream function (MSF) of the Eoi400 case (shading) and Eoi280 Eoi560 difference (contours every 2 Sv; black: positive; green: negative), taking the mean over the last 500 model years. Panel (b) is the same as (a) but for the E280 case and the Eoi280 E560 difference. Note the vertical stretching of the upper 1 km. (c) Oceanic meridional heat transport in the Eoi400 (solid lines; global: black; Indo-Pacific: red; Atlantic Ocean: blue), E280 (dashed lines, similar colouring), and Eoi560 (dotted line, only global). (d) Total meridional heat transport, induced by the top-of-model (TOM) net radiative fluxes (Eoi400: solid black; E280: dotted grey) and the corresponding atmospheric heat transport (Eoi400: solid red; E280: dotted blue). The difference in atmospheric heat transport with respect to the Eoi400 is shown for the E280 (dashed pink) and Eoi560 (dash-dotted pink), both of which are magnified 10-fold.


The differences in MSF pattern and strength between the Eoi400 and E280 cases are only partly reflected in the corresponding oceanic heat transports (OHTs), which are shown in Fig. 7c. Apart from a slight weakening of poleward heat transport at southern low latitudes (mostly in the Indo-Pacific sector), the Eoi400 OHT is very similar to the E280. The warmer Eoi560 case has an overall reduced poleward oceanic heat transport, but the relative change is small (<0.1 PW). The total meridional heat transport (MHT, shown in Fig. 7d) is slightly weaker in the Eoi400 compared to the E280. This is a direct response to the reduced meridional temperature gradient seen in the Eoi400, demanding a lower poleward heat transport on both hemispheres. As there are only minor changes in the oceanic component, the reduced total MHT must be mainly accounted for by the atmosphere. This is reflected by the atmospheric MHT difference in E280 with respect to Eoi400, having mostly the same sign as the total MHT. The only exception is seen at southern low latitudes, where a stronger southward OHT is compensated by a northward atmospheric MHT in the E280. The Eoi560 atmospheric MHT shows a small net positive offset with respect to the Eoi400 that may be explained by some remnant warming trend. In agreement with a further reduction of the meridional temperature gradient, the warmer Eoi560 case thus has a further reduced poleward MHT, although this transport small ( 0.1 W m−2) and is mostly accounted for by the ocean.

4.3.2 Atlantic Ocean

The Atlantic meridional overturning circulation (AMOC; Fig. 8a) clearly encompasses the entire deep northern overturning cell in the global MSF north of 30 S (Fig. 7a). The AMOC mean state is overall stronger in our Pliocene versus pre-industrial simulations, but its temporal variability is also considerably higher in all of the Pliocene cases (22.7 ± 1.6 Sv versus 18.3 ± 0.7 Sv; see also Fig. S5 in the Supplement). The Eoi560 AMOC undergoes a series of especially large (>10 Sv) intensity swings, although this is likely caused by the thermal adjustment of the deep ocean, and it is therefore mostly a transient model feature. These large swings stopped after  1000 model years, but the AMOC variability remains the largest of any of our cases (σ=2.4 Sv). The E280,P AMOC also exhibits a larger variability (σ=1.0 Sv) compared to the other pre-industrial cases, but it is still smaller than any of the Pliocene ones. A more detailed discussion on AMOC variability is provided in Sect. 4.6.2

Figure 8(a) Mean AMOC stream function over the last 500 years of the Eoi400 (shading) and E280 simulation (contours highlighted on the colour bar). (b) Difference in AMOC stream function with respect to the E280 for the Eoi400 (shading) and the E280,P (contours). (c) Eoi400 E280 (100-year mean) potential density difference, averaged over the upper 100 m and relative to the global average change. Contours show the annual mean maximum mixed-layer depth for the E280 (solid black), E280,P (solid blue), and Eoi400 (dotted black) and the 0.15 sea ice fraction (yellow; solid: E280; dashed: Eoi400). (d) Freshwater budget over the Arctic Ocean for the (left to right) E280, E280,P, and Eoi400 cases, including the net surface freshwater flux (SFWF) and transports across the Bering Strait, Canadian Archipelago (CAA), and Greenland–Scotland Ridge (GSR). The mean over the last 500 model years is taken, with the standard deviation of yearly averages indicated by the error bars.

The difference in AMOC streamfunction between Eoi400 and E280 is most pronounced between 1 and 3 km depth (Fig. 8b). The stronger and deeper AMOC cell can be partly explained by the altered vertical mixing parameters in our Pliocene model set-up, as the E280,P shows a reduced but very similar difference pattern with respect to the Eoi400. This is consistent with the stronger AMOC in the POP1 type vertical mixing simulations of Chandan and Richard Peltier (2017). In general, we find that the altered vertical mixing scheme in our Pliocene simulations does have an impact on the strength and behaviour of the AMOC, but not to the extent that is seen in any of the Pliocene cases. It is therefore likely that the altered boundary conditions and resulting circulation changes have a considerable impact on the AMOC strength and behaviour. A stronger mid-Pliocene AMOC was found consistently in the PlioMIP2 ensemble by Zhang et al. (2021), who also found a link with higher North Atlantic SSTs but no clear relation with Atlantic OHT. The latter is investigated in more detail by Weiffenbach et al. (2022), who decompose the OHT into the contributions from the overturning and gyre circulations, respectively. The PlioMIP2 boundary conditions applied here include the closure of several high-latitude oceanic gateways. This results in lighter upper-ocean waters across the Arctic Ocean in the Eoi400, but denser waters across most of the North Atlantic Ocean (Fig. 8c). The integrated net surface freshwater flux over the Arctic Ocean more than doubles in the Eoi400 case compared to the E280 one (Fig. 8d). This more than compensates for the missing transport from the Bering Strait closure. Lacking connectivity through the Canadian Archipelago demands for an equally large net southward freshwater transport across the Greenland–Scotland Ridge. The outflow of those light, low-salinity Arctic waters pushes the deep-water formation zone towards the south across the Labrador Sea into warmer and deeper waters. Increased salinity in the Labrador Sea and thus higher potential densities are tied to the closure of the Canadian Archipelago, as shown by the negative component in the E280 Arctic Ocean freshwater balance (i.e. net salt transport from the Labrador Sea). Such a southward shift of the deep-water formation zone is not present in the E280,P case.

4.4 Model–proxy comparison

As shown by Haywood et al. (2020), our Eoi400 case performs well when comparing the annual mean SSTs to the available PlioMIP2 time-specific proxy records (Foley and Dowsett2019; McClymont et al.2020). The zonally averaged, annual mean SST from our different simulations is shown in light of these proxies in Fig. 9a. Considering site-specific rather than zonally averaged SSTs from the Eoi400 simulations, we find a root-mean-square error (RMSE) of 2.7 C and a mean absolute error (MAE) of 2.2 C (using the combined U37K and Mg/Ca proxy record of McClymont et al.2020). Moreover, there is no significant warm or cold model bias across the SST or latitude range covered by the proxies. Relatively warm North Atlantic SSTs in the model are well reflected by the proxy record, consistent with the other CCSM and CESM models within PlioMIP2 as shown by de Nooijer et al. (2020). They also find that the CCSM and CESM model family best captures mid-Pliocene Arctic SST proxies, which is likely related to low sea ice cover.

Figure 9(a) Proxy-based annual mean SST reconstructions (yellow markers; diamonds: U37K; squares: Mg/Ca) versus model-derived zonally averaged SST for Eoi400 (black, including zonal variation in grey), E280 (blue), E560 (dotted green), Eoi280 (dash-dotted pink), and Eoi560 (dashed orange). The inset shows a point-wise comparison, including a simple linear regression. (b) Pliocene proxy–pre-industrial SST (ERSST; late 19th century reanalysis) and modelled zonal average SST difference from E280 (colour scheme is the same as in a; dash-dotted purple: E280,P). (c) Pliocene–pre-industrial SST difference over the Atlantic Ocean: proxy–pre-industrial (coloured markers) and model-derived Eoi400 E280 (shading).

Looking at the SST difference between Pliocene and pre-industrial conditions, rather than absolute values, reveals some more discrepancies between proxy records and our Eoi400 simulation (Fig. 9b). Most of the high positive temperature anomalies at northern midlatitudes and high latitudes remain well captured by the model (as do the much smaller differences at low latitudes). However, other sites at both northern and southern midlatitudes show a much poorer agreement between model-based and proxy-indicated Pliocene–pre-industrial SST differences. The lack of proxy data across most of the Southern Ocean makes it impossible to assess whether the large SST differences between our Pliocene and pre-industrial simulations are realistic. Despite the visually poor agreement compared to absolute mid-Pliocene SSTs, the RMSE and MAE are only slightly higher for the Pliocene–pre-industrial SST difference.

Surprisingly large Pliocene–pre-industrial SST differences (considering a limited CO2 effect) were found in the PlioMIP1, especially over the North Atlantic Ocean, but these are not well reflected in the model ensemble (Haywood et al.2013b). This seems to be greatly improved now, as shown by the large positive anomalies over much of the midlatitude to high-latitude North Atlantic Ocean in Fig. 9c. The largest discrepancies between proxy-based and model-derived Pliocene–pre-industrial SST differences are seen in the Mediterranean Sea and off the coasts of North America and southern Africa. These discrepancies can be largely explained by coastal upwelling and boundary currents, which are poorly resolved in our simulations because of the limited horizontal resolution (e.g. McClymont et al.2020; Li et al.2019). Although absolute Pliocene SSTs seem to agree better, they are already much too warm at some of these locations in our E280 case compared to reanalysis data.

Near-surface air temperatures reveal the same warming pattern between the Eoi400 and E280 case as is seen in the SSTs over low latitudes and midlatitudes (see also Fig. S9 in the Supplement). The highest temperature anomalies shift poleward over the Arctic Ocean and the Antarctic coastal region due to sea ice melt and reduced ice sheet cover. Our coolest Pliocene simulation (Eoi280) is similar to (or even warmer than) the E560 case at all latitudes. This indicates the importance of the implemented mid-Pliocene model boundary conditions and long-term impacts such as the partial loss of ice sheets.

4.5 Sensitivity to Pliocene boundary conditions versus atmospheric pCO2

Using our set of pre-industrial and Pliocene sensitivity simulations (i.e. E280, E560, Eoi280, and Eoi560), we can make an assessment of the effects of the applied Pliocene boundary conditions (Eoin En) versus external radiative forcing (X560 X280), as well as any state dependency in the model response. The difference in annual mean near-surface air temperature is shown in Fig. 10. Figure 10a and b show the modelled temperature response to a shift from pre-industrial to Pliocene boundary conditions (excluding CO2); Fig. 10c and d show the response to a doubling of atmospheric pCO2. Consequently, comparing Fig. 10a and c to Fig. 10b and d shows the state dependency of the effect of boundary conditions to a different pCO2 baseline (Fig. 10a and b) and vice versa (Fig. 10c and d).

Some generic temperature differences between the model cases can be easily identified in all of these comparisons. Stronger contrasts are seen over land than over the ocean, as a result of the thermal capacity and much larger potential for latent heat fluxes over the ocean. Polar amplification of the temperature change also seems universal across our simulations, caused mainly by ice-albedo feedbacks and reversed lapse rate feedbacks (i.e. positive at high latitudes, negative at low latitudes; see also Fig. S7 in the Supplement). As noted before, the ECS of our simulated Pliocene climate is very similar to that of the pre-industrial one (i.e. 3.2 C versus 3.17 C per CO2 doubling). It is therefore not surprising that the overall response to a doubling of atmospheric pCO2 is comparable, with mainly some regional differences as a result of the boundary conditions and local feedback mechanisms.

Changes in precipitation patterns between the four sensitivity simulations show large differences between the response to atmospheric pCO2 and other model boundary conditions. In line with the findings of Han et al. (2021), the effect of a CO2 doubling is mainly that wet regions get wetter and dry regions get drier, with this effect primarily being seen in the tropics. Differences in precipitation between Pliocene and pre-industrial cases are generally much more pronounced and widespread but are also enhanced between the high- and low-CO2 scenarios. Some of these differences are directly related to the boundary conditions (mainly topography and ice sheets) and their temperature response (e.g. sea ice, monsoons). In addition, large-scale precipitation changes include a westward and poleward shift in the tropics, resulting in wetter tropical Indian and West Pacific ocean conditions in our mid-Pliocene simulations. Han et al. (2021) show that these dynamic changes in precipitation are related to the meridional and zonal circulation patterns (e.g. Walker circulation influences the Indian–Pacific moisture exchange). Moreover, our CCSM4-Utr simulations are found to exhibit the largest asymmetry in hemispheric energy flux within the PlioMIP2 ensemble, explaining the significant shifts of the Pacific ITCZ and SPCZ. Further north, we see a poleward migration of North Atlantic storm tracks and the dynamic response of a prevailing North Pacific anticyclone. These shifts in precipitation correspond well with what is seen in Fig. 4d.

Figure 10Annual mean near-surface air temperature difference between our four sensitivity simulations: (a) Eoi280 E280, (b) Eoi560 E560, (c) E560 E280, and (d) Eoi560 Eoi280. Contours show the difference in mean annual precipitation at 100, 200, and 500 mm, with blue indicating a positive value and yellow indicating a negative value.

We also consider the zonally averaged temperature differences between the four sensitivity cases, along with the corresponding energy balance model (EBM) analysis and the contribution from its different components (Fig. 11). All of the zonally averaged temperature responses estimated using the EBM are in near-perfect agreement with the actual surface temperature differences in the model. Simply adding the contributions from albedo, emissivity, and MHT reproduces the full temperature response well, showing that non-linear effects between them are small. Both these findings are in line with the results of Hill et al. (2014) for PlioMIP1. Albedo plays a crucial role in the response to mid-Pliocene boundary conditions, especially at the surface (i.e. vegetation, ice sheets, and sea ice). The effect of planetary albedo is about a third lower compared to the surface, showing that the latter is partly compensated for by the shortwave contribution from clouds, aerosols, and water vapour. The emissivity of longwave radiation shows a more similar response between the effect of boundary conditions versus CO2 doubling, but it is still stronger in the former. About half of the temperature contribution through emissivity after a CO2 doubling is related to the direct radiative forcing, leaving the other half as being mainly the result of water vapour and lapse rate feedbacks. These feedbacks, together with the surface albedo, are thus the main drivers behind polar amplification of the temperature difference between all of the different cases considered here. The effect is partly mitigated by meridional heat transport and cloud feedbacks, but these have little impact on the global scale. This is again consistent with the findings of Hill et al. (2014), with a larger effect of the mid-Pliocene boundary conditions in PlioMIP2 versus PlioMIP1.

Figure 11The same as in Fig. 10 but showing zonally averaged temperature differences and contributions from different components within the energy balance model (EBM). The temperature difference between our simulations (total GCM; solid grey) is compared to the one predicted by the EBM (total EBM; dashed dark grey). Components of the EBM include: planetary albedo (solid blue), emissivity (solid red), meridional heat transport (MHT, dashed green), shortwave cloud forcing (negative, solid purple) and longwave (dashed dark red) cloud forcing. The sum of the contributions from albedo, emissivity, and MHT is also shown (Sum EBM; dashed black). Globally averaged values are shown for each component, with bracketed values added for surface albedo and emissivity excluding the effect of CO2 doubling.


Most of the effects of the Pliocene model boundary conditions on temperature (Fig. 10a, b) agree well with those seen earlier between our E280 and Eoi400 cases (Fig. 4b). The warming effect of removing or lowering ice sheets at high latitudes and the cooling effect of introducing lakes or removing desert at low latitudes are independent of the atmospheric CO2 level. Over the Arctic Ocean, reduced sea ice in the Pliocene simulations has an increased impact on temperatures above the surface towards lower pCO2 (Fig. 10a versus b). The opposite is seen over parts of the Southern Ocean, where sea ice cover is still relatively large in the E560 (see also Fig. S10 in the Supplement). Temperature changes that can be related to altered circulation patterns in the Pliocene versus pre-industrial conditions (e.g. monsoons, the South Pacific Convergence Zone, midlatitude ridging, storm tracks) are robust between different CO2 levels. The temperature response to a doubling of atmospheric pCO2 is very similar at low latitudes and midlatitudes between our pre-industrial and Pliocene simulations (Fig. 10c, d). Differences in sea ice cover (and to a lesser extent land ice and surface properties) amplify high-latitude warming in response to a CO2 doubling; this effect is most prominent over the Arctic in the pre-industrial cases and over the Southern Ocean in the Pliocene ones.

The near-surface air temperature differences between our four sensitivity cases (Fig. 10) are reflected in the SST (see Fig. S10 in the Supplement). The most prominent exception to this consistency is found over high-latitude surface waters, which are seen to change in temperature much less than the air above. This can be explained by the decoupling effect of sea ice, which tends to dampen the SST differences (but with large seasonal differences). Pliocene–pre-industrial SST differences over the Pacific Ocean are consistent between CO2 levels and are noticeably different from the effect of a pCO2 doubling. These patterns suggest a shift in the background state of both the ENSO and Pacific Decadal Oscillation (PDO) in our Pliocene simulations (see also Sect. 4.6). An even more prominent and consistent SST response to the Pliocene boundary conditions is found over the North Atlantic Ocean (also seen in Fig. 5b). This is likely the combined result of a stronger Pliocene AMOC and the isolation of the Labrador Sea through closure of the Canadian Archipelago. Large changes in SST between our Pliocene and pre-industrial simulations can be linked to the upwelling of relatively warm deep waters south of the ACC and is reflected in the air temperatures above (Fig. 10).

4.6 Modes of internal variability

4.6.1 El Niño–Southern Oscillation

In our simulated Pliocene climate both the mean state and several modes of variability are different from what is seen in the pre-industrial reference. As it plays a primary role in the tropics, we first look at the behaviour of the El Niño–Southern Oscillation (ENSO). In Fig. 12, time series of the Niño 1+2 and Niño 3.4 indices are shown for the E280 and Eoi400 case. The amplitude of both ENSO indices is greatly reduced in our Eoi400 versus E280 case (σ1+2: −54 %; σ3.4: −68 %). This is the largest such reduction seen among the PlioMIP2 ensemble (Fig. 2a of Oldeman et al.2021; OB21). Furthermore, the occurrence of strong and long-lasting (> 1 C; > 1 year) El Niño and La Niña events completely disappears in our mid-Pliocene simulation. Both of these findings are consistent among all of our mid-Pliocene (i.e. Eoi280, Eoi400, Eoi560) and pre-industrial (i.e. E280, E280,P, E560, E1120) cases (not shown).

Figure 12(a) ENSO time series for the E280 (blue) and Eoi400 (red) using monthly SST anomaly fields and a 5-month running mean. Note the different scaling. Temperature intervals are kept the same for visual comparison. (b) Multi-taper power spectrum, using 200 years of monthly data, including the 90 %, 95 %, and 99 % confidence levels. (c) Corresponding ENSO SST patterns, taking the mean over the 10% highest (shading) and 10% lowest (contours) monthly Niño 3.4 index values.

Spectral analysis of the Niño 3.4 indices also shows large differences between the ENSO behaviour in the E280 and Eoi400 case. The modelled pre-industrial ENSO variability is characterised by a broad spectral peak at periods of 3–10 years, with three statistically significant (at 99 % confidence) peaks and a dominant period of around 6 years. A similar broad peak is seen for the mid-Pliocene ENSO but with significant variability only at 4 years (99 %) and 9 years (95 %). Moreover, there is an increase of significant variability in the Eoi400 case at shorter periods of around 2 years or less. This agrees well with the predominantly weak and rather high frequent behaviour of the modelled Pliocene ENSO seen in the Niño 3.4 time series. A similar drop in spectral power of the mid-Pliocene ENSO is seen consistently in the PlioMIP2 ensemble, but this is more focussed around the 3–5-year period (Figs. 3 and 4 of OB21).

As would be expected from both of the Niño time series shown for each case, the amplitude of SST anomalies is seen to be much smaller throughout the tropical Pacific Ocean in our Eoi400 case versus the E280 case. Looking at the mean SST pattern, corresponding to the 10 % average highest and lowest monthly Niño 3.4 values, allows us to assess both the spatial distribution and strength of ENSO variability at the same time. In addition to an overall highly reduced amplitude, the associated region of SST variability expands westward and poleward across the central and western Pacific Ocean in our Pliocene simulation. This is consistent with the pattern shifts shown by the majority PlioMIP2 models (Fig. 5 of OB21). The mean difference in SST between the Eoi400 and E280 (Fig. 5b) over the main region of ENSO variability is smaller than the average in the tropics. The background state of the tropical Pacific Ocean is therefore not like El Niño but is rather slightly like La Niña. Only few models within the PlioMIP2 ensemble show such a strong La Niña-like pattern in the mean temperature of the Eoi400, but these are also the models with the largest reduction in ENSO variability (see also Figs. 9–11 of OB21).

4.6.2 Interannual to multi-decadal SST variability

Using 500 years of annual mean SSTs and AMOC strength, we assess the dominant modes of variability in the North Pacific and North Atlantic oceans (Fig. 13). We consider the dominant eigenmode of SST anomalies in the North Pacific Ocean to represent the PMV, and we consider the one that best correlates with average North Atlantic SSTs as the AMV. The explained variance by the PMV mode is seen to decrease in the Eoi400 versus E280 case, while that of the AMV greatly increases. As noted earlier, the AMOC variability is considerably larger in the Eoi400 case compared to the E280 one. Correlations between the PMV, AMV, AMOC, and annual Niño 3.4 indices can be found in Fig. S12 in the Supplement.

Figure 13Time series (upper panels; thick lines using a 20-year smoothing window), multi-tapered power spectra (middle panels), and spatial correlation patterns (lower panels) of North Pacific SST (left), North Atlantic SST (middle), and AMOC strength (right) variability for 500 years of annual SST data from our E280 (blue) and Eoi400 (red) simulations. The percentage of variance explained by the respective principal component is shown in black alongside the PMV and AMV time series. Spatial SST patterns show the temporal correlation between the SST at each location and the different variability indices. Boxed numbers also indicate the total correlation between each of these indices (i.e. time series in the upper panels.)

Spectral analysis shows that the E280 PMV has a broad peak at 10–20 years and therefore strongly resembles the PDO. This broad peak is seen to shift towards longer periods for the Eoi400 PMV, which is now accompanied by a sharper peak at just under 10 years. Looking at the patterns of SST variability associated with the PMV, an almost global footprint in the E280 is reduced to mostly a North Pacific signature in the Eoi400. This is likely related to the strongly decreased ENSO variability in the Eoi400 versus E280 case, as both the low-latitude SST patterns and the spectral properties at periods < 10 years show a high resemblance between PMV and ENSO in both simulations. This is supported by significant correlations (E280: 0.48; Eoi400: 0.34) between the PMV and yearly Niño 3.4 indices. The role of the AMV on global SST variability is much more limited compared to PMV but is again at its strongest and is more widespread in the E280 case. Meanwhile, the Eoi400 has a much higher percentage of explained variance and is strongly connected to AMOC variability. The spatial pattern of SST anomalies associated with the AMOC strength shows little resemblance to either the PMV or AMV in the E280 case, but both of their time series are significantly (yet weakly) correlated. In contrast, the AMOC strength strongly correlates with the AMV for the Eoi400, and their associated spatial patterns show a large resemblance. It is thus clear that the stronger and more variable AMOC in our Pliocene simulations has a profound impact on not only the mean North Atlantic SSTs but also their variability. Additionally, it seems that the different modes of SST variability are less connected between ocean basins in the Eoi400 versus E280 case.

5 Conclusions

We have completed a set of simulations using the CCSM4 and CESM1.0.5 within the PlioMIP2, including E280, E560, Eoi280, Eoi400, and Eoi560 cases (and E280,P and E1120 as additional sensitivity experiments). Our simulations show a warm mPWP climate that is 4.7 C warmer compared to our pre-industrial control in terms of globally averaged near-surface air temperature. With an estimated ECS of 3.1 C per CO2 doubling, which is consistent throughout the different pre-industrial and Pliocene cases, much of the mPWP signal is thus related to the effect of the applied PRISM4 boundary conditions. The effect of those boundary conditions is mainly to increase temperatures by about 3 C on average globally and 2.4 C in the ocean. The vertical distribution of oceanic temperatures is altered in our Pliocene simulations by the use of enhanced background vertical diffusivity. Our E280,P simulation can explain most of this differential warming pattern between pre-industrial and Pliocene cases. The enhanced vertical mixing acts to warm the deep ocean at the expense of the upper ocean, without affecting the globally averaged surface temperature. Most of the other results are thus not significantly influenced by the altered vertical diffusivity, with only regional differences seen at the surface that become negligible throughout the rest of the atmosphere.

The combined effect of increased atmospheric pCO2 and altered boundary conditions not only make our Eoi400 case warmer than the E280 one but also induce substantial polar amplification. This is mainly the result of latitudinally dependent surface albedo (ice and vegetation), water vapour, and lapse rate feedbacks. Cloud feedbacks are small in the global average, with the shortwave component acting to mitigate high-latitude warmth. The temperature differences are most pronounced where ice sheets are removed (including a large elevation effect) and in high-latitude oceans, where sea ice cover is greatly reduced in all of the Pliocene cases. The shortwave radiative feedbacks resulting from the mid-Pliocene boundary conditions are much stronger compared to those from a CO2 doubling. Furthermore, the latitudinal dependence in all of the considered radiative components is stronger, helping to explain the relatively warm high-latitude regions in our Pliocene simulations. The Eoi400 climate shows an increased precipitation rate with respect to the E280 one. Both equatorial and high-latitude regions are generally wetter in the Pliocene, while the subtropics become drier (with a distinct poleward shift of the storm tracks). Particularly high precipitation amounts are seen over most of the Indian Ocean and it surroundings, owing to enhanced North African and South (and Southeast) Asian monsoons in the Pliocene.

The warmer Eoi400 climate with a reduced meridional temperature gradient compared to the pre-industrial reference agrees well with the available proxies. A strong Pliocene warming signal over the North Atlantic Ocean is seen in both the proxy record and our simulations, which can be linked to an enhanced AMOC. This stronger AMOC, however, is not linked to an overall increase in oceanic meridional heat transport. Moreover, the total (atmosphere plus ocean) top-of-model-induced meridional heat transport is slightly reduced in the warmer simulations, in line with the reduced meridional temperature gradient. The stronger and more variable AMOC in our Pliocene versus pre-industrial simulations can in part be explained by the altered vertical diffusivity parameter. However, it is mostly the result of the applied mid-Pliocene boundary conditions, in particular the closure of several Arctic gateways.

In addition to differences in the mean state between our Eoi400 and E280 cases, there are clear shifts in the different modes of variability studied here: ENSO, PMV, and AMV. The ENSO amplitude is greatly reduced in our Pliocene simulations and characterised by shorter periodicity compared to the pre-industrial reference. The corresponding spatial pattern is also spread out across much of the tropical Pacific Ocean. Closely related to ENSO is the PMV in our E280 case, both having a distinct fingerprint on global SST variability on various timescales ranging from annual to multi-decadal. This teleconnection is lost in the Eoi400 case, with the PMV influence being mostly confined to just the North Pacific Ocean. Meanwhile, the AMV shows a strong connection to AMOC variability in our Pliocene simulations. Their mutual influence seems to be the dominant source of SST variability in the Eoi400, as opposed to ENSO or PMV in the E280.

As most of the PRISM4 boundary conditions, which are applied here as an external forcing to the model simulation, are in fact the result of long-term feedbacks (i.e. ice melt and vegetation changes), the Eoi400 can serve as a good analogue for future climatic changes. Our simulations show not only a strong warming compared to the pre-industrial reference but also considerable regional changes and shifts in the dominant modes of variability.

Data availability

PlioMIP2 model data, including those of the simulations presented here, can be downloaded from the server located at the School of Earth and Environment of the University of Leeds. Please contact Alan Haywood ( for access. The last 100 model years of our E280, E560, Eoi280, Eoi400, and Eoi560 are available within the dataset. PlioMIP2 data from CESM2, EC-Earth3.3, NorESM1-F, IPSLCM6A, and GISS2.1G can be obtained through the Earth System Grid Federation (ESGF) (, last access: 14 October 2021, Haywood et al.2020).


The supplement related to this article is available online at:

Author contributions

MAK, MLJB, and ASvdH set up the model simulations, MAK and MLJB managed the simulations and post-processing of the data. MLJB performed the analyses and set up the manuscript. All authors contributed to the final shape and contents of the manuscript.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “PlioMIP Phase 2: experimental design, implementation and scientific results”. It is not associated with a conference.


Simulations were performed at the SURFsara Dutch national computing facilities and were sponsored by NWO-EW (Netherlands Organisation for Scientific Research, Exact Sciences) (project nos. 17189 and 2020.022).

Financial support

This work was carried out within a programme of the Netherlands Earth System Science Centre (NESSC) and financially supported by the Ministry of Education, Culture and Science (OCW, grant no. 024.002.001).

Review statement

This paper was edited by Alan Haywood and reviewed by Dan Lunt and one anonymous referee.


Baatsen, M., von der Heydt, A. S., Huber, M., Kliphuis, M. A., Bijl, P. K., Sluijs, A., and Dijkstra, H. A.: The middle to late Eocene greenhouse climate modelled using the CESM 1.0.5, Clim. Past, 16, 2573–2597,, 2020. a, b, c, d, e, f, g

Bellucci, A., Gualdi, S., and Navarra, A.: The double-ITCZ syndrome in coupled general circulation models: The role of large-scale vertical circulation regimes, J. Clim., 23, 1127–1145,, 2010. a

Bitz, C. M., Shell, K. M., Gent, P. R., Bailey, D. A., Danabasoglu, G., Armour, K. C., Holland, M. M., and Kiehl, J. T.: Climate sensitivity of the community climate system model, version 4, J. Clim., 25, 3053–3070,, 2012. a

Blackmon, M., Boville, B., Bryan, F., Dickinson, R., Gent, P., Kiehl, J., Moritz, R., Randall, D., Shukla, J., Solomon, S., Bonan, G., Doney, S., Fung, I., Hack, J., Hunke, E., Hurrell, J., Kutzbach, J., Meehl, J., Otto-Bliesner, B., Saravanan, R., Schneider, E. K., Sloan, L., Spall, M., Taylor, K., Tribbia, J., and Washington, W.: The Community Climate System Model, B. Am. Meteorol. Soc., 82, 2357–2376,<2357:TCCSM>2.3.CO;2, 2001. a

Bloch-Johnson, J., Pierrehumbert, R. T., and Abbot, D. S.: Feedback temperature dependence determines the risk of high warming, Geophys. Res. Lett., 42, 4973–4980,, 2015. a

Burke, K. D., Williams, J. W., Chandler, M. A., Haywood, A. M., Lunt, D. J., and Otto-Bliesner, B. L.: Pliocene and Eocene provide best analogs for near-future climates, P. Natl. Acad. Sci. USA, 115, 13288–13293,, 2018. a

Caballero, R. and Huber, M.: State-dependent climate sensitivity in past warm climates and its implications for future climate projections, P. Natl. Acad. Sci. USA, 110, 14162–14167,, 2013. a

Chandan, D. and 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. a, b

Danabasoglu, G., Ferrari, R., and McWilliams, J. C.: Sensitivity of an ocean general circulation model to a parameterization of near-surface eddy fluxes, J. Clim., 21, 1192–1208,, 2008. a

Danabasoglu, G., Bates, S. C., Briegleb, B. P., Jayne, S. R., Jochum, M., Large, W. G., Peacock, S., and Yeager, S. G.: The CCSM4 ocean component, J. Clim., 25, 1361–1389,, 2012. a

de Nooijer, W., Zhang, Q., Li, Q., Zhang, Q., Li, X., Zhang, Z., Guo, C., Nisancioglu, K. H., Haywood, A. M., Tindall, J. C., Hunter, S. J., Dowsett, H. J., Stepanek, C., Lohmann, G., Otto-Bliesner, B. L., Feng, R., Sohl, L. E., Chandler, M. A., Tan, N., Contoux, C., Ramstein, G., Baatsen, M. L. J., von der Heydt, A. S., Chandan, D., Peltier, W. R., Abe-Ouchi, A., Chan, W.-L., Kamae, Y., and Brierley, C. M.: Evaluation of Arctic warming in mid-Pliocene climate simulations, Clim. Past, 16, 2325–2341,, 2020. a

Deser, C., Phillips, A. S., Tomas, R. A., Okumura, Y. M., Alexander, M. A., Capotondi, A., Scott, J. D., Kwon, Y. O., and Ohba, M.: ENSO and pacific decadal variability in the community climate system model version 4, J. Clim., 25, 2622–2651,, 2012. a

Dowsett, H., Robinson, M., Haywood, A., Salzmann, U., Hill, D., Sohl, L., Chandler, M., Williams, M., Foley, K., and Stoll, D.: The PRISM3D paeoenvironmental reconstruction, Stratigraphy, 7, 123–139, 2010. a

Dowsett, H., Dolan, A., Rowley, D., Moucha, R., Forte, A. M., Mitrovica, J. X., Pound, M., Salzmann, U., Robinson, M., Chandler, M., Foley, K., and Haywood, A.: The PRISM4 (mid-Piacenzian) paleoenvironmental reconstruction, Clim. Past, 12, 1519–1538,, 2016. a, b

Etminan, M., Myhre, G., Highwood, E. J., and Shine, K. P.: Radiative forcing of carbon dioxide, methane, and nitrous oxide: A significant revision of the methane radiative forcing, Geophys. Res. Lett., 43, 12614–12623,, 2016. a

Farnsworth, A., Lunt, D. J., Brien, C. L. O., and Foster, G. L.: Climate Sensitivity on Geological Timescales Controlled by Nonlinear Feedbacks and Ocean Circulation, Geophys. Res. Lett., 46, 9880–9889,, 2019. a

Fetterer, F., Knowles, K., Meier, W., Savaoie, M., and Windnagel, K.: Sea ice index, version 3, National Snow and Ice Data Center [data set],, 2017. a

Foley, K. and Dowsett, H.: Community sourced mid-Piacenzian sea surface temperature (SST) data: U.S. Geological Survey data release, Florence Bascom Geoscience Center [data set],, 2019. a

Gent, P. R. and McWilliams, J. C.: Isopycnal Mixing in Ocean Circulation Models, J. Phys. Oceanogr., 20, 150–155,<0150:IMIOCM>2.0.CO;2, 1990. a

Gent, P. R., Danabasoglu, G., Donner, L. J., Holland, M. M., Hunke, E. C., Jayne, S. R., Lawrence, D. M., Neale, R. B., Rasch, P. J., Vertenstein, M., Worley, P. H., Yang, Z. L., and Zhang, M.: The community climate system model version 4, J. Clim., 24, 4973–4991,, 2011. a

Ghil, M., Allen, M. R., Dettinger, M. D., Ide, K., Kondrashov, D., Mann, M. E., Robertson, A. W., Saunders, A., Tian, Y., Varadi, F., and Yiou, P.: Advanced spectral methods for climatic time series, Rev. Geophys., 40, 3-1–3-41,, 2002. a

Gregory, J. M., Ingram, W. J., Palmer, M. A., Jones, G. S., Stott, P. A., Thorpe, R. B., Lowe, J. A., Johns, T. C., and Williams, K. D.: A new method for diagnosing radiative forcing and climate sensitivity, Geophys. Res. Lett., 31, 2–5,, 2004. a, b, c

Han, Z., Zhang, Q., Li, Q., Feng, R., Haywood, A. M., Tindall, J. C., Hunter, S. J., Otto-Bliesner, B. L., Brady, E. C., Rosenbloom, N., Zhang, Z., Li, X., Guo, C., Nisancioglu, K. H., Stepanek, C., Lohmann, G., Sohl, L. E., Chandler, M. A., Tan, N., Ramstein, G., Baatsen, M. L. J., von der Heydt, A. S., Chandan, D., Peltier, W. R., Williams, C. J. R., Lunt, D. J., Cheng, J., Wen, Q., and Burls, N. J.: Evaluating the large-scale hydrological cycle response within the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2) ensemble, Clim. Past, 17, 2537–2558,, 2021. a, b, c

Haywood, A. M., Dolan, A. M., Pickering, S. J., Dowsett, H. J., McClymont, E. L., Prescott, C. L., Salzmann, U., Hill, D. J., Hunter, S. J., Lunt, D. J., Pope, J. O., and Valdes, P. J.: On the identification of a pliocene time slice for data-model comparison, Philosophical Transactions of the Royal Society A: Mathematical, Phys. Eng. Sci., 371, 20120515,, 2013a. a

Haywood, A. M., Hill, D. J., Dolan, A. M., Otto-Bliesner, B. L., Bragg, F., Chan, W.-L., Chandler, M. A., Contoux, C., Dowsett, H. J., Jost, A., Kamae, Y., Lohmann, G., Lunt, D. J., Abe-Ouchi, A., Pickering, S. J., Ramstein, G., Rosenbloom, N. A., Salzmann, U., Sohl, L., Stepanek, C., Ueda, H., Yan, Q., and Zhang, Z.: Large-scale features of Pliocene climate: results from the Pliocene Model Intercomparison Project, Clim. Past, 9, 191–209,, 2013b. a, b

Haywood, A. M., Dowsett, H. J., Dolan, A. M., Rowley, D., Abe-Ouchi, A., Otto-Bliesner, B., Chandler, M. A., Hunter, S. J., Lunt, D. J., Pound, M., and Salzmann, U.: The Pliocene Model Intercomparison Project (PlioMIP) Phase 2: scientific objectives and experimental design, Clim. Past, 12, 663–675,, 2016. a, b, c, d, e

Haywood, A. M., Tindall, J. C., Dowsett, H. J., Dolan, A. M., Foley, K. M., Hunter, S. J., Hill, D. J., Chan, W.-L., Abe-Ouchi, A., Stepanek, C., Lohmann, G., Chandan, D., Peltier, W. R., Tan, N., Contoux, C., Ramstein, G., Li, X., Zhang, Z., Guo, C., Nisancioglu, K. H., Zhang, Q., Li, Q., Kamae, Y., Chandler, M. A., Sohl, L. E., Otto-Bliesner, B. L., Feng, R., Brady, E. C., von der Heydt, A. S., Baatsen, M. L. J., and Lunt, D. J.: The Pliocene Model Intercomparison Project Phase 2: large-scale climate features and climate sensitivity, Clim. Past, 16, 2095–2123,, 2020 (data available at: a, b, c, d

Heinemann, M., Jungclaus, J. H., and Marotzke, J.: Warm Paleocene/Eocene climate as simulated in ECHAM5/MPI-OM, Clim. Past, 5, 785–802,, 2009. a

Herold, N., Buzan, J., Seton, M., Goldner, A., Green, J. A. M., Müller, R. D., Markwick, P., and Huber, M.: A suite of early Eocene (55 Ma) climate model boundary conditions, Geosci. Model Dev., 7, 2077–2090,, 2014. a

Hill, D. J., Haywood, A. M., Lunt, D. J., Hunter, S. J., Bragg, F. J., Contoux, C., Stepanek, C., Sohl, L., Rosenbloom, N. A., Chan, W.-L., Kamae, Y., Zhang, Z., Abe-Ouchi, A., Chandler, M. A., Jost, A., Lohmann, G., Otto-Bliesner, B. L., Ramstein, G., and Ueda, H.: Evaluating the dominant components of warming in Pliocene climate simulations, Clim. Past, 10, 79–90,, 2014. a, b, c

Huang, B., Thorne, P. W., Banzon, V. F., Boyer, T., Chepurin, G., Lawrimore, J. H., Menne, M. J., Smith, T. M., Vose, R. S., and Zhang, H. M.: Extended reconstructed Sea surface temperature, Version 5 (ERSSTv5): Upgrades, validations, and intercomparisons, J. Clim., 30, 8179–8205,, 2017. a

Hunke, E. C. and Lipscomb, W. H.: The Los Alamos sea ice model, documentation and software., Tech. Rep. LA-CC-06-012, version 4, 2008. a

Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., Lamarque, J. F., Large, W. G., Lawrence, D., Lindsay, K., Lipscomb, W. H., Long, M. C., Mahowald, N., Marsh, D. R., Neale, R. B., Rasch, P., Vavrus, S., Vertenstein, M., Bader, D., Collins, W. D., Hack, J. J., Kiehl, J., and Marshall, S.: The community earth system model: A framework for collaborative research, B. Am. Meteorol. Soc., 94, 1339–1360,, 2013. a

Kiehl, J. T., Shields, C. A., Hack, J. J., and Collins, W. D.: The climate sensitivity of the Community Climate System Model version 3 (CCSM3), J. Clim., 19, 2584–2596,, 2006. a

Knutti, R. and Rugenstein, M. A.: Feedbacks, climate sensitivity and the limits of linear models, Philosophical Transactions of the Royal Society A: Mathematical, Phys. Eng. Sci., 373, 20150146,, 2015. a

Large, W. G., McWilliams, J. C., and Doney, S. C.: Oceanic Vertical Mixing - a Review and a Model with a Nonlocal Boundary-Layer Parameterization, Rev. Geophys., 32, 363–403,, 1994. a

Lawrence, D. M., Oleson, K. W., Flanner, M. G., Thornton, P. E., Swenson, S. C., Lawrence, P. J., Zeng, X., Yang, Z.-L., Levis, S., Sakaguchi, K., Bonan, G. B., and Slater, A. G.: Parameterization improvements and functional and structural advances in Version 4 of the Community Land Model, J. Adv. Model. Earth Syst., 3, M03001,, 2011. a

Li, Z., Luo, Y., Arnold, N., and Tziperman, E.: Reductions in Strong Upwelling-Favorable Wind Events in the Pliocene, Paleoceanogr. Paleoclimatol., 34, 1931–1944,, 2019. a

Loptson, C. A., Lunt, D. J., and Francis, J. E.: Investigating vegetation–climate feedbacks during the early Eocene, Clim. Past, 10, 419–436,, 2014. a

Lunt, D. J., Bragg, F., Chan, W.-L., Hutchinson, D. K., Ladant, J.-B., Morozova, P., Niezgodzki, I., Steinig, S., Zhang, Z., Zhu, J., Abe-Ouchi, A., Anagnostou, E., de Boer, A. M., Coxall, H. K., Donnadieu, Y., Foster, G., Inglis, G. N., Knorr, G., Langebroek, P. M., Lear, C. H., Lohmann, G., Poulsen, C. J., Sepulchre, P., Tierney, J. E., Valdes, P. J., Volodin, E. M., Dunkley Jones, T., Hollis, C. J., Huber, M., and Otto-Bliesner, B. L.: DeepMIP: model intercomparison of early Eocene climatic optimum (EECO) large-scale climate features and comparison with proxy data, Clim. Past, 17, 203–227,, 2021. a, b, c

McClymont, E. L., Ford, H. L., Ho, S. L., Tindall, J. C., Haywood, A. M., Alonso-Garcia, M., Bailey, I., Berke, M. A., Littler, K., Patterson, M. O., Petrick, B., Peterse, F., Ravelo, A. C., Risebrobakken, B., De Schepper, S., Swann, G. E. A., Thirumalai, K., Tierney, J. E., van der Weijst, C., White, S., Abe-Ouchi, A., Baatsen, M. L. J., Brady, E. C., Chan, W.-L., Chandan, D., Feng, R., Guo, C., von der Heydt, A. S., Hunter, S., Li, X., Lohmann, G., Nisancioglu, K. H., Otto-Bliesner, B. L., Peltier, W. R., Stepanek, C., and Zhang, Z.: Lessons from a high-CO2 world: an ocean view from   3 million years ago, Clim. Past, 16, 1599–1615,, 2020. a, b, c, d

Neale, R. B., Richter, J., Park, S., Lauritzen, P. H., Vavrus, S. J., Rasch, P. J., and Zhang, M.: The Mean Climate of the Community Atmosphere Model (CAM4) in Forced SST and Fully Coupled Experiments, J. Clim., 26, 5150–5168,, 2013. a

Oldeman, A. M., Baatsen, M. L. J., von der Heydt, A. S., Dijkstra, H. A., Tindall, J. C., Abe-Ouchi, A., Booth, A. R., Brady, E. C., Chan, W.-L., Chandan, D., Chandler, M. A., Contoux, C., Feng, R., Guo, C., Haywood, A. M., Hunter, S. J., Kamae, Y., Li, Q., Li, X., Lohmann, G., Lunt, D. J., Nisancioglu, K. H., Otto-Bliesner, B. L., Peltier, W. R., Pontes, G. M., Ramstein, G., Sohl, L. E., Stepanek, C., Tan, N., Zhang, Q., Zhang, Z., Wainer, I., and Williams, C. J. R.: Reduced El Niño variability in the mid-Pliocene according to the PlioMIP2 ensemble, Clim. Past, 17, 2427–2450,, 2021. a

Oleson, K. W., Lawrence, D. M., Gordon, B., Flanner, M. G., Kluzek, E., Peter, J., Levis, S., Swenson, S. C., Thornton, E., Dai, A., Decker, M., Dickinson, R., Feddema, J., Heald, C. L., Lamarque, J.-f., Niu, G.-y., Qian, T., Running, S., Sakaguchi, K., Slater, A., Stöckli, R., Wang, A., Yang, L., Zeng, X., and Zeng, X.: Technical Description of version 4.0 of the Community Land Model (CLM), NCAR Tech. Note, NCAR/TN-47, 2010. a

Peng, G., Meier, W. N., Scott, D. J., and Savoie, M. H.: A long-term and reproducible passive microwave sea ice concentration data record for climate studies and monitoring, Earth Syst. Sci. Data, 5, 311–318,, 2013. a, b

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

Smith, R. D. and McWilliams, J. C.: Anisotropic horizontal viscosity for ocean models, Ocean Model., 5, 129–156,, 2003. a

Smith, R. D., Jones, P., Briegleb, B., Bryan, F., Danabasoglu, G., Dennis, J., Dukowicz, J., Eden, C., Fox-Kemper, B., Gent, P., Hecht, M., Jayne, S., Jochum, M., Large, W., Lindsay, K., Maltrud, M., Norton, N., Peacock, S., Vertenstein, M., and Yeager, S.: The Parallel Ocean Program (POP) reference manual: Ocean component of the Community Climate System Model (CCSM), Los Alamos National Laboratory Tech. Rep. LAUR-10-01853, 141, 1–141, (last access: 24 February 2017), 2010. a

Steele, M., Morley, R., and Ermold, W.: PHC: A global ocean hydrography with a high-quality Arctic Ocean, J. Clim., 14, 2079–2087,<2079:PAGOHW>2.0.CO;2, 2001. a

Trenberth, K. E. and Shea, D. J.: Atlantic hurricanes and natural variability in 2005, Geophys. Res. Lett., 33, 1–4,, 2006. a

Weiffenbach, J. E., Baatsen, M. L. J., Dijkstra, H. A., et al.: in preparation, 2022. a

Zhang, Z., Li, X., Guo, C., Otterå, O. H., Nisancioglu, K. H., Tan, N., Contoux, C., Ramstein, G., Feng, R., Otto-Bliesner, B. L., Brady, E., Chandan, D., Peltier, W. R., Baatsen, M. L. J., von der Heydt, A. S., Weiffenbach, J. E., Stepanek, C., Lohmann, G., Zhang, Q., Li, Q., Chandler, M. A., Sohl, L. E., Haywood, A. M., Hunter, S. J., Tindall, J. C., Williams, C., Lunt, D. J., Chan, W.-L., and Abe-Ouchi, A.: Mid-Pliocene Atlantic Meridional Overturning Circulation simulated in PlioMIP2, Clim. Past, 17, 529–543,, 2021. a

Short summary
The Pliocene was a period during which atmospheric CO2 was similar to today (i.e. ~ 400 ppm). We present the results of model simulations carried out within the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2) using the CESM 1.0.5. We find a climate that is much warmer than today, with augmented polar warming, increased precipitation, and strongly reduced sea ice cover. In addition, several leading modes of variability in temperature show an altered behaviour.