Technical Note : Calculating state dependent equilibrium climate sensitivity from palaeodata

The evidence from both data and models indicate that specific equilibrium climate sensitivity S[X] — the global annual mean surface temperature change (∆Tg) as a response to a change in radiative forcing X (∆R[X]) — is state dependent. Such a state dependency implies that the best fit in the scatter plot of ∆Tg versus ∆R[X] is not a linear regression, but for instance a higher order polynomial. While for the conventional linear case the slope (gradient) of the regression is correctly interpreted as the 5 specific equilibrium climate sensitivity S[X], the interpretation is not straightforward in the non-linear case. We here elaborate how such a state dependent scatter plot needs to be interpreted, and provide a theoretical understanding how to calculate S[X] in the non-linear case.


Introduction
One prominent approach to calculate (specific) equilibrium climate sensitivity S [X] from palaeodata is to evaluate the regression of scatter plots, in which global mean surface temperature change (∆T g ) has been plotted against radiative forcing changes (∆R [X] ), since following its definition, S [X] is easily obtained from the slope of a linear regression line, that needs to pass through the origin to avoid any biases.
Passing through the origin (∆T g = 0 K; ∆R [X] = 0 W m 2 ) implies that no temperature change (with respect to a defined reference climate state) is detected for conditions without forcing anomalies.Usually the pre-industrial climate state serves as reference.Here, X corresponds to the forcing processes considered, typically changes in CO 2 (sometimes also including the other greenhouse gases CH 4 and N 2 O), potentially corrected for some slow feedbacks such as planetary albedo changes causes by variations in land ice (LI), vegetation (VE) or dust (aerosols (AE)).This nomenclature above follows the suggestion of a recent review on this topic (PALAEOSENS-Project Members, 2012).
1 In the review of the PALAEOSENS group in 2012 a quantitative expression of S [X] based on Equation 1 was already included, but only for individual data points, or whole time series, and a state dependent character of S [X] was already detected, but could not be quantified in greater detail.Since then climate sensitivity from palaeodata has continued to be analysed by some regression analysis in the scatter plot of ∆T g versus ∆R [X] (e.g.von der Heydt et al., 2014;Martínez-Botí et al., 2015;Köhler et al., 2015).
The analysis of the problem is straightforward, if linear regression methods are applied, which implies that S [X] is a general property of the climate system.However, results point more and more in the direction, that climate sensitivity is state dependent, implying that for the ∆T g -∆R [X] -scatter plot non-linear regressions (for example higher order polynomials) are describing the data more appropriate than a simple linear fit (von der Heydt et al., 2014;Köhler et al., 2015).In such cases the quantification of S [X] becomes more intricate.
The aim of this technical note is to briefly highlight potential pitfalls when analysing S [X] from a data set, that indicates a state dependent behaviour, and offer understanding and solutions.These simple, but essential thoughts were not contained in our most recent contribution (Köhler et al., 2015).The key point of Köhler et al. (2015) was to show that given the available data climate sensitivity during the last 2.1 Myr is state dependent.One of the possible ways to quantify this state dependency is to provide a probability density function or PDF (approach I, contained in subsection 2.1).Here we generalise different methods to quantify climate sensitivity in case it is state dependent.Hence this technical note offers a more general understanding of state dependent climate sensitivity based on palaeodata.

Explaining different methods
We will in the following briefly expand on different ways how data can be analysed.Approaches are compared for the most simplistic case of two data points only (Fig. 1), but we will for illustrative purposes also apply the different methods to the ice core data of the last 800 kyr already investigated in detail in Köhler et al. (2015).Finally, we also also include data of the last 2.1 Myr based on only CO 2 proxy data from the Hönisch-lab (Hönisch et al., 2009) in our conclusions.Please refer to Köhler et al. (2015) for the origin of the data and details of the applied statistical methods by which the data have been analysed.Using a Monte-Carlo approach and F -test statistics a 3rd-order polynomial has been identified to best fit the scattered data of ∆T g against ∆R [CO2,LI] (Fig. 2c).

