Regional and global climate for the mid-Pliocene using the University of Toronto version of CCSM 4 and PlioMIP 2 boundary conditions

The Pliocene Model Intercomparison Project, Phase 2 (PlioMIP2) is an international collaboration to simulate the climate of the mid-Pliocene interglacial, marine isotope stage KM5c (3.205 Mya), using a wide selection of climate models with the objective of understanding the nature of the warming that is known to have occurred during the broader mid-Pliocene warm period. PlioMIP2 builds upon the successes of PlioMIP by shifting focus onto a specific interglacial and by using a revised set of geographic and orbital boundary conditions. In this paper, 5 we present the details of the mid-Pliocene simulations that we have performed with a slightly modified version of the Community Climate System Model version 4 (CCSM4), and the ‘enhanced’ variant of the PlioMIP2 boundary conditions, and discuss the simulated climatology through comparisons to our control simulations, and to proxy reconstructions of the climate of the mid-Pliocene. With the new boundary conditions, the University of Toronto version of the CCSM4 model simulates a mid-Pliocene which is more than twice as warm as that with the boundary 10 conditions used for PlioMIP Phase 1. The warming is more enhanced near the high-latitudes which is where most of the changes to the PlioMIP2 boundary conditions have been made. The elevated warming in the high-latitudes leads to a better match of the simulated climatology to proxy based reconstructions than possible with the previous version of the boundary conditions.


Introduction
The mid-Pliocene warm period, 3.3-3 million years ago, was the most recent time period during which the global temperature was higher than present for an interval of time longer than any of the Pleistocene interglacials.The prevalence of widespread warming during this time has been inferred from proxy-based sea surface temperature (SST) reconstructions from a number of widely distributed deepsea sedimentary cores (Robinson et al., 2008;Lawrence et al., 2009;Dowsett et al., 2010;Fedorov et al., 2012).The PRISM (Pliocene Research, Interpretation and Synoptic Mapping) version 3 (PRISM3; Dowsett et al., 2010) compilation of mid-Pliocene SSTs contains data from 100 sites.In contrast, surface air temperature (SAT) estimates are more difficult to reconstruct owing to the limited availability of land-based proxies that are useful for palaeothermometry.However, some records available from high latitudes in the Northern Hemisphere provide an important perspective on the magnitude of warming in that region.Rybczynski et al. (2013) found mid-Pliocene plant fossils in peat deposits on Ellesmere Island (north of the Arctic Circle in the Queen Elizabeth archipelago), which led them to estimate that the local mean annual SAT during the mid-Pliocene was 18.3 ± 4.1 • C warmer than present, while summer temperature hovered around 14 • C. Brigham-Grette et al. (2013) studied lacustrine records from the arctic lake El'gygytgyn in NE arctic Russia on the basis of which they inferred that the summer temperature there was ∼ 8 • C warmer than today.
The warming during the mid-Pliocene is expected to have resulted in a much higher eustatic sea level (ESL) owing to the melting of large amounts of polar ice over the extended Published by Copernicus Publications on behalf of the European Geosciences Union.D. Chandan and W. R. Peltier: PlioMIP2 and CCSM4 period of time during which warmer conditions prevailed.One of the earliest estimates for the mid-Pliocene ESL was provided by Dowsett and Cronin (1990), who derived a local estimate based on measurements of present-day elevation above sea level for mid-Pliocene marine deposits along the Orangeburg Scarp, near the border between South Carolina and North Carolina.After correcting the observed elevation for the regional long-term tectonic uplift rate, they arrived at a +35±18 m estimate of the mid-Pliocene ESL (the plus sign denotes a sea level higher than present).This value was used to perform the first global reconstruction of the mid-Pliocene climate: PRISM1 (Dowsett et al., 1994).
Based on an analysis of palaeoclimatic data from highlatitude deep-sea sediments, Kennett andHodell (1993, 1995) argued for a mid-Pliocene sea level of no more than +25 m and "probably significantly less than this for most of this interval".Their analysis led to a downward revision of the sea level estimate for PRISM2 (Dowsett et al., 1999) to +25 m.Subsequent versions of the PRISM reconstruction, namely PRISM3 (Dowsett et al., 2010) and PRISM4 (Dowsett et al., 2016), also include a +25 m mid-Pliocene ESL.A number of other estimates for the ESL all fall in the range of +10 to +30 m (Wardlaw and Quinn, 1991;Krantz, 1991;Dwyer and Chandler, 2009;Naish and Wilson, 2009;Rowley et al., 2013).Miller et al. (2012) performed a multi-proxy analysis using backstripped records from Virginia, New Zealand, and the Enewetak Atoll benthic foraminiferal δ 18 O constraints and Mg / Ca − δ 18 O estimates to argue that the mid-Pliocene sea level was 22±10 m higher than present within a margin of 2 standard deviations.More recently, Winnick and Caves (2015) argued for a reexamination of the transfer function used to convert mid-Pliocene benthic δ 18 O measurements into sea level in order to correct for the changing δ 18 O content of Antarctic ice.On the basis of the revised transfer function, they estimate that the mid-Pliocene ESL was only 9-13.5 m higher than present.
Although estimates for the mid-Pliocene ESL show considerable spread, all of them strongly suggest a substantial or total loss of the most vulnerable ice sheets: the Greenland Ice Sheet (GIS) and the West Antarctic Ice Sheet (WAIS).The GIS contains ∼ 7 m of ESL equivalent ice (Alley et al., 2005) and is known to have been greatly diminished under the influence of mid-Pliocene warmth (Lunt et al., 2008;Contoux et al., 2015).The WAIS, which presently contains ∼ 5 m of ESL equivalent ice (Fretwell et al., 2012), is mostly grounded below sea level and is therefore dynamically unstable and vulnerable to a runaway collapse (Weertman, 1974).The combined collapse of these two ice sheets could contribute up to 12 m to the mid-Pliocene ESL.The largest of the present-day ice sheets, the East Antarctic Ice Sheet (EAIS), is considered to be relatively stable owing to its bedrock lying mostly above sea level.However, the EAIS does contain a number of large subglacial basins; the largest are the Wilkes and Aurora basins, which are grounded below sea level and therefore susceptible to collapse similar to the WAIS.Conse-quently, mid-Pliocene ESL estimates above 12 m imply substantial amounts of additional melting in the most vulnerable sectors of the EAIS.
It would not be particularly surprising if higher warming and the accompanying melting of polar ice took place in a world with much higher radiative forcing from the presence of atmospheric greenhouse gases in higher concentrations.However, proxy reconstructions for the mid-Pliocene atmosphere do not show any significant increase in greenhouse gas concentrations compared to the present day.Various reconstructions for the atmospheric CO 2 concentration, such as those based on alkenone data (Pagani et al., 2010;Seki et al., 2010;Badger et al., 2013;Zhang et al., 2013a), δ 13 C measurements in marine organic matter (Raymo et al., 1996), δ 11 B measurements (Seki et al., 2010;Bartoli et al., 2011;Martínez-Botí et al., 2015), foraminiferal B / Ca ratios (Tripati et al., 2009), and stomatal density in fossilized leaves (Kürschner et al., 1996), show that during the mid-Pliocene CO 2 varied between ∼ 200 and ∼ 450 ppmv with most measurements reported in the range of 300-400 ppmv.A value of 400 ppmv is often used as a boundary condition in atmosphere-only and coupled climate models of the mid-Pliocene.Even if CO 2 concentrations during the mid-Pliocene were only equal to those characteristic of the present day, the significant melting of polar ice might still be reasonably expected due to the extended period of time during which this trace gas concentration was maintained.
Very little is known about the concentration of the very potent greenhouse gas methane, although it is expected to be at least as high as the modern concentration, if not higher, mainly due to its release from the thawing Arctic permafrost and the breakdown of methane clathrates, both of which are expected consequences of a warm Arctic.In the absence of proxy estimates for methane concentration in the mid-Pliocene atmosphere, the modern-day concentration value is usually used (for example in the context of the PlioMIP program) in climate models.The uncertainty in the methane concentration is partly compensated for in climate models by choosing a higher concentration of CO 2 .
The inference of warm temperatures and reduced land ice during a time period when the atmosphere was characterized by greenhouse gas concentrations not much different than the present, and significantly lower than projections for the future, is the primary reason that the mid-Pliocene has received considerable interest in recent decades.Since the future of our climate will be "Pliocene-like", it is of vital importance that we understand the magnitude and distribution of the mid-Pliocene warming, the mechanisms for the warming, and the consequences of the warming.
The Pliocene Model Intercomparison Project (PlioMIP; Haywood et al., 2011) was organized with the aim of understanding the magnitude and mechanisms of the degree of warmth that proxies indicate by using a diverse selection of climate models.Phase 1 of the project (henceforth simply PlioMIP) was concerned with simulating the averaged cli-mate of all warm interglacials in a 300 000-year "time slab" during the mid-Pliocene.Results from the eight models participating in PlioMIP found significant differences between the temperature reconstructions produced by the models and between the models and the data.While the models successfully simulated an overall warmer climate with some amplification of warming in high latitudes, they fell far short of capturing the magnitude of warming at high latitudes that is indicated by proxy-based inferences (Dowsett et al., 2013).
The outcomes of PlioMIP have helped in developing the revised strategy being employed in the context of PlioMIP2.The focus is now on a single warm interglacial, namely marine isotope stage KM5c (3.205 Mya).This approach is expected to improve the resolution of reconstructed temperature data and thereby improve confidence in the data.A new set of boundary conditions is available for modeling groups to use: the PRISM4 boundary conditions (Dowsett et al., 2016).In this paper we will employ these new boundary conditions in an attempt to better understand the expected characteristics of the mid-Pliocene climate.

