On the mechanisms of warming the mid-Pliocene and the inference of a hierarchy of climate sensitivities with relevance to the understanding of climate futures

We present results from our investigation into the physical mechanisms through which the mid-Pliocene, with ana atmospheric pCO2 of only ∼400 ppmv, could have supported the same magnitude of global warmth as that which has been projected for the climate at the end of the 21st century when pCO2 is expected to be three times higher. These mechanisms exploreallow us to understand the warming in terms of changes to the radiative properties of the surface, the clouds, greenhouse gases and changes to the meridional heat transport. We find that two-thirds of 5 the warming pervasive during the mid-Pliocene, compared to the pre-industrial, could be attributed to the reduction in the planetary emissivity owing to the higher concentrations of the greenhouse gases CO2 and water vapour, and the remaining one-third to the reduction in planetary albedo. We also find that changes to the orography and the pCO2 are the leading causes of the warming with each contributing in roughly equal parts to a total of 87 % of the warming and changes to the polar ice sheets responsible for the remaining warming. 10 Furthermore, we provide a mid-Pliocene perspective on ongoing efforts to understand the climate system’s sensitivity at various timescales and using multiple lines of evidence. The similarities in the boundary conditions between the mid-Pliocene and the present day, together with the globally elevated temperatures, make the mid-Pliocene an ideal palaeo time period from which to derive inferences of climate sensitivity and assess the impacts of various timescale-dependent feedback processes. We assess a hierarchy of climate sensitivities of increasing complexity in 15 order to explore the response of the climate over a very large range of timescales. The various sensitivities that we calculate provide insight on not only how the climate responds to a given forcing over a short timescale, but also on intermediate and very-long timescales. The latter category includes the impact of the feedback from the glacial isostatic adjustment of the Earth’s surface in response to the melting of the polar ice sheets.The picture that emerges is as follows: on the short timescale, owing to the influence of fast feedback processes, the climate sensitivity is 20 3.25 ◦C per doubling of CO2; sensitivity increases to 4.16 ◦C per doubling of CO2 on an intermediate timescale as the ice-albedo feedback becomes active and then sensitivity further increases to 7.0 ◦C per doubling of CO2 on long timescales due to the feedback from the glacial isostatic adjustment of the Earth’s surface in response to the melting of the polar ice sheets. Finally, once the slow feedbacks have stabilized, the sensitivity of the system drops to 3.35 ◦C per doubling of CO2. Our inference of the intermediate timescale climate sensitivity suggests that the projected 25


Introduction
The most recent assessment report by the Intergovernmental Panel on Climate Change (IPCC, 2013) expresses a high level of confidence that the global surface temperature at the end of the 21st century would be greater than 2 • C warmer than the 1850-1900 reference level in the two most extreme warming scenarios: RCP6 and RCP8.5.These levels of warming are driven by an increase in the radiative forcing to 6 W m −2 in RCP6 (Fujino et al., 2006;Hijioka et al., 2008) and 8.5 W m −2 in RCP8.5 (Riahi et al., 2007(Riahi et al., , 2011)), corresponding to atmospheric carbon dioxide partial pressures of (pCO 2 ) ∼ 850 and ∼ 1370 ppmv, respectively.
The mid-Pliocene time period (∼ 3-3.3 Mya) is also expected to have been roughly 2-3 • C warmer than the prein-Published by Copernicus Publications on behalf of the European Geosciences Union.dustrial (PI; Haywood and Valdes, 2004;Haywood et al., 2013;Kamae et al., 2016) and up to 3.8 • C warmer (Chandan and Peltier, 2017) using the latest reconstruction for the mid-Pliocene boundary conditions.Atmosphere-only modeling driven with proxy-inferred sea surface temperatures also yield a mid-Pliocene that is 2-3 • C warmer than the PI (Haywood et al., 2013).These levels of warming are nearly identical to those forecasted with the RCP6 and RCP8.5 scenarios, yet in contrast to those scenarios the mid-Pliocene warming took place under substantially less radiative forcing.Proxy inferences for pCO 2 yield values in the range of 300-400 ppmv (Seki et al., 2010;Martínez-Botí et al., 2015;Pagani et al., 2010;Badger et al., 2013), with some estimates pushing the lower bound to 200 ppmv (Bartoli et al., 2011;Tripati et al., 2009).A pCO 2 of 400 ppmv is typically used in coordinated efforts for the modeling of the mid-Pliocene, such as in the Pliocene Model Intercomparison Project phase 1 (PlioMIP; Haywood et al., 2011) and phase 2 (PlioMIP2; Haywood et al., 2016).See Chandan and Peltier (2017) for a recent review of proxy-based inferences of various climatological attributes of the mid-Pliocene.
There have been other hothouses during the recent Cenozoic, but the mid-Pliocene stands out among them as a unique "natural analogue" for the future, warmer climate.The recent interglacials were warmer than the mid-Pliocene, but they were ephemeral, and the climate quickly returned to a glacial state.In contrast, the mid-Pliocene was a much longer period of sustained warmth.More importantly, with regards to the future climate, the primary forcing during the interglacials was the increase in Northern Hemisphere summer insolation, whereas the primary forcing during the mid-Pliocene was the elevated pCO 2 .The mid-Pliocene also has the advantage of having had a continental configuration nearly identical to the present.As one travels further back in time to other hothouse worlds such as the middle Miocene climate optimum, the Oligocene warming, or the mid-Eocene climate optimum, the surface boundary conditions begin to increasingly depart from those characteristic of the present day.
Although the mid-Pliocene was in a state of quasiequilibrium, its pervasive warmth was set against the backdrop of a world that was undergoing a profound transformation that began in the mid-Eocene and ended at the onset of the 100 kyr Pleistocene glacial-interglacial cycles at ∼ 1 Mya.During this period, the global mean surface temperature was dropping (Hansen et al., 2013;Zachos et al., 2008), the atmospheric pCO 2 was declining (Zhang et al., 2013), ice sheets were beginning to form and grow in size in the polar regions (Zachos et al., 2001), and as a result the global mean sea level was falling (Haq et al., 1987;Hansen et al., 2013).Therefore, in Earth's history, the climate advanced from a hotter state to the (relatively cooler) mid-Pliocene; i.e., the climate evolved to the mid-Pliocene on a cooling trajectory.
But, for the purposes of understanding the future climate, we are actually interested in the transition from the present to a warmer climate.Therefore, it would appear to be worthwhile to examine this late Cenozoic transition in reverse.If we applied the changes known to have occurred in the climate system over the past 3 Myr in reverse using an appropriate climate model, we would reasonably expect to recover a "mid-Pliocene-like" state.The resulting quasi-equilibrium, mid-Pliocene-like state would be a function of the forcing and the internal feedbacks that are active on this warming trajectory.Therefore, this final state will invariably contain some information (an imprint) regarding the nature of the feedbacks and the climate system's response to those feedbacks and to the forcing itself.
The feedbacks that would operate in this expository scenario would also operate during the contemporary warming that is being driven by an increase in atmospheric pCO 2 .Some of these feedbacks are already well underway, and some will "come online" at different timescales in the future (see Sect. 4.4).If it were possible to extract some insight into the sensitivities of the Earth system to these feedbacks on the basis of the model's warm approach to the mid-Pliocene, then we could provide an independent constraint to the ongoing efforts to quantify, using a diversity of methods, the multiple sensitivities that characterize the potential of the climate to warm in the future.We therefore draw attention to the fact that the significance of the mid-Pliocene for understanding the future of the climate rests not only in its "warm state" persona but also in the fact that as a climate simulation approaches a mid-Pliocene-like state on a warming trajectory, the simulation will evolve under the influence of the same mechanisms, forcings, and feedbacks that the future of our climate will evolve.
Therefore, in several ways the mid-Pliocene is in an arguably privileged position to inform us about the warm climate of the future (a caveat to drawing parallels between the mid-Pliocene and the future is the role of aerosols and their indirect effects on the climate).In recognition of this advantage, in this paper we seek to conduct two fundamental analyses into the mid-Pliocene climate: first, we want to gain a deeper understanding into the origin of the warming within the mid-Pliocene simulation of Chandan and Peltier (2017) that was driven using the latest inference of the mid-Pliocene boundary conditions.It was found that with these new boundary conditions, it becomes possible for a climate model (specifically, our UofT-CCSM4 model) to, for the first time, reasonably accurately capture features of the proxy-inferred enhanced warming in the high latitudes during the mid-Pliocene.However, even with the new boundary conditions, our model shares a shortcoming with models using the older PlioMIP boundary conditions (Dowsett et al., 2013), which is that it is unable to simulate the proxyinferred negative SST anomalies at (i) the locations of tropical sites characterized by strong wind-driven oceanic equatorial divergence and (ii) at Mediterranean and Atlantic coastal sites.Therefore, understanding why the new boundary conditions have improved model reconstructions at some regions, namely the middle to high latitudes, but not at other regions is manifestly of interest to understanding the mid-Pliocene climate, and it can also guide strategies to refine boundary conditions so as to further reduce data-model disagreement.Furthermore, considering the similarities between the mid-Pliocene warmth and the forecasted warming discussed above, the insights gained from this investigation offer an independent constraint concerning where, for what reasons, and of what magnitude warming can be expected in the future.We have designed our investigation into this topic around three key questions.
1. What are the impacts of the differences in boundary conditions on the mid-Pliocene warming?(Sect.4.2) 2. What are the mechanisms sustaining mid-Pliocene warmth?(Sect. 4.3) 3. How do the modifications to the boundary conditions themselves change the effectiveness of these physical mechanisms?(Sect.4.3) The second theme of this paper is a mid-Pliocene inference of a hierarchy of climate sensitivity parameters that characterize the climate system's response to various forcings and feedbacks on various timescales.We build upon the work of PALAEOSENS Project Members ( 2012) and show how we can deduce these sensitivities using the suite of PlioMIP2 simulations.
Naturally, the quality of this effort can only be as good as the quality of the simulated mid-Pliocene because if we are not able to simulate a reasonable mid-Pliocene, then what possible grounds could we have to argue that the inferences of any sensitivity parameters are grounded in some reality?In this regard we are well positioned because, as mentioned earlier, our mid-Pliocene simulation is in very good agreement with proxy reconstructions of the climate of that time, notwithstanding the issues that continue to remain in the equatorial regions between approximately 45 • N and 45 • S.
This paper is organized as follows: in Sect. 2 we first introduce our numerical model and describe the model configuration employed and the process of the initialization of our numerical experiments.In Sect. 3 we provide a detailed introduction to and theoretical formulation of the analytical methods extensively employed in the interpretation of our simulation outcomes.The results are discussed in Sect. 4 and the conclusions of this study are summarized in Sect. 5.

The UofT version of CCSM4
The Coupled Climate System Model version 4 (CCSM4; Gent et al., 2011) is a numerical climate model that is comprised of four major sub-models: the atmosphere model, which is the Community Atmosphere Model version 4 (CAM4; Neale et al., 2013), the ocean model, which is the Parallel Ocean Program, version 2 (POP2; Smith et al., 2010), the land model, which is the Community Land Model version 4 (CLM4; Lawrence et al., 2012), and the sea ice model, which is the Community Ice Code version 4 (CICE4; Hunke and Lipscomb, 2008).These four components interact in a fully coupled manner with one another through a flux coupler without flux adjustment.The reader is referred to Chandan and Peltier (2017) for a summary of the improvements of these components over those that formed the basis of the previous CCSM3 model (Collins et al., 2006).
In Chandan and Peltier (2017) we introduced and employed a configuration of CCSM4 that was modified slightly from its default configuration in order to remove the impacts of two ocean model parameterizations that we argued to be highly tuned to modern-day conditions and therefore whose application to simulating paleo-oceans is questionable (see also Griffiths andPeltier, 2008, 2009;Peltier and Vettoretti, 2014).This configuration, which we referred to as the "University of Toronto version of CCSM4" (UofT-CCSM4) to distinguish it from its default configuration (NCAR CCSM4), has both the overflow parameterization and the tidal-mixing scheme disabled and the vertical profile of diapycnal diffusivity fixed to that used in the ocean component (POP1) of the CCSM3 model.Additionally, in order to be able to compare our mid-Pliocene simulations to the control simulations without any ambiguity associated with these changes, we also ran our PI and modern-day controls using the same model configuration.The new simulations presented in this paper are also run using this University of Toronto configuration.All simulations have been performed with the 1 • × 1 • resolution of the model.