Approach I: The point-wise approach
The easiest and most robust estimate of S [CO2,LI] can be obtained, when S [CO2,LI] is calculated individually for each time step t i out of ∆T g (t i ) and ∆R [CO2,LI] (t i ) (Fig. 2d).Taking the uncertainties of the individual data points into account a probability density functions (PDF) of S [CO2,LI] is calculated straightforward (Fig. 2f), from which the median, the most likely and the spread (uncertainty distribution) within S [CO2,LI] can be obtained.If the underlying data set of the PDF is split into various subsets (here distinguishing data for two different radiative forcing domains, Fig. 2e) a first, rough quantification of the state dependency of S [CO2,LI] is generated, and has already been obtained in Köhler et al. (2015).One known problem of this approach is that for small disturbances in the radiative forcing (∆R [CO2,LI] close to zero) one might obtain in the pointwise calculations of S [CO2,LI] unrealistically high and low values.Such values might be caused by dating uncertainties of the underlying palaeorecords or transient effects (de Boer et al., 2012).In our analysis in Köhler et al. (2015) we found 20 outliers (from 394 data points contained in Fig. 2) that did not match in the realistic range of S [CO2,LI] between 0 and 3 K W −1 m 2 , and they have been rejected from further analysis.Furthermore, from those data with ∆R [CO2,LI] close to zero, that have not been rejected, calculated values of S [CO2,LI] seemed to follow a different pattern than the rest of the data (Fig 2e).Again, we understand these anomalies to be probably based on dating uncertainties, non-negligible influence of transient climate response and the problem that the ratio from two small numbers might easily contain a large error.
A state dependent equation of the specific equilibrium climate sensitivity S [X] might be obtained from analysing the palaeodata in greater detail.To do so, a function has to be found that relates the temperature perturbations ∆T g to radiation perturba- . For reasons of simplicity both variables are described in the following by ∆T and ∆R.For non-linear description of ∆T as a function of ∆R a higher order polynomial is the most obvious choice (but other equations are possible): Climate sensitivity can then be calculated as: This approach is called point-wise (pw) since the derived equation in the S − ∆R data space agrees with the individual data points (Figs. 2e,3b).The reference climate has to be chosen such that a = 0, to ensure finite climate sensitivity at ∆R = 0.In the above case, climate sensitivity is constant (i.e.not state dependent) if the higher order terms are zero: c = d = 0. Otherwise, climate sensitivity is function of the radiation perturbation and therefore state dependent.This approach has been applied for the ice core data and is contained in Fig. 2e.

Approach II: Using local slopes (piece-wise linear analysis)
In the constant case (c = d = 0), climate sensitivity can also be found by taking the local slope of the function, therefore called However, in the state dependent case: Now S local (Equation 5) is evidently not equal to the point-wise-calculated climate sensitivity S pw (Equation 3).For illustrative purposes we have included some realisation of S based on local slopes in Fig. 3b.Clearly, they disagree from results obtained with the point-wise approach.Indeed, the suggested equations do not meet the scattered data of The condition for state dependency, however, remains the same: the higher order terms have to be non-zero.For the local slope case, this means that the slope is non-constant.

Combining data-based approaches and model results
Climate models usually perturb a reference climate {R 0 ; T 0 }, and end up with a new climate {R 1 ; T 1 }.They consider climate sensitivity in the following way: One might argue, that only radiative forcing anomalies are of interest, and not the absolute radiative forcings R 0 and R 1 , so the denominator in Equation 6 should be ∆R 1 .For the sake of generalisation we keep Equation 6as is, but the reader will see below (Equation 9) how relevant this formation might be.
Equation 2, the higher order polynomial fit of temperature versus forcing, is also valid, if based on absolute values in T and R, so one might fit in a model-output T versus R scatter plot: Using this approach would imply (using Eq. 7, here reduced for simplicity to a 2nd order polynomial): Now, if the reference climate (T 0 , R 0 ) happens to be the pre-industrial reference climate (which is not always the case), we with b = (b + 2cR 0 ).Equation 9 is equal to approach I, the point-wise quantification of climate sensitivity (S pw ), but the parameter b of the non-linear regression is depending on the reference climate, in detail the radiative forcing R 0 .
More generally S model can be obtained from S local , by the following equation: The benefit of using this local slope approach is that, in contrast to the point-wise approach, the sensitivity is not related to an arbitrarily chosen reference climate (often taken to be pre-industrial climate); it represents the slope of temperature as a function of radiation perturbations at each point.S local is therefore more readily comparable to climate model results.In using this approach though, one has to realise that in the state dependent case S local is not equal to the point-wise climate sensitivity S pw .