Model description
For the purposes of the analyses to be discussed herein, the Community Climate System Model version 4 (CCSM4; Gent et al., 2011) will be employed to simulate the mid-Pliocene.CCSM4 is a coupled climate model and consists of four major components that interact with each other through a flux coupler without flux adjustment.These components are briefly described in this section.We use a slightly modified version (described in Sects.2.1.1 and 3) of the default configuration of the model and refer to it as the "University of Toronto version of CCSM4" to distinguish it from the publicly released version of CCSM4 that has been run in default configuration.In this paper "CCSM4" refers to the Toronto version when it is mentioned in the context of one of our simulations and to the default version otherwise.CCSM4 has been used in many palaeoclimate studies, and it is one of the models previously employed in PlioMIP (Rosenbloom et al., 2013).

Atmosphere
The atmospheric component in CCSM4 is the Community Atmosphere Model version 4 (CAM4; Neale et al., 2013), which is the sixth-generation atmospheric general circulation model developed at the National Center for Atmospheric Research (NCAR).In this version, the default dynamical core has changed from a spectral core used in CAM3 to a finite volume (FV) dynamical core.The FV core improves tracer transport capabilities over the spectral core and also reduces the surface zonal wind bias that was present in CAM3.In coupled mode, this bias was responsible for an intensified Antarctic Circumpolar Current (ACC) and excessive transport through the Drake Passage (Large and Danabasoglu, 2006).CAM4 also includes changes to the parameterization of deep convection, which significantly improves the Madden-Julian Oscillation (Subramanian et al., 2011), and in coupled CCSM4 mode, the EL Niño-Southern Oscillation (ENSO) variability (Deser et al., 2012).All our simulations are performed on an FV grid with ∼ 1 • of horizontal resolution (192 grid cells in the latitude and 288 grid cells in the longitude) with 26 levels in the vertical.

Ocean
The ocean component in CCSM4 is the Parallel Ocean Program version 2 (POP2; Smith et al., 2010), which was developed at the Los Alamos National Laboratory.It includes significant improvements over version 1, such as the elimination of the cold bias in the equatorial Pacific cold tongue, an increase in the sharpness of the Pacific thermocline, and a reduction in the sea surface temperature (SST) and sea surface salinity (SSS) bias along the path of the North Atlantic Drift.
POP2 also comes with a new parameterization for the overflow of density-driven currents over oceanic ridges, such as the overflow of the Nordic Sea waters over the Denmark Strait and the Faroe Bank Channel, and the overflow of dense waters over the bottom topography in the Ross Sea and Weddell Sea regions of the Antarctic (Danabasoglu et al., 2010).These overflows are thought to be important in the formation of deep bottom waters (Briegleb et al., 2009).The new parameterization improves the penetration depth of the North Atlantic Deep Water (NADW), but introduces new biases that are worse in coupled simulations compared to uncoupled simulations (Danabasoglu et al., 2010).
An important caveat to be aware of when simulating palaeoclimate is that this parameterization is tuned to reproduce, to the greatest extent possible, the observations of the overturning circulation for the present day (such as those obtained by the RAPID array; Cunningham et al., 2007).As such, its impact when included in the simulation of a time period when the ocean bathymetry would have been significantly different remains unclear.Owing to concern regarding the introduction of biases stemming from this parameterization applied under Pliocene conditions when bathymetric depth is significantly modified, we have chosen to disable the parameterization in our Pliocene simulations.In addition, we have also chosen to disable this parameterization for the control simulations so that when the Pliocene simulations are compared to the controls, the differences arising from the choice of the PlioMIP2 boundary conditions can be assessed correctly, unencumbered by any other differences between the control and the Pliocene simulations.
In our simulations, the ocean model runs on the displaced pole grid in which the poles of the grid are positioned over Greenland and Antarctica, and therefore the ocean does not contain any grid singularities.The grid consists of 60 vertical levels with 320 × 384 cells at each level.This amounts to a nominal resolution of ∼ 1 • , although the dimensions of the grid cells are not uniform throughout the oceans.

Land
The land component in CCSM4 is the Community Land Model version 4 (CLM4; Lawrence et al., 2012).One of the chief deficiencies of CLM version 3 was its poor representation of the hydrological cycle.This is why version 4 includes a considerably revised hydrology with changes to the surface runoff scheme, the groundwater scheme, and soil sub-model hydrology, as well as revised parameterizations for canopy integration, canopy interception, and soil evaporation and frozen soil.These changes have led to notable improvements in soil water content, increased transpiration and photosynthesis, and better simulation of the annual cycle of land water storage.Another significant improvement in CLM4 over CLM3 is the addition of a fully prognostic carbon-nitrogen biogeochemical model.This sub-model is able to affect the climate through its control of the seasonal and annual vegetation phenology.The land model runs at the same resolution as the atmosphere model, CAM4, whereas the river transport model, which is included with the land model, runs at a resolution of 0.5 • .

Sea ice
The sea ice component in CCSM4 is based on the Community Ice Code version 4 (CICE4; Hunke and Lipscomb, 2006) developed at the Los Alamos National Laboratory.The main improvement in CICE4 over its predecessor, the Community Sea Ice Model version 5 (CSIM5; Briegleb et al., 2004), is in the shortwave radiative transfer scheme (Briegleb and Light, 2007).It now uses the microscopic optical properties of snow, ice, and external absorbers (such as black carbon and dust) to explicitly calculate the multiple scattering of solar radiation in sea ice using a delta-Eddington approximation to infer macroscopic optical properties, such as albedo and transmission.As a result, CCSM4 has a much better representation of sea ice albedo compared to CCSM3-CSIM5, in which the albedo was parameterized using the bulk properties (thickness and temperature) of sea ice and snow.The new radiative transfer scheme has also allowed for the inclusion of a simple melt pond parameterization based on the surface meltwater flux and inclusion of the effects of aerosol deposition and cycling.Holland et al. (2012) have found that the effects of direct radiative forcing by melt ponds and aerosols on the Arctic sea ice is much higher than on the Antarctic sea ice.The sea ice component operates on the same grid as the ocean component.Recent analyses of the importance of these parameterization schemes in the sea ice module of the NCAR coupled model have been presented in the context of analyses of the "Snowball Earth" phenomenon in the papers of Yang et al. (2012a, b).

