The influence of ice sheets on the climate during the past 38 million years

Since the inception of the Antarctic ice sheet at the Eocene-Oligocene Transition (∼34 Myr ago), land ice has played a crucial role in Earth’s climate. Through the ice-albedo and surface-height-temperature feedbacks, land ice variability strengthens atmospheric temperature changes induced by orbital and CO2 variations. Quantification of these feedbacks on long time scales has hitherto scarcely been undertaken. In this study, we use a zonally averaged energy balance climate model bi-directionally coupled to a one-dimensional ice sheet model. The relative simplicity of these models allows us to perform 5 integrations over the past 38 Myr in a fully transient fashion, using a benthic oxygen isotope record as forcing to inversely simulate CO2. Output of the model are mutually consistent records of CO2, temperature, ice volume-equivalent sea level and benthic δO. Firstly, we investigate the relation between global temperature and CO2, which changes once the model run has experienced high CO2 concentrations. Secondly, we study the influence of ice sheets on the evolution of global temperature and polar amplification by comparing runs with ice sheet-climate interaction switched on and off. We find that ice volume 10 variability has a strong enhancing effect on atmospheric temperature changes, particularly in the regions where the ice sheets are located. As a result, polar amplification in the Northern Hemisphere decreases towards warmer climates as there is little land ice left to melt. Conversely, decay of the Antarctic ice sheet increases polar amplification in the Southern Hemisphere in the high-CO2 regime. Our results also show that in cooler climates than the pre-industrial, the ice-albedo feedback predominates the surface-height-temperature feedback, while in warmer climates they are more equal in strength. 15


Introduction
The most abundant information source with the highest resolution on Cenozoic global climate change are stacked benthic oxygen isotope (δ 18 O) records (Lisiecki and Raymo, 2005;Zachos et al., 2008;Cramer et al., 2009), which have been well-studied and statistically analysed (e.g.Mudelsee et al., 2014).The benthic δ 18 O signal is known to be comprised of two factors (e.g.Chappell and Shackleton, 1986): 1) the deep-sea temperature, and 2) the volume of land ice on Earth.An additional independent record of either one is therefore required to separate the signal into its individual constituents.Deep-sea temperature records can be reconstructed based on the Mg/Ca proxy (Lear et al., 2000;Sosdian and Rosenthal, 2009;Elderfield et al., 1 Clim. Past Discuss., doi:10.5194/cp-2016-109, 2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.2012), but a global average deep-sea temperature is hard to obtain.Sea level records are also available, but are subject to the same problem of inferring a global mean (Miller et al., 2005;Kominz et al., 2008;Rohling et al., 2014).Studies using sea level records face the additional challenge of converting local sea level to ice volume, which is not straightforward mainly because of dynamic topography (Mitrovica and Milne, 2003;Kendall et al., 2005).Alternatively, calculation of benthic δ 18 O can be incorporated in coupled ice sheet-climate models, by using parameterisations of the contribution of deep-sea temperature (Duplessy et al., 2002), and the isotopic content of ice sheets (Cuffey, 2000).Hitherto, studies using this approach have mostly focused on relatively short time intervals surrounding important climatic events, such as the Eocene-Oligocene Transition (33.9 Myr ago; Tigchelaar et al. (2011); Ladant et al. (2014); Wilson et al. (2013)), the Mid-Miocene Climatic Optimum (13.9 Myr ago; Langebroek et al. (2010); Gasson et al. (2016)), and the Pliocene-Pleistocene Transition (2.6 Myr ago, Willeit et al. (2015)), using models of varying complexity.Computer power remains a limiting factor for these coupled ice sheet-climate studies.
Nevertheless, a study using a model of reduced complexity simulated the past 3 million years (Berger et al., 1999).However, they did not include the Southern Hemisphere and the Antarctic ice sheet in their model.Here, we will use a coupling between an energy balance global climate model and a one-dimensional ice sheet model of all major ice sheets, to perform transient simulations over the past 38 Myr.
Our model approach builds on the inverse routine to derive atmospheric temperature from benthic δ 18 O, that was introduced by Oerlemans (2004).This methodology was consequently developed further to force stand-alone three-dimensional ice sheet models over the past 1 Myr (Bintanja et al., 2005;De Boer et al., 2013), the past 3 Myr (Bintanja and Van de Wal, 2008), and the past 5 Myr (De Boer et al., 2014).With this inverse routine, the past 40 Myr were simulated (De Boer et al., 2010) and further analysed (De Boer et al., 2012), using a one-dimensional ice sheet model that calculates all land ice on Earth.In Stap et al. (2014), this ice sheet model was coupled to a zonally averaged energy balance climate model (Bintanja, 1997), and run over the past 800 kyr forced by the EPICA Dome C ice-core CO 2 record (EPICA community members, 2004).Inclusion of a climate model added CO 2 to global temperature and sea level as an integrated component of the simulated system.In addition, it rendered the possibility of investigating ice sheet-climate interactions, specifically the ice-albedo and surface-height-temperature feedbacks.Furthermore, instead of annual mean and globally uniform temperature perturbations to present-day climate, seasonal meridional temperature distributions were used to force the different ice sheets.In a subsequent study, the inverse routine was transformed to yield CO 2 concentrations using the benthic δ 18 O as input, making CO 2 a prognostic variable (Stap et al., 2016a).The resulting values were used to force the coupled model over the past 5 Myr.In Stap et al. (2016b), the model was run over the period 38 to 10 Myr ago, and the influence of Antarctic topographic changes on the simulated CO 2 was investigated.
Here, we will first explore a hysteresis that occurs in our coupled model, which necessitates us to reconfigure our setup.
Using a new reference simulation, we will attempt to quantify the influence of ice sheets on climate change during the past 38 million years, specifically on global temperature perturbations and polar amplification.To this end, we compare runs of the climate model where ice sheet-climate interaction is switched off to the new reference run.We find that ice sheets modify the Earth System Sensitivity more strongly in cooler climates.Furthermore, they have a large regional impact, which causes CO 2induced climate change to be heavily dependent on latitude.Moreover, polar amplification is determined by the background climate state.