Discussion
When using PDFs approach I, each data point is weighted with equal weight.This is different from approaches in which regression functions are applied.For example, in linear regressions which would be applied for the state independent case, the frequently applied regression method of ordinary least squares (OLS) minimises the sum of squared residuals.For the pre-conditioned case (regression lines in the ∆T g − ∆R [X] -space have to pass through the origin, thus a = 0) only the slope of the regression line is determined, that is also equivalent to S [X] .We briefly illustrate the fact that in such linear regressions data points further away from the origin get a higher weight by having a closer look on the most simplistic case with only 2 data points, i = 1, 2 with point i: ∆T i , ∆R i , S i = ∆T i /∆R i OLS calculates a mean slope or S after from which one can easily determine that the squares in the equation lead to a higher weight for data points further away from the origin.
One interesting effect of the application of regression statistics (being either linear, or non-linear) is the fact that the uncertainties in the parameters values to this higher order polynomial are very small, compared to the more general uncertainties (width of the PDF) in approach I.

Conclusions
We find the following conclusions: 1.In case no state dependency in S [X] is found the data in the ∆T g − ∆R [X] can be analysed by linear regressions to derive the slope that is after S . Be aware that approaches that combine point-wise-derived values of S [X] in a probability density function to a more general number treat all individual data points with the same weight, while data far away from the origin get a higher weight in linear regression analysis using OLS. 3. Taken at face value, the more sophisticated quantification of S [CO2,LI] as a function of radiative forcing perturbation ∆R [CO2,LI] obtained here from the data sets described in Köhler et al. (2015) leads to numbers, which are by a factor of about 2 to 2.7 higher during climate conditions representing interglacials of the Pleistocene (last 2.1 Myr) than during full glacials during this period (Fig. 4).