Design of the numerical experiments
This section describes the boundary conditions and the initial conditions used in our simulations along with the methodology through which they have been integrated into the CCSM4 model.The simulations are referred to using the nomenclature first employed by Lunt et al. (2012) and subsequently adopted with modifications for PlioMIP2 by Haywood et al. (2016).In this notation, simulations are referred to by an abbreviated form Ex c , where c is the concentration of atmospheric CO 2 in ppmv and x represents boundary conditions that have been changed from the preindustrial (PI) such that x can be absent (for cases in which no boundary conditions have been modified) or it can be "o" for a change in orography and/or "i" for a change in land ice configuration.
Our control and mid-Pliocene simulations were initially integrated for several hundred model years with the diapycnal diffusivity in the ocean, κ, fixed to a constant depth-independent background value of 0.16 cm 2 s −1 ; this is the modern-day pelagic value.Subsequent integrations of the simulations were continued with κ fixed to the depthdependent profile that was used in the ocean component, POP1, of the CCSM3 model (Collins et al., 2006).This profile, which we will henceforth refer to simply as the "POP1 profile", is shown in the Supplement Fig. S1 and is characterized by a depth-dependent hyperbolic tangent function that gives an upper ocean κ of 0.16 cm 2 s −1 and a deep ocean κ of 1.0 cm 2 s −1 with the smooth transition between these asymptotic values centered at a depth of 1000 m; this is the approximate depth of the main thermocline.A simulation performed with the POP1 profile will be denoted by appending a P to the name derived using the nomenclature mentioned earlier and will be referred to as the "POP1 variant".
The motivation behind the choice of these two diapycnal mixing profiles in the context of palaeoclimate simulations, instead of the more complex and spatially varying profile that is part of CCSM4, has been previously articulated by Peltier and Vettoretti (2014) (henceforth, PV14).Essentially, the spatially varying component of the CCSM4 mixing field includes a significant contribution from the turbulent mixing that is generated by the flow of the barotropic tide over rough bottom topography and tuned in the CCSM4 model for the modern-day tidal regime (determined in part by the modern-day bathymetry and coastlines).During the mid-Pliocene there were important changes to the coastlines and bathymetry that would have affected regions that are hot spots of tidal mixing in the present day, such as the Maritime Continent, the eastern Pacific Ocean, the Labrador Sea, and the Denmark Strait and Faroe Bank Channel region in the North Atlantic.In PV14 with a focus on the Dansgaard-Oeschger oscillation phenomenon, the CCSM4 diapycnal turbulent diffusivity scheme was similarly disabled based on the explicit demonstration in the papers of Griffiths andPeltier (2008, 2009) that the barotropic tidal regime of the glacial ocean was dramatically different from the modern.Additionally, the elevations of the mid-ocean ridges could have differed from the present day by as much as 500 m (see Fig. 7 in Dowsett et al., 2016, comparing PRISM4 bathymetry to modern) due to the influence of dynamic topography (Moucha et al., 2008), which would have also had an important impact on the generation of turbulent mixing.It is for these reasons that we have opted to employ simplified models of the diapycnal diffusivity for the mid-Pliocene ocean and also for the oceans in the control simulations to facilitate comparison with the mid-Pliocene results without introducing ambiguity associated with the choice of different mixing schemes.Furthermore, the two different choices of the diapycnal diffusivity will allow us to assess the degree of impact that this parameter has on the strength of the meridional overturning circulation under mid-Pliocene boundary conditions.Given that one of the defining largescale characteristics of the global oceans is the vertical variation in the turbulent diapycnal diffusivity by an order of magnitude, from low values above the thermocline to high values closer to the rough ocean floor (Waterhouse et al., 2014), our POP1 variant simulations will constitute our most accurate PlioMIP2 simulations.It is these simulations that should be used by other groups for intercomparison.

Control experiments
We have simulated two PlioMIP2 control experiments, E 280 and E 400 (and their POP1 variants E 280 P and E 400 P), with atmospheric CO 2 concentrations of 280 and 400 ppmv, respectively.Henceforth, we will refer to the former as the PI control and to the latter as the modern control (which resembles modern only inasmuch as the atmospheric CO 2 is close to modern; all other trace gases are identical to the PI control and no urban or agricultural land units are included in the land model.)Both E 280 and E 400 are started as CCSM "hybrid runs" from an existing 3500-year PI control simulation (called "cesmpifv1mts"), which was run with an atmospheric CO 2 concentration of 280.4 ppmv and initialized in a similar manner to the glacial simulations discussed in Vettoretti and Peltier (2013), and PV14, in which modern-day temperature and salinity were assumed (Levitus and Boyer, 1994) and the ocean and the atmosphere were assumed to be at rest.The model was then run continuously for 3500 years with the overflow parameterization and the tidal mixing parameterization turned off throughout the simulation.
Cesmpifv1mts is unique to this paper but very closely related to the PI control employed for the purpose of the Dansgaard-Oeschger analyses discussed in PV14 and further in Vettoretti andPeltier (2015, 2016).Because the PI control employed in PV14 was created by branching from an NCAR PI control simulation at year 863 (simulation b40.1850.track1.1deg.006run in the default CCSM4 configuration) and run for an additional 1200 years with the overflow parameterization and the tidal mixing parameteri-zation both turned off, we were concerned that this PI control could have inherited memory of branching from the standard configuration of CCSM4.We therefore considered it wise to produce cesmpifv1mts as an entirely new PI control.In retrospect this PI control is found to be very close to that employed in PV14.Time series for the evolution of the globally averaged ocean temperature, SST, salinity, and Arctic sea ice area and volume in cesmpifv1mts are shown in the Supplement Sect. 2. In the last 500 years of this simulation, the global ocean temperature was drifting at a rate of only 0.01 • C century −1 .Another slight difference between cesmpifv1mts and the PI control in PV14 (a technical difference essentially irrelevant to the simulated climate) is that the modern-day topography and bathymetry for cesmpifv1mts was generated in-house from ETOPO1.The PI control in PV14, however, used the default NCAR generated modernday boundary conditions based on ETOPO2v2 because that simulation was branched from an NCAR PI control.
Therefore our two PlioMIP2 control experiments, E 280 and E 400 , have a common ancestor (cesmpifv1mts) and they have identical orography, bathymetry, land ice, and river directions.The vegetation, soil, and wetland and lake distribution in the land model in these two controls is the same as employed in the 1850 configuration in CCSM4 (component set label B_1850_CN).The concentrations of the atmospheric trace gases and orbital parameters are set to those prescribed for PlioMIP2 and listed in Table 1.The Supplement Sect. 3 presents a comparison of the PI SST and sea ice concentration between the POP1 variant of our control (E 280 P) and the HadISST dataset (Rayner et al., 2003).

Pliocene experiments
We report on two PlioMIP2 mid-Pliocene experiments, namely the Core simulation Eoi 400 and the Tier 1 simulation Eoi 450 .In keeping with the notation discussed above, both of these simulations include the modified orography (and bathymetry) and land ice specified by the PLioMIP2 "enhanced" boundary conditions described in Dowsett et al. (2016).The POP1 variants of these simulations have also been performed.Additional PlioMIP2 mid-Pliocene experiments and experiments (not included in the PlioMIP2 experiment list) exploring the sensitivity of the climatology to variations in PlioMIP2 boundary conditions were also performed, but the results of these simulations will not be discussed in the present paper.In order to keep our discussion concise and straightforward, we will focus only on the 400 ppmv Pliocene simulations (Eoi 400 and Eoi 400 P) and refer to results from the 450 ppmv simulations only in selected discussions.The following subsections describe the process through which the PlioMIP2 enhanced boundary conditions have been integrated into the CCSM4 model.

Orography and the atmosphere model
The orography, or above-sea-level topography that includes the topography of ice sheets, appears as a boundary condition primarily in the atmospheric component of the coupled model where this field is ingested by the model in the form of the surface geopotential height field.Another application of orography is in the generation of the global dataset of river direction vectors, which is used by the land model and discussed further below.Following the anomaly method (Haywood et al., 2016), we have first computed the anomaly between the PlioMIP2 mid-Pliocene orography and the PlioMIP2 modern orography and then superimposed this anomaly, shown in Fig. 1, on our local modern orography.This method leads to subtle differences between the resulting land-sea mask (LSM) and the LSM prescribed in PlioMIP2.Therefore in an effort to maintain the integrity of the PlioMIP2 LSM, we conditionally overwrite the orography of any grid cell in which the LSM differs from that described by PlioMIP2.

Bathymetry and the ocean model
The vertical levels in the POP2 ocean model are referred to as "KMT levels" with values in the range of 1-60.Prescribing the mid-Pliocene bathymetry thus requires generating a twodimensional array containing the KMT levels of each ocean grid cell.An initial KMT field is easily obtained by using the POP2 model's ancillary scripts and the PlioMIP2 mid-Pliocene bathymetry.However, this raises a number of issues that need to be addressed.These issues fall into two categories: (i) conflict between the LSM on the displaced dipole ocean grid and the PlioMIP2 LSM, particularly along straits and regions that contain archipelagos or complex coastlines, and (ii) issues pertaining to the requirements of the ocean model, such as the minimum KMT level and the minimum width of narrow straits in terms of the number of grid cells required to permit flow through them.
We have created a suite of graphical tools to enable the targeted editing of boundary condition data products for CCSM4, such as editing the KMT levels.Figure 2 shows a screenshot from one of the tools, KMTEditor, which is seen here zoomed over a region of the North Atlantic.The thick black rectangle is a cursor over a pixel, and the pixels with thin black borders indicate cells with values that have been edited.The ocean model is initialized from a state of rest with modern-day temperature and salinity derived from Levitus and Boyer (1994).No adjustment was made to the global salinity for the +25 m mid-Pliocene sea level because any applicable adjustment would be negligible.In comparison, for the case of the Last Glacial Maximum, at which time the globally averaged sea level was lower than present by approximately 120 m, coupled climate simulations increase the ocean salinity by 1 PSU compared to modern to account for the reduced volume of the oceans (see Vettoretti and Peltier, 2013;Peltier and Vettoretti, 2014).