Methodology
We use the same simplified coupled ice sheet-climate model setup as Stap et al. (2016b) (Fig. 1).The climate component is a zonally averaged energy balance climate model (based on North, 1975;Bintanja, 1997) (1) Here, α is surface albedo, and Q local radiation obtained from Laskar et al. (2004).Ice-sheet dependent tuning factors C abl determine the threshold for which ablation starts (listed in Stap et al., 2014).The ice sheet model (De Boer et al., 2010) calculates the surface height change and ice sheet extent of five ice sheets on Earth: Greenland, North America, Eurasia, East Antarctica and West Antarctica.This information is used to update the land ice fraction and surface height profile in the climate model for the next time step.Exchange of variables takes place every 500 model years.The coupled model is forced with insolation data (Laskar et al., 2004) and an inverse routine, which yields CO 2 concentrations from the difference between the modeled benthic δ 18 O value and an observed value an output-timestep later (Stap et al., 2016a).The radiative forcing anomaly with respect to present day is multiplied by a factor 1.3 to account for non-CO 2 greenhouse gases.The result of the model consists of mutually consistent records of benthic δ 18 O, atmospheric CO 2 , temperature, and ice-volume equivalent sea level.
First, we extend the run over the period 38 to 10 Myr ago forced with the stacked benthic δ 18 O of Zachos et al. (2008) described in Stap et al. (2016b), to include the past 10 Myr (this run is called '38 Myr').We compare the results of the past 5 Myr to the reference run of Stap et al. (2016a), that simulated the past 5 Myr using the record of Lisiecki and Raymo (2005) as forcing (this run is called '5 Myr').The 38 Myr run shows CO 2 values that are much lower over the past 5 Myr than the results of the old reference run, and therefore also far below the EPICA Dome C record (EPICA community members, 2004).
To regain agreement with the EPICA Dome C record, we define a new reference run in this study (new REF).In this new reference run we increase the cloud optical thickness parameter τ cl from 3.11 to 3.41.We opt to alter this parameter because it was already used as a tuning parameter in the original climate model (Bintanja, 1997), and in the ice sheet-climate model coupling (Stap et al., 2014).Increasing τ cl will lower the temperatures calculated by the climate model, such that for the same benthic δ 18 O higher simulated CO 2 levels are obtained, in better agreement with the EPICA Dome C record.However, this will also raise the threshold CO 2 level for the inception of the East Antarctic ice sheet (EAIS).By increasing the ablation threshold parameter C abl in the insolation-temperature-melt calculation (Eq.( 1)) of the EAIS (from −30 to −10), this ice sheet glaciates at lower temperatures and therefore at lower CO 2 concentrations.This compensates the unintended CO 2 threshold increase.In this run, the Zachos et al. (2008) record is once more used to force the model over the past 38 Myr after a 2 Myr spin-up, similarly as in Stap et al. (2016b).Furthermore, to study the effect of ice sheet-climate interactions, we perform three to 0.69 K (Fig. 3b, blue and cyan lines).Keeping the ocean overturning strength fixed at PD also leads to a small reduction; the difference becomes 0.73 K (Fig. 3b, dark green and light green lines).The combined effect of uncoupled ice and ocean overturning strength is still not sufficient to eliminate the hysteresis (Fig. 3b, red and orange lines).Even when in addition sea ice and snow cover are kept constant, a small hysteresis is present (not shown).This means that the hysteresis is inherent to the core of the climate model: the parameterisation of vertical and horizontal energy transfer in the ocean and atmosphere.The factors mentioned above act to enhance this hysteresis.
Over the past 800 kyr, the 5-Myr run was calibrated to the EPICA Dome C ice core proxy-CO 2 record (EPICA community members, 2004) and shows negligible bias (-3.9 ppm), whereas with the same parameter setting the 38-Myr run shows much lower values than this proxy record with a -47.7 ppm bias (Fig. 2a, mind that here we show 1-kyr values instead of 40-kyr averages).2c).In the time between these events, CO 2 in the new reference run is modestly higher, up to 100 ppm.This is because the deep-sea temperatures are lower at the same CO 2 , and therefore contribute less to the δ 18 O anomaly with respect to present day.Compensating for the lower deep-sea temperatures, higher CO 2 increases the δ 18 O anomaly, by increasing both deep-sea temperature and the contribution of ice volume, hence raising sea levels (not shown).After the MMCO, when the EAIS has stabilised to near-PD size, the new reference simulation shows higher CO 2 values.Over the past million years, the new reference simulation (Fig. 2a, red line) agrees much better with the 5-Myr run (Fig. 2a, green line) and with the EPICA proxy-CO 2 record (Fig. 2a, cyan line); the bias with respect to this proxy-record is reduced to 13.6 ppm.
Even after re-calibration, the simulated new reference CO 2 remains lower than in the 5-Myr run during the Pliocene and early   4, orange dots).The warmest anomaly is only increased by 21 % (factor 1.21) when ice is coupled, by 9 % when only albedo is coupled, and by 3 % when only surface height is coupled.This means the surface-height-temperature feedback becomes relatively more important in warmer climates.Since decreased ESS at higher CO 2 is still present in the run where ice volume is kept at PD level (Fig. 4, blue dots), reduced ice volume variability is not its only determining factor.
The influence of ice sheets on the climate is strongest in the region where they are situated, leading to increased polar amplification.This is demonstrated by the relation between global temperature and Northern Hemispheric (40 to 80 • N, Fig. 5a), or Antarctic temperature (60 to 90 • S, Fig. 5b).In the Northern Hemisphere, the minimum local temperature with respect to PI is -2.0 K in the uncoupled case, and -9.5 K in the run with fully coupled land ice.When only the albedo or surface height changes are coupled, the Northern Hemispheric temperature anomaly reaches -6.4 K and -2.8 K low points respectively.Conversely, the amount of land ice lost in warmer climates is relatively small, as only the Greenland ice sheet (∼ 7 m.s.l.e.) is left to melt.
Consequently, the Northern Hemispheric temperature is then not affected much by not including land ice changes.The remaining polar amplification in the Northern Hemisphere is hence mostly caused by other factors, such as sea ice and snow cover variability.In the Southern Hemisphere, the lowest temperature is similar for the coupled and uncoupled simulations, although it is achieved at a higher global temperature in the uncoupled case.These Southern Hemispheric temperatures are similarly low because the Antarctic ice sheet grows relatively little in size towards colder-than-PI conditions (see also Stap et al., 2014).
When Antarctica is allowed to melt in warm climates, however, the local temperature increase becomes much stronger: 11.6 instead of 5.9 K with respect to PI.In these conditions, we find that coupling albedo changes leads to a maximum Antarctic temperature anomaly of only 7.0 K (Fig. 5b, black dots).When only surface height changes are coupled, this anomaly reaches 7.4 K (Fig. 5b, orange dots).This result implies that albedo changes are relatively less important in Antarctica than in the high latitudes of the Northern Hemisphere.The reason is that the Antarctic continent remains snow covered throughout most of the year when the land ice retreats, which reduces the albedo change (see also Stap et al., 2016a).Since temperature changes are strongest in the Southern Hemisphere in warmer-than-PI climates, this explains the increased relative importance of the surface-height-temperature feedback on ESS in these climates.
6 Finally, we compare the relation between global temperature and logarithmic CO 2 in three model runs with uncoupled ice (Fig. 6).In one run the ice sheets are kept at PD condition as before (now called PD ice, blue dots), in another one we use the LGM condition (LGM ice, black dots), and in the last one all ice is removed (no ice, red dots).Naturally, the more ice is present on Earth, the colder the climate becomes, so the LGM ice run is colder than the PD ice run, which in turn is colder than the no ice run.The difference between the PD ice and the no ice run is fairly uniform over the whole CO 2 range.The difference between the LGM ice and the no ice run, however, is larger in cold climates than in warm climates as it shrinks from ∼2.8 to ∼1.6 K.This is explained by the extra land ice in the LGM ice run cooling the climate and increasing the area on Earth covered by snow and sea ice.As a result of this area increase, the land surface has a higher albedo, which cools the climate further.In cold climates this effect is stronger because the snow-and sea ice-covered area grows more towards the equator, where there is more incoming solar radiation.Consequently, the albedo increase is more effective as it leads to absorption of more energy, and thus to a stronger temperature decrease.