Numerical experiments
The analyses presented in this paper involve several experiments that have been proposed for the PlioMIP2 program (see Haywood et al., 2016, Table 3).These simulations are referred to using the nomenclature first employed by Lunt et al. (2012b) and subsequently adopted with modifications for PlioMIP2.In this notation simulations are referred to by the form Ex c where c is the concentration of atmospheric CO 2 and x represents boundary conditions that have been changed from the PI such that x can be absent (for the case in which no boundary conditions have been modified) or it can be either or both "o" for a change in orography and "i" for a change in ice sheets (Table 1.)The "o" experiments also include changes to the bathymetry and the land surface.Therefore, in our paper whenever we refer to orography-related changes we mean the changes that arise due to the collectivity of modifications made to the orography, bathymetry, and the land surface boundary conditions.In all our simulations, vegetation is kept fixed to that applicable to the experiment case (mid-Pliocene or modern) and the concentration of trace gases and aerosols are kept fixed to PI values (Table 1).The PlioMIP2 simulations can be categorized as control simulations, mid-Pliocene simulations, or sensitivity simulations.

Control simulations
Two PlioMIP2 experiments, the Core simulation E 280 and the Tier 2 simulation E 400 , constitute our set of control simulations.These simulations were introduced in Chandan and Peltier (2017) in which they were referred to as E 280 P and E 400 P, respectively.Both simulations are configured with modern-day topography, bathymetry, land-sea mask (LSM), vegetation, and ice sheets.The atmospheric CO 2 concentration in the former simulation is set to the PI value of 280 ppmv and as such we will refer to that simulation as the "PI control", whereas the latter simulation has a modern-daylike CO 2 concentration of 400 ppmv and we will therefore refer to it as the "modern control."This simulation is "modern" only insofar as the CO 2 concentration is concerned; it does not have modern concentrations of other trace greenhouse gases or modern land units such as urban areas or agricultural land units.

Mid-Pliocene simulations
In this paper by "mid-Pliocene simulation" we mean a simulation that has full mid-Pliocene surface boundary conditions, i.e., the Eoi case with mid-Pliocene orography, vegetation, LSM, lakes, and ice sheets, and has an atmospheric CO 2 concentration in the range of most mid-Pliocene reconstructions (350-450 ppmv).Therefore, we do not refer to the PlioMIP2 Eoi 280 simulation, which has mid-Pliocene surface boundary conditions but a 280 ppmv atmospheric pCO 2 concentration, as a mid-Pliocene simulation and instead refer to it as a "sensitivity simulation" (next section) designed to explore the response of the climate to PI-like low atmospheric pCO 2 under mid-Pliocene surface boundary conditions.The boundary conditions we use are from the latest mid-Pliocene reconstruction from the PRISM group, PRISM4 (Dowsett et al., 2016), which has been adopted as the mandatory boundary conditions set for the PlioMIP2 program.We use the "enhanced" variant of the PRISM4-PlioMIP2 boundary conditions.The data set provides reconstructions for the mid-Pliocene orography, bathymetry, ice sheet distribution, lakes, vegetation, and soil texture and color.We have simulated two mid-Pliocene simulations, the Core simulation Eoi 400 and the Tier 1 simulation Eoi 450 , both of which were described in Chandan and Peltier (2017) wherein they were referred to as Eoi 400 P and Eoi 450 P, respectively.

Sensitivity simulations
We refer to those PlioMIP2 simulations that are neither entirely mid-Pliocene nor entirely present day (or PI) as sensitivity simulations.These have been proposed to investigate the impact on the climate system from changes to each of the boundary conditions that differentiate the mid-Pliocene from the PI.There are three configurations that constitute this set of sensitivity experiments: Ei, Eo, and Eoi.The Eo and Ei sensitivity experiments have two variants with 280 and 400 ppmv atmospheric pCO 2 , while the Eoi sensitivity experiment has 280 ppmv atmospheric pCO 2 (as discussed above, we categorize the other Eoi variants with higher CO 2 concentrations as mid-Pliocene simulations).We therefore report on five new PlioMIP2 simulations in this paper: Ei 280 , Ei 400 , Eo 280 , Eo 400 , and Eoi 280 .
Figure 1 shows the topography for each of these three configurations and the orographic anomaly of these configurations with respect to the PlioMIP2 modern orography.These anomalies are added to the local-modern orography to produce the atmospheric boundary conditions for each configuration in accordance with the anomaly method that we have previously employed (Chandan and Peltier, 2017).In the Ei configuration only those regions covered by presentday ice sheets are modified to yield the distribution of the mid-Pliocene ice sheets.This change is clearly visible in the orography anomaly for this configuration (Fig. 1d) in which large reductions in orography over most of Greenland, West Antarctica, and the subglacial basins of East Antarctica are the result of the reduced mid-Pliocene ice cover.The mid-Pliocene Greenland ice sheet reconstruction in PlioMIP2 is based on the modeling results from the PLISMIP project (Dolan et al., 2012).The configuration of the Antarctic ice sheets is provided by the PRISM3D (Hill et al., 2007;Dowsett et al., 2010) reconstruction, which was also used in the previous phase of PlioMIP (Haywood et al., 2013).In the PRISM3D reconstruction there is some thickening of the ice sheet over the interior of East Antarctica compared to the present day.
For this configuration a complication arises over West Antarctica where the marine grounded ice sheets inextricably link the local ice sheet distribution with the local LSM.In this case replacing the present-day ice configuration with the mid-Pliocene configuration would also result in unwanted changes to the LSM in the region (such as the changes seen in the Eoi configuration).Therefore, in order to maintain the LSM integrity in the ice sheet sensitivity experiments, we have replaced all regions where there is ice in the present day but not in the mid-Pliocene with very low orography (Fig. 1a).Lastly, in this configuration all other boundary conditions, including vegetation and river routing in the land model, are kept identical to present day (Table 1).
In the Eo configuration an anomaly in the orography exists at all locations where there is mid-Pliocene land, except at the locations of present-day ice sheets (Fig. 1b, e).In addition, anomalies also exist over parts of West Antarctica where the replacement of present-day marine grounded ice sheets by ocean has led to large negative values of the anomaly.
The presence of marine grounded ice sheets in Antarctica again introduces complications; in this case, changing the mid-Pliocene LSM would unavoidably alter the local ice distribution, whereas the aim of this configuration is to leave the ice sheets unchanged from the present day.Therefore the integrity of present-day ice sheets cannot be entirely maintained for this configuration.In addition to the changes to the orography and the LSM, this configuration also implements full mid-Pliocene bathymetry, and consequently the ocean model sees the same mid-Pliocene bathymetry as that implemented in the mid-Pliocene simulations of Chandan and Peltier (2017).The vegetation distribution, soil type, lakes, and the river routing schemes are also changed to those used in our mid-Pliocene simulations.
The methodology for implementing the Ei and the Eo boundary conditions into the UofT-CCSM4 model is identical to that employed in Chandan and Peltier (2017) for the implementation of the full-Pliocene (Eoi) configuration and the reader is referred to that publication for the details of this process.The Eoi configuration includes the complete set of changes made to all the surface boundary conditions.The topography and the orographic anomaly with respect to PlioMIP2 modern conditions are shown in Fig. 1c and f.

Model initialization
Simulating five sensitivity experiments in addition to the control and the mid-Pliocene experiments already reported in Chandan and Peltier (2017) (Levitus and Boyer, 1994).The simulation cesmpifv1mts is a non-PlioMIP2 preexisting PI control simulation.
the expenditure of a significant amount of computational resources.It was therefore necessary to carefully consider the process through which we initialized the five sensitivity experiments to ensure that the requirements for computational resources could be kept reasonable.Because the ocean component of the climate takes the longest time to equilibrate, we decided to initialize the sensitivity experiments wherever possible from an intermediate state obtained during the evolution of one of the PlioMIP2 simulations of Chandan and Peltier (2017), whose ocean is expected to be close to the ocean in the sensitivity experiments.An intermediate state, as opposed to an equilibrated state, was used to ensure that we do not incur an overwhelming influence of the simulation from which the initial conditions were derived while retaining the benefits of a faster spin-up.This initialization process is illustrated in Fig. 2. The two Ei experiments were initialized from the two control simulations because it was hypothesized that the changes to the ice sheets would be limited in their ability to significantly modify the state of the global ocean.On the other hand, the changes to the orography (and bathymetry) were expected to be more important not only because these changes are more pervasive, but also because the modified bathymetry would interact with the ocean circulation directly.For this reason, we expect that the ocean in the mid-Pliocene simulation would be close to the equilibrium state that the oceans would reach in the Eo experiments.On the basis of the same arguments, Eoi 280 could also have been branched from our 400 ppmv mid-Pliocene experiment; however, since both the Eoi mid-Pliocene experiments reported in Chandan and Peltier (2017) were initialized from modernday conditions (Levitus and Boyer, 1994), we decided to do the same for Eoi 280 .

The factorization technique
There are four critical differences between the boundary conditions of the mid-Pliocene and the PI: changes to orography, changes to ice sheet configuration, change in the concentration of atmospheric pCO 2 , and changes to ocean bathymetry.We will focus on the impacts of the first three of these differences.The impacts on the mid-Pliocene climate owing to the difference in bathymetry are expected to be important, particularly since a different bathymetry can alter the locations of deep water formation at high latitudes and modify the strength of the overturning circulation.However, we choose not to address this issue directly here not only because no sensitivity experiment for the ocean bathymetry is specified in PlioMIP2, but also because in using the factorization methodology discussed below, the treatment of ocean bathymetry as a separate factor would require eight additional experiments.Other differences between the mid-Pliocene and the modern day, which are mostly confined to the land and include changes to the soil type, vegetation distribution, and the distribution of wetlands and lakes, are expected to contribute much less to the mid-Pliocene warmth (Lunt et al., 2012b).
To understand the impact upon some observable of the climate due to a perturbation made to a specific boundary condition, one can simply difference the value of that observable between two simulations that differ from one another only with regards to that particular perturbation.For instance, if one wants to study the effect on the SAT caused by changes made to the elevation of a certain mountain range, then studying the SAT difference between two simulations (one with the revised elevations and one without) would suffice.However, the situation is clearly more complex when several boundary conditions are modified in concert and one wishes to discriminate the impacts on the observable from each of the individual boundary condition changes.This situation is commonly encountered in paleoclimatology as the past usually differs from the control by way of more than one boundary condition, such as the case between the mid-Pliocene and the PI.
When the impacts of two or more changes are sought, the same method can be extended with the help of additional simulations.However, this linearized approach would not give any information regarding the nonlinear interactions between the modified boundary conditions.Stein and Alpert (1993) proposed a factorization method to separate the impacts of several "factors" and quantify the magnitude of the nonlinear interactions between those factors, also referred to as their "synergy".This method has seen extensive application in paleoclimatology, such as for studying the following: the effects of ocean-vegetation interaction (Wohlfahrt et al., 2004) and the effects of soil schemes (Stärz et al., 2016) on the climate of the mid-Holocene; the impacts of boundary conditions during the last interglacial (Loutre et al., 2014), contributions to the warming during the last interglacial (Crucifix and Loutre, 2002), and other interglacials over the past 800 000 years (Yin and Berger, 2012); vegetation dynamics during the last glacial maximum (LGM; Claussen et al., 2013;Jahn et al., 2005) and the impacts of soil schemes on the LGM climate (Stärz et al., 2016); and contributions to the mid-Pliocene warming (Lunt et al., 2012b) and the strengths of the East Asian monsoons during the mid-Pliocene (Zhang et al., 2015).Lunt et al. (2012b) argue that the formulation of Stein and Alpert (1993) suffers from the shortcoming that it leads to a nonunique solution.To overcome this limitation they proposed a reformulation of this technique in which the synergy terms are partitioned evenly between the contributions of each factor to yield a symmetric and unique factorization.The principal drawback of these factorization methods is that they require 2 N simulations, where N is the number of factors being studied.For our case in which three boundary conditions are changed, this method requires a total of eight simulations.Due to the power dependence, adding another category, such as ocean bathymetry, would require eight additional simulations, whereas two additional categories (five in total) would require 24 additional simulations.Clearly, this approach is only feasible for small values of N.
In the factorization technique to be employed in this paper, we seek a decomposition of the form where ξ is the total anomaly of some observable of the climate, ξ , between the control simulation and some other simulation of primary interest (in our case these are respectively the PI and the mid-Pliocene simulations).The quantities dξ CO 2 , dξ orog , and dξ ice , herein referred to as the "factorized components or terms", are respectively the contributions to the total anomaly due to changes made to the atmospheric pCO 2 , orography, and ice sheets.The mathematical expressions for these terms are (Lunt et al., 2012b;Haywood et al., 2016) dξ where the right-hand side involving our experiment IDs should be read to imply the manner in which the observable ξ is to be combined from those experiments.