Land model and river transport model
The land model is capable of dynamically predicting vegetation types as a response to evolving local climate.However, the simulations described here are run with the dynamical vegetation scheme turned off.Instead, we prescribe the mid-Pliocene vegetation reconstruction by Salzmann et al. (2008), which was also used in PlioMIP.This reconstruction is generated using the BIOME4 scheme (Prentice et al., 1992) and is available either as a 28-type "biome" dataset or a 9-type "mega-biome" dataset.Because of the uncertainties associated with vegetation reconstructions so far back in time, we use the mega-biome dataset instead of the biome dataset to minimize error.
The CLM4 model represents vegetation in terms of 17 plant functional types (PFTs) and is able to represent multiple PFTs and land units within a grid cell using a sub-grid hierarchy.In order to use the Salzmann et al. (2008) reconstruction, the mega-biome types have to be mapped into PFTs.We have generated one such mapping by projecting the modernday mega-biome types onto the modern-day PFTs.Such a map is a one-to-many map, but it is readily accommodated in CLM4 because of its ability to represent multiple PFTs per grid cell.The detailed properties of this mapping are described in the Supplement.
The soil color was kept fixed to modern (Lawrence and Chase, 2007).The direction vectors for river routing were generated from the gradient of the topography.These were inspected and edited manually using another graphical tool to ensure that all rivers reach the oceans (or an inland sea) and that there are no river loops that would prevent freshwater from reaching the oceans and consequently leading to salinity drift in the ocean.

Results and discussion
We have organized the results of our numerical experiments into seven individual subsections.In each section, the discussion is based on the climatology simulated for the mid-Pliocene in comparison to that of the control experiments.These climatologies represent the average over 30 model years (averaging years indicated in Table 2), which is the required averaging duration agreed upon at the PlioMIP2 workshop in Leeds in 2016.We mentioned earlier that the POP1 variant simulations will be our primary PlioMIP2 simulations in view of the more physically appropriate form of the diapycnal diffusivity profile employed for these simulations.Results for the constant kappa profile will be discussed only for comparison purposes when we wish to discuss the sensitivity of our results to variation in the diapycnal diffusivity.

Model evolution
The objective of PlioMIP2 is to simulate the equilibrium climate of a typical interglacial in the mid-Pliocene.Therefore, the climate simulations need to be run for a length of time that is sufficiently long to ensure that the various components of the coupled model system have come into equilibrium.In particular, this is in order to ensure that the ocean, which takes much longer than the atmosphere to equilibrate, has also reached an equilibrium state.Figures 3 and 4 show the evolution of the ocean temperature in our simulations, and Fig. 5 shows the evolution of the air temperature in the atmospheric layer closest to the surface.In these figures, a "fork" in a time series denotes the branching of the POP1 variant of that experiment (see Sect. 3) from the constant κ variant.In all cases, the original simulations were also allowed to evolve for a considerable amount of time.The total number of model years for which each of our simulations have been run and the top of the atmosphere (TOA) energy imbalance  over the climatology are listed in Table 2.All simulations have a very low TOA energy imbalance, indicating that the models are very close to statistical equilibrium states.
It is observed that while the SSTs (Fig. 3a) in all models have come into equilibrium (Table 3), the deep ocean is continuing to warm (Fig. 3d).It is also seen that the introduction of the POP1 profile for κ leads to a greater exchange of heat between the upper ocean and the deeper ocean as expected on physical grounds.This increases the rate of warming in the deeper ocean and decreases the rate of warming for the upper and the middle ocean.The ocean as a whole continues to take up heat in all models (Fig. 4), although the trends are quite small (Table 3).The ocean in the PI control shows the slowest rate of warming (0.03 • C century −1 ), even after the introduction of the POP1 profile for κ.This is due to the fact that this control run was initialized from an existing and well-equilibrated (Supplement Fig. S3) 3500year control and integrated for a further 1700 years, which has given the deep ocean sufficient time to come into equilibrium.The ocean in the modern control is warming at a slightly faster rate of 0.06 • C century −1 .The oceans in the Eoi 400 P and Eoi 450 P mid-Pliocene simulations are warming at rates of 0.05 • C century −1 and 0.05 • C century −1 , respectively at the end of over 2500 years of integration.

Surface air temperature
The mean annual surface air temperature (MASAT) anomalies (i.e., anomalies computed using a 2 m air temperature) in the Eoi 400 P mid-Pliocene simulation compared to both of our control simulations are shown in Fig. 6.A polar stereoscopic projection is employed for these graphics because the MASAT anomalies exhibit distinct spatial features only in the high latitudes.Compared to the PI, the mid-Pliocene temperatures at both poles are much higher (Fig. 6a-b).The pattern in the Northern Hemisphere bears resemblance to the polar amplification that is observed in contemporary measurements of surface temperatures.However, the mid-Pliocene differs from PI (and modern) in the land-sea mask, ice sheet configuration, and vegetation.Because of these changes, the amplification near the poles during the mid-Pliocene is expected to include feedback mechanisms that are not active under present-day conditions.The significant warming of the mid-Pliocene SATs over the high-latitude North Atlantic, the Arctic, and the Labrador Sea compared to PI, as shown in Fig. 6a, is also partly due to the significant difference in the sea ice concentration between the mid-Pliocene and the PI control (Sect.4.7).
The MASAT anomaly compared to modern is shown in Fig. 6c-d.Because both the mid-Pliocene and the modern experiments have the same atmospheric CO 2 concentration, this anomaly is a good approximation that excludes the warming signal that would arise from differences in atmospheric CO 2 , such as was the case for the MASAT anomaly compared to PI conditions.Therefore the MASAT anomaly with respect to modern conditions represents the impacts of other changes related to the mid-Pliocene boundary conditions and the feedbacks associated with those changes.Naturally, the amplitude of warming is expected to be reduced compared to what would be obtained in comparison to PI conditions.The zonal averages of both MASAT anomalies (solid lines in Fig. 8) demonstrate that the anomaly compared to modern is about 1 • C lower than the anomaly compared to PI throughout the tropics and the subtropics.However, at high northern latitudes, we begin to see a difference between the two anomalies.In the Southern Hemisphere the difference between the anomalies as compared to both control simulations continues to remain at a level of approximately 1 • C, peaking at ∼ 2 • C around 70 • S. In the Northern Hemisphere the difference between the two anomalies diverges rapidly.
The larger anomaly with respect to the PI appears to be due to the 120 ppmv CO 2 difference between the mid-Pliocene and PI and to the dramatic difference between the Northern Hemisphere LSM.
The sensitivity of the simulated SAT anomalies to the choice of the ocean diapycnal diffusivity is assessed by repeating this analysis with the set of simulations performed with constant diapycnal diffusivity.Figure 8 shows the MASAT anomalies (dashed lines) in the Eoi 400 simulation with respect to the PI control (E 280 in this case) and the mod- ern control (E 400 ) alongside the corresponding anomalies obtained (and discussed above) within the set of simulations performed with the POP1 diffusivity profile.This comparison shows an asymmetric response of the anomalies between the two hemispheres.The anomalies obtained with constant κ simulations (dashed lines) are appreciably reduced in the Southern Hemisphere compared to the anomalies obtained with POP1 type simulations (solid lines), while the opposite is true for the Northern Hemisphere.There are no appreciable differences throughout the tropical and extra-tropical latitudes.We believe this analysis provides a useful first estimate for the community regarding the magnitudes of changes that could be expected and the regions where those changes could occur due to changes in κ.This would assist in better understanding the differences between results from the various models participating in PlioMIP2.Furthermore, the differences in surface air temperature at high latitudes (originating from the choice of κ) are certainly going to be an important consideration for any future study with the goal of simulating the response of high-latitude mid-Pliocene ice sheets using ice sheet coupled models forced by simulated climate as boundary conditions.
Table 4 compares the MASAT simulated in the Eoi 400 P experiment to the PRISM3 interval (3.3-3.0Mya) MASATs, which were used in the original PlioMIP program to com- pare the simulated temperatures to proxy-based inferences of terrestrial temperatures (Salzmann et al., 2013).A caveat that must be kept in mind when comparing our simulated results to this compilation is that the PRISM3 proxy estimate  5).The 400 ppmv mid-Pliocene simulation is 3.8 • C warmer than the PI and 1.8 • C warmer than the modern control, while the 450 ppmv mid-Pliocene simulation is 4.3 • C warmer than PI and 2.3 • C warmer than modern control.For the set of simulations performed with a constant diapycnal diffusivity, the 400 ppmv mid-Pliocene simulation is 3.5 • C warmer than the PI and 1.5 • C warmer than the modern control.Therefore, irrespective of the choice of the ocean diapycnal diffusivity, the new PlioMIP2 boundary conditions simulate a mid-Pliocene with a globally averaged MASAT anomaly with respect to PI control larger than the magnitude of the anomaly predicted by every model exercised previously in the context of the original PlioMIP program (Haywood et al., 2013) and double the anomaly that Rosenbloom et al. (2013) found (1.86 • C) using CCSM4.
Our simulated anomaly is also much higher than that found by Kamae et al. (2016) (2.4 • C) using the PlioMIP2 boundary conditions and the MRI-CGGM2.3coupled climate model.A likely explanation for why their anomaly is lower than ours could be their choice of the "standard" boundary conditions, which do not require changes to the land-sea mask in the model.A more likely explanation of the difference, however, is simply that the relatively short integration length of 500 years would not have been sufficient to enable the ocean in their model to reach a state of quasiequilibrium.
The Earth system sensitivity (ESS) calculated from our mid-Pliocene simulations is provided in Table 5.The ESS inferred on the basis of our 400 ppmv mid-Pliocene experiment is 7.4 • C per doubling of CO 2 .Although this number is fairly high, it agrees well with the 7.1±1.0-9.7±1.3 range of estimates that Pagani et al. ( 2010) inferred using proxy-based methods.Our estimate for the ESS is double that obtained using CCSM4 for PlioMIP of 3.51 • C and significantly higher than the PlioMIP multi-model mean of 5.01 • C.
Both of our mid-Pliocene simulations have the greatest warming occurring during the months of JJA and the least warming during the months of DJF (Table 6).The warming during JJA is more than 1 • C higher than during DJF for both mid-Pliocene models and compared to both of our controls.In fact, while the JJA-DJF temperature difference in our control simulations is ∼ 3.5 • C, the difference is ∼ 4.7 • C in the mid-Pliocene simulations.This represents an increase in seasonality in the mid-Pliocene over that under either PI or modern conditions.
In Fig. 7 we show the Northern Hemisphere SAT anomalies in Eoi 400 P for JJA and DJF, which captures this seasonality.The first row shows the anomalies relative to the PI and the second row shows the anomalies relative to modern.The increase in seasonality in the mid-Pliocene is readily apparent, especially in comparison to modern conditions such that the Pliocene winter (summer) is colder (warmer) than modern over land.The significant and widespread cooling of the Northern Hemisphere in the mid-Pliocene winter compared to modern is not unlike the temperature trend that has been observed in recent decades (Sun et al., 2016;Cohen et al., 2013;Overland et al., 2011).The trends for recent decades show that while the Arctic has been warming during the winter, large portions of Eurasia and North America have experienced a cooling trend.A large body of work has been produced in an attempt to understand this (see a recent review www.clim-past.net/13/919/2017/Clim.Past, 13, 919-942, 2017 by Vihma, 2014).One of the most promising explanations for the observed cooling over land has come in the form of the dynamical connection that has been proposed between the Arctic surface warming, as a result of the September-October-November (SON) Arctic sea ice loss, and the corresponding equatorward shift in the jet stream, which leads to a greater intrusion of cold arctic air over the midlatitude landmasses (Overland et al., 2016;Barnes and Screen, 2015;Deser et al., 2015;Cohen et al., 2014).An equatorward shift in the jet stream has been shown to occur in a dry dynamical core model from an imposed surface warming in the Arctic (mimicking the warming from sea ice loss) by Butler et al. (2010, also see references therein).More recently, Deser et al. (2015) have applied the CCSM4 model to the analysis of the atmospheric response to Arctic surface warming from sea ice loss by imposing an artificial longwave forcing on the sea ice component of the model to simulate sea ice loss.They found that in response to the SON Arctic sea ice loss, the CCSM4 model simulates an equatorward shift in the jet stream.
Similarly, we find that the largest sea ice loss in the 400 ppmv mid-Pliocene simulation compared to the modern control occurs in SON (Supplement Fig. S9).The reduced sea ice coverage in our mid-Pliocene simulation leads to an increased heat flux into the atmosphere during the winter season (when the sea-air temperature contrast is largest), which leads to the warming of the atmospheric column above the Arctic more than over the surrounding midlatitudes.This in turn leads to a larger increase in geopotential height over the Arctic compared to the midlatitudes and drives the jet stream equatorward in accordance with the thermal wind relation, which brings the polar front closer to the midlatitudes resulting in the cooling that is seen in Fig. 7c. Figure 9 shows the boreal winter zonal-mean zonal wind anomaly between the 400 ppmv mid-Pliocene simulation and the modern con- trol, showing this equatorward shift in the zonal wind and an equatorward shift in the upper stratospheric winds.

