On the state dependency of the equilibrium climate sensitivity during the last 5 million years

It is still an open question how equilibrium warming in response to increasing radiative forcing – the specific equilibrium climate sensitivity S – depends on background climate. We here present palaeodata-based evidence on the state dependency of S, by using CO2 proxy data together with a 3-D ice-sheet-model-based reconstruction of land ice albedo over the last 5 million years (Myr). We find that the land ice albedo forcing depends non-linearly on the background climate, while any non-linearity of CO2 radiative forcing depends on the CO2 data set used. This nonlinearity has not, so far, been accounted for in similar approaches due to previously more simplistic approximations, in which land ice albedo radiative forcing was a linear function of sea level change. The latitudinal dependency of icesheet area changes is important for the non-linearity between land ice albedo and sea level. In our set-up, in which the radiative forcing of CO2 and of the land ice albedo (LI) is combined, we find a state dependence in the calculated specific equilibrium climate sensitivity, S[CO2,LI], for most of the Pleistocene (last 2.1 Myr). During Pleistocene intermediate glaciated climates and interglacial periods, S[CO2,LI] is on average ∼ 45% larger than during Pleistocene full glacial conditions. In the Pliocene part of our analysis (2.6–5 MyrBP) the CO2 data uncertainties prevent a well-supported calculation for S[CO2,LI], but our analysis suggests that during times without a large land ice area in the Northern Hemisphere (e.g. before 2.82 MyrBP), the specific equilibrium climate sensitivity, S[CO2,LI], was smaller than during interglacials of the Pleistocene. We thus find support for a previously proposed state change in the climate system with the widespread appearance of northern hemispheric ice sheets. This study points for the first time to a so far overlooked non-linearity in the land ice albedo radiative forcing, which is important for similar palaeodata-based approaches to calculate climate sensitivity. However, the implications of this study for a suggested warming under CO2 doubling are not yet entirely clear since the details of necessary corrections for other slow feedbacks are not fully known and the uncertainties that exist in the ice-sheet simulations and global temperature reconstructions are large.


Introduction
One measure to describe the potential anthropogenic impact on climate is the equilibrium global annual mean surface air temperature rise caused by the radiative forcing of a doubling of atmospheric CO 2 concentration.While this quantity, called equilibrium climate sensitivity (ECS), can be calculated from climate models (e.g.Vial et al., 2013), it is also important for model validation to make estimates based on palaeodata.This is especially relevant since some important feedbacks of the climate system are not incorporated into all models.For example, when coupling a climate model Published by Copernicus Publications on behalf of the European Geosciences Union.
interactively to a model of stratospheric chemistry, including ozone, the calculated transient warming on a 100-year timescale differs by 20 % from results without such an interactive coupling (Nowack et al., 2015).
The ultimate cause for orbital-scale climate change is latitudinal and seasonal variation in the incoming solar radiations (Milankovitch, 1941;Laskar et al., 2004), which are then amplified by various feedbacks in the climate system (Hays et al., 1976).So far, seasonality in incoming solar radiation is not resolved in our approach.
A major restriction of any geological-data-based estimate of climate sensitivity is that there was no period in Earth's history during which the atmospheric CO 2 concentration and global temperature changed as rapidly as today.Therefore, in all these data-based approaches (including our study here), ECS defined as global equilibrium temperature rise in response to a doubling of atmospheric CO 2 can only be roughly estimated.Such data-based studies are nevertheless important to find any specific pattern of how global temperature changed with respect to a given variation in the radiative forcing.Our approach focuses on the contribution of various climate feedbacks to the reconstructed global temperature changes (PALAEOSENS-Project Members, 2012).When using palaeodata to calculate climate sensitivity, one has to correct for slow feedbacks, whose impacts on climate are incorporated into the temperature reconstructions.Slow feedbacks are of interest in a more distant future (Zeebe, 2013) but are not yet considered in climate simulations using fully coupled climate models underlying the fifth assessment report of the IPCC (Stocker et al., 2013).More generally, from palaeodata the specific equilibrium climate sensitivity, S [X] , is calculated, which is, in line with the proposed nomenclature of PALAEOSENS-Project Members (2012), the ratio of the equilibrium global (g) surface temperature change ( T g ) over the specific radiative forcing ( R) of the processes X; hence S . In this concept "slow feedbacks" are treated the same as the radiative forcing for practical reasons.The division into "forcing" and "feedback" is based on the timescale of the process.PALAEOSENS-Project Members (2012) found that a century is a welljustified timescale that might distinguish fast feedbacks from slow forcings.All relevant processes that are not considered in the forcing term X will nevertheless impact on climate change as feedbacks and are contained in the calculated climate sensitivity.This has to be kept in mind for comparing model-based and data-based approaches, and it makes their comparison difficult, since in model-based results only those processes implemented in the model have an impact on calculated temperature change.
In practical terms, the palaeodata that are typically available for the calculation of S are the radiative forcing of CO 2 and surface albedo changes caused by land ice (LI) sheets.Thus, S [CO 2 ,LI] can be calculated containing the radiative forcing of two processes, which are most important on glacial-interglacial timescales of the late Pleistocene (Köhler et al., 2010).The whole approach, therefore, relies on the simplification that the climate response of the CO 2 radiative forcing and the surface albedo radiative forcing are similar.We are aware that such a simplification might not be possible for every radiative forcing, since Shindell (2014) showed that the per unit radiative forcing of well-mixed greenhouse gases (e.g.CO 2 or CH 4 ) leads to a different climate response than that of aerosols or ozone.However, we are not aware that a difference in the response has been shown for radiative forcing from surface albedo changes ( R [LI] ) and CO 2 ( R [CO 2 ] ).Hence we combine them linearly.
Both model-based (e.g.Crucifix, 2006;Hargreaves et al., 2007;Yoshimori et al., 2011;Yin and Berger, 2012;Caballero and Huber, 2013;Goldner et al., 2013;Kutzbach et al., 2013;Meraner et al., 2013) and palaeodata-based (PALAEOSENS-Project Members, 2012;von der Heydt et al., 2014) approaches have already indicated that S varies for different background climates; see also a recent review of Knutti and Rugenstein (2015) on the limits of linear models to constrain climate sensitivity.The majority of simulation studies shows a rise in climate sensitivity for a warmer background climate.One of the exceptions, based on analysis for mainly colder than present climates (Kutzbach et al., 2013), found the opposite (a rise in climate sensitivity for colder climate) with various versions of the Community Climate System Model (CCSM), which points to disagreements that still exist between models.However, Caballero and Huber (2013), using the same model, found rising climate sensitivity for warmer climates, as did the majority of studies.
The state-dependent character of S based on palaeodata was only recently investigated more systematically in von der Heydt et al. (2014).It was found that the strength of some of the fast feedbacks depends on the background climate state.This is in agreement with other model-based approaches, which proposed a state dependency of water vapour (Meraner et al., 2013) or clouds (Crucifix, 2006;Hargreaves et al., 2007).Distinguishing different climate regimes in palaeodata covering the last 800 000 years (0.8 Myr), the time for which ice core records exist, von der Heydt et al. ( 2014) revealed a ∼ 36 % larger S [CO 2 ,LI] for "warm" background climates when compared to "cold" climates.However, a limitation in this analysis was that average "warmer" climates were still colder than in the present day and interglacial periods were largely undersampled.A recent investigation (Martínez-Botí et al., 2015) found that S [CO 2 ,LI] for the late Pleistocene and the Plio-Pleistocene transition similarly suggest that no state dependency in the specific equilibrium climate sensitivity is observed in their proxy data.
Here, we consider changes in S [CO 2 ,LI] over the last 5 Myr.We go beyond previous studies in various ways.First, we increase the amount and spread of the underlying data, which offers the possibility to calculate S [CO 2 ,LI] based on palaeodata covering the Pleistocene and most of the Pliocene.The latter is the comparatively warm epoch between ∼ 2.6 and 5.3 Myr BP that has been suggested as a palaeoanalogue for the future (Haywood et al., 2010).Second, we calculate the radiative forcing of the land ice albedo from a detailed spatial analysis of simulated land ice distribution obtained with 3-D ice-sheet models.Our approach therefore enhances the embedded complexity of the underlying physical climate system with respect to previous studies.Third, polar amplification was previously assumed to be constant over time (e.g.van de Wal et al., 2011).However, climate models (Haywood et al., 2013) indicate that during the Pliocene, when less ice was present on the Northern Hemisphere, the temperature perturbations were more uniformly spread over all latitudes.We incorporate this changing polar amplification into our global temperature record.Fourth, we explicitly analyse for the first time whether the relationship between temperature change and radiative forcing is better described by a linear or a nonlinear function.If the applied statistics inform us that the T g -R relationship contains a non-linearity, then the specific equilibrium climate sensitivity is state-dependent.Any knowledge on a state dependency of S is important for the interpretation of palaeodata and for the projection of long-term future climate change.