One-dimensional energy balance model
Here we discuss the 1-D EBM that we will use to explore several mechanisms through which the climate system adjusts its energy balance and determine to what extent each of them is responsible for the mid-Pliocene warmth.These mechanisms affect the zonal-mean energy balance by changing the Earth's radiative properties and its meridional heat transport.This approach was formulated by Heinemann et al. (2009) and first applied to understanding the Paleocene-Eocene climate.Subsequently, it has been applied to the Eocene (Lunt et al., 2012a;Loptson et al., 2014), the middle Miocene (Krapp and Jungclaus, 2011), the mid-Pliocene (Hill et al., 2014), and the snowball Earth (Voigt et al., 2011;Liu et al., 2013Liu et al., , 2017)).We begin with a review of the model as formulated by Heinemann et al. (2009) and show how the impacts of these mechanisms may be estimated using standard outputs from any GCM.Then we show how we can extend the 1-D EBM to further understand how the individual changes to the boundary conditions themselves impact each of these mechanisms.
The 1-D EBM (henceforth, just EBM) is written as where φ is the latitude, σ is the Stefan-Boltzmann constant, and SW ↓ TOA (φ), α(φ), H (φ), (φ), and τ (φ) are respectively the zonal-mean downwelling shortwave flux at the top of the atmosphere (TOA), the zonal-mean planetary albedo, the divergence of the zonal-mean meridional heat transport, the zonal-mean planetary emissivity, and the zonal-mean surface temperature.Henceforth, all the physical quantities are to be understood as zonal means and therefore for readability we stop using the adjective "zonal-mean" and drop the explicit dependence of these quantities on φ.The planetary albedo and the planetary emissivity can be easily calculated from the The surface temperature τ in the EBM is determined by the parameters α, σ , and H so that by defining τ ≡ T (α, , H ) the solution to the EBM as formulated in Eq. (3) becomes Now, the perturbations in surface temperature in response to a change in either α, σ , or H can be expressed as where the primes are used to designate modified values.If we assume that the changes in the parameters are small, then the total change in the temperature between the two simulations, T = T (α, , H ) − T (α , , H ), can be approximated as Equation ( 8) breaks down the total warming into contributions from three mechanisms: changes to the planetary albedo, changes to the planetary emissivity, and changes to the meridional heat transport.Since planetary albedo includes both the surface albedo and the cloud albedo we can further partition dT alb to split the effects of both by expressing it to first order as where dT swc is the temperature change induced due to changes in the shortwave forcing by the clouds and dT salb is the temperature change induced due to changes to the surface albedo.The latter can be diagnosed from the GCM's clear-sky fluxes (denoted by subscript "cs") as: Similarly, the planetary emissivity includes cloud emissivity and the emissivity of the atmospheric greenhouse gases (including water vapor), and a change to either or both can contribute to a change in the planetary emissivity.Thereby, dT emiss can also be partitioned to first order as where dT gg is the effect on temperature from a change in emissivity arising from a change in the concentration of atmospheric greenhouse gases and dT lwc is the effect on temperature due to a change in the longwave forcing by the clouds.Using GCM clear-sky fluxes, dT gg can be computed as So far we have reviewed the EBM methodology as it was formulated in Heinemann et al. (2009) for application to the Eocene.Equations ( 3) and ( 6)-( 12) allow us to quantify the contributions to the mid-Pliocene minus PI temperature anomaly from changes to the climate system's radiative parameters, the concentration of atmospheric greenhouse gases, and the meridional heat flux.However, they do not help us understand which of the modifications made to the various boundary conditions has contributed to changes to the climate's physical regime.An important motivation behind the PlioMIP2 sensitivity simulations is to isolate the impact of the various changes to the boundary conditions on the climate, in particular the surface temperature.Quantifying how the boundary condition changes themselves affect the mechanisms of warming is important to developing that understanding.We therefore need to be able to extend the EBM framework to satisfy this objective.
We achieve this by incorporating concepts from the factorization methodology to extend the Heinemann et al. (2009) formulation.We begin by using Eq. ( 1) to break each of the terms on the right-hand side of Eq. ( 8) into contributions arising from changes to orography, pCO 2 , and ice sheets as dT emiss = dT CO 2 ,emiss + dT orog,emiss + dT ice,emiss We illustrate the methodology through which the nine quantities on the right side of Eq. ( 13) can be evaluated with the help of examples.Analogous to Eq. ( 7), a quantity, for example dT CO 2 ,emiss , can be expressed as the temperature difference: where d CO 2 is the perturbation to the planetary emissivity caused by a change to the pCO 2 .This perturbation can be estimated by factorizing the PI to mid-Pliocene planetary emissivity change, , using Eq. ( 1) as where each term on the right-hand side is evaluated using Eq. ( 2) for the case ξ = .Similarly, the quantity dT orog,H , for example, can be written as where the change in the divergence of the meridional heat transport due to a change in orography, dH orog , is obtained by factorizing H .Likewise, all nine quantities in Eq. ( 13) can be evaluated.Additionally, the factorization of the clearsky components of α, , and H can be used in conjunction with the method outlined earlier in this section to disentangle the impact of clouds in each of the nine quantities.

Model evolution
The analyses presented in this paper are based on the climatologies computed over 30-year intervals.Table 2 lists the model years over which each simulation's climatology has been computed.These climatologies are computed towards the end of appropriately long numerical integrations in order to ensure that each of the simulations is in a state of quasi-equilibrium.We therefore begin the discussion of our results by first discussing the evolution and the quasiequilibrium characteristics of the five sensitivity experiments that are original to this paper.
Figure 3 shows the evolution of the 2 m surface air temperatures (SATs) and the TOA energy imbalances in the new simulations.The time series for the evolution of the SATs demonstrate that all models have reached the desired state of quasi-equilibrium.The TOA energy imbalance time series show that towards the end of their integrations, all models have either reached a steady state of energy imbalance (Eo 280 and Ei 280 ) or are progressing on a very slow declining trend.In all cases the climatological TOA energy imbalances are quite small (Table 2).Both these diagnostics strongly suggest that all models have reached quasi-equilibrium.Since the planetary TOA energy imbalance is controlled by the oceans, as they are the largest reservoirs of heat in the cli-mate system, it is instructive to examine the mean state of the oceans throughout the water column to further assess the reasonableness of these quasi-equilibria.Figure 4 shows the evolution of the mean ocean temperatures at all depths as Hovmöller diagrams for each of the sensitivity simulations.We see that in Ei 280 the temperature of the upper ocean (0-550 m) has nearly stabilized, while the ocean below continues to warm at a very slow rate of ∼ 0.04 • C century −1 (see Table S1 in the Supplement).A very similar situation is observed for Ei 400 but with a slightly greater rate of warming of the middle (550-1850 m) and the deep oceans (below 1850 m), which is clearly visible in the Hovmöller diagram.The SST and the upper ocean temperatures, however, have stabilized and do not exhibit any trend.This simulation has settled to a TOA imbalance of 0.16 W m −2 , which even after almost 1500 years of integration remains stubbornly identical to the climatological TOA imbalance of 0.17 W m −2 in the E 400 simulation of Chandan and Peltier (2017) from which it was branched (see Sect. 2.2.3).
As mentioned previously, the oceans in the Eo experiments were initialized from the same restart files produced at an intermediate time during the integration of the 400 ppmv pCO 2 mid-Pliocene experiment (Fig. 2).The ocean in this initial condition was previously evolving under the influence of a 400 ppmv pCO 2 atmosphere and is therefore much too warm for a 280 ppmv pCO 2 atmosphere.Thereupon, Eo 280 evolves by first eliminating all the excess heat contained in its oceans (Fig. 3b) from the model at the TOA to reach quasi-equilibrium after just 700 model years and remains in quasi-equilibrium for the subsequent 200 years of integration.Some of the excess heat also permeates downward to warm the deep oceans.The TOA energy imbalance, at only 0.03 W m −2 , is the best among all our experiments and the oceans from the surface down to a depth of ∼ 1900 m show zero temperature trend.The Eo 400 simulation, on the other hand, does not show much variation in its TOA energy imbalance during the course of its evolution.The absence of any large initial transient in the TOA energy imbalance suggests that the initial ocean state was a very close match to the Eo 400 case.The TOA energy imbalance starts at ∼ 0.15 W m −2 ,  then slowly increases to just over 0.20 W m −2 during which the deep ocean warms up (Fig. 3d), and then the TOA imbalance begins a very gradual decrease before settling into the climatological value of 0.07 W m −2 .The SST and the upper ocean show a trend of ∼ 0.03 • C century −1 , while the deepest ocean is warming at twice that rate.Among the five new experiments in this paper, only Eoi 280 was started from modern-day ocean temperatures (Levitus and Boyer, 1994) to maintain conformity to the Eoi experiments in Chandan and Peltier (2017).Therefore, while the Ei and the Eo experiments benefited from the spin-up already performed for their antecedent experiments, Eoi 280 had to undergo full spin-up.In hindsight, starting from modern conditions was a good choice as it appears that the modern-day initial state was fairly close to the state that the model wanted to evolve towards, at least in the upper ocean.The TOA energy intake during the first ∼ 1000 years was invested almost exclusively in warming the deep oceans as seen in Fig. 3e.After this warming, the TOA energy imbalance stabilized and has remained essentially fixed.