Precipitation
The anomaly in the mid-Pliocene annual precipitation compared to our two control simulations is shown in Fig. 10.We see the presence of a strong double ITCZ, which is a recognized problem with CAM4 (Gent et al., 2011).The anomaly is larger with respect to the PI than it is with respect to modern conditions.This is expected as the Pliocene atmosphere is much warmer than in the PI owing to greater atmospheric CO 2 , whereas the atmosphere in the modern con- trol has the same CO 2 concentration as the Pliocene.Despite the differences in the magnitude of the anomalies, the broad features remain the same.The precipitation increases over mountain belts, such as the Andes, the Himalayas, and the Tibetan Plateau, due to the increased mid-Pliocene orography of these regions (Fig. 1) and the higher moisture content of the air rising above the mountain belts.

Ocean temperature
Figure 11 shows the 400 ppmv mid-Pliocene SST and the SST anomalies compared to the PI control and the modern control.Our simulation is characterized by a fairly extensive expansion of the warm pool in the mid-Pliocene, the existence of which has been previously suggested on the basis of proxy-based reconstructions of SSTs (Brierley et al., 2009;Dowsett et al., 2012;Fedorov et al., 2012).The blue and red contour lines in Fig. 11a show the extent of 30 • C and 31 • C waters of the equatorial warm pool.Such warm waters are not present in the PI control simulation in which the temperature throughout the equatorial warm pool is 1.5-2 • C lower than the mid-Pliocene (Fig. 11b).Our decision to keep the ocean configuration identical between the mid-Pliocene and the control simulations, as well as our decision to perform experiments with two choices of κ, has enabled us to determine (Supplement Sect. 4)   1931, 1960 , and PI control  E 280 P 5131, 5160 for two oceanic basins.On the left are the zonal means over the Atlantic and Arctic basins, and the zonal mean over the Southern, Indian, and Pacific oceans is shown on the right.The data points are the PRISM3 estimates for SSTs (Dowsett et al., 2010), which have been categorized into three confidence categories (Dowsett et al., 2012).The shaded region highlights the range of the simulated mid-Pliocene temperature in each basin.The denotes the range of model years over which the climatology was computed.
pool is solely due to the choice of PlioMIP2 boundary conditions.Warm waters are also present at high latitudes in the North Atlantic and in the Southern Ocean, where the SST anomalies show amplified warming compared to the rest of the ocean.The related shift in the east-west temperature gradient across the equatorial Pacific is expected to have a (perhaps significant) impact on the ENSO process, an impact that will be discussed in detail elsewhere.
In comparing the mid-Pliocene SST to the modern control SST (Fig. 11c), a large region of negative anomalies is seen off the west coast of North America, implying that the ocean surface in the modern control is warmer than in the mid-Pliocene.This is because of an increase in the mid-Pliocene surface wind stress compared to the present day.In this region, the wind stress is responsible for forcing the prominent present-day California Current, which drives the coastal upwelling that brings colder waters from below, thereby reducing the ocean surface temperature.The increased wind stress under mid-Pliocene conditions is responsible for greater coastal upwelling, the reduction of SSTs (compared to modern), and the formation of the "cold tongue" that extends from the coast.
Another way in which the mid-Pliocene ocean can be compared to our two controls is through the meridional profile of zonal mean SST, which is shown in Fig. 12 for the Atlantic and Arctic basins and the Southern, Indian, and Pacific basins.The mid-Pliocene ocean is warmer than the PI by at least 2.5-3.5 • C and warmer than the modern by at least 1-2 • C. The largest anomalies occur in the 45-65 • N region of the Atlantic where most of the NADW forms.The data points in Fig. 12 are estimates of the mid-Pliocene SST from the PRISM3 dataset categorized into three confidence levels (Dowsett et al., 2012): very high, high, and medium.
Before we discuss the comparison to PRISM3 it needs to be noted that the PRISM3 reconstruction was generated for the original PlioMIP program in which the aim was to simulate the average climate of the warm intervals in a 300 000-year time slab from ∼ 3.3 to 3 millions years ago; therefore, the PRISM3 boundary conditions and the SST reconstructions are an average over that time period.The PRISM3 dataset is therefore not strictly applicable to PlioMIP2 in which the focus is on a time slice centered on the single interglacial peak at MIS KM5c.However, since the revised dataset that will eventually be applicable to PlioMIP2 is not yet available, we will compare our results to PRISM3.
We find that our simulated meridional SST profile is in rather good agreement with the limited number of data points that are available.The shaded region shows the range of the simulated mid-Pliocene temperatures over the specific ocean basin.Disagreement is seen near the equator in the Indo-Pacific basin where the SST estimates are ∼ 2 • C lower than our simulated SSTs.In this region, proxy estimates in fact match better to the modern control than to the mid-Pliocene.However, this is in precise agreement with what O'Brien et al. ( 2014) have recently argued regarding the inability of Mg / Ca and alkenone-based proxies (which have been used to reconstruct some of the equatorial region SSTs in the PRISM3 dataset) to capture the significantly warmer temperatures in this region.Specifically, they argue that the insensitivity of alkenone proxies to temperatures > 29 • C and the dependence of the Mg / Ca calibration on seawater Mg / Ca ratios cause both proxies to fail when it comes to recording the warmer water temperatures during the mid-Pliocene.This has led to speculations concerning the possible operation of "thermostat-like mechanisms" that might have limited the warm pool temperatures during the Pliocene.Using the TEX 86 proxy and a revised calibration of the Mg / Ca proxy, O'Brien et al. (2014) argue that the equatorial warm pool temperatures were about 2 • C warmer than the present day, a suggestion with which our simulated results agree.Another region of potential disagreement between the model prediction and proxy reconstruction is the high-latitude North Atlantic.
We have also compared the simulated mid-Pliocene SST anomalies with respect to PI to the large compilation of SST anomaly estimates from the PRISM3 program.Figure 13 shows the PRISM3 estimate for the mid-Pliocene SST anomalies (compared to PI) categorized into the three confidence levels mentioned above.There are roughly 100 sites in the PRISM3 dataset distributed over all ocean basins that have been arranged as a function of latitude in the figure .The data indicate that the mid-Pliocene ocean was on average a few degrees warmer than the present day and that the warming was particularly pronounced in the high latitudes of the Northern Hemisphere.The SST anomalies at PRISM3 sites obtained from our mid-Pliocene simulations Eoi 400 P and Eoi 450 P are shown in Fig. 13b as blue and green dots, respectively.Both simulations are able to capture the highlatitude warming in the Northern Hemisphere, except at the locations of the most northern data points, the reliability of which has recently been called into question in the Pliocene community.Our mid-Pliocene ocean is warmer than the data suggest in the Southern Ocean, but the differences are not extreme.
Also shown in Fig. 13a are the mid-Pliocene SST anomalies obtained with CCSM4 (black dots; Rosenbloom et al., 2013) using the PRISM3 boundary conditions (Sohl et al., 2009) for the original PlioMIP program and the multimodel mean (MMM) from PlioMIP (black squares).With the PRISM3 boundary conditions, the original CCSM4 simulation had difficulty simulating the Northern Hemisphere high-latitude warming.The response of the model is very flat across all latitudes.The MMM from PlioMIP similarly did not suggest any enhanced warming in the high latitudes.However, with the PRISM4 boundary conditions the Toronto version of the CCSM4 model simulates the amplification near the poles very well when the model is run to near statistical equilibrium (Fig. 13b).Recently, using sensitivity studies conducted with the CCSM4 model, Otto-Bliesner et al. ( 2017) have shown that the closure of the Bering Strait and the Canadian Arctic Archipelago, both of which are the major differences between the revised PRISM4 and the older PRISM3 paleoenvironmental reconstructions in the Northern Hemisphere, leads to greater warming over the North Atlantic.The authors have suggested that this is due to the inhibition of the transport of fresher waters from the Pacific to the North Atlantic, leading to a stronger Atlantic Meridional Overturning Circulation and larger northward oceanic heat transport.