Discussion and conclusions
We  et al., 2008).This allows us to simulate periods further back in time, for which CO 2 data are highly uncertain (Beerling and Royer, 2011), and when the Antarctic ice sheet was very dynamic.This constitutes a wider simulated range of climates.For a comparison between the results of our model and the LLN-2D model, the reader is referred to Stap et al. (2014).Our coupled model also represents an improvement upon the model setup of De Boer et al. ( 2010), who used the same ice-sheet model in stand-alone form to simulate the past 40 Myr (De Boer et al., 2010, 2012).The inclusion of a climate model enables us to simulate, and force the different ice sheets with, seasonal meridional temperature distributions instead of globally uniform perturbations to present-day climate with a fixed seasonal cycle.Nonetheless, we recognise that our coupled ice sheet-model is relatively simple.However, more sophisticated models, such as GCMs coupled to three-dimensional thermodynamic ice models, are as of yet not suitable to perform multi-million year integrations, because of limited computer power.Our results should therefore be seen as a first step in the direction of simulating these time spans with coupled ice sheet-climate models, thereby identifying interesting phenomena and potential obstacles.When results of these more sophisticated models are achieved, they can be compared to ours to see which features appear in the full hierarchy of models and which are specific to more comprehensive models including more physics.
In our model, the results for CO 2 concentrations lower than roughly 450 ppm depend on the transient evolution of CO 2 .
When during the run the model has previously experienced high CO 2 values, temperatures are higher than when this is not the case.This hysteresis is persistent even in runs without any change in albedo due to snow-, sea ice-or permanent land ice-cover  (2016a).It remains debatable which simulation is the most veracious over the past 5 Myr.On the one hand, the long simulation carries a longer memory, which would be closer to the state of the actual climate system.On the other hand, it is uncertain how accurately our climate model simulates very warm climates; the climate model is designed and tested for PD and LGM climates (Bintanja and Oerlemans, 1996).This argument favours the shorter 5-Myr run as the more trustworthy result.Regardless, it remains unknown whether the hysteresis is an artefact of our model, or is also exhibited by other models.We therefore suggest that in the future, climate models should be tested for this behaviour by confronting them with high CO 2 values before simulating cooler climates. of Gasson et al. (2012).A sigmoidal shape is also apparent in the modeled relation between logarithmic CO 2 and sea level, coherent with the results from Foster and Rohling (2013).The threshold CO 2 value for glaciation of the Southern Hemisphere is distinctly higher than for the Northern Hemisphere in our model, which was also found by DeConto et al. (2008).In our model, this is a consequence of lower forcing temperatures, as the Antarctic ice sheet initiates at higher latitudes.
When ice sheets are kept at PD level, the relation between logarithmic CO 2 and global temperature shows a declining slope.
This may be compared to the hybrid data-model results for climate sensitivity of Köhler et al. (2015), as well as to the modeled climate sensitivity of Friedrich (2015).Köhler et al. (2015) investigated the relation between the radiative forcing of proxy-data CO 2 (EPICA community members, 2004;Hönisch et al., 2009) -which is linearly related to logarithmic CO 2 (Myhre et al., 1998) -and modeled global temperature (De Boer et al., 2014, scaled).Friedrich (2015) forced the intermediate complexity climate model LOVECLIM over the past 800 kyr using the EPICA Dome C ice-core record (EPICA community members, 2004) and an Northern Hemispheric ice sheet reconstruction (Ganopolski and Calov, 2011).The resulting climate sensitivity of these studies is opposite to ours, as they show increased climate sensitivity at higher CO 2 concentrations.These studies, however, consider a smaller range of CO 2 .Furthermore, their approach is different, as they do take into account ice volume variations, but compensate for their effect by adding their radiative forcing to the forcing induced by CO 2 variations (see PALAEOSENS Project Members, 2012).Implicitly assumed in their approach is that these radiative forcings have the same effect on temperature, which may not generally be the case (Yoshimori et al., 2011).Our findings are in agreement with Ritz et al. (2011), who used a two-dimensional energy balance climate model that showed an increase of climate sensitivity from 3.0 K per CO 2 doubling at PI conditions to 4.3 K at LGM conditions.Ice volume changes enhance the modeled effect of CO 2 on temperature, via the ice-albedo and the surface-height-temperature feedbacks.As the growth and decay of the Antarctic ice sheet is spread over a larger range of logarithmic CO 2 , this effect is stronger on the low CO 2 branch than on the high CO 2 branch.This leads to an increased curve in the ln(CO 2 )-temperature relation when land ice is allowed to vary.We do not find Earth System Sensitivity (ESS) to be constant, opposite to the implicit assumption in Van de Wal et al. (2011).In fact, in our model ESS is stronger at lower CO 2 , similar to the findings of Hansen et al. (2013), who performed CO 2 doubling experiments using the simplified atmosphere-ocean model of Russell et al. (1995).
However, their ESS is consistently higher than what our model calculates.They eventually also find increased sensitivity again at high CO 2 (2480 to 9920 ppm), which is outside the range we simulate during our time span.
Our results further clearly show a non-linear relation between temperatures on both hemispheres, caused by asymmetric ice sheet-climate interaction.In cooler climates, the Northern Hemispheric ice sheets change in size, causing large fluctuations in the temperature on this hemisphere.The ice-albedo feedback is much stronger than the surface-height-temperature feedback in these conditions (see also Stap et al., 2014).This is reflected in the Northern Hemispheric (40 to 80 • N), Antarctic (60 to 90 • S), and global temperature profiles.Oppositely, in warmer climates the Antarctic ice sheet is more dynamic, so that temperature changes more strongly on the Southern Hemisphere.Now, the surface-height-temperature feedback becomes relatively more important.When the ice sheets are kept constant, the temperature perturbations are more uniformly distributed over the globe.
The different response of the northern and southern high latitudes to CO 2 changes challenges the approach of De Boer et al.
(2010) and De Boer et al. ( 2012), who reconstructed a single high-latitude temperature anomaly.Furthermore, their record cannot readily be translated to global conditions by a constant factor (as is done in e.g Martínez-Botí et al., 2015), because the conversion depends on the prevailing climate state.This also holds for the interpretation of proxy data from high northern and southern latitudes, both marine and ice cores.
Future work on polar amplification and ESS on long timescales may involve using more sophisticated models.When results of these models become available, a comparison to this work can be made, to see what the added complexity invokes in the simulated system.Alternatively, the simplified model used here can be extended with more aspects of the climate system, moving towards a full Earth System Model.Important aspects hitherto neglected in our model are the effects of dynamic vegetation (e.g.Knorr et al., 2011;Liakka et al., 2014), dynamic topography (Mitrovica and Milne, 2003;Kendall et al., 2005), and land surface changes (e.g.tectonics, sedimentation) (Stap et al., 2016b;Wilson et al., 2012;Gasson et al., 2015).Ultimately, this model could also be coupled to a carbon cycle model, e.g.BICYCLE (Köhler and Fischer, 2004), in order to simulate climate using only insolation data as input.
with 5 • latitudinal resolution, including a simple ocean model mimicking meridional seawater circulation with varying strength based on the density difference between the polar and equatorial waters.The climate model provides the temperature input (T ) to the mass balance module of a onedimensional ice sheet model that uses the Shallow Ice Approximation.The ablation (M ) in the mass balance parameterisation is calculated by an insolation-temperature melt equation: model tests over the past 38 Myr, where we use the same CO 2 input as the new reference run, but keep the ice sheet extent and Clim.Past Discuss., doi:10.5194/cp-2016-109,2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.surfaceheight constant at present-day (PD ice/ice uncoupled) or Last Glacial Maximum (LGM ice) level, or remove all ice in the climate model completely (no ice).In two last tests, we keep the surface height constant at present-day (PD) level, but still vary ice sheet extent (height uncoupled) or oppositely keep the ice sheet extent at PD level and only vary the surface height (albedo uncoupled).This allows us to distinguish the effects of the surface-height-temperature feedback and the ice-albedo feedback.We show 40-kyr running averages of all variables, because shorter scale variability is too strong in our model during the period 38 to 10 Myr ago (seeStap et al., 2016b).To further explore the difference we find between the results of the 5 Myr and the 38 Myr run, we additionally conduct four pairs of experiments with the model in forward mode.In forward mode, we do not use the inverse routine, but force the model by a-priori designed CO 2 scenarios.We run the model fully coupled, or alternatively by keeping ice sheets (PD ice), or ocean overturning strength (PD OT), or both (PD ice + OT) constant at their present-day levels.Starting from a CO 2 concentration of 450 ppm, we force the model by changing the CO 2 input in steps of 50 ppm every 50 kyr.In one set of experiments (named 'up'), the CO 2 is first raised from 450 ppm to 1200 ppm, then lowered to 150 ppm, and increased again to 600 ppm.In the other set (named 'down'), the CO 2 is initially dropped from 450 to 150 ppm, then raised to 1200 ppm, and ultimately decreased again to 300 ppm.Insolation is kept at PD level in all these equilibrium experiments, and the parameters of the coupled model are set to their old reference values (C abl = −30 for the EAIS, and τ cl = 3.11)(Stap et al., 2014(Stap et al., , 2016b)).3Results: hysteresis and new referenceWhen we compare the final 5 Myr of our 38-Myr simulation to our 5-Myr simulation, we notice that the 38-Myr simulation shows much lower CO 2 concentrations (Fig.2b, green and blue lines).These contradicting results cannot be explained by the use of different forcing records -Zachos et al. (2008) for 38 Myr as opposed toLisiecki and Raymo (2005) for 5 Myr -as these show similar values during this time (0.02 ‰ average difference, not shown).Instead, this incongruence is caused by model hysteresis, to wit a dependency of the simulated global temperature on the evolution of CO 2 during the run as well as the prevailing CO 2 level.This becomes apparent when we study the relation between CO 2 and global temperature in the forward runs with stepwise changing CO 2 described in Sect.2, performed with the fully coupled model (Fig.3a).Starting at 450 ppm CO 2 , the 'up' and 'down' runs show the same initial global temperature.However, in the 'down' run, where the CO 2 progresses stepwise downward first and then upward, the global temperatures at low (< 450 ppm) CO 2 values are approximately a degree lower than those in the 'up' run, where CO 2 is first raised and then lowered.When the 'down' run is integrated over another CO 2 cycle, it shows the same global temperatures as the up run (not shown).This means that once the coupled model has experienced high CO 2 values during its run, as is the case in the 38-Myr run but not in the 5-Myr run, the climates at lower CO 2 are warmer.This has important consequences for the simulated CO 2 concentrations as they have to decline further to obtain similare temperature values, which is what happens in the transient 38-Myr simulation forced by the inverse routine.This behavior is a form of hysteresis as results depend on previous conditions of the model.The question now arises what the cause of this hysteresis is.The global temperature difference between the 'up' and 'down' run is 0.94 K at 150 ppm CO 2 .When the ice sheet model is uncoupled, and the climate model is directly forced by PD ice sheets, this reduces Clim.Past Discuss., doi:10.5194/cp-2016-109,2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.
Pleistocene (Fig.2b, green line), as a consequence of the hysteresis.Although it is more variable than the reconstruction based on a constant Earth System Sensitivity (ESS) byVan de Wal et al. (2011) (Fig. 2b, black line), the long-term means are now similar.4Results: ice sheet-climate interactionIn this section, we investigate the influence of ice sheet-climate interaction on polar amplification, and on the Earth System Sensitivity (ESS).The ESS is defined as the global temperature response to a radiative forcing caused by changing CO 2 , tak-Clim.Past Discuss., doi:10.5194/cp-2016-109,2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.inginto account all climate feedbacks(PALAEOSENS Project Members, 2012).This radiative forcing by CO 2 is proportional to the logarithmic change of CO 2(Myhre et al., 1998).In Fig.4(red dots), we therefore show the relation between global temperature anomalies from pre-industrial (PI) and the logarithm of CO 2 divided by a reference PI value of 280 ppm in our new reference run.Evidently, this relation is not constant, as in warm climates the global temperature increase for a given CO 2 increase is less strong.The slope of a least squares linear regression shows a value of 10.6 K for ln(CO 2 /CO Clim.Past Discuss., doi:10.5194/cp-2016-109,2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.
have presented mutually consistent transient simulations of CO 2 , temperature, ice volume and benthic δ 18 O over the past 38 million years.They were obtained using a coupling between a zonally averaged energy balance climate model and a onedimensional ice sheet model.A similar coupled ice sheet-climate model of reduced complexity was introduced by Gallée et al. (1992): the LLN-2D model.However, in contrast to the LLN-2D model we include the Southern Hemisphere and the Antarctic ice sheet.Furthermore, we use an inverse routine that yields atmospheric CO 2 from an observed benthic δ 18 O record (Zachos Clim.Past Discuss., doi:10.5194/cp-2016-109,2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.and without changes in ocean overturning strength.However, these factors do enhance it.The hysteresis leads to a mismatch in runs over the past 38 Myr between the CO 2 concentrations of the past 800 kyr and the EPICA Dome C proxy-record.Since this proxy record served as a main target for our model performance (see Stap et al., 2016a), we have re-calibrated the model configuration.We have achieved more satisfactory results by increasing the cloud optical thickness parameter, which decreases temperatures in the climate model.Therefore, simulated CO 2 concentrations are higher, in order to obtain the same benthic δ 18 O values.However, even after this re-calibration, simulated CO 2 values are still lower than in the 5-Myr run of Stap et al.
The relation between temperature and ice volume in our model can roughly be divided into three regimes (see also De Boer et al. (2010) and Van de Wal et al. (2011)): 1) at low CO 2 values, strong ice volume variability due to dynamic Northern Hemispheric ice sheets, 2) at intermediate CO 2 values, weaker variability, 3) at high CO 2 values, strong variability due to a dynamic Antarctic ice sheet.This constitutes a sigmoidal temperature-sea level relation, similar to the data-analysis results Clim.Past Discuss., doi:10.5194/cp-2016-109,2016 Manuscript under review for journal Clim.Past Published: 8 November 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Schematic overview of the coupling of the zonally averaged energy balance climate model and the one-dimensional ice sheet model.