Surface warming: impacts of boundary conditions
In order to understand the impacts of the changes to the orography, the ice sheets, and the pCO 2 on the mean annual surface air temperature (MASAT), we begin by first reviewing the attributes of the MASAT anomalies obtained by means of simple differencing between the sensitivity experiments and the control experiments.The MASAT anomalies of the Ei and the Eo sensitivity experiments with respect to control experiments are shown in Fig. 5.In discussing these anomalies, our intention is to focus upon the impacts of changes to the surface boundary conditions (orography and ice sheets), and therefore the control experiments for the computation of these anomalies are chosen such that both the control and the experiment from which we are differencing the control have the same pCO 2 .This ensures that the anomalies of interest are not overwhelmed by the anomaly from the difference in pCO 2 .The effect of changing the pCO 2 is discussed subsequently.
Two MASAT anomalies that stem from the changes to the ice sheets, Ei 280 − E 280 and Ei 400 − E 400 , are shown in Fig. 5a-b.Both anomalies bear close resemblance to one another in the Southern Hemisphere (SH), indicating that changing the distribution of the ice sheets over the Antarctic continent (Sect.2.2.3) yields the same outcome regardless of the background value of pCO 2 .The strong warming seen in the Australian and the French sectors of Antarctica is due to the local reduction in orography resulting from the removal of some of the marine grounded ice sheets within the Aurora and the Wilkes subglacial basins and the accompanying decrease in ice sheet elevation in the surrounding regions (see orography anomalies in Fig. 1d).The cooling over the interior of Antarctica is also related to changes to the ice sheet elevation; in this case, a thicker East Antarc-tic Ice Sheet in the PRISM4 reconstruction.However, it is seen that the anomalies in the Northern Hemisphere (NH) do depend on the pCO 2 background.An extensive region of > 5 • C warming south of Greenland and extending from the Labrador Sea through the GIN Seas and into the Barents Sea, which is present in the anomaly of the Ei 280 experiment, is absent from the anomaly of the Ei 400 experiment.The anomaly over Greenland itself nevertheless remains identical between the two experiments and is orographic in origin as a result of the removal of much of the Greenland ice sheet.
Two MASAT anomalies that stem from the changes to the orography, Eo 280 − E 280 and Eo 400 − E 400 , are shown in Fig. 5c-d.In contrast to the ice-sheet-induced anomalies, these anomalies are larger in magnitude, spatially much more extensive, and almost entirely positive.In the SH the largest anomalies are found clustered in the coastal regions of Antarctica where present-day marine grounded sections of ice sheets have been removed to make way for an open ocean during the mid-Pliocene (see orography anomalies in Fig. 1e).This includes all of East Antarctica and small coastal pockets in West Antarctica, such as those at the locations of the present-day Amery Ice Shelf and the Wilkes and Aurora subglacial basins.In these locations the MASAT anomalies result from both (i) the reduction in orography and (ii) the change in the surface boundary conditions from ice to water.In the NH the warming over the Queen Elizabeth Archipelago and the Arctic coast of Siberia is directly associated with the increased fraction of land and the changes to local vegetation.Once again we see that there is a broad region of extensive warming around Greenland that is present in the anomaly of the 280 ppmv experiment, but absent or at least greatly diminished in the anomaly of the 400 ppmv experiment.
In addition to the two ice-sheet-related anomalies and the two orography-related anomalies that are shown in Fig. 5, there are still other combinations of PlioMIP2 simulations that could also have been used to assess the impacts of changes to these two BCs (Sect.3.1).For example, two additional differences, Eoi 280 − Eo 280 and Eoi 400 − Eo 400 , also in principle capture the impact of the same changes made to the ice sheets.Similarly, there are a total of four combinations that explore the impact of changes to the orography.These additional anomalies are not included in Fig. 5 because we wanted to focus upon just those differences that would be the most natural choices for an investigator wanting to test for the effects of changing the boundary conditions.Furthermore, by contrasting only two such differences we are able to satisfactorily make our point regarding the dependence of the results upon the background state.In the next section we will revisit the impacts of these boundary conditions using the factorization technique, which employs our complete set of experiments.
here are the 280 ppmv pCO 2 experiments.It is readily apparent that the background, which in the case of these anomalies is characterized by the ice sheet configuration and the orography, again yields considerable influence on the anomaly.The broad, elevated warming around Greenland, seen in the difference E 400 − E 280 and that is related to the local reduction in sea ice is absent in the other differences.This is because the extensive region of thin, fragile sea ice that is present in E 280 and whose melting in the E 400 experiment leads to this warming, is absent in the Ei, Eo, and Eoi cases (because of the warming introduced by these modifications) and is no longer available to be thawed by the warming from the increased pCO 2 .Curiously, the ≥ 5 • C warming over the Arctic present in the difference E 400 − E 280 and in the difference between the Ei simulations is absent from the other two differences.Essentially, this indicates that the impact of changing pCO 2 is less in the full Eoi case than when changing pCO 2 separately in Ei and Eo and adding them together.
The discussions thus far highlight the dependence of the response of the climate system to a perturbation in its boundary conditions on the stationary background state of the system.Both columns in Fig. 5 are intended to represent the Table 3. Annual and seasonal global mean SAT anomalies between the mid-Pliocene and the PI and the contributions to these anomalies from changes to the boundary conditions.The residual is the quantity T − (dT CO 2 + dT orog + dT ice ) and is seen to be very small in all cases.All values are in units of impacts of the same changes to the boundary conditions, but the inferred impacts differ significantly between them.This is the motivation for the application of the factorization methodology previously described and to which we now turn.
Climatological MASAT anomalies from a change to pCO 2 shown for all four combinations through which it can be deduced within our set of simulations.The 1st and the 99th percentiles of the anomalies are shown above each panel.The climatology years are noted in Table 2.

Insights from factorization
The total MASAT anomaly between the mid-Pliocene and the PI, T , and its decomposition using Eq. ( 2) into contributions originating from changes to pCO 2 , (dT CO 2 ), orography, (dT orog ), and ice sheets, (dT ice ), is shown in Fig. 7.The zonal means of these anomalies are shown in Fig. 8 and their global means are listed in Table 3.The difference T is the same difference previously discussed in Chandan and Peltier (2017).It is readily seen that the factorization has exposed distinct features between the anomalies of the sensitivity simulations shown in Figs. 5 and 6, which provides confidence in the methodology we are applying.
The warming from the 120 ppmv increase in pCO 2 between the PI and the mid-Pliocene is nearly uniform in the zonal direction (Fig. 7b) and increases in magnitude towards the poles where dT CO 2 reaches as high as 5 • C (Fig. 8).This polar amplification is reminiscent of the pattern of temperature increase that the climate system is experiencing at present.Away from the poles, dT CO 2 is nearly uniform with a warming ∼ 1 • C. The two lobes of locally enhanced warming on either side of the Antarctic peninsula are due to the enhancement of the warming from the reduction of sea ice.The global means in Table 3 show that dT CO 2 is the largest con-tributor to the mid-Pliocene warming.This conclusion was also reached by Lunt et al. (2012b) using the much older PRISM2 paleoenvironmental reconstruction (Dowsett et al., 1999).
The orographic contribution to the mid-Pliocene warming (Fig. 7c) is largely confined to the coastlines of Antarctica and to high latitudes in the NH where there have been significant changes to the orography and the LSM.In fact, the largest dT orog values exist exactly over those regions where the LSM has been changed.This makes sense as the change from ice-permafrost to (warmer) ocean represents a drastic change to the lower boundary conditions of the atmosphere.Surrounding these are expansive regions of reduced but still relatively high temperature anomalies over the oceans, which are due to changes to the sea ice concentration (the anomaly is seen to closely follow the contours of sea ice reduction).The elevated warming over the Tibetan Plateau present in the total anomaly, T , is seen to be picked up by the factorization process as an almost exclusively orography-related change.In the zonal mean, dT orog is seen to have an amplitude as high as dT CO 2 .In the midlatitudes of the SH the zonal mean dT orog reaches a global minimum, which is due to the fact that there is very little land in those latitudes and www.clim-past.net/14/825/2018/Clim.Past, 14, 825-856, 2018 MASAT anomaly between the mid-Pliocene and the PI ( T ) and its factorization into contributions from changes in atmospheric pCO 2 (dT CO 2 ), orography (dT orog ), and ice sheet configuration (dT ice ) using Eq. ( 2).The 1st and the 99th percentiles of the quantities plotted are shown above each panel.The factorization has consigned dT ice almost entirely to the high latitudes.This is seen both in the spatial pattern of the anomaly and in the zonal mean.In the latter, it is seen that dT ice has almost no influence throughout the tropics and the midlatitudes.This is again a positive outcome because we would be rightfully suspicious of the factorization if it had assigned any effect on the temperature anomaly in regions far removed from the locations of land ice removal.
These discussions show that the results from the factorization of the MASAT are physically comprehensible, which is certainly a requirement for us to be confident in the outcome of the factorization.A final check on the validity of these results is the assessment of how faithfully the total temperature anomaly and its factorized components satisfy Eq. ( 1).This is assessed with the help of the residual T − (dT CO 2 + dT orog + dT ice ), which is shown in Fig. S1 in the Supplement.The residual shows that the factorization has done a remarkably good job throughout the globe.A ribbon of large-valued residuals that extends from the Labrador Sea to the Barents Sea is the only notable exception.A possible reason for why factorization has performed poorly in this region is presented in the observation that the ribbon closely follows the margins of sea ice.As elaborated further in Sect.4.5 below, the complex thermodynamics of sea ice at its margins makes it difficult for the factorization method to disentangle the effects of the various changes to boundary conditions on the sea ice concentration.Since the sea ice is involved in thermal exchange with the air in its vicinity, the complications owing to the thermodynamics of sea ice are also communicated to the surface air, leading to the failure of the SAT factorization in those regions.This interaction with the sea ice should, however, also lead to the presence of large residuals off the coast of Antarctica.But this is not found to be the case.It is likely that this asymmetry between the two hemispheres is due to the additional impacts of northward heat transport and deep water formation.
Supplement Figs.S2-S5 show the factorization of the SAT anomalies for each of the seasons.These figures along with the global means of the factorized component for each season (Table 3) show that the warming from the elevated pCO 2 remains remarkably steady across all seasons.The warming from changes to the ice sheets is also seen to be relatively steady.The orographic contribution to the warming, on the other hand, is characterized by extremely large-amplitude differences between the seasons such that its impact during JJA is twice that during DJF.It is the variation in this component that is seen to reduce the total warming in winter and that leads to the increase in seasonality (compared to the PI) that was previously reported in Chandan and Peltier (2017).The global mean residuals of the seasonal factorizations are given in Table 3 and these are also found to be very small.