Meridional overturning circulation
In Fig. 14 we show the evolution of the Atlantic Meridional Overturning Circulation (AMOC) maximum in our simulations.The AMOC in all our simulations appears to have reached equilibrium.The mean state of the AMOC maximum for each simulation over the 30-year climatology is listed in Table 2.It is seen that across all simulations the AMOC is ∼ 2 Sv stronger in the variants with the POP1 diffusivity profile compared to the variants that employ only the background pelagic value of the diapycnal diffusivity.The dependence of the AMOC on the nature of vertical mixing is expected, and therefore it is something to keep in mind when comparing AMOCs from different models with their own vertical mixing schemes (Zhang et al., 2013b).Additionally, the AMOC shows more variability in simulations with the POP1 mixing profile.
The strengths of the AMOC in the PI control are 20 and 21.5 Sv for simulations E 280 and E 280 P, respectively, and in the modern control simulations E 400 and E 400 P the AMOC strengths are 21.9 and 24.2 Sv.Our climatological estimates of the PI AMOC are in satisfactory agreement with the 17.2 Sv estimated over a short time span (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) by Mc-Carthy et al. (2015) using measurements obtained from the RAPID monitoring array.
Several different proxies, on the basis of which it has been argued that one may infer the strength of NADW or Antarctic bottom water (AABW) formation, have been invoked to argue for a more vigorous mid-Pliocene AMOC, sometimes called the "super conveyor", compared to the present day.These include arguments from comparisons of mid-Pliocene benthic δ 13 C values in ocean basins to modern-day values (Billups et al., 1997;Ravelo and Anderson, 2000;Raymo et al., 1996), measurements of Nd and Pb isotope composition recorded in ferromanganese crusts and nodules (Frank et al., 2002), oceanic carbonate dissolution history (Frenz et al., 2006), and reconstructions of past marine ice sheet extents in the Ross Ice Shelf regions (McKay et al., 2012).However, similar to our own findings that we will report below, coupled climate models have not been able to reproduce such invigorated AMOC in the mid-Pliocene (Haywood and Valdes, 2004;Zhang et al., 2013b).For the end of the century, which like the mid-Pliocene would correspond to a warm climate, multi-model projections show a very high likelihood of a reduction in the strength of the NADW cell (Collins et al., 2013).
The inability of models in the PlioMIP project to simulate an energized AMOC (Zhang et al., 2013b) coincides with the additional failure of these models to simulate the Northern Hemisphere SST amplification that proxy records suggest (Fig. 13a).It has therefore been previously argued (Haywood et al., 2016;Hu et al., 2015) that the closure of the Bering Strait (in PlioMIP2) could lead to larger AMOC and consequently greater oceanic heat transport to higher latitudes in the North Atlantic.If this were to happen, it would not only reconcile the model predictions with SST proxies, but also with the proxies that suggest an intensified AMOC.Recently, Hu et al. (2015) reported on the effects of the Bering Strait closure on AMOC strength and the meridional heat transport using CCSM3 and CCSM2 under boundary conditions for the present, 15 thousand years before present, and 112 thousand years before present.They found that under all these conditions, i.e., regardless of the background climate state and for both models, the closure of the Bering Strait resulted in a strengthening of the AMOC by ∼ 2-3 Sv.
We find that the AMOCs in our 400 and 450 ppmv mid-Pliocene simulations are almost identical and their strengths are only 2-3 Sv higher than in the PI control (Fig. 14).This represents an increase in the mid-Pliocene AMOC by just 10 % over the PI AMOC and therefore does not lend support to the idea of a significantly intensified AMOC dur-ing the mid-Pliocene.This increase is comparable to the increase that has been estimated by Hu et al. (2015) to result from the closure of the Bering Strait and would therefore lead to speculation that the stronger AMOC in our mid-Pliocene simulations is the result of the Bering Strait closure.However, this argument is complicated by the fact that the AMOC in the 400 ppmv modern control, which like the PI control is characterized by an open Bering Strait, is stronger than that in the mid-Pliocene simulations.Additionally, our 400 ppmv Pliocene AMOC is weaker than that simulated with the CCSM4 model in PlioMIP (which had an open Bering Strait; Zhang et al., 2013b), although this could be due to the much shorter model run (550 years) on the basis of which the PlioMIP CCSM4 diagnostic was computed.Therefore, it is presently not possible to conclude that the marginal strengthening of the AMOC seen in our mid-Pliocene simulations is due to the closing of the Bering Strait.