Figure 4 .Figure 5 .
Figure 4. Relation between the logarithm of CO2 divided by the PI value of 280 ppm, and global temperature anomalies with respect to PI (of the reference run), for the reference simulation (red dots), the simulation with uncoupled ice (blue dots) and the simulation with only surface height (orange dots) or albedo (black dots) coupled.

Figure 6 .
Figure 6.Relation between the logarithm of CO2 divided by the PI value of 280 ppm, and global temperature anomalies with respect to PI (of the reference run), for the simulation with ice kept at PD level (blue dots), at LGM level (black dots) and the simulation with all land ice removed (red dots).
2,ref )values below 0 (CO 2 < 280 ppm: coldest climates), and 3.7 K for values above 0.69 (CO 2 > 560 ppm: warmest climates), a reduction of 65 %.In the ice uncoupled run, the slope reduces by only 46 % from 5.6 K to 3.0 K going from the coldest to the warmest ln(CO 2 /CO 2,ref ) regime (Fig.4, blue dots).In this case, the standard error of an linear regression through all data points is reduced by 58 % with respect to the fully coupled run, from 0.0050 K to 0.0021 K.The fact that the relation between ln(CO 2 /CO 2,ref ) and global temperature is better approximated by a linear fit when land ice is uncoupled means the log(CO 2 )-T relation is more linear.The coldest anomaly is amplified by 79 % (factor 1.79) if land ice changes are incorporated, by 50 % if only albedo is coupled (Fig.4, black dots), and by 4 % if only surface height is coupled (Fig.