Mechanisms of warming
We next focus on understanding the mechanisms that lead to the observed warming in our 400 ppmv pCO 2 mid-Pliocene simulation compared to the PI.Since the orbital parameters and the solar constant are the same between the mid-Pliocene and the PI simulations, the TOA shortwave flux and its latitudinal distribution are also identical between the two.Therefore the warming has to do with how much of the incident energy is retained by the planet and how it is transported in the meridional direction to balance energy deficits.
In Fig. 9 we show the zonal-mean MASAT anomaly between the PI and the 400 ppmv pCO 2 mid-Pliocene experiment and its dissociation using the EBM methodology discussed in Sect.3.2 into contributions arising from changes to the Earth's radiative properties and the MHT.The curve labeled "GCM" is the MASAT anomaly obtained directly from the UofT-CCSM4 model and is identical to the curve identified as T in Fig. 8.The dashed gray line labeled "EBM" is the zonal-mean MASAT anomaly predicted using the EBM solution given in Eq. ( 6).This line is difficult to see in the figure because it virtually overlaps the GCM anomaly.The .Factorization of the zonal-mean mid-Pliocene minus PI temperature anomaly computed using the 1-D EBM (dashed gray) into contributions arising from changes to cloud albedo (red), cloud emissivity (purple), surface albedo (blue), greenhouse gas concentrations (including water vapor; green), and the total meridional heat transport (yellow).The solid gray line is the mid-Pliocene minus PI temperature anomaly computed directly from the GCM output.The dashed black line is the sum of the aforementioned individual contributions and is seen to almost perfectly overlap the total anomaly computed from both the EBM and the GCM.The abscissa is linear in the sine of latitude.
impressive agreement between the EBM and the GCM results is further reinforced through their meridional averages, which are 3.82 and 3.80 • C for the EBM and the GCM, respectively.
The other curves in Fig. 9 show the EBM-inferred contributions to the MASAT anomaly from changes to cloud albedo (Eq.9), surface albedo (Eq.10), cloud emissivity (Eq.11), emissivity from greenhouse gases (Eq.12), and the MHT (Eq.7).The warming from each of these mechanisms should add up to the EBM-inferred temperature anomaly if our methodology is soundly formulated.The dashed black curve in Fig. 9 is the sum of these quantities and it is indeed nearly indistinguishable from the EBM-deduced (and the GCM-deduced) profiles of temperature anomaly.The meridional averages (planetary averages) of these mechanisms given in Table 5 sum to 3.87 • C, which is only marginally greater than either the EBM or the GCM planetary average.
The dominant mechanism of warming during the mid-Pliocene is found to be changes to the planetary emissivwww.clim-past.net/14/825/2018/Clim.Past, 14, 825-856, 2018 The planetary emissivity reduction is almost entirely driven by the reduction in the atmosphere's emissivity owing to the increase in the concentration of the greenhouse gases (dT gg ) CO 2 and water vapor.The impact from changes to cloud emissivity (dT lwc ) also leads to a net warming, but it is negligible in comparison to the greenhouse effect.However, this is not a "textbook warming" resulting from a higher imposed pCO 2 and the concomitant change to the atmospheric water vapor content; it is seen that this effect, dT CO 2 ,gg , on its own would only contribute 1.25 • C, or roughly 54 % (Table 5) of the total warming that we have inferred in our analysis to be due to increased greenhouse gases.The other 46 % of the greenhouse effect comes from the increase in atmospheric water vapor as a result of the changes to the orography and the ice sheets, which contribute 0.79 and 0.26 • C to the greenhouse effect, respectively.This increase in water vapor is a direct result of the warming introduced by changes to these boundary conditions and the subsequent effect on the atmosphere's vapor-carrying capacity as required by the Clausius-Clapeyron relationship.
Water vapor is therefore an extremely potent greenhouse gas during the mid-Pliocene.Figure S6 contrasts the zonal mean of the vertically integrated total precipitable water content between the PI and the mid-Pliocene.The total precipitable water in the mid-Pliocene atmosphere is 1.48×10 16 kg, which is an increase of ∼ 22 % over the PI amount of 1.21 × 10 16 kg.In this regard, the mid-Pliocene greenhouse effect is a textbook greenhouse warming enhanced by the long-timescale feedbacks from changes to orography and ice sheets (Sect.4.4).Considering that the changes to the ice sheets, and in turn most of the changes to the orography, are ultimately due to the presence of high pCO 2 , the water vapor greenhouse effect from the changes to these two boundary conditions is itself an indirect effect of the high pCO 2 .
Warming owing to changes to the planetary albedo, dT alb , (shown in Fig. 10a) contributes 1.47 • C to the mid-Pliocene warmth.dT alb results from the net effect of cooling driven by changes to the cloud albedo (dT swc , −0.65 • C) and warming driven by changes to the surface albedo (dT salb , 2.12 • C).Cloud albedo changes are seen to cool the planet everywhere except in the tropics and the midlatitudes of the SH (Fig. 10c).The cooling is most pronounced at high latitudes where temperature anomalies as low as −10 • C are observed.Despite the presence of these large anomalies at high latitudes, the meridionally integrated dT swc only comes out to a modest −0.65 • C.This is because the high latitudes occupy a very small fraction of the Earth's surface area and the global average is thus influenced more by the modest-magnitude anomalies in the tropics and the midlatitudes that together account for most of the surface area.In these regions, the see-saw pattern in dT swc between the two hemispheres leads to cancelations that further reduce the net impact of changes to the cloud albedo.
It has been argued by Burls and Fedorov (2014) and Fedorov et al. ( 2015) that a reduction in the meridional gradient of cloud albedo during the mid-Pliocene could have been responsible for sustaining the reduced zonal and meridional SST gradients that have been proposed by some investigators (Wara et al., 2005;Ravelo et al., 2006;Fedorov et al., 2013) to have existed during that time (see Yang et al., 2016, for a recent review of several hypotheses that have been proposed to explain such a reduction in gradients).Burls and Fedorov (2014) were able to show, very impressively, the control that the meridional gradient of cloud albedo exerts on tropical SST gradients by explicitly modifying the cloud albedo gradient within the NCAR version of CCSM4.However, how such changes to cloud albedo could have manifested on their own is not clear.We find that no significant change to the meridional cloud albedo gradient between our mid-Pliocene and PI experiments arises naturally within our experimental setup; the small change in the tropic to extratropic cloud albedo that is seen (Fig. S7c) is an order of magnitude smaller than what would be required based on the analysis of Burls and Fedorov (2014) to simulate mid-Pliocene-like reduced SST gradients.Figure S7b for the zonal-mean anomaly of the surface albedo shows that throughout the range of latitudes, the mid-Pliocene surface albedo is either less than or equal to the PI surface albedo.Consequently, the surface albedo contribution to the MASAT anomaly is either positive or identically zero (Fig. 9).Table 5 shows that, globally, surface albedo change is the second most dominant (after the greenhouse effect) of the five mechanisms of warming.Regionally, its magnitude peaks at 60-70 • S and north of 80 • N. It makes perfect sense for the largest reduction in surface albedo (and the greatest temperature increase) to occur in these latitudes since the direct changes in the removal of land ice and the accompanying changes in the reduction in sea ice and snowfall on land have drastically reduced the surface albedo in these regions.
The largest contributor to surface albedo change (Fig. S7b) and correspondingly the largest contributor to surfacealbedo-induced warming (Fig. 10b) is orography.Besides the obvious impacts near the poles, orography has also been implicated by the factorization process to be the sole cause for the modest surface-albedo-induced warming at 15 and 30 • N.This warming at 15 • N is due to the albedo reduction from the vegetation changes that are made in conjunction with the orography changes.Specifically, it is related to the increase in grassland and dry shrubland at the southern margins of the Sahara Desert, which leads to a decrease in surface albedo (Fig. S8).The warming at 30 • N is related to the orographic changes over the Tibetan Plateau (Fig. S8) and the associated changes in ground snow cover.It is encouraging to see that these orography-related changes are correctly identified by the factorization process.The impacts of ice sheets and pCO 2 on surface albedo are confined to the polar regions, as there are no mechanisms (in the absence of a dynamic vegetation model) for either of them to affect the surface albedo away from the poles (CO 2 can also be an influence through changes to the snowfall rate over the tops of mountain ranges, but this is likely to be a small effect compared to the influence of dynamic vegetation).
The MASAT anomalies at high latitudes are seen to result from a delicate balance between individual contributions with very large magnitudes and opposing signs.In particular, the partial cancellation of the two effects with the largest amplitudes -surface albedo changes and cloud albedo changes -is very important in determining the MASAT anomaly.Since clouds are heavily parameterized in climate models, there is considerable uncertainty within these models with regards to the effects of clouds.In this case, what this uncertainty means is that any change to the cloud component would clearly impact the effect of planetary albedo and its partition into the effects coming from surface albedo changes and cloud albedo changes.The practical implication of this is that our ability to model the mid-Pliocene climate accurately depends in large part on the ability of our climate models to properly simulate clouds.This point is further underscored by the recent analysis of Sagoo and Storelvmo (2017), which demonstrated the profound impact that the microphysics scheme employed in the description of mixedphase clouds has on the climate of both high-dust time periods such as the Last Glacial Maximum and low-dust time periods such as the mid-Pliocene.
The redistribution of heat in the meridional direction cannot affect the total energy balance of the planet and therefore cannot lead to net cooling or warming, as it simply moves heat from areas of local excess to areas of local deficit.Accordingly, the planetary average of the heat transport should theoretically be zero.Table 5 shows that in our analysis the meridionally integrated dT H (−0.03 • C) does come out sufficiently close to zero.Although heat transport does not affect the net planetary energy balance, it does affect the local energy balance.Figure 9 shows that the local impact on tem-perature from a change to the MHT is non-negligible only in the high latitudes of both hemispheres.In these regions the EBM analysis shows that there is a cooling from a reduction in the MHT during the mid-Pliocene.