Meridional heat transport
The atmosphere and the ocean are together responsible for transporting the excess heat that accumulates near the equator to the high latitudes, where this heat can be radiated to space in the form of longwave radiation, thereby helping to maintain an equilibrium climate.Under present-day conditions the maximum transport of heat poleward is ∼ 5.5 PW that peaks at 30-40 • latitude in each hemisphere (Trenberth and Carson, 2001).The atmosphere dominates the heat transport poleward of the subtropics, and the peak transport of heat is ∼ 5 PW at 40 • latitude.By comparison, the ocean carries much less heat and dominates only in the deep tropics.The maximum heat transported by the ocean is just under 2 PW.These characteristics of the present-day meridional heat transport are well represented in our control simulations (Fig. 15).
Figure 15 shows that the total meridional heat transport in our Eoi 400 P mid-Pliocene simulation is lower than both the PI and modern controls.The reduction in the transport of heat is seen in both the atmosphere and the ocean.The atmospheric heat transport (AHT) in all our simulations (mid-Pliocene and controls) is essentially identical throughout the tropics, the subtropics, and the Southern Hemisphere.The only notable difference arises in the mid-to-high latitudes of the Northern Hemisphere where the mid-Pliocene AHT is lower than both controls (Fig. 15a, c).A very small difference in AHT between the mid-Pliocene and the control is also noticed close to 65 • S, which might be due to the substantial differences in topography and grounded ice sheets between the two simulations at this latitude.
The reduction in the oceanic heat transport (OHT) in the mid-Pliocene is primarily due to the reduction in OHT in the Indo-Pacific basin (Fig. 15b, d).It is likely that this reduction can be attributed to the closing of the Bering Strait.Hu et al. (2015) have shown that closing the Bering Strait can lead to a decrease in the northward heat transport in the North Pa-  Although a comparison between our mid-Pliocene simulations and our control simulations is not the same as comparing two simulations with identical boundary conditions save for the differences concerning the Bering Strait, as Hu et al. (2015) have done, our analysis suggests that their results regarding the impact of the Bering Strait closure on the North Pacific OHT is a robust response by the climate system that persists even when there are other differences in boundary conditions.
We find that the southward OHT in the South Pacific during the mid-Pliocene is also reduced compared to the PI.This, however, is in contrast to the consistent increase in the southward OHT in the South Pacific that Hu et al. (2015) found as a consequence of the Bering Strait closure.This suggests to us that the difference in the atmospheric radiative forcing and the geographical changes between the mid-Pliocene and the PI control have had an impact on the southward OHT in the South Pacific in addition to what would be expected from the analysis of Hu et al. (2015).However, we have to be cautious concerning this comparison, as the simulated differences could also be due to the differences between the ocean components in these versions of CCSM.Both CCSM2 and CCSM3 use POP1 as their ocean model, whereas CCSM4 uses the POP2 ocean model.
The OHT in the Atlantic basin during the mid-Pliocene is only marginally higher, by ∼ 3 %, than the PI.This in- crease is substantially less than the ∼ 10 % increase in the strength of the AMOC.The mid-Pliocene Atlantic OHT is also nearly identical to the modern control despite geographical differences between the two cases.In fact, in the latitude range of 50-70 • N, the Atlantic in the modern control transports more heat than the Atlantic in the mid-Pliocene.This further increases the difficulty in assessing whatever impact the closed Bering Strait might have on the mid-Pliocene climate through the reorganization of oceanic heat transport, and points to the possibility of increasing meridional heat transport in the Atlantic sector without the need for closing the Bering Strait.
Finally, we note the presence of what appears to be a local Bjerknes compensation around 65 • S latitude in our 400 ppmv mid-Pliocene simulation (Fig. 15a, c).In the vicinity of this latitude, the mid-Pliocene AHT is reduced compared to both the controls.It is then left up to the ocean to transport the excess heat southward of this latitude, and consequently the OHT becomes greater than in the controls.This compensation ensures that there is no anomalous change in the net meridional heat transport near this latitude.

Sea ice
The mid-Pliocene sea ice is considerably reduced during both boreal and austral winters compared to PI (Fig. 16).In the Northern Hemisphere the greatest loss in sea ice occurs (besides Hudson Bay, which is expected not to have been present during the mid-Pliocene) in the Labrador Sea and the Green-land and Norwegian seas in the Atlantic sector and in the Barents and Kara seas in the Arctic sector.The reduction in Northern Hemisphere sea ice is particularly pronounced in the summer months to such an extent that the Arctic can be considered ice free.
Sea ice is uniformly reduced along the coastlines of Antarctica during the austral winter.The largest reduction in sea ice is seen in the Weddell Sea and off the coast of Queen Maud Land and Wilkes Land.The sea ice is most concentrated in the vicinity of the mid-Pliocene Antarctic archipelago and in the region presently occupied by the Filchner-Ronne ice shelf.The presence of the archipelago has allowed for the ice that today forms in the Bellingshausen and Amundsen seas to move closer to the pole and therefore leads to a poleward retreat of the ice margin.In austral summer the only concentration of sea ice is in the archipelago, while the rest of the coastline of Antarctica is largely ice free.

Conclusions
In this paper we have described the implementation of the revised boundary conditions for the mid-Pliocene epoch in the CCMS4 coupled climate model and employed the new structure in the reconstruction of mid-Pliocene climate conditions.We have performed two mid-Pliocene experiments, the core experiment denoted Eoi 400 and the Tier 1 experiment denoted Eoi 450 , along with the core control experiments E 280 and the Tier 2 control experiment E 400 .In ad-dition, we have two versions of these simulations that are differentiated by the ocean's vertical mixing profile.The first version has a constant (pelagic) value of diapycnal diffusivity throughout the ocean.The second version has a mixing profile fixed to that used by the ocean component POP1 of the CCSM3 model.The discussions and analysis in this paper are based on the climatology that is simulated by the second version of the experiments.
We find that the PRISM4 boundary conditions mandated in PlioMIP2 lead to greater warming in the mid-Pliocene, in particular to enhanced warming at high latitudes compared to that inferred using the same CCSM4 model and the PRISM3 boundary conditions from the original PlioMIP.The simulated 400 ppmv mid-Pliocene climate has a global MASAT that is 3.8 • C higher than the PI control and 1.8 • C higher than the modern-day control.These anomalies are larger than the anomalies predicted by every model previously exercised in PlioMIP and more than double the anomaly obtained with CCSM4 in this context (Haywood et al., 2013), demonstrating that the changes to the boundary conditions have had a considerable impact on the climate.We expect that it is also important to run such integrations to statistical equilibrium.In addition, we find that the globally averaged temperature difference between the seasons JJA and DJF has increased during the mid-Pliocene compared to both the PI and the modern controls.While the JJA-DJF temperature difference in both of our controls is ∼ 3.5 • C, the difference increases to ∼ 4.7 • C in the mid-Pliocene.
The mid-Pliocene ocean that we have simulated is characterized by (i) a fairly expansive warm pool where the temperatures are 1.5-2 • C warmer than the PI and (ii) elevated levels of warming at high latitudes in the Southern Ocean and the North Atlantic.The SST anomalies with respect to PI agree rather well with the proxy-inferred SST anomalies compiled for the PRISM3 reconstruction.Both the 400 and 450 ppmv mid-Pliocene simulations are able to capture the mid-to-high latitude warming that is seen in the PRISM3 dataset.The agreement between the results of our simulations with the new boundary conditions and the PRISM3 dataset is much better than what was possible with any of the models in the original PlioMIP (Dowsett et al., 2013).The caveat to our present result is that the PRISM3 SST reconstruction was generated for the original PlioMIP program in which the aim was to simulate the average climate of the warm intervals in a 300 000-year time slab from ∼ 3.3 to 3 million years ago; therefore, the PRISM3 boundary conditions and the SST reconstructions are an average over that time period.The PRISM3 dataset is therefore not strictly applicable to PlioMIP2 in which the focus is on a time slice centered on the single interglacial peak at MIS KM5c.However, the revised dataset that will be applicable to PlioMIP2 is not yet available.When the dataset has been made available, we will revisit this data-model comparison.
The larger warming that is seen in both SAT and SST in our simulations at high latitudes is accompanied by a de-crease in the net meridional heat transport compared to that in either of our controls.The reduced meridional heat transport is the result of a reduction in both the atmospheric and oceanic transports.Partitioning of the OHT into contributions from the major oceanic basins shows that the northward transport of heat is greatly reduced in the Indo-Pacific basin.In the Atlantic basin, the meridional transport of heat is only marginally increased (∼ 3 %) compared to the PI control.Compared to the modern control, it is either identical or, in the high northern latitudes, lower than the control.This suggests that the amplified warming at the high latitudes in the mid-Pliocene, inferred from proxies and supported by our simulations, could have more to do with the local positive feedback processes activated by the changes in geography, ice sheets, and vegetation than with the increased northward transport of heat.
Lastly, we note that our simulations do not support the case for a mid-Pliocene AMOC that is substantially stronger than in the PI control.The existence of a stronger AMOC has been argued from various proxy-based inferences (Billups et al., 1997;Ravelo and Anderson, 2000;Raymo et al., 1996;Frank et al., 2002;Frenz et al., 2006;McKay et al., 2012) and has been considered as a possible remedy to the inability of previous climate models to simulate an amplified highlatitude mid-Pliocene warming.Our mid-Pliocene AMOC is only ∼ 10 % stronger than the PI AMOC and slightly weaker than the AMOC in the modern control (Table 2; Fig. 14).It has been argued that the closed Bering Strait in the latest boundary conditions should lead to a stronger AMOC (Haywood et al., 2016;Hu et al., 2015).Although our mid-Pliocene AMOC is indeed stronger than the PI AMOC, the fact that it is also weaker than the modern control AMOC makes it difficult to support the idea that a closed Bering Strait would necessarily lead to a stronger AMOC.
Data availability.All the PlioMIP2/PRISM4 boundary condition data are available from the USGS PlioMIP2 web page: http: //geology.er.usgs.gov/egpsc/prism/7.2_pliomip2_data.html.The results of our mid-Pliocene and control simulations will be released in due time to the data repository for the PlioMIP2 program that has been set up at the University of Leeds (Haywood et al., 2016).
Competing interests.The authors declare that they have no conflict of interest.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.
www.clim-past.net/13/919/2017/Clim.Past, 13, 919-942, 2017 Acknowledgements.Computations were performed on the TCS supercomputer at the SciNet HPC Consortium.SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, the Ontario Research Fund -Research Excellence, and the University of Toronto.We are grateful to Guido Vettoretti for guiding DC through the intricate process of implementing boundary conditions for palaeoclimate simulations with CESM.We are also grateful to NCAR for organizing the annual CESM tutorials, one of which was attended by DC with partial support from NCAR.We would also like to thank Bette Otto-Bliesner for funding a short visit for DC to NCAR, which allowed him to understand more about the CESM model and helped him to implement the mid-Pliocene boundary conditions.DC is also very grateful to the Centre for Global Change Science at the University of Toronto, which funded multiple trips to conferences and workshops related to the work described in the paper.We thank two anonymous reviewers and the editor for their constructive comments, which have significantly improved our paper.The research of WRP at the University of Toronto is funded by NSERC Discovery Grant A9627.
Edited by: Alan Haywood Reviewed by: two anonymous referees