4.
A previous approach based on first-order local slopes (von der Heydt et al., 2014) already suggested higher S [CO2,LI] during interglacials than during glacials for data of the last 800 kyr, though only by 40%, while other approaches (Martínez-Botí et al., 2015) did not considered state dependency within the data set covering the last 800 kyr.These previous results agree more with the results we obtain for full glacial conditions.This might be the case, because linear regressions might not be forced through the origin (parameter a was not necessarily zero).When comparing results based on PDFs with the other approaches (as done in Figure 9 in Köhler et al. (2015)) the difference might also be explained because in linear regressions data points further away from the origin get a higher weight.The reason of this cold bias is not necessarily caused by the fact that most data are available for cold climates, as data may be binned (von der Heydt et al., 2014;Köhler et al., 2015).
5. Interglacial climate conditions (∆R = 0 W m −2 , Fig. 4) have a specific equilibrium climate sensitivity S [CO2,LI] between 2.0 K W −1 m 2 (data of the last 2.1 Myr based on Hönisch-lab CO 2 -proxies) and 2.7 K W −1 m 2 (ice core CO 2 data of the last 800 kyr).They are clearly at the upper end of ranges of S [CO2,LI] reported so far from various data sets of the last 65 Myr, e.g. in PALAEOSENS-Project Members (2012), Figure 3, the 95% probability in S [CO2,LI] ranged from 0.48 to 1.91 K W −1 m 2 .Transfering these results based on data of the last Pleistocene to the near future remains difficult, since these palaeodata cover mainly conditions with negative radiative forcing anomaly, while for the future positive radiative forcing anomalies related to a rise in CO 2 are of interest, which obstructs a comparison as we showed that S [CO2,LI] depends on ∆R [CO2,LI] .Furthermore, S [CO2,LI] needs to be corrected by factors related to fast and slow feedbacks to derive S a , the actual or Charney climate sensitivity that is also constrained with climate models (see Equation 2 in PALAEOSENS-Project Members, 2012).However, such corrections and a detailed calculation of S a are not readily available for ∆R [CO2,LI] = 0 W m −2 and a detailed calculation of S a is beyond the scope of this technical note.
Nevertheless does the palaeodata-based analysis suggest that the equilibrium climate sensitivity for present-day based on palaeodata is more at the high end with respect to reported values in the IPCC AR5 report (e.g.Thematic Focus Element 6 in Stocker et al., 2013).-7 individual slopes ( T/ R) average of individual slopes local slopes of 2nd order polynom   LI] caused by CO2 (Bereiter et al., 2015), corrected by land ice albedo feedback (de Boer et al., 2014;Köhler et al., 2015).including the different fits to the data for both approaches I (S pw , Equation 3) and the uncorrected approach II (S slope , Equation 5).Results of approach II agree with results of approach I when they are corrected following Equation 11 Clim.Past Discuss., doi:10.5194/cp-2016-23,2016 Manuscript under review for journal Clim.Past Published: 15 March 2016 c Author(s) 2016.CC-BY 3.0 License.
Clim.Past Discuss., doi:10.5194/cp-2016-23,2016 Manuscript under review for journal Clim.Past Published: 15 March 2016 c Author(s) 2016.CC-BY 3.0 License.with ∆R i = R i − R 0 .This equation is generally valid, and agrees in case of the reference climate being the pre-industrial periods with S pw (approach I).In other words, point-wise climate sensitivity is a measure of the mean local slope climate sensitivity over the radiation perturbation interval.
Clim.Past Discuss., doi:10.5194/cp-2016-23,2016 Manuscript under review for journal Clim.Past Published: 15 March 2016 c Author(s) 2016.CC-BY 3.0 License.2. If state dependency of (specific) equilibrium climate sensitivity is found one has to be careful when quantifying equilibrium climate sensitivity S [X] .The results based on local slopes are not directly comparable to estimates based on models or point-wise calculations, they can however be transferred into each other: the results based on local slopes transfer in results based on models or on the point-wise analysis by calculating the integral over the radiation perturbation interval following Eq.9.

Figure 1 .Figure 2 .
Figure 1.Conceptualised figure how different regressions and slopes might be calculated for two data points only in the ∆Tg-∆R [CO 2 ,LI] data space, regressions are partly (a = 0) forced through the origin.According to the definition S [CO 2 ,LI] = ∆Tg/∆R [CO 2 ,LI] , S [CO 2 ,LI] corresponds to the slope of a straight line passing through the origin.The average of individual slopes (light green) is obtained by unweighted average (comparable to analysis within PDFs), it differs significantly from the linear regression line (dark green, broken line) because regressions are weighting individual results differently.

Figure 2 .
Figure 2. (Previous page) Data of the last 0.8 Myr analysed for specific equilibrium climate sensitivity S [CO 2 ,LI] .In all subplots data are split into data from "warm" (grey) and "cold" (blue) periods distinguished by ∆R [CO 2 ,LI] = −3.5 W m −2 .(a) Time series of radiative forcing

Figure 4 .
Figure 4. Revised calculation of S [CO 2 ,LI] , not only based on PDF distributions, but following the quantification of the point-wise approach I (S pw ).Results based on the ice core data of the last 0.8 Myr and based on the Hönisch data over the last 2.1 Myr (Hönisch et al., 2009) are compared with previous estimates of von der Heydt et al. (2014) (vdHeydt2014) and Martínez-Botí et al. (2015) (M2015).
The data of ∆T g and ∆R [CO2,LI] based on ice core CO 2 and own model-based deconvolution of benthic δ 18 O into ice