A mid-Pliocene perspective on climate sensitivity
Climate sensitivity has been a cornerstone concept in climate science for decades owing to its popularity as a means of quantifying the response of the climate to anthropogenic greenhouse gases.However, in the most general form, it is simply defined as the change in the globally averaged MASAT in response to an imposed radiative forcing.Several specific forms of climate sensitivity are possible within this general definition, subject to which processes are treated as forcings and as internal feedbacks.Much confusion has been created in the literature in recent years due to the introduction of several different forms of climate sensitivity together with inadequate differentiation between the forcings and the feedbacks (PALAEOSENS Project Members, 2012).As Previdi et al. (2013) note in a recent review of the subject: "It is important to distinguish between these forms (of climate sensitivity) in order to avoid confusion and to reconcile results from different studies that employ alternative sensitivity definitions." Two measures of climate sensitivity appear in the climate science literature more often than others.The first is the equilibrium climate sensitivity (ECS), more commonly known as the Charney sensitivity (CS) after Jule Gregory Charney who was the chair of the 1979 National Research Council "study group on carbon dioxide and climate" (National Research Council, 1979), which is widely recognized as the first modern assessment of the impact of greenhouse gases on the warming of the lower atmosphere.The CS is defined as the change in surface temperature per W m −2 of applied radiative forcing, once all the fast feedback processes such as changes to water vapor, sea ice, aerosols, and clouds have come to equilibrium.Since these processes equilibrate very quickly, with even the slowest among these, sea ice, taking only a few decades to reach equilibrium, the typical timescale over which CS is considered to be applicable is 100 years.For this reason CS has become important in discussions of the end-of-century climate, particularly within the context of the work that has been coordinated by the IPCC.
The second measure of climate sensitivity that is widely employed, especially in the paleoclimate community, is the Earth system sensitivity (ESS), which extends the CS by also including internal feedbacks from slow (long-timescale) processes such as changes to land-based ice sheets, vegetation, ocean circulation, and the carbon cycle and is therefore important to addressing the question of the state of the climate over tens of thousands of years.However, in practice, only some of these slow feedbacks are included in calculations for the ESS, and others are either omitted or included as forcings.To make matters more complicated, proxy-based estimates of ESS are also affected by orographic changes (as the orography in the past was different from today) that should ideally be considered as a separate forcing, but owing to the difficulty in quantifying the magnitude of forcing from the orographic changes, the orographic effect is unfortunately considered to be an internal feedback in response to an external forcing.A further difficulty regarding the orographic changes, as we discuss below, is that some orographic changes are in response to the melting of polar ice sheets and the subsequent adjustment of the surface of the solid Earth, and these may in fact be reasonably considered to constitute slow internal feedbacks as opposed to externally applied forcing.
We will therefore seek herein to infer a hierarchy of climate sensitivities of increasing complexity in order to explore the response of the climate over a very large range of timescales and in response to various choices of forcings and feedbacks.The sensitivities that we calculate provide insight not only on how the climate responds to a given forcing over a short timescale, but also on intermediate and very long timescales.We begin with a discussion of the CS that is relevant on the shortest timescales (∼ 100 years) and has therefore become important to discussions of the climate at the end of the 21st century.Then we assess the sensitivity on an intermediate timescale on which the feedback from ice sheets starts to become important.On the longest timescale we consider the ESS for which we expand upon the work previously presented in Chandan and Peltier (2017).
Before proceeding with the discussion of these sensitivities, it is important to make a comment regarding the notation that we will be using in the remainder of this paper.PALAEOSENS Project Members (2012) (henceforth, simply PP) have proposed a nomenclature with the aim of building a common vernacular among researchers to unambiguously discuss climate sensitivities and to avoid the confusion that Previdi et al. (2013) have cautioned us about.We will therefore employ the PP notation, and we begin with a brief overview of the fundamentals of this notation.The climate sensitivity, S, in its general form is expressed as where R is some imposed radiative perturbation and T is the ensuing globally averaged change in MASAT.Furthermore, for some abstract process P , we label the radiative perturbation from that process as R [P ] and the temperature response due to that perturbation and from all associated fast feedbacks as T [P ] .Finally, the "specific sensitivity" (a term introduced by PP) of the climate system to the process P is defined as We now begin with the discussion of the shortest-timescale sensitivity, CS.Since the only radiative forcing is that from the difference in the pCO 2 between the mid-Pliocene and the present day, R in Eq. ( 17) is just the specific forcing The change in temperature in response to this forcing, T [CO 2 ] , is just that due to both the imposed radiative forcing and feedbacks from all fast-timescale processes.Accordingly, CS (which PP denotes as S a where the superscript stands for "actuo", meaning present day) is given as It is thus seen that CS is just the specific sensitivity to the change to pCO 2 , i.e., S a = S [CO 2 ] .The radiative forcing is easily calculated (Myhre et al., 1998) as In contrast, it is not at all straightforward to compute T [CO 2 ] .To be sure, given our objective of contributing a mid-Pliocene perspective on CS, T [CO 2 ] must be derived in some way from our set of simulations that incorporates the effect of the mid-Pliocene boundary conditions.This is one reason why we cannot use for T [CO 2 ] the temperature difference between the control simulations E 400 and E 280 , as neither of these simulations communicates any information regarding the unique mid-Pliocene surface boundary conditions.Furthermore, a general argument as to why we cannot simply use the temperature difference between any two of our simulations, differing only in the CO 2 concentration (e.g., Eoi 400 and Eoi 280 ), is that the results above show that there is no reason to prefer any one such difference over any of the other differences.This is where the factorization methodology becomes useful because it tries to circumvent the subtlety surrounding the dependence of the response of a system on the background to tease out the specific response from changing CO 2 .
Factorization of the MASAT anomaly comes to our rescue and provides a way (Eq.2) to estimate T [CO 2 ] , whose value was previously calculated to be 1.67 • C (Table 3).With this we infer the CS of the UofT-CCSM4 model to be 0.88 • K W −1 m 2 or, equivalently, 3.25 • K per doubling of pCO 2 .This is remarkably close to the CS estimate of 3.2 • K per doubling of pCO 2 that has been reported by Bitz et al. (2012) for the closely related NCAR-CCSM4 model.However, there is an important difference between their study and ours; their result is a modern-day inference of the CS as their simulations are performed against a background defined by modern-day orography, LSM, ice sheet distribution, bathymetry, and vegetation.This is different from our case since the temperature difference we use is the average (Eq.2) of all the possible ways in which the temperature difference between a 280 and 400 ppmv atmosphere can be calculated with differing background states.
Charney sensitivity is a reasonable measure of the climate system's sensitivity only on a timescale of ∼ 100 years because it does not incorporate the effects of the slower feedback processes.Among the various slow feedback processes that operate in the climate system, one which holds particular relevance for the future is the feedback from the large-scale wasting of polar ice sheets.This process operates on timescales of a few hundred to a thousand years, which is not fast enough to qualify for the purposes of CS but is still considerably faster than the slow feedback processes such as changes to ocean circulation and the carbon cycle and changes to orography induced by glacial isostatic adjustment (GIA; for recent discussions of the global process of isostatic adjustment see Peltier et al., 2015, andRoy andPeltier, 2017).State-of-the-art numerical modeling of the Antarctic ice sheet, forced with atmospheric boundary conditions derived from regional climate model simulations using projected emission scenarios, suggests that the entire West Antarctic Ice Sheet could collapse in just 250 years (Pollard and DeConto, 2016).The feedback from land ice is therefore relevant to the understanding of our climate a few hundred years into the future and indispensable to discussions of lasting climate change mitigation strategies.We can quantify the effect of this feedback with a new intermediate-timescale sensitivity that is constructed by extending CS to include the land ice feedback (Fig. 11b).We denote this new sensitivity as S i and write it as where the total temperature change is now expressed as the sum of the temperature change due to the higher pCO 2 , T [CO 2 ] , and the temperature change as a result of the feedback from the melting of the ice sheets, T [ice] .The factorization of the MASAT anomaly gives T [ice] = 0.47 • C (Table 3) and we therefore infer that S i = 1.12 • K W −1 m 2 or, equivalently, 4.16 • K per doubling of pCO 2 .This is ∼ 30 % higher than CS.
On the longest timescales the potential for warming in the climate system is quantified by its ESS.Chandan and Peltier (2017) arrived at a mid-Pliocene inference for the ESS of ∼ 7 • C per doubling of CO 2 (or ∼ 2 K W −1 m 2 ).Their calculation assumed that the only forcing upon the mid-Pliocene climate was from the change in the concentration of the atmospheric CO 2 and that the responses from changes to the ice sheets and orography could be treated as internal feedbacks operating over very long timescales.The treatment of the ice sheet response as an internal feedback is a reasonable assumption since land ice extent and volume respond to the direct radiative perturbation introduced by the greenhouse gases.However, as we will discuss below, for the case of orography this assumption is only partially accurate.
Almost all of the changes to the LSM during the mid-Pliocene occurred due to the reduction in land ice volume and the ensuing increase in global sea level, which can be considered as a slow internal feedback to a higher pCO 2 .In addition, the component of the mid-Pliocene orography that is due to the isostatic rebound of the surface of the solid Earth (orographic change when the surface is above water, bathymetric change when the surface is below water) following the removal of the ice sheet loads can also be considered as an internal response.Because both of these processes are controlled by the GIA process we refer to these two components of orography collectively as the GIA component of orography.Figure 11d illustrates the feedbacks from the interaction between land ice and the GIA component of orography included in the ESS definition of Chandan and Peltier (2017).Mathematically, this definition is expressed as But two additional components of the mid-Pliocene orography, (i) that due to mantle-convection-induced deformation on the surface, usually referred to as dynamic topography (see Forte et al., 1993, for the modern description), and (ii) the filling over northern Canada of the erosional features Clim.Past, 14, 825-856, 2018 www.clim-past.net/14/825/2018/created by numerous Pleistocene glaciations are unrelated to the forcing from CO 2 and should actually be treated as external forcings.However, in the calculations of Chandan and Peltier (2017) these were also treated as internal feedbacks since there is no entirely accurate means available at present to isolate them from the GIA component of orographic forcing.To the extent that the contribution from these components is not disproportionately large compared to the contribution from the GIA components, the estimate from Chandan and Peltier ( 2017) is a good representation of the ESS.Furthermore, this choice of the definition of ESS allowed Chandan and Peltier (2017) to make a straightforward comparison to proxy-derived ESS estimates (Pagani et al., 2010) and to the estimates from the original PlioMIP program (Haywood et al., 2013).However, if the assumption of the dominance of the GIA components of the orography anomaly over the non-GIA components is incorrect, then the inference by Chandan and Peltier (2017) could be an overestimation.In this case, we need to consider orography as a forcing rather than as an internal feedback.A different sensitivity measure (Fig. 11c) is then necessary, which is the specific sensitivity to changes in pCO 2 and orography (denoted S p by PP): where R [OR] is the orographic forcing and T is the total temperature change between the PI and the mid-Pliocene, which is 3.80 • C (Chandan and Peltier, 2017).The pCO 2 forcing R [CO 2 ] is again inferred using Eq.(20).To compute R [OR] , following PP we assume that the efficacy (Hansen et al., 2005), e, of orographic forcing is the same as the efficacy of pCO 2 forcing; i.e., we assume . By eliminating e from these two equations and substituting the result in Eq. ( 23) we can convert the problem of inferring the orographic radiative forcing into one of inferring the orographically induced temperature change: Factorization of the MASAT anomaly (Eq.2) yields or, equivalently, 3.84 • K per doubling of pCO 2 .A few points need to be made regarding this calculation: firstly, the efficacy of the orographic forcing is almost certainly not the same as that from pCO 2 (Hansen et al., 2005).Therefore, a more general expression for the forcing in the denominator of Eq. ( 24) would be where f is the ratio of the orography efficacy to the pCO 2 efficacy.Therefore, the inference of S [CO 2 , OR] is subject to the uncertainty in the magnitude of f .For example, if f = 0.9, i.e., if the orography-induced radiative imbalance were only 90 % as effective as a pCO 2 -induced radiative imbalance of the same magnitude, then S [CO 2 , OR] = 1.09 • K W −1 m 2 or, equivalently, 4.03 • K per doubling of pCO 2 .On the other hand, if f = 1.10, then S [CO 2 , OR] = 0.99 • K W −1 m 2 or, equivalently, 3.66 • K per doubling of pCO 2 .Indeed, the value of f could be well beyond this ±10 % bound.For instance, Modak et al. (2016) recently studied the efficacy of the solar forcing and found it to be ∼ 80 % relative to the CO 2 forcing in the NCAR CAM5 model, which is a successor to the CAM4 atmospheric model in the presently employed UofT-CCSM4 coupled climate model.Secondly, regardless of the efficacy of the orographic forcing (i.e., regardless of f ), the radiative forcing in S [CO 2 , OR] is really just the forcing R [CO 2 ] scaled by a factor, which in our case is ∼ 1.922, and therefore this total radiative forcing can be solved for the atmospheric pCO 2 that would on its own result in the same magnitude of forcing and comes out to 555 ppmv, or nearly double the PI concentration.It is well known from experiments with GCMs in which pCO 2 is varied that the efficacies of various feedback processes in the climate, and thereby the climate sensitivity itself, depend on the magnitude of forcing (Hansen et al., 2005;Colman and McAvaney, 2009) and generally climate sensitivity reduces with increased forcing.It is partly for this reason and partly because of the uncertainty in f that S [CO 2 , OR] is less than the ESS.
To this point we have considered two forcings, pCO 2 and orography, and have either ignored land ice altogether or treated it as an internal feedback.We now consider the climate sensitivity S [CO 2 , OR, LI] , in which land ice is also treated as a forcing (Fig. 11e).Since the Greenland ice sheet (GIS) and the West Antarctic Ice Sheet (WAIS) were absent during the mid-Pliocene (Jansen et al., 2000;Kleiven et al., 2002;De Schepper et al., 2014), the positive ice-albedo feedback involving these ice sheets was clearly not operational (only once the ice sheets have begun forming does this positive feedback become active).However, in the absence of these ice sheets, local warming will be sustained through the presence of low albedo on the Earth's surface and a low orography.This warming in turn makes the inception of these ice sheets difficult.In this manner, the ice sheet configuration acts as a forcing.This consideration is of particular significance for the GIS as it has been shown on the basis of a twodimensional continent and ocean-resolving energy balance analysis (North et al., 1983;Hyde et al., 1990) that within specific ranges of parameters, such as insolation and CO 2 , Greenland is ideally suited for the inception of an ice sheet based solely on its location and the difference in heat capacity between its landmass and the surrounding ocean.Consequently, as the global drawdown of CO 2 continued during the mid-Pliocene and the CO 2 level started to become favorable for the inception of the GIS, the preexisting forcing that existed in the absence of the GIS would have acted as a barrier that the system would have first needed to overcome.Similarly for Antarctica, DeConto et al. (2008) show that the CO 2 threshold for glaciation is ∼ 750 ppmv, but in spite of the mid-Pliocene CO 2 level remaining much less than this threshold value there was only intermittent glaciation over West Antarctica and large subglacial basins of East Antarctica (De Schepper et al., 2014).It could again be argued that the absence of the ice sheets in these parts of Antarctica and the ensuing warmth they helped sustain made it difficult for the ice sheets to develop.
Therefore, during the mid-Pliocene before the inception of the GIS, the WAIS, and extension of the EAIS into the subglacial basins of East Antarctica (and assuming relative stability of the interior sections of the EAIS) the response of the climate system to changes in forcing would likely have been characterized by S [CO 2 , OR, LI] .Similarly, if the future radiative forcing could be maintained at a level comparable to the present day (equivalent to the mid-Pliocene), then in the very distant future, long after the feedback processes that influence the ESS have stabilized, the sensitivity of the climate system would be characterized by S [CO 2 , OR, LI] , which is given as R [LI] is the forcing from land ice that can be obtained by appealing to the same efficacy argument that was used to calculate R [OR] earlier.We get S [CO 2 , OR, LI] = 0.9 • K W −1 m 2 or, equivalently, 3.35 • K per doubling of pCO 2 .

Sea ice
We include a brief discussion of sea ice in our results because it is a major determinant of the planetary albedo, the reduction of which was found to be the second largest contributor to mid-Pliocene warming.The climatological sea ice concentration anomalies of the five sensitivity experiments and the mid-Pliocene experiment with respect to control experiments are shown in Fig. 12.The anomalies are shown on a form of Lambert azimuthal equal area projection to show both poles in one figure.The analysis of these figures suggests that the orography changes induce a greater reduction in sea ice than the changes to the ice sheets.
We again turn to factorization to see if we can better understand the contributions of each of the boundary condition changes.Figure S9 shows the outcome of applying the factorization methodology to the sea ice concentration anomaly between the mid-Pliocene and the PI.It is seen that the increase in pCO 2 and the changes made to the orography together account for most of the loss of sea ice.The largest orographycontributed reductions to sea ice are confined to the Hudson Bay, the Queen Elizabeth Archipelago, and the coastlines of the Arctic Ocean.This is a consequence of the fact that in the factorization process we have explicitly ascribed all the anomaly at any location where the LSM has changed to orography.However, unlike the MASAT, the residuals from the factorization of sea ice anomaly are very large (Fig. S10).This greatly diminishes any prospects of gaining meaningful insights into sea ice changes using the factorization methodology.The large residuals are confined to the margins of the sea ice where the sea ice is thinnest, least concentrated, and most prone to melting.The warming induced from all three boundary condition changes will compete to melt the ice.The thermodynamics of sea ice is a notoriously nonlinear process with many factors influencing the melting and does not simply scale with the magnitude of warming.Along with the fact that there is only so much ice to melt, disentangling the effects from the three processes understandably becomes challenging for the factorization process.
The annual means of the sea ice area and volume in each of our simulations are provided in Table 2.The Arctic sea ice is found to be much more susceptible than the Antarctic sea ice to the climatic changes induced as a result of making a change to any of the boundary conditions.The percentage of reduction in volume is always much greater than the percentage of reduction in sea ice area.This points to thinning of the sea ice in response to the perturbations.As a result of the thinning from the changes, the average thicknesses of the NH and the SH sea ice drop from ∼ 4.10 and ∼ 1.48 m, respectively, in our PI control to ∼ 0.75 and ∼ 0.88 m, respectively, in the 400 ppmv pCO 2 mid-Pliocene experiment.
A careful analysis of the data in Table 2 also reveals a very peculiar dependence of the sea ice response to changes to surface boundary conditions on the pCO 2 of the system.In any given experimental configuration of the climate model, it is found that the NH sea ice responds more strongly to the changes when the system itself is in the low pCO 2 state than when it is in the high pCO 2 state; i.e., for any given response the relative change is greater in the low pCO 2 case.The opposite situation is found for the SH sea ice, which exhibits a greater response to the boundary conditions when in the high pCO 2 state than in the low pCO 2 state.