Figure 1 .
Figure 1.The topography anomaly field, which when superimposed on our local modern-day topography, results in the PRISM4 mid-Pliocene topography.

Figure 2 .
Figure 2. A screenshot of one of our graphical tools for editing the ocean model grid.The dark black square on the left side plot is the cursor, which can be used to navigate to individual cells in order to change their values.The cells outlined by a thin dark border denote cells with values that have been edited.

Figure 3 .
Figure 3. Evolution of the ocean temperature at different depths.When a new curve "forks" from an existing curve, it represents the branching of the P version of that simulation in which the ocean model κ has been fixed to the POP1 type profile.In all cases, the original simulation (with κ = 0.16 throughout the ocean) was also continued further in time.(a) The evolution of the sea surface temperature and the volume-averaged temperatures in the upper ocean at 0-550 m of depth (b), the middle ocean at 550-1850 m of depth (c), and the lower ocean below 1850 m (d).

Figure 4 .
Figure 4. Similar to Fig. 3 but showing the evolution of the globally averaged ocean temperature.

Figure 5 .
Figure5.Evolution of the globally integrated air temperature of the atmospheric model layer closest to the surface.When a new curve "forks" from an existing curve, it represents the branching of the P version of that simulation in which the ocean model κ has been fixed to the POP1 type profile.In all cases, the original simulation (with κ = 0.16 throughout the ocean) was also continued further in time.

Figure 6 .
Figure 6.Climatological surface air temperature anomalies for the Eoi 400 P mid-Pliocene simulation 2691, 2720 compared to the preindustrial control E 280 P 5131, 5160 (a, b) and compared to the modern control E 400 P1931, 1960 (c, d).The denotes the range of model years over which the climatology was computed.Anomaly is defined as Pliocene minus control.

Figure 7 .Figure 8 .
Figure 7. Similar to Fig. 6 but showing the mean seasonal surface air temperature anomalies: (a) and (c) are for DJF, (b) and (d) are for JJA.

Figure 9 .
Figure 9.The boreal winter (DJF) climatological zonal-mean zonal wind anomaly between the 400 ppmv mid-Pliocene simulation (Eoi 400 P) 2691, 2720 and the modern control (E 400 P) 1931, 1960 , showing the prominent equatorward shift in the zonal winds over the midlatitudes.The denotes the range of model years over which the climatology was computed.Anomaly is defined as Pliocene minus control.

Figure 10 .
Figure 10.The annual precipitation anomaly between the Eoi 400 P mid-Pliocene simulation 2691, 2720 and (a) the PI control E 280 P 5131, 5160 and (b) the modern control E 400 P 1931, 1960 .Anomaly is defined as Pliocene minus control.The denotes the range of model years over which the climatology was computed.

Figure 11 .
Figure 11.(a) The mid-Pliocene SST from the Eoi 400 P simulation 2691, 2720 .The blue and the red contours are for 30 • C and 31 • C, respectively.The SST anomaly with (b) PI control E 280 P 5131, 5160 and (c) modern control E 400 P 1931, 1960 .Anomaly is defined as Pliocene minus control.The denotes the range of model years over which the climatology was computed.

Figure 12 .
Figure12.The zonal mean SSTs for the 400 ppmv mid-Pliocene Eoi 400 P 2691, 2720 , modern control E400 P 1931, 1960 , and PI control  E 280 P 5131, 5160 for two oceanic basins.On the left are the zonal means over the Atlantic and Arctic basins, and the zonal mean over the Southern, Indian, and Pacific oceans is shown on the right.The data points are the PRISM3 estimates for SSTs(Dowsett et al., 2010), which have been categorized into three confidence categories(Dowsett et al., 2012).The shaded region highlights the range of the simulated mid-Pliocene temperature in each basin.The denotes the range of model years over which the climatology was computed.

Figure 13 .
Figure 13.Data-model comparison of the mid-Pliocene SST anomaly compared to preindustrial.The data in both (a) and (b) are the PRISM3 estimates (Dowsett et al., 2010), which have been categorized into three confidence categories (Dowsett et al., 2012).(a) Comparison of the PRISM3 estimates to the multi-model mean from PlioMIP1 (transparent black square markers) and to the anomalies obtained by Rosenbloom et al. (2013) for PlioMIP1 (black dots) using the same CCSM4 model that we use in this study.(b) Comparison of the PRISM3 estimates to the anomalies obtained from our two mid-Pliocene simulations, Eoi 400 P 2691, 2720 and Eoi 450 P 2621, 2650 , shown as blue and green dots, respectively.The denotes the range of model years over which the climatology was computed.

Figure 14 .
Figure 14.Evolution of the maximum strength of the Atlantic Meridional Overturning Circulation.The results shown have been filteredwith a 20-year running mean to remove high-frequency variability.For each simulation, the curve that starts at model year zero is for the version that uses a vertically constant background diapycnal diffusivity, while the curve with identical color that starts midway is for the version that uses the CCSM3/POP1 type vertical profile of diapycnal diffusivity.

Figure 15 .
Figure 15.Total meridional heat transport and its decomposition.The top (bottom) row compares the Eoi 400 P mid-Pliocene simulation 2691, 2720 to the PI E 280 P 5131, 5160 (modern E 400 P 1931, 1960 ) control.The left column shows the total meridional heat transport and its decomposition into the atmospheric and oceanic components.The right column decomposes the oceanic components into the transport in the Atlantic and the Indo-Pacific basins.The denotes the range of model years over which the climatology was computed.

Figure 16 .
Figure 16.Surface area covered by sea ice for the mid-Pliocene 2691, 2720 (a-d) and PI control 5131, 5160 (e-h).The climatological mean for the austral winter months JJA is presented in the first two columns, and the boreal winter months DJF are in the last two columns.The denotes the range of model years over which the climatology was computed.

Table 1 .
Configuration common to all experiments described in this paper.

Table 2 .
Model details and diagnostics.

Table 3 .
Ocean temperature trends globally and for the different oceanic depths considered in Fig.3.All values are in units of • C century −1 .

Table 4 .
Proxy reconstructed mid-Pliocene mean annual surface air temperatures (MASATs) during the PRISM3 interval (3.3-3.0Mya)fromSalzmannetal. (2013)compared to that simulated with our 400 ppmv mid-Pliocene simulation Eoi 400 P. Subject to this caveat, it is observed that except for the two Siberian sites (Chara Basin and Lake Baikal) our simulated MASATs are in very good agreement with proxy inferences.The globally averaged MASAT is 16.8 • C for model Eoi 400 P and 17.3 • C for model Eoi 450 P (see Table

Table 6 .
Mean seasonal surface air temperatures and anomalies.