Methods
We calculate the radiative forcing of CO 2 and land ice albedo, R [CO 2 ,LI] , by applying the same energy balance model as used before for the late Pleistocene (Köhler et al., 2010).This approach uses CO 2 data from ice cores as well as from proxies from three different labs published for the last 5 Myr and calculates changes in surface albedo from zonally averaged changes in land ice area.The latter are here based on results from 3-D ice-sheet model simulations (de Boer et al., 2014) that deconvolved the benthic δ 18 O stack LR04 (Lisiecki and Raymo, 2005) into its temperature and sea level (ice volume) component.The time series of global temperature change, T g , over the last 5 Myr used here is also based on this deconvolution.The reconstructed records of ice volume and temperature changes are therefore mutually consistent.A state dependency in S [CO 2 ,LI] is then supported by the data if a non-linear function (higher-order polynomial) gives a statistically better fit to the scattered data of T g versus R [CO 2 ,LI] than a linear fit.
2.1 Ice-sheet models, changes in surface albedo, and radiative forcing, ∆R [LI] Using an inverse modelling approach and the 3-D ice-sheet model ANICE (de Boer et al., 2014), the benthic δ 18 O stack LR04 (Lisiecki and Raymo, 2005) is deconvolved in deepocean temperature, ice-volume-based sea-level variations, and a representation of the four main ice sheets in Antarctica, Greenland, Eurasia, and North America.The spatial resolution (grid cell size) for the Antarctic, Eurasian, and North American ice sheets is 40 km × 40 km, while Greenland is simulated by cells of 20 km × 20 km.In the vertical dimension, velocities and temperature are calculated for 15 layers.In ANICE, shallow ice and shallow shelf approximations are used.With respect to the full Stokes 3-D description that completely describes the temporal and spatial evolution of an ice body, some higher-order stress terms are therefore neglected in ANICE in order to allow for long transient runs.
A detailed description of the model is found in de Boer et al. (2013).This approach combines palaeodata and mass conservation for δ 18 O with physical knowledge on ice-sheet growth and decay.It therefore includes a realistic estimate of both volume and surface area of the major ice sheets.The calculated change in deep-ocean temperature is in this ice-sheetcentred approach connected with temperature anomalies over land in the Northern Hemisphere (NH) high-latitude band (40-85 • N, T NH ), in which the Greenland, Eurasian, and North American ice sheets grow.Temporal resolution of all simulation results from the 3-D ice-sheet models is 2 kyr.
From these results, published previously (de Boer et al., 2014), the latitudinal distribution of land ice area in latitudinal bands i of 5 • ( A LI (i)) is calculated (Fig. 1b), which leads to changes in the land-ice-sheet-based radiative forcing, R [LI] , with respect to preindustrial times.R [LI] (i) for every latitudinal band (Fig. 1c) is calculated from local surface insolation (I S (i)), changes in ice-sheet area ( A LI (i)), and surface albedo anomalies ( α), normalized to a global impact (by division by the Earth's surface area A Earth , and integrated thereafter.For the calculation of I S (i), the annual mean insolation at the top of the atmosphere (TOA) at each latitude, I TOA (i), (Fig. 1a) is reduced by absorption a and reflection α A within the atmopshere (I S (i The values of the parameters a = 0.2 and α A = 0.212 are here held constant at their present values derived in Köhler et al. (2010).This approach to calculate R [LI] is based on surface albedo anomalies ( α), implying that latitudes that are always ice-free contribute nothing to R [LI] .It is assumed that ice sheets cover land when growing; thus, local surface albedo α rises as assumed previously (Köhler et al., 2010) from 0.2 to 0.75.For calculating I TOA (i) (Fig. 1a), which varies due to orbital configurations (Laskar et al., 2004) early 2003 as published by the SORCE (Solar Radiation and Climate Experiment) project (http://lasp.colorado.edu/home/sorce) (Kopp and Lean, 2011).Changes in solar energy output are not considered but are -based on present knowledge (Roth and Joos, 2013) -smaller than 1 W m −2 during the last 10 kyr, and, following our earlier approach (Köhler et al., 2010), presumably smaller than 0.2 %.
For validation of the ANICE ice-sheet model, we compare the spatially and temporally variable results in R [LI] obtained for Termination I (the last 20 kyr) with those based on the land ice sheet distribution of Peltier (2004).This paper describes an approach called ICE-5G, in which data on sea level change, which include the contribution from glacial isostatic adjustment, are used to obtain a physically consistent picture that also considers viscoelastic modelling of the mantle of Earth in order to derive how the land ice sheet distribution during the last deglaciation might have looked like.For this comparison the ICE-5G data are treated similarly to those from ANICE, e.g.only data every 2 kyr are considered and averaged on latitudinal bands of 5 • .The spatial distribution of land ice in the most recent version of ICE-6G (Peltier et al., 2015) are similar to ICE-5G, and therefore no significant difference to ICE-6G is expected and the comparison to that version is omitted.

Global temperature change, ∆T g
In the ANICE model (de Boer et al., 2014) the temperature anomaly of the deep ocean (deconvolved from the benthic δ 18 O stack) is coupled to the NH temperature change, T NH , by a fixed ratio that has been derived in a series of transient climate runs.A more extensive analysis of this parameterization is presented in de Boer et al. (2010).
We calculate global surface temperature change, T g , from these ANICE-based NH temperature anomalies, T NH , using a polar amplification factor, f pa , which itself depends on climate (Fig. 2).
The value of the polar amplification factor, f pa , was constrained for certain times from simulation results of the PMIP3@LGM PlioMIP@mPWP T g1 : f pa linear T g2 : f pa step function T g3 : f pa linear matching MIP T g @LGM and @mPWP T g1 : f pa linear T g2 : f pa step function T g3 : f pa linear matching MIP T g @LGM and @mPWP  (Haywood et al., 2013), respectively.In our standard application, T g1 (black line), f pa is calculated as a linear function depending on northern hemispheric temperature change, T NH (insert), inter-and extrapolated between these two PMIP3 and PlioMIP-based values.Alternatively ( T g2 , orange line), f pa varies as a step function with high values for the Pleistocene (periods with Northern Hemisphere land ice sheets) and low values for the Pliocene (periods mainly without NH land ice sheets) with the step between both values occurring at 2.82 Myr BP, when our results indicate large changes in NH land ice.In T g3 (blue line), f pa varied freely to meet T g reconstructed for LGM by PMIP3 (−4.6 K) and for the mPWP by PlioMIP (+2.7 K).See methods for further details.Panel (b): NH temperature change, T NH , as deconvolved from the benthic δ 18 O stack LR04 (Lisiecki and Raymo, 2005) by applying a 3-D ice-sheet model in an inverse mode (de Boer et al., 2014).Uncertainty in T NH (grey) is the 1σ error calculated from eight different model realizations (de Boer et al., 2014).Panel (c): global surface temperature change, T g , as used here based on T g = T NH • f −1 pa .Results for T g based on all three approaches for f pa are given (same colour code as in subfigure a).Symbols show T g ± 1σ as derived within PlioMIP (mPWP, red circle) and PMIP3-CMIP5 (LGM, blue square).Red vertical lines mark the time period of the mPWP.(PlioMIP).For the Last Glacial Maximum (LGM, about 20 kyr BP), f pa was determined from PMIP3-CMIP5 (Braconnot et al., 2012) to be around 2.7 ± 0.3; for the mid-Pliocene Warm Period (mPWP, about 3.2 Myr BP), f pa was determined to be around 1.6 ± 0.1 based on PlioMIP results (Haywood et al., 2013).In our standard set-up (calculating T g1 ), we linearly inter-and extrapolate f pa as function of T NH based on these two anchor values for all background climates found during the last 5 Myr (insert in Fig. 2a).Climate models already suggest that polar amplification is not constant, but how it is changing over time is not entirely clear (Masson-Delmotte et al., 2006;Abe-Ouchi et al., 2007;Hargreaves et al., 2007;Yoshimori et al., 2009;Singarayer and Valdes, 2010).We therefore calculate an alternative global temperature change, T g2 , in which we assume polar amplification, f pa , to be a step function, with f pa = 1.6 (the mPWP value) taken for times with extensive northern hemispheric land ice (according to our results before 2.82 Myr BP) and with f pa = 2.7 (the LGM value) thereafter.This choice is motivated by investigations with a coupled ice-sheet-climate model, from which northern hemispheric land ice was identified to be the main controlling factor for the polar amplification (Stap et al., 2014).
At the LGM, T g was, based on the eight PMIP3 models contributing to this estimate in f pa , −4.6 ± 0.8 K, so slight colder but showing considerable overlap with the most recent LGM estimate of T g = −4.0 ± 0.8 K (Annan and Hargreaves, 2013).If we take into consideration that the MARGO (Multiproxy Approach for the Reconstruction of the Glacial Ocean surface) sea surface temperature (SST) data underlying this LGM temperature estimate (Annan and Hargreaves, 2013) are potentially biased towards tropical SSTs that are too warm (Schmidt et al., 2014), the PMIP3 results are a good representation of LGM climate.For both choices of f pa (varying linearly as a function of T NH or as step function over time), the global temperature change at the LGM obtained in our reconstruction is T g = −5.7 ± 0.6 K, so slightly colder than other approaches but overlapping with the PMIP3-based results within the uncertainties.
As a third alternative ( T g3 ), we constrain the global temperature changes by the values from PMIP3 for the LGM (−4.6 K) and from PlioMIP for the mPWP (+2.7 K) and vary f pa freely.In this case, f pa rises to ∼ 3.3 during glacial maxima of the Pleistocene and to ∼ 1.0 during the Pliocene.Our approach based on the T NH reconstruction is not able to meet all four constraints given by PMIP3 and PlioMIP ( T g , f pa for both the LGM and the mPWP) at the same time.This mainly illustrates that the approach used in de Boer et al. ( 2014), although coherently solving for temperature and ice volume underestimates polar temperature change prior to the onset of the NH glacial inception, since it only calculates ice volume and deep-water temperature change from benthic δ 18 O.
Throughout the following, our analysis is based on the temperature time series T g1 .However, we repeat our analysis with the alternative temperature time series to investigate the robustness of our approach to the selected time series.As can been seen in the results our main conclusions and functional dependencies are independent from the choice of T g and are also supported if based on either T g2 or T g3 (see Sect. 3.2).
The modelled surface-air temperature change, T NH , has already been compared (de Boer et al., 2014) with three independent proxy-based records of sea surface temperature (SST) change in the North Atlantic (Lawrence et al., 2009), the equatorial Pacific (Herbert et al., 2010), and the Southern Ocean (Martínez-Garcia et al., 2010) which cover at least the last 3.5 Myr.The main features of the simulated temperature change and the data-based SST reconstruction agree: the overall cooling trend from about 3.5 to 1 Myr ago is found in the simulation results and in all SST records and so is the strong glacial-interglacial (100 kyr) variability thereafter.

Radiative forcing of CO
Several labs have developed different proxy-based approaches to reconstruct atmospheric CO 2 for times before 0.8 Myr, for which no CO 2 data from ice cores exist yet.Since we are interested in how CO 2 might have changed over the last 5 Myr, and in its relationship to global climate, we only consider longer time series for our analysis.Thus, some approaches, e.g. based on stomata, with only a few data points during the last 5 Myr are not considered (see Beerling and Royer, 2011).The considered CO 2 data are, in detail, as follows (Fig. 3): a. Ice core CO 2 data were compiled by Bereiter et al. (2015) into a stacked ice core CO 2 record covering the last 0.8 Myr, including a revision of the CO 2 data from the lowest part of the EPICA Dome C ice core.Originally, the stack as published (Bereiter et al., 2015) con-  (Hönisch et al., 2009;Foster, 2008;Martínez-Botí et al., 2015) or alkenones (Pagani et al., 2010;Zhang et al., 2013) from the three labs, Hönisch, Foster, and Pagani.Panel (b): zoom-in on the time period covered by ice core data (last 0.8 Myr.) tained data points before the year 1750 CE, the beginning of the industrialization, but it was here resampled to the 2 kyr time step of the ice-sheet simulation results by averaging available data points and reducing the sample size to n = 394.The stack contains data from the ice cores at Law Dome (Rubino et al., 2013;MacFarling-Meure et al., 2006; 0-2 kyr BP), EPICA Dome C (Monnin et al., 2001(Monnin et al., , 2004;;Schneider et al., 2013;Siegenthaler et al., 2005;Bereiter et al., 2015;2-11 kyr BP, 104-155 kyr BP, 393-806 kyr BP), West Antarctic Ice Sheet Divide (Marcott et al., 2014;11-22 kyr BP), Siple Dome (Ahn and Brook, 2014;22-40 kyr BP), Talos Dome (Bereiter et al., 2012;40-60 kyr BP), EPICA Dronning Maud Land (Bereiter et al., 2012;60-104 kyr BP), and Vostok (Petit et al., 1999;155-393 (Foster, 2008;Martínez-Botí et al., 2015) are available for the last 3.3 Myr (n = 105), obtained via δ 11 B from ODP site 999 in the Caribbean Sea.CO 2 purely based on G. ruber δ 11 B was reconstructed for the last glacial cycle (Foster, 2008) and for about 0.8 Myr during the Plio-Pleistocene transition (Martínez-Botí et al., 2015).We take both these data sets, using identical calibration as plotted previously (Martínez-Botí et al., 2015).The overlap of the data with the ice core CO 2 is reasonable, with a tendency to overestimate the maximum anomalies in CO 2 (by more than +50 ppmv during warm previous interglacials and by −25 ppmv during the LGM; Fig. 3b).
d. CO 2 reconstructions based on alkenone from the Pagani lab (Pagani et al., 2010;Zhang et al., 2013) (n = 153) cover the whole 5 Myr and are derived from different marine sediment cores.Site 925 is contained in both publications, although with different uncertainties.
From site 925 we use the extended and most recent CO 2 data of Zhang et al. (2013) containing 50 data points over the last 5 Myr, 18 points more than initially published.Data from the sites 806, 925, and 1012 are different from the ice core CO 2 reference during the last 0.8 Myr by +50 to +100 ppmv, while data from site 882 have no overlapping data points with the ice cores.It is not straightforward to correct these CO 2 data from the Pagani lab that are different from the ice core CO 2 .Therefore, we refrain from applying any corrections but keep these offsets in mind for our interpretation.
Other CO 2 data based on B/Ca (Tripati et al., 2009) are not considered here, since critical issues concerning their calibration have been raised (Allen et al., 2012).A second δ 11 Bbased record of the Hönisch lab (Bartoli et al., 2011) from G. sacculifer obtained from ODP site 999 is not used for further analysis because δ 11 B was measured on other samples than those proxies that are necessary to determine the related climate state (e.g.δ 18 O).Thus, a clear identification of whether glacial or interglacial conditions were prevailing for individual data points was difficult.Furthermore, these calculated CO 2 values (Bartoli et al., 2011) have very high uncertainties: 1σ is 3 times larger than in the original Hönisch lab data set (Hönisch et al., 2009).These CO 2 data of Bartoli et al. (2011) disagree with the most recent data from the Foster lab (Martínez-Botí et al., 2015), especially before the onset of northern hemispheric glaciation around 2.8 Myr ago.Another CO 2 time series form the Foster lab (Seki et al., 2010), based on a mixture of both alkenones and δ 11 B approaches covering the last 5 Myr, is not considered here, since the applied size correction for the alkenone pro-ducers was subsequently been shown to be incorrect (Badger et al., 2013).
Radiative forcing based on CO 2 is calculated us- • ln(CO 2 /CO 2,0 ), with CO 2,0 = 278 ppmv being the preindustrial reference value (Myhre et al., 1998).The specific equilibrium climate sensitivity for a forcing X is defined as . In an analysis of S [X] when calculated for every point in time for the last 0.8 Myr based on ice core data, PALAEOSENS-Project Members (2012) revealed the range of possible values, which fluctuated widely, not following a simple functionality even when analysed as moving averages.This study also clarified that S [X] based on small disturbances in T g or R [X] is due to dating uncertainties prone to unrealistically high or low values.Only when data are analysed in a scatter plot does a non-linear functionality between T g and R [X] , and therefore a state dependency of S [X] , emerge as a signal out of the noisy data (von der Heydt et al., 2014).
Here, T g is approximated as a function of R [X] by fitting a non-linear function (a polynomial up to the third order, . The individual contribution of land ice albedo and CO 2 to a state dependency of S [CO 2 ,LI] can be investigated by analysing both S [CO 2 ] and S [CO 2 ,LI] .If the best fit follows a linear function, e.g. for state-independent behaviour of S [X] , its values might be determined from the slope of the regression line in the T g -R [X] space.However, note that here a necessary condition for the calculation of S [X] over the whole range of R [X] , but not for the analysis of any state dependency, is that any fitting function crosses the origin with R [CO 2 ,LI] = 0 W m −2 and T g = 0 K, implying for the fitting parameters that a = 0.This is also in line with the general concept that without any change in the external forcing, no change in global mean temperature should appear.Since the data sets have apparent offsets from the origin, we first investigate which order of the polynomial best fits the data by allowing parameter a to vary from 0. For the calculation of mean values of S [CO 2 ,LI] , we then analyse the LI] is first calculated individually for every data point and then stacked for different background conditions (described by R [CO 2 ,LI] ).In doing so, we circumvent the problem which appeared in the T g -R [X] space that the regression function needs to meet the origin.Some of the individual values of S [CO 2 ,LI] are still unrealistically high or low; therefore, values in S [CO 2 ,LI] outside the plausible range of 0-3 K W −1 m 2 are rejected from further analysis.
The scattered data of S [CO 2 ,LI] as a function of R [CO 2 ,LI] are then compiled in a probability density function (PDF), in which we also consider the given uncertainties of the individual data points.For the calculation of the PDFs, we distinguish between a few different climate states, when supported by the data.For the time being, the data coverage is too sparse and uncertainties are too large to calculate any state-dependent values of S [CO 2 ,LI] in greater detail.
The fitting routines (Press et al., 1992) use the method of general linear least squares.Here, a function χ 2 = n i y is minimized, which calculates the sum of squares of the offsets of the fit from the n data points normalized by the average variance σ 2 y .Since established numerical methods for calculating a non-linear fit through data cannot consider uncertainties in x, we base our regression analysis on a Monte Carlo approach.Data points are randomly picked from the Gaussian distribution described by the given 1σ standard deviation of each data point in both directions x ( R [X] ) and y ( T g ).We generated 5000 of these data sets, calculated their individual non-linear fits, and further analysed results based on averages of the regression parameters.The Monte Carlo approach converges if the number of replicates exceeds 1000, e.g.variations in the mean of the parameters are at least 1 order of magnitude smaller than the uncertainties connected with the averaging of the results.We used the χ 2 of the fitting routines in F tests to investigate whether a higher-order polynomial would describe the scattered data in the T g -R [X] parameter space better than a lower-order polynomial, and we use the higher-order polynomial only if it significantly better describes the data at the 1 % level (p value of F test: p ≤ 0.01).

Uncertainty estimates
As previously described in detail (Köhler et al., 2010), standard error propagation is used to calculate uncertainties in T and R. For R [LI] , changes in surface albedo are assumed to have a 1σ uncertainty of 0.1.Simulated changes in land ice area have a relative uncertainty of 10 % in the various simulation scenarios performed in de Boer et al. (2014).The different approaches in reconstructing CO 2 all have different uncertainties, as plotted in Fig. 3. Ice core CO 2 has a 1σ uncertainty of 2 ppmv, while those based on other proxies have individual errors connected with the data points that are of the order of 20-50 ppmv.Radiative forcing based on CO • ln(CO 2 /CO 2,0 ), has in addition to the uncertainty in CO 2 itself also another 10 % 1σ uncertainty (Forster et al., 2007).The uncertainty in the incoming insolation is restricted to variations in the solar constant known to be of the order of 0.2 %.Annual mean global surface temperature, T g , is solely based on the polar amplification factor, f pa , and T NH .Uncertainty in T NH is estimated based on eight different model realizations of the deconvolution of benthic δ 18 O into sea level and temperature (de Boer et al., 2014).Based on the analysis of the PMIP3 and PlioMIP results, the polar amplification factor has a relative uncertainty of 10 % (see Fig. 2a).
These uncertainties used in an error propagation lead to the σ T g , σ R [CO 2 ] , and σ R [CO 2 ,LI] of the individual data points and are used to constrain the Monte Carlo statistics.The stated uncertainties of the parameters of the polynomials fitting the scattered T -R data given in Table 1 and used to plot and calculate S [CO 2 ,LI] are derived from averaging results of the Monte Carlo approach.Note that higher-order polynomials give more constrains on the parameters and therefore lead to smaller uncertainties.

Individual radiative forcing contributions from land ice albedo and CO 2
We calculate a resulting radiative forcing of CO 2 , R [CO 2 ] , that spans a range of −2.8W m −2 to +2.5 W m −2 when compared to the forcing of preindustrial conditions, e.g. when R [CO 2 ] = 0 W m −2 (Fig. 4b).The uncertainty in R [CO 2 ] depends on the proxy.It is about 10 % in ice cores, and generally less than 0.5 W m −2 for other proxies, with the exception of some individual points from the Pagani lab with uncertainties of around 1 W m −2 .
In contrast to these rather uncertain and patchy results, the ice-sheet simulations lead to a continuous time series of surface albedo changes and R [LI] ranging between −4 W m −2 during the ice ages of the late Pleistocene and +1 W m −2 during the interglacials of the Pliocene (Fig. 4c).During warmer than preindustrial climate, R [LI] is thus rather small and between 4.2 and 3.0 Myr ago only slightly higher than R [orbit] , the radiative forcing due to global annual mean insolation changes caused by variations in the orbital parameters of the solar system (Laskar et al., 2004) (Fig. 4c).
Reconstructed R [LI] for the last 20 kyr agrees reasonably well with an alternative based on the ICE-5G ice-sheet reconstruction of Peltier (2004) (Fig. 5).Changes in land ice fraction in the northern high latitudes around 15 kyr ago are more abrupt around 70 • N in ICE-5G than in ANICE (Fig. 5b, e).The northward retreat of the southern edge of the NH ice sheets happens later in ICE-5G than in ANICE.In combination, both effects lead to only small differences in the spatial and temporal distribution of the radiative forcing, R [LI] , when based on either ANICE or ICE-5G (Fig. 5c and f).
The ice-albedo forcing, R [LI] , has a non-linear relationship to sea level change (Fig. 6a), which is caused by the use of sophisticated 3-D ice-sheet models.Hence, other approaches which approximate R [LI] directly from sea level (Hansen et al., 2008;Martínez-Botí et al., 2015) or from simpler 1-D ice-sheet models or calculate R [LI] from global land ice area changes without considering latitudinal dependency (Köhler et al., 2010;von der Heydt et al., 2014) lack an important non-linearity of the climate system.This nonlinearity in the R [LI] sea level relationship is also weakly ).CO 2 data are from ice cores (Bereiter et al., 2015) and based on δ 11 B (Hönisch lab (Hönisch et al., 2009), Foster lab (Foster, 2008;Martínez-Botí et al., 2015)) and on alkenones (Pagani lab (Pagani et al., 2010;Zhang et al., 2013)).Panel contained in results based on ICE-5G for Termination I (Fig. 6a).However, when plotting identical time steps of Termination I from results based on ANICE, the non-linearity is not yet persistent.This implies that a larger pool of results from various different climates needs to be averaged in or- der to obtain a statistically robust functional relationship between R [LI] and sea level (as done in this study).
The combined forcing, R [CO 2 ,LI] , can only be obtained for the data points for which CO 2 data exist (Fig. 4d).The combined forcing ranges from −6 to −7 W m −2 during the LGM to, in general, positive values during the Pliocene with a maximum of +3 W m −2 .Between 5.0 and 2.7 Myr ago (most of the Pliocene), the ice-sheet area and R [LI] , apart from two exceptions around 3.3 Myr and after 2.8 Myr ago (Fig. 4c), are smaller than today suggesting warmer temperatures throughout.Proxy data suggest that CO 2 and R [CO 2 ] were mostly higher in the Pliocene than during preindustrial times.

Detecting any state dependency in S [CO 2 ,LI]
As explained in detail in Sect.2.4, S [CO 2 ,LI] can be considered state-dependent if the scattered data of T g against R [CO 2 ,LI] are better described by a non-linear rather than a linear fit.The plots for the different CO 2 approaches reveal proxy-specific results (Fig. 7).Ice core data (r 2 = 0.72) are best described by a third-order polynomial and the Hönisch data (r 2 = 0.79) by a second-order polynomial, while for the Foster (r 2 = 0.61) and the Pagani (r 2 = 0.45) data, a secondorder fit is not statistically significantly better than a linear fit (Table 1).
The fit through the Hönisch data agrees more with the fit through the ice core CO 2 data than with the fit through the other CO 2 -proxy-based approaches; however, the Hönisch data set only extends 2.1 Myr back in time and contains no CO 2 data in the warm Pliocene.Thus, our finding of a state dependency in climate sensitivity obtained from the ice core data and covering predominately colder than present periods -for which a first indication was published in von der Heydt et al. ( 2014) -is extended to the last 2.1 Myr.However, we can still not extrapolate this finding to the warmer than present climates of the last 5 Myr since the ice core and the Hönisch data do not cover these periods and the Foster and the Pagani data do not suggest a similar relationship.These findings remain qualitatively the same if our analyses are based on the alternative global temperature changes T g2 and T g3 (Table 1 n: number of data points in data set.χ 2 : weighted sum of squares following either a linear fit (first-order) or a non-linear fit (second-order polynomial); for some data sets (labelled: a ), there are also secondor third-order polynomials.F : F ratio for F test to determine where the higher-order fit describes the data better than the lower-order fit (first-vs.second-order polynomial or second-vs.third-order polynomial).p: p value of the F test.L: significance level of F test (/: not significant (p > 0.01); * : significant at 1 % level (0.001 < p ≤ 0.01); * * : significant at 0.1 % level (p ≤ 0.001)).When analysing the contribution from land ice albedo ( R [LI] ) and CO 2 radiative forcing ( R [CO 2 ] ) separately, we find a similar non-linearity in the T g -R [CO 2 ] scatter plot only in the CO 2 data from ice cores (Fig. 7a).The relationship between temperature and radiative forcing of CO 2 is best described by a linear function in the Hönisch and the Pagani data sets (Fig. 7c and g, Table 1) and in data from the Foster lab even by a second-order polynomial with inverse slope leading to a decline in S [CO 2 ] for rising R [CO 2 ] (Fig. 7e).This inverse slope obtained for the Foster data between T g and R [CO 2 ] is the only case in which a detected non-linearity partly depends on the use of the temperature change time series, e.g. the relationship is linear when based on T g3 (Table 1).Furthermore, this inverse slope might be caused by the under-representation of data for negative radiative forcing.However, when data in the T g -R [X] scatter plots are binned in x or y direction to overcome any uneven distribution of data, no change in the significance of the non-linearities is observed.The data scatter is large and regression coefficients between R [CO 2 ] and T g for Foster (r 2 = 0.42) and Pagani (r 2 = 0.03) are small.This large scatter and weak quality of the fit in the Pagani data are probably caused by some difficulties in the alkenone-based proxy, e.g. in the size dependency, and by variations in the degree of passive vs. active uptake of CO 2 by the alkenone-producing coccolithophorids (Bolton and Stoll, 2013).Furthermore, van de Wal et al. (2011) already showed that the relationship of CO 2 to temperature change during the last 20 Myr is opposite in sign for alkenone-based CO 2 than for other approaches.
The ice-albedo forcing, R [LI] , in our simulation results based on 3-D ice-sheet models (de Boer et al., 2014) has a specific relationship to global temperature change.Here both a step function or a linear change in the polar amplification factor, f pa , lead to similar results (Fig. 6b).However, when overly simplified approaches to calculate R [LI] are applied (e.g. based on 1-D ice-sheet models (de Boer et al., 2010), related to sea level (Hansen et al., 2008;Martínez-Botí et al., 2015), or based on global land ice area changes without considering their latitudinal changes in detail (Köhler et al., 2010;PALAEOSENS-Project Members, 2012;von der Heydt et al., 2014)), the T g -R [LI] relationship is more linear.The range of R [LI] proposed for the same range of T g is then reduced by 30 % (Fig. 6b and c).R [LI] is effected by ice-sheet area rather than ice-sheet volume.Threedimensional ice-sheet models include this effect better than calculations based on 1-D ice-sheet models or calculations directly from sea level variations.This non-linearity between ice volume (or sea level) and ice area is supported by data and theory from the scaling of glaciers (Bahr, 1997;Bahr et al., 2015).In addition, latitudinal variation in land ice distribution affects the radiative forcing, R [LI] , in a non-linear way (Fig. 1) and thereby likely contributes to a state dependency in the equilibrium climate sensitivity, S [CO 2 ,LI] .
To verify the robustness of our findings with respect to the uncertainties attached to all data points, we performed an additional sensitivity study by artificially reducing the uncertainties in T g (σ T g ) and R [CO 2 ,LI] (σ R ) by a factor of 2 or 10.For both reduction factors, we find statistically the  (b, d, f, h): radiative forcing of CO 2 and land ice albedo ( R [CO 2 ,LI] ).Lines show average best fits (first-, second-, or thirdorder polynomials) to 5000 Monte Carlo realizations of the data (details in Table 1).Subfigures differ by the CO 2 data they are based on: (a, b) -ice cores (Bereiter et al., 2015) ; (c, d)δ 11 B from the Hönisch lab (Hönisch et al., 2009); (e, f)δ 11 B from the Foster lab (Foster, 2008;Martínez-Botí et al., 2015); (g, h) -alkenones from the Pagani lab (Pagani et al., 2010;Zhang et al., 2013).Each row contains information on the number of data points n; each subfigure shows the mean uncertainty of the fit by dividing χ 2 (the weighted sum of squares from the regression analysis) by n and the correlation coefficient r 2 .Uncertainties show 1σ .same non-linearities in the T g -R [CO 2 ,LI] -scattered data than with the original uncertainties in all four CO 2 data sets (non-linearity between T g and R [CO 2 ,LI] in the data sets based on CO 2 in ice cores and from the Hönisch lab; a linear relationship between both variables only exists if based on Foster or Pagani lab CO 2 data; Table 2).Our proposed state dependency of S [CO 2 ,LI] is therefore independent of the assumed uncertainties.Any calculated value of S [CO 2 ,LI] nevertheless depends in detail on the assumed uncertainties in the underlying data.
Since a first detection of any state dependency in S [CO 2 ,LI] has already been performed for the ice core CO 2 data in von der Heydt et al. (2014), it is of interest to investigate which of our improvements with respect to this earlier analysis are most important.We therefore performed a further sensitivity study in which one or all of the three times series T g , R [CO 2 ] , and R [LI] were identical to the approach of von der Heydt et al. ( 2014).However, since in this earlier study all data have been resampled to 100 yr, we have to preprocess these data sets prior to Monte Carlo statistics to 2-kyr averages to match the temporal resolution of the 3-D ice-sheet models used here.In this additional analysis (Table 2), we find that even when all three data sets are substituted with those used in von der Heydt et al. ( 2014), we find a non-linearity in the T g -R [CO 2 ,LI] scatter plot that points to a state dependency in S [CO 2 ,LI] .However, the r 2 is then 10% smaller than in our results indicating a weaker correlation between temperature change and radiative forcing, and a second-order polynomial is sufficient to fit the data, while in our best guess these ice-core-based CO 2 data are best described by a third-order polynomial.If data are binned before analysis, similarly as in von der Heydt et al. ( 2014), we find a non-linearity in the scattered data only for the data sets used in this study or when CO 2 is substituted by the previous time series but not when the previous versions of R [LI] or T g are used.In these binned data both our improved time series of T g and R [LI] are necessary to generate this non-linearity indicating a state dependency in S [CO 2 ,LI] .The analysis of both studies are still different in detail (higher-order polynomial versus piecewise linear regressions), and therefore the absence of any non-linearity in the binned data when all three time series have been substituted by those from the previous study are not contradictory to our stated non-linearity.
In model-based approaches the final radiative forcing R including all feedbacks from an obtained temperature change leads to a different nomenclature, in which temperature change is the independent variable, typically plotted on the x axis (e.g.Bloch-Johnson et al., 2015).Our approach differs from those studies since feedbacks are not contained in R (but in S), which we only understand as the forcing terms.Therefore, in our study R is the independent variable that determines the background condition of the climate system.
Table 2. Sensitivity analyses: (1) investigating the importance of the uncertainties on the regression results by artificially reducing both σ T g and σ R by a factor of 2 or 10; (2) investigating the importance of the three variables T g , CO 2 , and R [LI] with respect to the previous analysis of the ice-core-based CO 2 data of von der Heydt et al. ( 2014) (cited here as vdH2014).Here, all data are resampled to 2 kyr, while in vdH2014 data are resampled to 100 yrs and binned in T g before any regression analysis.The coefficients a, b, c, and d describe the linear or higher-order polynomial that best fits the data (F test statistics based on 5000 Monte Carlo-generated realizations of the scattered T g -R [CO 2 ,LI] data).The data are randomly picked from the entire Gaussian distribution described by the 1σ of the given uncertainties in both T g and R [CO 2 ,LI] .The parameter values of fitted polynomials are given as mean ±1σ uncertainty from the different Monte Carlo realizations.In all scenarios summarized here, T g vs. R [CO 2 ,LI] with T g = T g1 was investigated.n: number of data points in data set.χ 2 : weighted sum of squares following either a linear fit (first order) or a non-linear fit (second-order polynomial); for some data sets (labelled: a ), there are also second-or third-order polynomials.F : F ratio for F test to determine where the higher-order fit describes the data better than the lower-order fit (first-vs.second-order polynomial or second-vs.third-order polynomial).p: p value of the F test.L: significance level of F test (/: not significant (p > 0.01); * : significant at 1 % level (0.001 < p ≤ 0.01); * * : significant at 0.1 % level (p ≤ 0.001)).r 2 : correlation coefficient of the fit.a, b, c, d: derived coefficients of fitted polynomial y(x) = a + bx + cx 2 + dx 3 .

Calculating the specific equilibrium climate sensitivity, S [CO 2 ,LI]
The non-linear regression of the T g -R based on the pointwise results (Fig. 8).For both the Pagani and the Foster data sets, the slopes of the linear regression lines in T g -R [CO 2 ,LI] might in principle be used to calculate S [CO 2 ,LI] .However, both data sets have a rather large offset in the y direction ( T g ) (y interception is far away from the origin) that might bias these results.These offsets are nearly identical when calculations are based on the alternative global temperature changes T g2 or T g3 (Table 1).Note that S [CO 2 ,LI] as calculated for each data point in Fig. 8 also contains 20 and 11 outliers in the ice core and the Hönisch data sets, respectively, that do not fall into the most P. Köhler et al.: State dependency of the equilibrium climate sensitivity plausible range of 0.0-3.0K W −1 m 2 .These outliers are typically generated when dividing smaller anomalies in T g and R [CO 2 ,LI] during interglacials, when already small uncertainties generate a large change in the ratio that defined S [CO 2 ,LI] .They are omitted from further analysis.
S [CO 2 ,LI] based on the ice core and the Hönisch lab data rarely falls below 0.8 K W −1 m 2 (Fig. 8).In both data sets, we distinguish between "cold" from "warm" conditions using the threshold of R [CO 2 ,LI] = −3.5 W m −2 to make our results comparable to the piecewise linear analysis of "warm" and "cold" periods in von der Heydt et al. (2014).For the ice core data of the last 0.8 Myr, the S [CO 2 ,LI] is not normally distributed but has a long tail towards higher values (Fig. 8c).However, this long tail is partially caused by data points with R [CO 2 ,LI] not far from 0 W m −2 , which are prone to high uncertainties.Only conditions during "cold" periods, representing glacial maxima, have a nearly Gaussian distribution in S [CO 2 ,LI] , with a mean value of 1.05 +0.23  −0.21 K W −1 m 2 .For "warm" periods the PDF is skewed, with S [CO 2 ,LI] = 1.56 +0.60  −0.44 K W −1 m 2 .Results based on the Hönisch data covering the last 2.1 Myr are nearly identical with S [CO 2 ,LI] = 1.07 +0.29  −0.24 K W −1 m 2 ("cold") and S [CO 2 ,LI] = 1.51 +0.68  −0.55 K W −1 m 2 ("warm").Both data sets thus consistently suggest that during Pleistocene warm periods, S [CO 2 ,LI] was about ∼ 45 % larger than during Pleistocene cold periods.
In a piecewise linear regression analysis of data covering the last 0.8 Myr, a state dependency in climate sensitivity was already detected (von der Heydt et al., 2014), including a rise in S [CO 2 ,LI] from 0.98 ± 0.27 K W −1 m 2 during "cold" periods to 1.34 ± 0.12 K W −1 m 2 during "warm" periods of the late Pleistocene.To allow a direct comparison with our study, we here cite results shown in the Supporting Information of von der Heydt et al. (2014) in which the global temperature anomaly was similar to our T g .Some important details, however, of our study and the previous study (von der Heydt et al., 2014) differ because (i) the assumed changes in temperature and land ice albedo are based on different time series and (ii) we here use CO 2 as resampled to the 2 kyr temporal spacing of the 3-D ice-sheet models while all data are resampled at 100 years time steps and binned before analysis in von der Heydt et al. (2014).Note, that we tested that data binning does not lead to large changes in our results and conclusions.Nevertheless, the calculated S [CO 2 ,LI] of the "cold" periods (von der Heydt et al., 2014) matches our glacial values of S [CO 2 ,LI] derived from the ice cores within the uncertainties, but the values for the "warm" periods are smaller in the previous estimates of von der Heydt et al. (2014) than in our results (Fig. 9).This difference in the "warm" period for both studies is caused by the revised R [LI] , which mainly leads to differences with respect to previous studies for intermediate glaciated and interglacial climates.
The calculated PDFs of S [CO 2 ,LI] (Fig. 9) based on ice cores or Hönisch lab data are qualitatively the same if based on the alternative assumptions regarding polar amplification, which also include a case with a constant polar amplification during the Pleistocene.The mean values of the PDF of S [CO 2 ,LI] are then shifted by less than 0.15 K W −1 m 2 for "cold" periods and by less than 0.25 K W −1 m 2 for "warm" periods towards smaller values.
The 5 Myr long data sets from the Foster and the Pagani lab show no indication of state dependency.One might argue that these 5 Myr long time series should be split into times when large ice sheets in the NH were present and when they were absent because their presence is expected to have an influence on the climate and its sensitivity.According to our simulation results (Fig. 1b), large NH land ice first appeared around 2.82 Myr BP, which is also the time which has been suggested by Sarnthein (2013) for the onset of NH land ice and for which Martínez-Botí et al. (2015) found a pronounced decline in CO 2 .Note that the start of northern hemispheric glaciation in our 3-D ice-sheet simulations was gradual at first and intensified around 2.7 Myr ago (Fig. 1b), in agreement with other studies (Raymo, 1994;Haug et al., 2005).We tested the Foster lab data for any changes in the regression analysis if the data set was split into two time periods, one with and one without NH ice sheets.We found significantly different relationships between temperature change and radiative forcing for most of the Pleistocene than for either an ice-free NH Pliocene  or all available Pliocene data (Foster lab data 2.5-3.3Myr BP) (Fig. 10).For the Pleistocene, LI] data are in themselves non-linear (thus, S [CO 2 ,LI] is state-dependent), and for the Pliocene the relationship seems to be linear (thus, S [CO 2 ,LI] seems to be constant) over this period.However, the fit through T g -R [CO 2 ,LI] is of low quality (r 2 = 0.04 for 2.82-3.3Myr BP and r 2 = 0.23 for 2.5-3.3Myr BP,) which prevents us from calculating any quantitive values of S [CO 2 ,LI] based on these data.Remember, that in all regression analyses we consider the uncertainties in both x and y direction in all data points by the application of Monte Carlo statistics, something which also distinguishes our approach from Martínez-Botí et al. (2015) and possibly contributes to different results.
Nevertheless, our data compilation clearly points to a regime shift in the climate system with different climate sensitivities before and after 2.82 Myr BP.From the available proxy-based data indicating CO 2 around 400 ppmv during various times in the Pliocene, together with our simulated global temperature change of around 2 K and ice-sheet albedo forcing of about 0.5 W m −2 (Fig. 4), we can estimate that in the NH ice-free Pliocene, S [CO 2 ,LI] was around 1 K W −1 m 2 , in agreement with Martínez-Botí et al. (2015).This estimate is of a similar size as our results for the full glacial conditions of most of the Pleistocene, but it is smaller than during intermediate glaciated to interglacial conditions of the late Pleistocene.A possible reason could be that in the warm Pliocene, the sea ice albedo feedback might have been weaker or even absent (von der Heydt et al., 2014), but some studies (Stevens and Bony, 2013;Fedorov et al., 2013) also suggest that processes are missing in state-of-the-art climate models.A recent study (Kirtland Turner, 2014) concluded that at the onset of the northern hemispheric glaciation, a fundamental change in the interplay of the carbon cycle and the climate system occurred leading to a switch from in-phase glacial-interglacial changes in deep-ocean δ 18 O and δ 13 C to antiphase changes.If true such a change in the carbon cycleclimate system might also affect climate sensitivity.
A more direct calculation of the specific equilibrium climate sensitivity, S [CO 2 ,LI] , as a function of background climate state that goes beyond the PDFs provided so far is desirable, but with the available data and within the given theoretical and methodological framework, it is not yet possible.-Botí et al. (2015) recently analysed the ice core CO 2 and the new CO 2 data from the Foster lab around the end of the Pliocene separately, finding S [CO 2 ,LI] of 0.91 ± 0.10 and 1.01 ± 0.19 K W −1 m 2 , respectively.Both results are nearly indistinguishable within their uncertainties; thus, Martínez-Botí et al. (2015) concluded that S [CO 2 ,LI] is not state-dependent, since it did not change between the Pliocene and Pleistocene.However, since they based the radiative forcing of land ice albedo ( R [LI] ) on a linear function of sea level, they miss an important non-linearity of the climate system.We find that the large uncertainty in R [CO 2 ] might also be another reason for state independency in S [CO 2 ,LI] in the Foster lab data set.S [CO 2 ,LI] based on the ice core analysis of Martínez-Botí et al. (2015) is slightly smaller than our results based on the cold periods from the ice core data set (Fig. 9).This indicates that the information which is relevant to suggest any state dependency in S [CO 2 ,LI] is mainly contained in data covering the so-called "warm" climates of the Pleistocene.Thus, especially the land ice area distribution and R [LI] from intermediate glaciated states are important here.However, it should be emphasized that Martínez-Botí et al. (2015) never attempted to detect any state dependency in S [CO 2 ,LI] within either the Pleistocene or the Pliocene data sets.In searching for non-linearities in the scattered data of T g versus R [CO 2 ,LI] by statistical methods, we here go beyond their approach.

Martínez
Comparing data-based estimates of S [CO 2 ,LI] directly with climate model results (e.g.Lunt et al., 2010) is not straightforward, and it is not done in the following because in climate models only those processes considered explicitly as forcing will have an impact on calculated temperature change, while the data-based temperature reconstruction contains the effect of all processes (PALAEOSENS-Project Members, 2012).Furthermore, in Fedorov et al. (2013) climate simulation results have been discussed to understand which processes and mechanisms were responsible for the spatially very heterogeneous changes observed during the last 5 Myr, e.g. the increase in the polar amplification factor over time.Since the results of Fedorov et al. (2013) were unable to explain all observations, it was concluded that a combination of different dynamical feedbacks is underestimated in the climate models.We are not able to generate spatially explicit results.However, from our analysis we could conclude that the equilibrium climate sensitivity represented by S [CO 2 ,LI] was a function of background climate state and probably changed dramatically between conditions with and without Northern Hemisphere land ice.
The contribution of greenhouse gas radiative forcing and of seasonally and latitudinally variable incoming solar radiation to the simulated global temperature anomalies of the last eight interglacials have been analysed individually before (Yin and Berger, 2012).It was found that the greenhouse gas forcing was the main driver of the simulated temperature change with the incoming solar radiation amplifying or dampening its signal for all but one interglacial (Marine Isotope Stage (MIS) 7), with two interglacials (MIS 1 and MIS 19) having variations close to zero.Furthermore, they calculated the ECS (temperature rise for a doubling of CO 2 ) for the different interglacial background conditions and found ECS to decrease with increasing background temperature.A calculation of climate sensitivity for individual points in time has been performed before (PALAEOSENS-Project Members, 2012) but has been rejected due to large uncertainties, mainly during interglacials since in the definition of S, one then needs to calculate the ratio of two small numbers in T g and R [CO 2 ,LI] , which has typically a low signalto-noise ratio.At first glance, this might seem contrary to our finding of a larger climate sensitivity during late Pleistocene interglacials when compared to late Pleistocene full glacial conditions.However, as mentioned already in the previous paragraph, the comparison of (palaeo)data-based calculations of S with ECS calculated from climate models is not directly possible.Furthermore, in our approach we include changes in land ice sheet (albedo forcing or R [LI] ), while Yin and Berger (2012) kept ice sheets at present state.When investigating S [CO 2 ,LI] over the whole range of climate states (from full glacial conditions to a warm Pliocene with a (nearly) ice-free Northern Hemisphere resulting in a variable forcing term R [LI] ), we therefore probe a completely different climate regime, which is not directly comparable with results obtained from simulations of interglacials only.
-7 There exist some intrinsic uncertainties in our approach based on the underlying data sets, which are not included in the Monte Carlo statistic.For example, the global temperature anomaly in the LGM still disagrees in various approaches (Annan and Hargreaves, 2013;Schmittner et al., 2011;Schmidt et al., 2014), and Pliocene sea level and ice-sheet dynamics are still a matter of debate (Rohling et al., 2014;Dolan et al., 2015;Koenig et al., 2015;Rovere et al., 2014;de Boer et al., 2015).Taking these issues into account might lead to changes in our quantitative estimates but not necessarily to a revision of our main finding of state dependency in S [CO 2 ,LI] .In the light of the existing uncertainties, our findings must be supported by other modelling approaches to come to firm conclusions.Furthermore, our assumption that we can estimate equilibrium climate sensitivity from palaeodata implicitly assumes that these data represent predominately equilibrium climate states.This might be a simplification, but since filtering out data points in which temperature changed abruptly led to similar results (PALAEOSENS-Project Members, 2012), it should only have a minor effect on the conclusions.
To calculate in detail the effect of climate change on temperature, it would be important to also include other forcing agents, e.g.CH 4 , N 2 O, or aerosols.For the Pliocene, strong chemistry-climate feedbacks have been proposed (Unger and Yue, 2014), suggesting high ozone and aerosol levels and potentially high CH 4 values.This implies that the relationship of CO 2 to other forcing agents might have been different for the cold climates of the late Pleistocene than for the warm climates of the Pliocene.Therefore, assumptions on the influence of other slow feedbacks based on data of the late Pleistocene (Köhler et al., 2010) cannot be extrapolated to the Pliocene.Hence, we restrict our analysis of the Pliocene data to S [CO 2 ,LI] and again emphasize that an estimate of climate sensitivity for the present day, S a , from our palaeo sensitivity (PALAEOSENS-Project Members, 2012) is not straightforward, especially for these data.
For the Pleistocene data we might roughly approximate the implications of our findings for equilibrium temperature changes under CO 2 doubling, or ECS, by considering the so far neglected feedbacks (CH 4 , N 2 O, aerosols, or vegetation).However, we are aware that this is a simplification, since it was already shown that the per unit radiative forcing climate effect of well-mixed greenhouse gases and aerosols differs (Shindell, 2014).In palaeodata of the last 0.8 Myr, the equilibrium climate sensitivity considering all feedbacks was only about two thirds of S [CO 2 ,LI] (PALAEOSENS-Project Members, 2012).A CO 2 doubling would then lead to an equilibrium rise in global temperature of, on average, 2.5 K (68 % probability range: 2.0-3.5 K) or to, on average, 3.7 K (68 % probability range: 2.5-5.5 K) during Pleistocene full glacial climates ("cold") or Pleistocene "warm" climates (intermediate glaciated to interglacial conditions), respectively.Both average values of ECS are well within the range proposed by palaeodata and models so far (PALAEOSENS-Project Members, 2012;Stocker et al., 2013), but we especially emphasize the potential existence of a long tail of S [CO 2 ,LI] towards higher values.Such estimates of ECS are very uncertain due to the different effect of various forcings and are not yet possible for Pliocene climate states (see above).These long-term temperature change estimates for a doubling of CO 2 are mainly of interest for model validation.To be applicable to the not so distant future, these equilibrium estimates need to be corrected for oceanic heat uptake to calculate any transient temperature response (Zeebe, 2013).Whether climate in the future is more comparable to the climate states of interglacials of the late Pleistocene or to the warm Pliocene is difficult to say, although this has, according to our results, major implications for the expected equilibrium temperature rise.The Greenland ice sheet might completely disappear (Levermann et al., 2013) in the longterm for the projected future greenhouse gas emissions, but it might reduce its ice volume in the next 2000 years by less than 50 %.Another study (Loutre and Berger, 2000) suggests that the Greenland ice sheet might also disappear in the long run for atmospheric CO 2 concentrations between 200 and 300 pmmv.These studies suggest that for the coming millennia, the Earth might still contain a significant amount of P. Köhler et al.: State dependency of the equilibrium climate sensitivity northern hemispheric land ice, and thus the climate and the proposed climate sensitivity S [CO 2 ,LI] will probably be more comparable to the interglacials of the late Pleistocene.In the more distant future, the Northern Hemisphere may become free of land ice and its climate and climate sensitivity more comparable to the warm Pliocene.
When compared to the two most recent contributions to this topic (von der Heydt et al., 2014;Martínez-Botí et al., 2015), our study goes beyond them by four improvements that have been laid out in detail in the introduction.The most important of these improvements is the systematical detection of state dependency in S [CO 2 ,LI] using Monte Carlo statistics.However, we have been able to extend the finding of state dependency in S [CO 2 ,LI] from the ice core data of the last 800 kyr to the last 2.1 Myr.Furthermore, the improvements in the underlying time series of R [LI] have been important to obtain a data set in which the state dependency S [CO 2 ,LI] can be detected.The role of the T g time series seems at first glance to be of similar importance to that of R [LI] .However, state dependency in S [CO 2 ,LI] was also obtained for the alternative temperature time series T g2 or T g3 , and therefore a detailed knowledge of T g is of minor importance for our overall conclusions.

Conclusions
In conclusion, we find that the specific equilibrium climate sensitivity based on radiative forcing of CO 2 and land ice albedo, S [CO 2 ,LI] , is state-dependent if CO 2 data from ice cores or from the Hönisch lab, based on δ 11 B, are analysed.The state dependency arises from the non-linear relationship between changes in radiative forcing of land ice albedo, R [LI] , and changes in global temperature.Previous studies were not able to detect such a state dependency because land ice albedo forcing was not based on results from 3-D ice-sheet models, which contain much of this nonlinearity.So far, the state dependency of S [CO 2 ,LI] based on ice core CO 2 , which was derived from predominately glacial conditions of the late Pleistocene, can be extrapolated to the last 2.1 Myr.During intermediate glaciated and interglacial periods of most of the Pleistocene, S [CO 2 ,LI] was, on average, about ∼ 45 % higher (mean: 1.54 K W −1 m 2 ; 68 % probability range: 1.0-2.2K W −1 m 2 ) than during full glacial conditions of the Pleistocene (mean 1.06 K W −1 m 2 ; 68 % probability range: 0.8-1.4K W −1 m 2 ).Before 2.1 Myr BP the published CO 2 data are too sparse, depend on the applied methodology, and have uncertainties that are too large to come to a statistically well-supported conclusion on the value of S [CO 2 ,LI] .The data available so far suggest that the appearance of northern hemispheric land ice sheets changed the climate system and accordingly influenced climate sensitivity.In the Pliocene, S [CO 2 ,LI] was therefore probably smaller than during the interglacials of the Pleistocene.

Figure 1 .
Figure 1.Radiative forcing of land ice sheets averaged for latitudinal bands of 5 • .Panel (a): annual mean insolation at the top of the atmosphere, I TOA , based on orbital variations (Laskar et al., 2004).Panel (b): fraction of each latitudinal band of 5 • covered by land ice as simulated by the 3-D ice-sheet model ANICE (de Boer et al., 2014).Panel (c): calculated radiative forcing of land ice sheets, R [LI] , normalized to global-scale impact.

Figure 2 .
Figure 2. Calculating global surface temperature change, T g .Panel (a): polar amplification factor, f pa , the ratio between Northern Hemisphere (NH) land temperature change, T NH , and global temperature change, T g , as a function of time based on values for LGM (blue square) and mid-Pliocene warm period (mPWP) (red circle) derived from the model intercomparison projects (MIPs) PMIP3-CMIP5 and PlioMIP(Haywood et al., 2013), respectively.In our standard application, T g1 (black line), f pa is calculated as a linear function depending on northern hemispheric temperature change, T NH (insert), inter-and extrapolated between these two PMIP3 and PlioMIP-based values.Alternatively ( T g2 , orange line), f pa varies as a step function with high values for the Pleistocene (periods with Northern Hemisphere land ice sheets) and low values for the Pliocene (periods mainly without NH land ice sheets) with the step between both values occurring at 2.82 Myr BP, when our results indicate large changes in NH land ice.In T g3 (blue line), f pa varied freely to meet T g reconstructed for LGM by PMIP3 (−4.6 K) and for the mPWP by PlioMIP (+2.7 K).See methods for further details.Panel (b): NH temperature change, T NH , as deconvolved from the benthic δ 18 O stack LR04(Lisiecki and Raymo, 2005) by applying a 3-D ice-sheet model in an inverse mode(de Boer et al., 2014).Uncertainty in T NH (grey) is the 1σ error calculated from eight different model realizations(de Boer et al., 2014).Panel (c): global surface temperature change, T g , as used here based on T g = T NH • f −1 pa .Results for T g based on all three approaches for f pa are given (same colour code as in subfigure a).Symbols show T g ± 1σ as derived within PlioMIP (mPWP, red circle) and PMIP3-CMIP5 (LGM, blue square).Red vertical lines mark the time period of the mPWP.

2. 4
How to calculate the specific equilibrium climate sensitivity, S [CO 2 ,LI]

Figure 4 .
Figure 4. Changes in temperature and radiative forcing over the last 5 Myr.Panel (a): global mean surface temperature change, T g , calculated with the polar amplification factor, f pa , being a linear function of the Northern Hemisphere land temperature change, T NH .Marked are the mid-Pliocene warm period (mPWP) (red horizontal bar), global warming calculated within PlioMIP (red circle), and global cooling during the LGM derived from PMIP3 and CMIP5 (blue square).Panel (b): changes in radiative forcing based on atmospheric CO 2 ( R [CO 2 ]).CO 2 data are from ice cores(Bereiter et al., 2015) and based on δ 11 B (Hönisch lab(Hönisch et al., 2009), Foster lab(Foster, 2008;Martínez-Botí et al., 2015)) and on alkenones(Pagani lab (Pagani et al., 2010;Zhang et al., 2013)).Panel (c) shows radiative forcing of land ice, R [LI] , and, for comparison, global annual mean insolation changes due to orbital variation, R [orbit] .Panel (d): the sum of the radiative forcing changes due to CO 2 and land ice sheets ( R [CO 2 ,LI] ) whenever CO 2 data allow its calculation.Uncertainties show 1σ .
Figure 4. Changes in temperature and radiative forcing over the last 5 Myr.Panel (a): global mean surface temperature change, T g , calculated with the polar amplification factor, f pa , being a linear function of the Northern Hemisphere land temperature change, T NH .Marked are the mid-Pliocene warm period (mPWP) (red horizontal bar), global warming calculated within PlioMIP (red circle), and global cooling during the LGM derived from PMIP3 and CMIP5 (blue square).Panel (b): changes in radiative forcing based on atmospheric CO 2 ( R [CO 2 ]).CO 2 data are from ice cores(Bereiter et al., 2015) and based on δ 11 B (Hönisch lab(Hönisch et al., 2009), Foster lab(Foster, 2008;Martínez-Botí et al., 2015)) and on alkenones(Pagani lab (Pagani et al., 2010;Zhang et al., 2013)).Panel (c) shows radiative forcing of land ice, R [LI] , and, for comparison, global annual mean insolation changes due to orbital variation, R [orbit] .Panel (d): the sum of the radiative forcing changes due to CO 2 and land ice sheets ( R [CO 2 ,LI] ) whenever CO 2 data allow its calculation.Uncertainties show 1σ .

Figure 5 .
Figure 5. Comparing the calculation of radiative forcing of land ice sheets for the last 20 kyr for two different ice-sheet set-ups.Left: the 3-D ice-sheet model ANICE used in this study (de Boer et al., 2014); right: calculation based on 1 • × 1 • model output from ICE-5G (Peltier, 2004), results for radiative forcing of land ice sheets, R [LI] , then based on similar aggregation to latitudinal bands of 5 • as for ANICE.Panels (a, d): annual mean insolation at the top of the atmosphere, I TOA , based on orbital variations (Laskar et al., 2004).Panels (b, e): fraction of each latitudinal bands of 5 • covered by land ice as simulated by the 3-D ice-sheet models.Panels (c, f): calculated radiative forcing of land ice sheets, R [LI] , normalized to global-scale impact.

Figure 6 .
Figure 6.Details on land ice albedo forcing ( R [LI] ).Panel (a): scatter plot of sea level change (de Boer et al., 2014) against land ice albedo forcing (this study), R [LI] , based on the 3-D ice-sheet model ANICE.Data are approximated with a third-order non-linear fit.For comparison, a fit based on sea level change as applied in other applications(Hansen et al., 2008;Martínez-Botí et al., 2015) is shown as dashed line.Furthermore, for Termination I (T-I), results based on ANICE and on ICE-5G(Peltier, 2004) are compared.Panels (b, c): relationship between global surface temperature change, T g , and land ice albedo forcing, R [LI] , for different set-ups.Results plotted over the whole last 5 Myr (one data point every 2 kyr).Panel (b): standard set-up with T g = T g1 = T NH • f −1pa using a polar amplification, f pa , that varies linearly as a function of T NH .R [LI] as based on 3-D ice-sheet models as calculated in this study (see Fig.1c).Panel (c): set-up with a constant f pa = 2.7 as applied previously in van deWal et al. (2011).R [LI] is based on 1-D ice-sheet model results and is calculated from sea level change with 0.0308 W m −2 per metre of sea level change.Underlying 1-D ice-sheet model results of T NH and sea level have been published previously(de Boer et al., 2010) and used elsewhere(Martínez-Botí et al., 2015).

PFigure 7 .
Figure 7. Scatter plots of data of global temperature change, T g , against radiative forcing R [X] .T g is calculated with the polar amplification factor, f pa , being a linear function of T NH .Left column (a, c, e, g): radiative forcing of CO 2 ( R [CO 2 ] ).Right column(b, d, f, h): radiative forcing of CO 2 and land ice albedo ( R [CO 2 ,LI] ).Lines show average best fits (first-, second-, or thirdorder polynomials) to 5000 Monte Carlo realizations of the data (details in Table1).Subfigures differ by the CO 2 data they are based on: (a, b) -ice cores(Bereiter et al., 2015); (c, d)δ 11 B from the Hönisch lab(Hönisch et al., 2009); (e, f)δ 11 B from the Foster lab(Foster, 2008;Martínez-Botí et al., 2015); (g, h) -alkenones from the Pagani lab(Pagani et al., 2010;Zhang et al., 2013).Each row contains information on the number of data points n; each subfigure shows the mean uncertainty of the fit by dividing χ 2 (the weighted sum of squares from the regression analysis) by n and the correlation coefficient r 2 .Uncertainties show 1σ .

Figure 9 .
Figure 9. Probability density function of different approaches to calculate specific equilibrium climate sensitivity, S [CO 2 ,LI] .Results of this study are based on pointwise analysis of the ice core (last 0.8 Myr) and the Hönisch (last 2.1 Myr) data for "cold" periods ( R [CO 2 ,LI] < −3.5 W m −2 ) and "warm" periods ( R [CO 2 ,LI] > −3.5 W m −2 ).In von der Heydt et al. (2014), a similar split of the ice core data was performed.We show their results based on similar T g as obtained here published in the Supporting Information in von der Heydt et al. (2014).Martínez-Botí et al. (2015) calculated S [CO 2 ,LI] either for ice core data of the whole last 0.8 Myr or based on δ 11 B for 0.8 Myr of the Pliocene between 2.5-3.3Myr BP.Vertical lines and labels give the mean of the different results.
P.Köhler et al.:State dependency of the equilibrium climate sensitivity State dependency of the equilibrium climate sensitivity cated in the eastern equatorial Atlantic.The data go back to 2.1 Myr BP and agree favourably with the ice core CO 2 during the last 0.8 Myr.c. CO 2 data from the Foster lab kyr BP).b.CO 2 based on δ 11 B isotopes measured on planktonic shells of G. sacculifer from the Hönisch lab (Hönisch et al., 2009) (n = 52) is obtained from ODP668B lo-P.Köhler et al.:

Table 1 .
Fitting a linear or a non-linear function to the data: 5000 Monte Carlo-generated realizations of the scatteredT g -R [CO 2 ] or T g -R [CO 2 ,LI]were analysed.The data are randomly picked from the entire Gaussian distribution described by the 1σ of the given uncertainties in both T g and R [X] .The parameter values of fitted polynomials are given as mean ±1σ uncertainty from the different Monte Carlo realizations.Data sets differ in the underlying T g and CO 2 data.T g : either T g or polar amplification, f pa , are fixed at LGM and mPWP at values from PMIP3 and PlioMIP with different functionality for f pa (see methods for details).CO 2 data from ice cores and Hönisch, Foster, and Pagani labs.T g1 : based on a fixed polar amplification factor, f pa , at LGM and mPWP, with f pa being a linear function of T NH elsewhere analysing T g vs. R [CO 2 ] ).P.Köhler et al.:State dependency of the equilibrium climate sensitivity g2 : based on a fixed f pa that steps from its mPWP value to its LGM value at 2.82 Myr BP analysing T g vs. R [CO 2 ] g3 : fixed T g at LGM and mPWP, based on f pa being a linear function of T NH elsewhere analysing T g vs. R [CO 2 ] Sensitivity analysis 2: investigating the importance of T g , CO 2 , and R[LI]in the data set from ice cores with respect to vdH2014 Data at 2 kyr intervals (if available) [CO 2 ,LI] scatter plot revealed that both the ice core CO 2 and the Hönisch lab data contain a state dependency in S [CO 2 ,LI] .As explained in Sect.2.4, we analyse the mean and uncertainty in S [CO 2 ,LI] for both data sets from probability density functions for different background climate states represented by R [CO 2 ,LI] Figure 10.Best-guess 3.3 Myr scatter plot of global temperature change, T g , against the radiative forcing of CO 2 and land ice albedo ( R [CO 2 ,LI] ).The Hönisch lab (Hönisch et al., 2009) data for the last 2.1 Myr (most of the Pleistocene) and the Pliocene part of the Foster lab data (Martínez-Botí et al., 2015) -entire data (2.5-3.3Myr BP) and only the data for the almost land-ice-free Northern Hemisphere times (2.82-3.3Myr BP) -are compiled to illustrate how the functional dependency between T g and R [CO 2 ,LI] changed as function of background climate state.