Meridional heat transport
The impacts of the changes to the boundary conditions within PlioMIP2 upon the planetary meridional heat transport (MHT) are relevant not only for understanding the mid-Pliocene puzzle, but they are also imperative to continuing efforts to understand the influence on this immensely important component of the climate system under contemporary global warming.It was recently shown in Chandan and Peltier (2017) that with the PlioMIP2 boundary conditions there is a small but not insignificant reduction in the MHT compared to both the PI and the modern-day climates and in both the atmospheric and the oceanic components of the transport.Here, our objective is to understand how each of the individual changes to the boundary conditions have affected the MHT.In what follows, we adhere to the format  2. used earlier for the discussion of the SAT; first we consider MHT changes through selected comparisons, then we factorize the MHT anomaly and assess how reasonable our conclusions derived from the comparisons were.
Figure 13a-b compare the atmospheric, oceanic, and total MHTs in the Ei 280 and the Ei 400 simulations to the E 280 and the E 400 controls, respectively.It is readily evident that modifications to the polar ice sheets have had a fairly inconspicuous effect on the total MHT and the oceanic heat transport (OHT).Some differences are observable upon careful examination, but the differences are significantly smaller than the magnitudes of the heat transports (i.e., the relative deviations are small).It is also clear that there is very little dependence on the background state characterized by the atmospheric pCO 2 , except in the midlatitudes in the NH.The OHTs can be further decomposed into transports in the Atlantic and the Indo-Pacific basins.These decompositions for the Ei experiments are shown in Fig. 14a-b.Overall, the modifications to the ice sheets have resulted in an increase in OHT in the Atlantic basin in the low pCO 2 case and a decrease in the high pCO 2 case.At 30 • N, for example, the Atlantic in Ei 280 transports ∼ 6 % more heat than in E 280 , whereas the Atlantic in Ei 400 transports ∼ 5 % less heat than in E 400 .A more stark difference is seen around 50 • N where the Atlantic OHT in Ei 280 is ∼ 17 % larger than E 280 , while that in Ei 400 is ∼ 13 % less than the control.
The MHTs in the Eo experiments are likewise compared to their equivalent pCO 2 control experiments in Figs.13d-e and 14d-e.In the NH, the orography changes lead to a robust decrease in the total MHT and the magnitude of the decrease is significantly larger than any decrease in the MHT from changes to the ice sheets.For instance, compared to the PI, the largest decrease in the total MHT in Eo 280 is ∼ 5 %, while it is < 2 % in the case of Ei 280 .This decline in the northward total MHT is seen to be driven almost entirely by the reduction in the northward AHT.In contrast, the effect of orography on the total MHT and the AHT in the SH is very small.The dependence of the MHT on the background pCO 2 is seen to be small in the NH and almost nonexistent in the SH.The OHT in the Atlantic basin appears to be the only transport noticeably affected by the background pCO 2 (Fig. 14d-e) wherein we see that, analogous to the Ei case, in the midlatitudes of the NH there is an increase in the OHT compared to control in the low pCO 2 state and a decrease in the high pCO 2 state.Also, it should be noted that unlike the Ei case, orography changes are seen to lead to a large reduction in the northward OHT in the Indo-Pacific basin.
effects are combining nearly linearly, thereby enabling the larger-magnitude effects from the changes to the orography to dominate the much smaller impacts from changes to the ice sheets.The impacts from changing pCO 2 are presented for three configurations of the surface boundary conditions in the rightmost columns of Figs. 13 and 14.It is seen that regardless of the surface boundary conditions of the experiments, the impacts of changing pCO 2 are nearly identical in all cases.Three important observations can be made: firstly, changing pCO 2 has no impact whatsoever on the AHT in either hemisphere.Secondly, in addition to having no impact on the AHT in the NH, there is also no effect on the total MHT or the OHT in the NH.And lastly, the only impact of changing pCO 2 appears to be a small increase in the southward total MHT driven by the increase in OHT in the oceans of the SH.

Factorization results
The factorizations for the MHT, AHT, and OHT anomalies (mid-Pliocene minus PI) are shown in Fig. 15.The first feature to notice is the asymmetry between the two hemispheres in the total MHT anomaly: during the mid-Pliocene there is a net reduction in the heat transport in the NH and a net increase in the heat transport in the SH.In the midlatitudes of the NH, the net decrease is largely the consequence of the decrease in the AHT, whereas in the tropics and the subtropics where the OHT is strongest, it is the reduction in the OHT that drives the decline in the MHT.OHT is also largely responsible for the increase in the total MHT throughout the SH midlatitudes.
The AHT anomalies are largest in the midlatitudes, where AHT is the dominant mode of heat transport.Orography is seen to be the dominant driver of OHT anomalies, while the impacts from ice sheets and atmospheric pCO 2 are considerably smaller.The AHT anomaly is entirely negative in the NH and mostly positive in the SH.The factorization performs very well at all latitudes except for between 40 and 60 • N where the sum of all factorized components (dashed black line in Fig. 15b) underestimates the total anomaly.
In the mid-Pliocene climate the OHT changes sign at ∼ 5 • S latitude (Fig. 14); north of this latitude the heat transport is northward, and it is southward south of this latitude.The OHT anomalies appear to be antisymmetric about this latitude such that the northward transport is inhibited, while  2.
the southward transport is energized.The absolute magnitudes of the changes are nearly the same on either side.The factorization shows that the reduction in the northward transport is dominated by the orographic component, whose impact outside of the range of this OHT anomaly is negligible.In contrast, the positive OHT anomaly in the SH is dominated by the pCO 2 response, which is also much smaller outside of the range of this SH anomaly.

Discussion and conclusions
In this paper our focus has been on two critical issues, namely understanding the origins of the mid-Pliocene warmth that has been recently simulated by Chandan and Peltier (2017) using the latest inference of the mid-Pliocene boundary conditions and arriving at a mid-Pliocene inference regarding the sensitivity of the climate system on a range of timescales.We have organized our discussion in this section around these two topics.

Origins of mid-Pliocene warmth
Regarding the origin of the mid-Pliocene warmth, we have designed our investigation around the following three questions discussed in the Introduction as providing the primary motivation for this work.
5.1.1What are the impacts of the differences in boundary conditions on the mid-Pliocene warming?
With our suite of PlioMIP2 simulations, we have been able to perform the very first nonlinear factorization of the MASAT anomaly between the mid-Pliocene simulation forced with the most recent PRISM4 boundary conditions and the PI.A total of eight simulations were needed to employ the technique of Stein and Alpert (1993) and Lunt et al. (2012b) for three categories of modifications: orography, ice sheets, and atmospheric pCO 2 .We find that the factorization approach performs remarkably well in disentangling the effect of each of the three changes on the total warming.The residuals from the factorization are negligible in magnitude over > 95 % of the surface area of the Earth.The higher pCO 2 in the mid-Pliocene atmosphere (relative to the PI) is found to be the primary cause of warming, contributing ∼ 45 % to the total MASAT anomaly.The contribution from orography change is a very close second accounting for ∼ 42 % of the anomaly.The remaining 13 % of the temperature difference is associated with the reduction in polar ice sheets.In contrast, using the much older PRISM2 boundary conditions for the mid-Pliocene along with the UK Met Office coupled ocean-atmosphere GCM HadCM3 version 4.5 (Gordon et al., 2000), Lunt et al. (2012b) found the contributions from CO 2 , orography, and ice sheets to be 48, 31, and 21 %, respectively.This indicates that, save for the inter-model differences, the large changes to the orography between the newest PRISM4 reconstruction and the older PRISM2 reconstruction have had a significant impact on the warming.Since these orography changes are concentrated at high latitudes, the impacts of these changes at a re-gional scale are even higher than those at the global scale.The differences in the ice sheet contribution are also due to the differences in the polar ice sheets between PRISM2 and PRISM4.
From the factorization of the seasonal SAT anomalies it is found that the contributions from atmospheric pCO 2 and ice sheets are relatively stable through all seasons.In contrast, the contribution of orography is found to change by 100 % between DFJ and JJA.The reduction in the orographic warming during DJF is found to be the reason for the increased seasonality during the mid-Pliocene that was previously identified and discussed in Chandan and Peltier (2017).
In this paper we have extensively applied the factorization technique as formulated in Lunt et al. (2012b), and in doing so we have come to recognize some of its limitations.The performance is not identical across the range of climate observables.While the method is extremely adept at breaking the SAT anomaly into various components, for the case of sea ice concentration, factorization is seen to do a poor job as evidenced by the large residuals that it produces.Similarly, the method performs very well in breaking down the changes to albedo, emissivity, and total MHT and OHT, but it comes up short in factorizing the AHT.It is therefore very important that practitioners of the method consider the residuals to be sure that their conclusions are sound.

What are the mechanisms sustaining mid-Pliocene warmth?
To investigate the physical mechanisms through which the mid-Pliocene could warm compared to the PI, we have adopted a widely employed one-dimensional EBM formulation (Heinemann et al., 2009;Hill et al., 2014;Lunt et al., 2012a).Through this we have been able to disaggregate the zonal-mean temperature anomaly into contributions from changes to cloud albedo and surface albedo (together, planetary albedo), cloud emissivity and greenhouse gas emissivity (together, planetary emissivity), and changes to the meridional heat transport.
The principal mechanism of warming, responsible for 63 % or almost two-thirds of the total warming, is found to be the reduction in the planetary emissivity.This is mostly due to the reduction in the atmosphere's emissivity from an increase in the concentration of the atmospheric greenhouse gases CO 2 and water vapor.The effect of cloud emissivity is found to be negligible.
Changes to the planetary albedo accounts for the remaining one-third (37 %) of the total warming and it is found to result from a balance between the cooling effects of the clouds due to an increase in cloud albedo and warming due to a decrease in the surface albedo.Both these effects attain their largest values at high latitudes where the delicate balance between these opposing effects is pertinent to determining the magnitude of total warming.It is therefore clear that the ability of our climate models to characterize the impacts of clouds will be a large determiner of our ability to simulate mid-Pliocene warmth at high latitudes.
In general, our findings regarding the dominant mechanisms of warming and the spatial features of the various mechanisms of warming are consistent with the outcomes of the previous phase of PlioMIP as reported in Hill et al. (2014) and with Feng et al. (2017); however, the amplitudes of the mechanisms are now greatly amplified in the presence of the revised boundary conditions.Finally, we note that the 2/3-1/3 partition of the contribution to warming between the decrease in planetary emissivity and the decrease in planetary albedo was also found by Heinemann et al. (2009) for the case of another hothouse in the recent Earth's past: the late Paleocene-early Eocene.We have formulated a simple extension to the EBM model with the theoretical formulations of the factorization method to enable us to further evaluate the influence of each of the boundary condition differences on the effectiveness of the mechanisms of warming mentioned above.It is found that changes to pCO 2 , orography, and ice sheets respectively contribute 54, 34, and 12 % to the warming from increased greenhouse gases: CO 2 and water vapor.While the pCO 2 is prescribed within our simulations, the water vapor content is free to adjust to the atmosphere's carrying capacity according to the Clausius-Clapeyron relationship and is found to be a significant contributor to the warming.In fact, the total precipitable water in the mid-Pliocene is ∼ 22 % higher than the PI.While some of the increased water vapor is the result of the feedback from CO 2 -induced warming, a significant fraction of water vapor originates as feedback due to orography and ice sheet changes.Planetary albedo change, the second-to-leading mechanism of warming, is found to be dominated by changes to orography that contribute 54 % to the albedo-change-induced warming.Orography changes are also the highest-amplitude changes and although these dominate at the high latitudes of both hemispheres there is also notable warming in the midlatitudes of the NH owing to changes to the orography of the Tibetan Plateau and changes to vegetation in sub-Saharan Africa.The contributions from higher pCO 2 and ice sheet changes to the planetary-albedo-change-induced warming are found to be 33 and 13 %, respectively.The impacts of CO 2 and ice sheet changes on the surface-albedoinduced warming are found to be confined to the high latitudes.
Changes to the MHT do not change the planetary energy balance, but they do affect the local energy balance.The large-scale pattern of MHT change is characterized by a persistent reduction in MHT in the NH and an increase in the MHT in the midlatitudes of the SH.We find that the modi-fications to the orography almost entirely explain the MHT reduction in the NH, whereas the increase in MHT in the SH is primarily due to the elevated pCO 2 and supported in small part by the increase in MHT due to orography.The orography-driven reduction in the NH MHT is the result of a reduction in the OHT in the tropics and a reduction in the AHT in the midlatitudes.The impact of increased pCO 2 is only notable in the SH tropics and midlatitudes, where the increased pCO 2 is found to lead to an increase in the OHT.Its impact on the AHT is found to be entirely negligible.

Mid-Pliocene inference of climate sensitivity
We have attempted to infer the sensitivities of the climate system on various timescales through our suite of simulations.Our mid-Pliocene-derived inference of the CS (applicable on timescales of ∼ 100 years) of the UofT-CCSM4 model is S a = 3.25 • K per doubling of pCO 2 (Table 4).This is essentially indistinguishable from the estimate S a = 3.2 • K per doubling of pCO 2 that has been reported by Bitz et al. (2012) for the closely related NCAR-CCSM4 model.The accord between these two inferences of CS using very different lines of analysis, and despite the dependence of feedback mechanisms on the temperature (Colman and McAvaney, 2009), shows that the CCSM4 CS of 3.2 • K per doubling of pCO 2 is a very robust result.The difference between the UofT-CCSM4 and the NCAR-CCSM4 model with regards to parameterizations in the ocean model (Chandan and Peltier, 2017) obviously has no consequence for the short-timescale response of the climate as the ocean responds on a much longer timescale.
Following the initial, fast-timescale response of the climate to an applied radiative perturbation (which is characterized by the CS), the long-term evolution of the climate is governed by the slow feedbacks that start "coming online" at various timescales.On the timescale of a few hundred years after the initial radiative forcing is applied, which we refer to here as intermediate timescale, the albedo-driven positive feedback from the melting of the polar ice sheets (Pollard and DeConto, 2016) becomes important.We infer the climate system's sensitivity including this feedback, S i , to be 4.16 • K per doubling of pCO 2 , which is 30 % larger than the estimate for the CS.On a still longer timescale, the surface of the solid Earth begins to respond to the removal of the ice load and starts interacting with the climate through changes to the orography and the bathymetry.The climate system's long-term sensitivity, which includes this GIA response, is characterized by its ESS, which we have previously inferred in Chandan and Peltier (2017) to be ∼ 7 • K per doubling of pCO 2 .Finally, once the slow feedbacks that influence the ESS have stabilized, the sensitivity of the system drops to 3.35 • K per doubling of CO 2 .
Although the CS is the focus of much study, we wish to note the particularly large inference for S i , which one should be cognizant of while planning long-term climate mitiga-tion strategies.The significance of S i is best illustrated by comparing the potential warming implied by this parameter to predicted levels of warming at year 2300 CE obtained through modeling efforts employing the Extended Concentration Pathway (ECP; Meinshausen et al., 2011) scenarios.These are future emission scenarios developed as simple extensions of the Representative Concentration Pathway (RCP; van Vuuren et al., 2011) scenarios for years 2101-2300 CE for the purposes of assessing the response of the climate to anthropogenic emissions beyond the end of this century (Collins et al., 2013).Among the several extension scenarios, ECP4.5, in which the radiative forcing stabilizes at 4.5 W m −2 after year 2150 CE (although the forcing has already begun fluctuating about this value at least 100 years earlier), is of particular interest to us since in this scenario the pCO 2 stabilizes to a value roughly double (∼ 560 ppmv) that during the PI.Although it is abundantly clear by now that under no situation will the future emissions follow this optimistic scenario, this is nevertheless a convenient case for illustrating the impact of feedbacks.The radiative forcing in this scenario is a little over the 3.7 W m −2 that we would expect from a doubling of pCO 2 due to the presence of other greenhouse gases.
Intercomparison work performed for the IPCC AR5 (Collins et al., 2013;Zickfeld et al., 2013) using Earth system models of intermediate complexity (EMICs) has suggested that under the ECP4.5 scenario we are committed to a warming of 3.11 ± 0.60 • C by 2300 CE (compared to the 1850-1900 reference level).Even longer integrations out to year 3000 CE, performed with EMICs under the assumption that radiative forcing remains stationary at year 2300 CE values, suggest that the warming remains approximately the same.In contrast, our inference of S i would suggest that even if radiative forcings stabilized today at a pCO 2 level that is nearly identical to that prescribed in our mid-Pliocene experiment and ∼ 160 ppmv less than in the ECP4.5 scenario, then the warming by 2300 CE could still be ∼ 2.1 • C greater than the PI.Alternatively, in the case of emissions stabilizing to ∼ 560 ppmv in year 2150 CE as in the ECP4.5 scenario, our inference of S i suggests that we would be firmly committed to a warming of > 4 • C.This implies that the EMIC results, which do not take into account the intermediate-timescale feedbacks, could be underestimating the potential for warming on this timescale by ∼ 1 • C.
Similarly, the potential for long-term future warming suggested by our high inference of the ESS can be contextualized by comparing to results of even longer EMIC simulations.The utility of such a comparison is significantly less than that presented above for S i owing to the extremely unconstrained task of forecasting an emission scenario for the timescales over which ESS is meaningful.Zickfeld et al. (2013) considered an even longer extension of the ECP4.5 scenario in which the 2300 CE radiative forcing is kept stationary until 3000 CE after which the pCO 2 is allowed to evolve freely.They found that even out to year 4000 CE the EMICs do not predict any increase in surface temperature beyond the 2300 CE estimate.Our mid-Pliocene-informed estimate for the ESS would suggest, however, that there is potential for a greater magnitude of warming 2000 years into the future than suggested by EMICs.

Figure 1 .
Figure 1.The topography for the Ei, Eo, and Eoi experiment configurations (a) and the orographic anomaly of those configurations from PlioMIP2 modern (b).The land-based ice sheets in the top row are displayed on the purple-colored scale that is shown separately.
shortwave upwelling flux at the TOA, longwave upwelling flux at the TOA, and longwave upwelling flux at the surface of the Earth, respectively.The divergence of the meridional heat transport is similarly diagnosed as the sum of the net shortwave and the net longwave fluxes at the TOA:

Figure 3 .
Figure 3. Evolution of the (a) 2 m surface air temperature and (b) the top of the atmosphere energy imbalance for the sensitivity experiments.

Figure 4 .
Figure 4. Hovmöller diagrams for the evolution of ocean temperatures as a function of depth in the five sensitivity simulations.

Figure 5 .
Figure 5. Climatological MASAT anomalies of the Ei and the Eo sensitivity experiments with respect to controls.The control experiment is chosen such that the pCO 2 is identical between the sensitivity experiment and the control.The left column shows the anomalies for the 280 ppmv pCO 2 experiments and the right column shows the anomalies for the 400 ppmv pCO 2 experiments.The 1st and the 99th percentiles of the anomalies are shown above each panel.The climatology years are noted in Table2.

Figure 8 .
Figure 8. Zonal mean of the MASAT anomalies shown in Fig. 7.The abscissa is linear in the sine of latitude.
Figure9.Factorization of the zonal-mean mid-Pliocene minus PI temperature anomaly computed using the 1-D EBM (dashed gray) into contributions arising from changes to cloud albedo (red), cloud emissivity (purple), surface albedo (blue), greenhouse gas concentrations (including water vapor; green), and the total meridional heat transport (yellow).The solid gray line is the mid-Pliocene minus PI temperature anomaly computed directly from the GCM output.The dashed black line is the sum of the aforementioned individual contributions and is seen to almost perfectly overlap the total anomaly computed from both the EBM and the GCM.The abscissa is linear in the sine of latitude.

Figure 10 .
Figure 10.Factorization of the PI to mid-Pliocene change in zonal-mean temperature (red curve) due to change in (a) planetary albedo, (b) clear-sky or surface albedo, (c) cloud albedo, (d) planetary emissivity, (e) clear-sky emissivity or emissivity change due to change in greenhouse gas concentrations including water vapor, (f) cloud emissivity, and (g) meridional heat transport into contributions arising from changes made to the following boundary conditions: orography (purple), ice sheets (green), and atmospheric pCO 2 (blue).The abscissa is linear in the sine of latitude.

Figure 11 .
Figure 11.Schematic representations of various measures of climate sensitivity.Forcings are denoted in green, whereas the slow internal feedbacks are shown in orange.See main text for details.

Figure 12 .
Figure 12.Climatological sea ice concentration anomalies of the five sensitivity experiments and the 400 ppmv mid-Pliocene experiment with respect to controls.The control experiment is chosen such that the pCO 2 is identical between the experiment and the control.(a) The anomalies for 280 ppmv pCO 2 experiments and (b) the anomalies for 400 ppmv pCO 2 experiments.The climatology years are noted in Table2.

Figure 13 .
Figure 13.Various comparisons of the climatological total meridional heat transport and the individual transports by the atmosphere and the ocean.The climatology years are noted in Table2.

Figure 14 .
Figure 14.Various comparisons of the climatological total meridional heat transport by the ocean and the individual transports in the Atlantic and the Indo-Pacific basins.The total oceanic heat transport curves are the same curves that are labeled "Ocean" and shown in Fig. 13.The climatology years are noted in Table2.

Figure 15 .
Figure 15.Meridional heat transport anomalies (mid-Pliocene minus PI) and their factorization.(a) Total meridional heat transport, (b) atmospheric heat transport, and (c) oceanic heat transport.The dashed black curve in each figure is the sum of the factorized components.
5.1.3How do the modifications to the boundaryconditions themselves change the effectiveness of these physical mechanisms?

Table 1 .
Boundary conditions for the control, Pliocene, and sensitivity experiments.The orbital parameters common to all simulations are the solar constant S 0 = 1365 W m −2 , eccentricity e = 0.016724, obliquity = 23.446• , and the longitude of perihelion ω = 102.04• .
a Includes orography, bathymetry, LSM, and river routing.b Includes vegetation, soil, and lakes.

Eoi400 Eoi450 cesmpifv1mts E280 E400 3501 Eo400 Eoi280 Eo280 Ei400 Ei280 4795 4569 1000 Figure 2. Relationship
among our PlioMIP2 simulations.The simulations colored green are previously introduced in Chandan and Peltier (2017), while the blue-colored simulations are new to this paper.The numbers at the locations of branching indicate the model year from which one or more new simulations were branched.All simulations without an ancestor have their ocean components initialized from modern-day conditions

Table 4 .
Various measures of climate sensitivity calculated from our experiments.

Table 5 .
Estimates of the contributions to the mid-Pliocene minus PI temperature difference from various mechanisms.All values are in units of • C. The last row in the column indicating the total contribution of a mechanism is calculated from the decomposition shown in Fig.9and therefore may differ slightly from the sum of the contribution from ice sheets, orography, and CO 2 .The total row is not the sum of the rows above it.bdT swc + dT salb + dT lwc + dT gg + dT H ity. The reduction in emissivity impedes the mid-Pliocene atmosphere's ability to radiate energy to space and thereby allows it to retain a greater proportion of the energy introduced by the incident shortwave radiation.Planetaryemissivity-induced warming, dT emiss , is found to contribute up to 2.43 • C, or roughly 63 % of the total PI to mid-Pliocene temperature anomaly. a