Including the efficacy of land ice changes in deriving climate sensitivity from paleodata

The equilibrium climate sensitivity (ECS) of climate models is calculated as the equilibrium global mean surface air warming resulting from a simulated doubling of the atmospheric CO₂ concentration. In these simulations, long-term processes in the climate system, such as land ice changes, are not incorporated. Hence, climate sensitivity derived from paleodata has to be compensated for these processes, when comparing it to the ECS of climate models. Several recent studies found that the impact these long-term processes have on global temperature cannot be quantified directly through the global radiative forcing they induce. This renders the prevailing approach of deconvoluting paleotemperatures through a partitioning based on radiative forcings inaccurate. Here, we therefore implement an efficacy factor ε_[LI] that relates the impact of land ice changes on global temperature to that of CO₂ changes in our calculation of climate sensitivity from paleodata. We apply our refined approach to a proxy-inferred paleoclimate dataset, using ${\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}={\mathrm{0.45}}_{-\mathrm{0.20}}^{+\mathrm{0.34}}$ based on a multi-model assemblage of simulated relative influences of land ice changes on the Last Glacial Maximum temperature anomaly. The implemented ε_[LI] is smaller than unity, meaning that per unit of radiative, forcing the impact on global temperature is less strong for land ice changes than for CO₂ changes. Consequently, our obtained ECS estimate of 5.8±1.3 K, where the uncertainty reflects the implemented range in ε_[LI], is ∼50 % higher than when differences in efficacy are not considered.

1 Introduction

Equilibrium climate sensitivity (ECS) expresses the simulated equilibrated surface air temperature response to an instantaneous doubling of the atmospheric CO₂ concentration. The simulated effect of the applied CO₂ radiative forcing anomaly includes the Planck response, as well as the fast feedbacks such as those involving changes to snow, sea ice, lapse rate, clouds, and water vapour. ECS varies significantly between different state-of-the-art climate models; for instance, the CMIP5 ensemble shows a range of 1.9 to 4.4 K (Vial et al., 2013). Several ways have been put forward to constrain ECS, for example, through the use of paleoclimate data (e.g. Covey et al., 1996; Edwards et al., 2007), which is also the focus of this study. However, unlike results of models, temperature reconstructions based on paleoclimate proxy data always contain a mixed signal of all processes active in the climate system. Among these are long-term processes (or slow feedbacks) such as changes in vegetation, dust, and, arguably most importantly, land ice changes, which are kept constant in the climate model runs used to calculate ECS. Therefore, it is necessary to correct paleotemperature records for the influence of these processes, in order to make a meaningful comparison to ECS calculated by climate models.

In a coordinated community effort, the PALAEOSENS project proposed to relate the temperature response caused by these long-term processes to the globally averaged radiative forcing they induce (PALAEOSENS Project Members, 2012). Consequently, the paleotemperature record can be disentangled on the basis of the separate radiative forcings of these long-term processes (e.g. von der Heydt et al., 2014; Martínez-Botí et al., 2015; Köhler et al., 2015a, 2017b, 2018; Friedrich et al., 2016). If all processes are accounted for in this manner, the effect of CO₂ changes and the accompanying short-term feedbacks, as described by the ECS, can be estimated. However, several studies have shown that, depending on the type of radiative forcing, the same global-average radiative forcing can lead to different global temperature changes (e.g. Stuber et al., 2005; Hansen et al., 2005; Yoshimori et al., 2011). For instance, in a previous article (Stap et al., 2018a), we simulated the separate and combined effects of CO₂ changes and land ice changes on global surface air temperature using the intermediate-complexity climate model CLIMBER-2 and showed that the specific global temperature change per unit radiative forcing change depends on which process is involved. As a possible solution to this problem, Hansen et al. (2005) formulated the concept of “efficacy” factors, which express the impact of radiative forcing by a certain process in comparison to the effect of radiative forcing by CO₂ changes.

Based on the concept of Hansen et al. (2005), here we introduce an efficacy factor for radiative forcing by albedo changes due to land ice variability in our method of deriving climate sensitivity from paleodata. We first illustrate our refined approach by applying it to transient simulations over the past 5 Myr using CLIMBER-2 (Stap et al., 2018a), obtaining a quantification of the effect on global temperature of CO₂ changes and the accompanying short-term feedbacks from a simulation forced by both land ice and CO₂ changes. We compare this result to a simulation where CO₂ changes are the only operating long-term process. In this manner, we can assess the error resulting from using a constant efficacy factor. Thereafter, we refine a previous estimate of climate sensitivity based on a paleoclimate dataset of the past 800 kyr (Köhler et al., 2015a, 2018). In this dataset, the sole effect of CO₂ is not a priori known. We therefore investigate the influence of the introduced efficacy factor on the calculated climate sensitivity. To do so, we appraise the influence of land ice changes and the associated efficacy using a range that is given by different modelling efforts of the Last Glacial Maximum (LGM; ∼21 kyr ago) (Shakun, 2017). The climate sensitivity resulting from applying this range provides a quantification of the consequence of the uncertain efficacy of land ice changes.

2 Material and methods

In this section, we first summarize the approach to obtaining climate sensitivity from paleodata that has been used in numerous earlier studies (e.g. PALAEOSENS Project Members, 2012; von der Heydt et al., 2014; Martínez-Botí et al., 2015; Köhler et al., 2015a, 2017b, 2018; Friedrich et al., 2016). We then discuss our main refinement to that approach, which is the inclusion of the efficacy of land ice changes, and a further small refinement that is meant to unify the dependent variable in cross-plots of radiative forcing and global temperature anomalies.

2.1 Approach to obtain climate sensitivity from paleodata

In climate model simulations used to quantify ECS, fast feedbacks, i.e. processes in the climate system with timescales of less than ∼100 years, are accounted for. However, slower processes, such as those involving changes to ice sheets, vegetation, and dust, are commonly kept constant. The resulting response is also sometimes called Charney sensitivity (Charney et al., 1979). Following the notation of PALAEOSENS Project Members (2012), the ratio of the temperature change ($\mathrm{\Delta }{T}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$) to the radiative forcing due to the CO₂ change ($\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$) yields S^a (in K W⁻¹ m², and where “a” stands for actuo):

$$\begin{array}{}\text{(1)}& {\displaystyle }{S}^{\mathrm{a}}={\displaystyle \frac{\mathrm{\Delta }{T}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}.\end{array}$$

The subscript denotes that CO₂ is the only long-term process involved. Analogously, paleoclimate sensitivity (S^p) can be deduced from paleotemperature reconstructions and paleo-CO₂ records as

$$\begin{array}{}\text{(2)}& {\displaystyle }{S}^{\mathrm{p}}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}.\end{array}$$

In this case, the average global paleotemperature anomaly with respect to the pre-industrial period (PI) (ΔT_g) is, however, also affected by the long-term processes that are typically neglected in climate simulations. Therefore, a correction to the paleotemperature record is needed to obtain $\mathrm{\Delta }{T}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ from ΔT_g:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{\Delta }{T}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{\Delta }{T}_{\mathrm{g}}(\mathrm{1}-f),\end{array}$$

or equivalently S^a from S^p:

$$\begin{array}{}\text{(4)}& {\displaystyle }{S}^{\mathrm{a}}=S^{\mathrm{p}}(\mathrm{1}-f)={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}(\mathrm{1}-f).\end{array}$$

Here, f represents the effect of the slow feedbacks on paleotemperature (e.g. van de Wal et al., 2011). To obtain f, PALAEOSENS Project Members (2012) proposed an approach, which has subsequently been used in numerous studies aiming to constrain climate sensitivity from paleodata (e.g. von der Heydt et al., 2014; Martínez-Botí et al., 2015; Köhler et al., 2015a, 2017b, 2018; Friedrich et al., 2016) and paleoclimate modelling studies (e.g. PALAEOSENS Project Members, 2012; Friedrich et al., 2016; Chandan and Peltier, 2018). They suggested quantifying the influence of the long-term processes (X) by the radiative forcing change they induce (ΔR_[X]), relative to the total radiative forcing perturbation:

$$\begin{array}{}\text{(5)}& {\displaystyle }f={\displaystyle \frac{\mathrm{\Delta }{R}_{\left[X\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+\mathrm{\Delta }{R}_{\left[X\right]}}}=\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+\mathrm{\Delta }{R}_{\left[X\right]}}}.\end{array}$$

Combining Eqs. (4) and (5) and following the PALAEOSENS nomenclature, we can then derive the “specific” paleoclimate sensitivity $S_{[{\mathrm{CO}}_{\mathrm{2}},X]}$, where X represents the processes that are accounted for in the calculation of f:

$$\begin{array}{ll}{\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},X]}& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}(\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[X\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+\mathrm{\Delta }{R}_{\left[X\right]}}})\\ \text{(6)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+\mathrm{\Delta }{R}_{\left[X\right]}}}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{[{\mathrm{CO}}_{\mathrm{2}},X]}}}.\end{array}$$

If, for instance, only the most important slow feedback in the climate system, namely radiative forcing anomalies induced by albedo changes due to land ice (LI) variability, is taken into account, then one can correct S^p to derive the following specific climate sensitivity:

$$\begin{array}{}\text{(7)}& {\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]}}}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}}}.\end{array}$$

Using this approach, several studies performed a least-squares regression through scattered data from paleotemperature and radiative forcing records (Martínez-Botí et al., 2015; Friedrich et al., 2016; Köhler et al., 2015a, 2017b, 2018) relating ΔT_g to $\mathrm{\Delta }{R}_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ in a time-independent manner, from which $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ could be determined. In the course of those studies, a state dependency of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ as function of background climate has been deduced for those data which are best approximated by a non-linear function. Furthermore, the quantification of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ for those state-dependent cases has been formalized in Köhler et al. (2017b). A synthesis of estimates of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ from both colder- and warmer-than-present climates has been compiled by von der Heydt et al. (2016).

2.2 Refinement 1: taking the efficacy of land ice changes into account

The validity of the PALAEOSENS approach to calculating f is contingent on the notion that identical global-average radiative forcing changes lead to identical global temperature responses, regardless of the processes involved. However, it has been demonstrated that the horizontal and vertical distribution of the radiative forcing affects the resulting temperature response (e.g. Stuber et al., 2005; Hansen et al., 2005; Yoshimori et al., 2011; Stap et al., 2018a) because, for example, different fast feedbacks are triggered depending on the location of the forcing. To address this issue, Hansen et al. (2005) introduced the concept of efficacy factors, which we will explore further in this study. These factors (ε_[X]) relate the strength of the temperature response to radiative forcing caused by a certain process X (ΔT_[X]∕ΔR_[X]) to a similar ratio caused by CO₂ radiative forcing ($\mathrm{\Delta }{T}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}/\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$). This introduction of efficacy requires a reformulation of f as f_ε:

$$\begin{array}{ll}{\displaystyle }{f}_{\mathit{\epsilon }}& {\displaystyle }={\displaystyle \frac{{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }=\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}},\end{array}$$

and hence also of $S_{[{\mathrm{CO}}_{\mathrm{2}},X]}$ as $S_{[{\mathrm{CO}}_{\mathrm{2}},X]}^{\mathit{\epsilon }}$:

$$\begin{array}{}\text{(9)}& {\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},X]}^{\mathit{\epsilon }}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}}.\end{array}$$

In these reformulations, where in principle ε_[X] can take any value, we introduce the superscript ε. This serves to clearly distinguish these newly derived sensitivities from those of the PALAEOSENS approach, in which efficacy was not taken into account, implying that identical radiative forcing of different processes leads to identical temperature changes.

To calculate $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$, we constrain the efficacy factor for radiative forcing by land ice changes (ε_[LI]), using the following formulation, which is based on but slightly modified from Hansen et al. (2005):

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{\Delta }{T}_{\left[\mathrm{LI}\right]}}{\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]}}}={\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}{\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}-\mathrm{\Delta }{T}_{\left[\mathrm{LI}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}.\end{array}$$

This leads to

$$\begin{array}{}\text{(11)}& {\displaystyle }{\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}={\displaystyle \frac{\mathit{\omega }}{\mathrm{1}-\mathit{\omega }}}{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]}}},\end{array}$$

where ω represents the fractional relative influence of land ice changes on the global temperature change ($\mathit{\omega }=\mathrm{\Delta }{T}_{\left[\mathrm{LI}\right]}/\mathrm{\Delta }{T}_{\mathrm{g}}$). If ε_[LI] is assumed to be constant in time (see Sect. 3.2), it can be calculated using Eq. (11) from data of any time step in the record of $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ and ΔR_[LI], and consequently applied to the whole record (Fig. 1a, c). As before, with this ε_[LI] a quantification of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ can be obtained by performing a least-squares regression through scattered data from paleotemperature and radiative forcing records, now relating ΔT_g to $(\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]})$ in a time-independent manner.

Figure 1 Time series of radiative forcing anomalies (ΔR) caused by CO₂ (red) changes and land ice changes (blue) and global temperature anomalies (ΔT_g) with respect to PI from (a–b) the CLIMBER-2 model dataset (Stap et al., 2018a), with temperature data for experiment OIC in black and for experiment OC in green, and from (c–d) the proxy-inferred dataset (Köhler et al., 2015a), with solid lines for the whole dataset and dots for the data used in this study which exclude times with strong temperature–CO₂ divergence (see Sect. 4.1). Note the differing axis scales.

Note that apart from the formulation based on Hansen et al. (2005) followed here, other formulations of the efficacy factor are possible. For instance, one can define an alternative efficacy factor (ε_[LI],alt) such that it relates the effect of land ice changes on global temperature directly to the radiative forcing anomaly caused by CO₂ changes, leading to

$$\begin{array}{}\text{(12)}& {\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},X],\mathrm{alt}}^{\mathit{\epsilon }}={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[\mathrm{LI}\right],\mathrm{alt}}\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}.\end{array}$$

In this alternative case, the efficacy factor ε_[LI],alt relates to our original ε_[LI] as

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{\epsilon }}_{\left[\mathrm{LI}\right],\mathrm{alt}}={\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}.\end{array}$$

This implies that if ε_[LI] is indeed constant, any non-linearity in the relation between $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ and ΔR_[LI] would demand a more complex formulation of the alternative efficacy factor ε_[LI],alt (e.g. via a higher-order polynomial). Since we find such a non-linearity in our data (Fig. 2), using an F test to determine that a second-order polynomial is a significantly (p value < 0.0001) better fit to the data than a linear function, we refrain from following this alternative formulation further.

2.3 Refinement 2: unifying the dependent variable

To calculate $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$, previous studies have used cross-plots of global temperature anomalies and radiative forcing. The latter is caused by a combination of CO₂ and land ice changes, which is cumbersome if one wants to compare $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ to other specific paleoclimate sensitivities $S_{[{\mathrm{CO}}_{\mathrm{2}},X]}^{\mathit{\epsilon }}$, where more and/or different long-term processes are considered. Here, we therefore reformulate our quantification of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ to unify the dependent variable as $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$.

$$\begin{array}{ll}{\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},X]}^{\mathit{\epsilon }}& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{g}}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}}\\ \text{(14)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }{T}_{[-X]}^{\mathit{\epsilon }}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}\end{array}$$

Here, $\mathrm{\Delta }{T}_{[-X]}^{\mathit{\epsilon }}$ is the global temperature change (with respect to PI) stripped of the inferred influence of processes X, defined as

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{\Delta }{T}_{[-X]}^{\mathit{\epsilon }}:=\mathrm{\Delta }{T}_{\mathrm{g}}{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[X\right]}\mathrm{\Delta }{R}_{\left[X\right]}}}.\end{array}$$

Hence, for the calculation of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$, we use

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}:=\mathrm{\Delta }{T}_{\mathrm{g}}{\displaystyle \frac{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}+{\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]}}}.\end{array}$$

Now, we quantify $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ by performing a least-squares regression (regfunc) through scattered data from $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$ and $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$. We use the precondition that no change in CO₂ is related to no change in $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$, meaning the regression intersects the y axis at the origin ($(x,y)=(\mathrm{0},\mathrm{0})$). Following Köhler et al. (2017b), for any non-zero $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$, we calculate $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ as

$$\begin{array}{}\text{(17)}& {\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}{|}_{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}={\displaystyle \frac{\mathrm{regfunc}}{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}}{|}_{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}}.\end{array}$$

If $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{0}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, as is, for instance, the case for pre-industrial conditions, $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ is quantified as

$$\begin{array}{}\text{(18)}& {\displaystyle }{S}_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}{|}_{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{0}}={\displaystyle \frac{\mathit{\delta }\left(\mathrm{regfunc}\right)}{\mathit{\delta }\left(\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}\right)}}{|}_{\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{0}}.\end{array}$$

Equations (17) and (18) yield a quantification of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$, which can be compared to the value obtained for $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ using the PALAEOSENS approach that does not consider efficacy differences (equivalent to using ε_[LI]=1) (Köhler et al., 2018).

To obtain S^a, one needs to multiply $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ by a conversion factor $\mathit{\varphi }=\mathrm{0.64}\pm \mathrm{0.07}$ (1σ uncertainty) that accounts for the influence of other long-term processes, namely vegetation, aerosol, and non-CO₂ greenhouse gas changes (PALAEOSENS Project Members, 2012). Note that this multiplication by ϕ ignores any possible state dependencies in ϕ and assumes unit efficacy for processes other than land ice changes. Because a comprehensive analysis of the efficacy and state dependency of these other processes is beyond the scope of this study, it is a source of uncertainty to be investigated in future research. Finally, we obtain the equivalent ECS by multiplying S^a by $\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}=\mathrm{3.71}\pm \mathrm{0.37}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (1σ uncertainty), the radiative forcing perturbation representing a CO₂ doubling (Myhre et al., 1998).

3 Illustration of the approach using model simulations

In this section, we illustrate our refined approach, which considers efficacy differences, by applying it to transient simulations over the past 5 Myr using CLIMBER-2 (Stap et al., 2018a). We obtain a quantification of the effect on global temperature of CO₂ changes and the accompanying short-term feedbacks from a simulation forced by both land ice and CO₂ changes. We compare this result to a simulation where CO₂ changes are the only operating long-term process. By doing so, we assess the error resulting from using a constant efficacy factor.

3.1 CLIMBER-2 model simulations

Currently, long (∼10⁵ to ∼10⁶ years) integrations of state-of-the-art climate models, such as general circulation models and Earth system models, are not yet not feasible due to limited computer power. This gap can be filled by using models of reduced complexity (Claussen et al., 2002; Stap et al., 2017). Using the intermediate-complexity climate model CLIMBER-2 (Petoukhov et al., 2000; Ganopolski et al., 2001), climate simulations over the past 5 Myr were performed and analysed in Stap et al. (2018a). CLIMBER-2 combines a 2.5-dimensional statistical–dynamical atmosphere model, with a three-basin zonally averaged ocean model (Stocker et al., 1992) and a model that calculates dynamic vegetation cover based on the temperature and precipitation (Brovkin et al., 1997). The simulations could be forced by solar insolation changes due to orbital (O) variations (Laskar et al., 2004), by land ice (I) changes in both hemispheres (based on de Boer et al., 2013), and by CO₂ (C) changes (based on van de Wal et al., 2011). In the reference experiment (OIC) all these factors are varied, while in other model integrations the land ice (experiment OC) or the CO₂ concentration (experiment OI) is kept fixed at PI level. The synergy of land ice and CO₂ changes is negligibly small, meaning their induced temperature changes add approximately linearly when both forcings are applied. Furthermore, the influence of orbital variations is also very small, so that experiment OC approximately yields the sole effect of CO₂ changes on global temperature (ΔT_[OC]). As in Stap et al. (2018a), we use the simple energy balance model of Köhler et al. (2010) to analyse the applied radiative forcing of land ice albedo and CO₂ changes and simulated global temperature changes, after averaging to 1000-year temporal resolution (Fig. 1a, b).

3.2 Analysis

First, we analyse experiment OC, which will serve as a target for our approach as deployed later in this section. We use a least-squares regression through scattered data of $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ and ΔT_[OC] to fit a second-order polynomial (Fig. 3a). Using a higher-order polynomial rather than a linear function allows us to capture the state dependency of paleoclimate sensitivity. Fitting even higher-order polynomials leads to negligible coefficients for the higher powers and is not pursued further. From the fit, we calculate a specific paleoclimate sensitivity $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ of 0.74 K W⁻¹ m² for PI conditions ($\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{0}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) using Eq. (18). Note that, in this case, $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ is equal to $S_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}^{\mathit{\epsilon }}$, $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$, and $S_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ as there are no land ice changes and therefore also no efficacy differences. The fit further shows decreasing $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ for rising $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$.

Figure 2 The relation between radiative forcing anomalies caused by CO₂ changes ($\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$) and land ice changes (ΔR_[LI]) from the whole proxy-inferred dataset (Köhler et al., 2015a) (pink dots). The red line represents a second-order polynomial least-squares regression through the scattered data.

Now, we apply our approach to the results of experiment OIC, in which both CO₂ and land ice cover vary over time, with the aim of deducing the sole effect of CO₂ changes on global temperature. We calculate the efficacy of land ice changes for the LGM (21 kyr ago) from experiment OI, in which the CO₂ concentration is kept constant. We obtain $\mathit{\omega }=\mathrm{\Delta }{T}_{\left[\mathrm{LI}\right]}/\mathrm{\Delta }{T}_{\mathrm{g}}=\mathrm{\Delta }{T}_{\left[\mathrm{OI}\right]}/\mathrm{\Delta }{T}_{\left[\mathrm{OIC}\right]}=\mathrm{0.54}$. Consequently, we find ε_[LI]=0.58 from Eq. (11) and apply this value to the whole record of $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ and ΔR_[LI]. In this manner, we calculate $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$ using Eq. (16). In principle, ε_[LI] can be obtained using data from any time step of the record, preferably when the radiative forcing anomalies are large to prevent outliers resulting from divisions by small numbers. For example, using the results from all glacial marine isotope stages of the past 810 kyr (Marine Isotope Stage (MIS) 2, 6, 8, 10, 12, 14, 16, 18, and 20), instead of only the LGM, leads to a mean (±1σ) ε_[LI] of 0.56±0.09.

We then fit a second-order polynomial to the scattered data of the thusly obtained $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$ from the results of experiment OIC and $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ (Fig. 3b, c). Between $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=-\mathrm{0.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{0.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, outliers resulted from division by small numbers (not shown in Fig. 3b). To remove these outliers, we first calculate the root mean square error (RMSE) between the fit and the data in the remainder of the domain. We then exclude all values from the range $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=-\mathrm{0.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ to $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{0.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, where the fit differs from the data by more than 3 times the RMSE, and perform the regression again. This yields an $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ of 0.72 K W⁻¹ m² for PI (Fig. 3b) in the LGM-only case and ${\mathrm{0.73}}_{-\mathrm{0.05}}^{+\mathrm{0.06}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ in the case where all glacial periods are used (Fig. 3c). This supports our approach since it is only slightly lower than the $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ of 0.74 K W⁻¹ m² obtained from experiment OC, which it should approximate. The relationship between $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$ and $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ (Fig. 3b) is less linear than that between ΔT_[OC] and $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ (Fig. 3a); hence, the state dependency of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ is enhanced. However, the difference between the $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ obtained from both experiments remains smaller than 0.07 K W⁻¹ m² throughout the entire 5 Myr interval in the LGM-only case, indicating that a constant efficacy is an acceptable assumption which only introduces a negligible additional uncertainty. However, the possible time dependency of efficacy could be investigated more rigorously in future research using more sophisticated climate models.

Figure 3 Temperature anomalies with respect to PI over the last 5 Myr from CLIMBER-2 (Stap et al., 2018a) against imposed radiative forcing of CO₂. (a) Simulation with fixed PI land ice distribution (experiment OC) (ΔT_[OC]). (b) Calculated global temperature perturbations from experiment OIC stripped of the inferred influence of land ice ($\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$) using Eq. (16) with ε_[LI]=0.58. Here, ε_[LI] is obtained from matching climate sensitivity with the target value at the LGM. (c) Same as in panel (b) but using ε_[LI]=0.47 (cyan dots), ε_[LI]=0.56 (pink dots), and ε_[LI]=0.65 (yellow dots), Here, ε_[LI] is obtained from the mean (±1σ) of matching climate sensitivity with the target value at all glacial marine isotope stages of the past 810 kyr (MIS 2, 6, 8, 10, 12, 14, 16, 18, and 20). (d) Same as in panel (b) but using ε_[LI]=1, which is equivalent to the PALAEOSENS approach where efficacy differences were not considered. The red lines – and in panel (c) also the orange and blue lines – represent second-order polynomial least-squares regressions through the scattered data.

The PALAEOSENS approach that does not consider efficacy differences (ε_[LI]=1) yields a PI $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ of 0.54 K W⁻¹ m² (Fig. 3d). This is clearly much more off-target than the results of our approach, signifying the importance of considering efficacy.

4 Application to proxy-inferred paleoclimate data

In this section, we compare our refined approach to calculating $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ incorporating efficacy to our previous quantification of $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ (Köhler et al., 2018) by reanalysing the same paleoclimate dataset (introduced in Köhler et al., 2015a). Other than for climate model simulations, in proxy-based datasets the influence of land ice changes on global temperature perturbations cannot be directly obtained and hence is a priori unknown. We therefore base the value of ε_[LI] we implement here on a multi-model assemblage of simulated relative influences of land ice changes on the LGM temperature anomaly (Shakun, 2017).

4.1 Proxy-inferred paleoclimate dataset

The dataset to be investigated contains reconstructions of ΔT_g, $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$, and ΔR_[LI] for the past 800 kyr. Although the dataset covers the past 5 Myr, here we focus only on the past 800 kyr (Fig. 1c, d) because over this period $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ is constrained by high-fidelity measurements of CO₂ within ice cores, whereas Pliocene and Early Pleistocene CO₂ levels are still heavily debated (e.g. Badger et al., 2013; Martínez-Botí et al., 2015; Willeit et al., 2015; Stap et al., 2016, 2017; Chalk et al., 2017; Dyez et al., 2018). Radiative forcing by CO₂ is obtained from Antarctic ice core data compiled by Bereiter et al. (2015), using $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=\mathrm{5.35}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\cdot \ln ({\mathrm{CO}}_{\mathrm{2}}/(\mathrm{278}\mathrm{ppm}\left)\right)$ (Myhre et al., 1998). The revised formulation for $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ from Etminan et al. (2016) leads to very similar results with less than 0.01 W m⁻² differences between the approaches for typical late Pleistocene CO₂ values (Köhler et al., 2017a). Radiative forcing caused by land ice albedo changes, as well as the global surface air temperature record (ΔT_g), are based on results of the 3-D ice-sheet model ANICE (de Boer et al., 2014) forced by northern hemispheric temperature anomalies with respect to a reference PI climate. The ANICE results are here considered to be proxy-inferred because, unlike climate models, ANICE is not constrained by climatic boundary conditions such as insolation and greenhouse gases. The temperature anomalies follow directly from a benthic δ¹⁸O stack (Lisiecki and Raymo, 2005) using an inverse technique. Nevertheless, the results are model-dependent and therefore subject to uncertainty. ANICE provided geographically specific land ice distributions and hence radiative forcing due to albedo changes with respect to PI in both hemispheres. In Köhler et al. (2015a), the northern hemispheric (NH) temperature anomalies (ΔT_NH) are translated into global temperature perturbations (ΔT_g1 in Köhler et al. (2015a)) using polar amplification factors ($f_{\mathrm{PA}}=\mathrm{\Delta }{T}_{\mathrm{NH}}/\mathrm{\Delta }{T}_{\mathrm{g}}$) as follows: at the LGM, f_PA=2.7 is taken from the average of PMIP3 model data (Braconnot et al., 2012), while at the mid-Pliocene Warm Period (mPWP, about 3.2 Myr ago), f_PA=1.6 is calculated from the average of PlioMIP results (Haywood et al., 2013). At all other times, f_PA is linearly varied as a function of the NH temperature. In Appendix A, we investigate the influence of the chosen polar amplification factor on our results. The temporal resolution of the dataset is 2000 years.

Analysing this dataset, Köhler et al. (2018) found a temperature–CO₂ divergence appearing mainly during, or in connection with, periods of decreasing obliquity related to land ice growth or sea level fall. For these periods, a significantly different $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ was obtained than for the remainder of the time frame. However, in the future, we expect sea level to rise; hence, these intervals of strong temperature–CO₂ divergence should not be considered for the interpretation of paleodata in the context of future warming, e.g. by using paleodata to constrain ECS. In the following analysis, we therefore exclude these times with strong temperature–CO₂ divergence, leaving 217 data points as indicated in Fig. 1c, d.

Figure 4 The global temperature perturbations stripped of the inferred influence of land ice ($\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$) calculated using Eq. (16) against $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ from the proxy-inferred paleoclimate dataset (Köhler et al., 2015a), using (a) ε_[LI]=0.79 (maroon dots), ε_[LI]=0.45 (cyan dots), and ε_[LI]=0.25 (green dots). Here, ε_[LI] is obtained by converting the multi-model assemblage of simulated relative influences of land ice changes on the LGM temperature anomaly (0.46±0.14) (Shakun, 2017). (b) Same as in panel (a) but using ε_[LI]=1 (grey dots), which is equivalent to the PALAEOSENS approach. The brown, blue, dark green (a), and black lines (b) represent second-order polynomial least-squares regressions through the data.

4.2 Analysis

Shakun (2017) compiled model-based estimates of the relative impact of land ice changes on the LGM temperature anomaly (ω in Eq. 11) using an ensemble of 12 climate models and estimated ω to be 0.46±0.14 (mean ±1σ; full range 0.20–0.68). Applying these values, in combination with the LGM values (taken here as the mean of the data at 20 and 22 kyr ago) $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=-\mathrm{2.04}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and $\mathrm{\Delta }{R}_{\left[\mathrm{LI}\right]}=-\mathrm{3.88}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, yields ${\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}={\mathrm{0.45}}_{-\mathrm{0.20}}^{+\mathrm{0.34}}$. Implementing this range for ε_[LI] in Eq. (16), we calculate $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$ over the whole 800 kyr period. Fitting second-order polynomials by least-squares regression to the scattered data of $\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$ and $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$, we infer a PI $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ of ${\mathrm{2.45}}_{-\mathrm{0.56}}^{+\mathrm{0.53}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ (Fig. 4a). The substantial uncertainty given here only reflects the 1σ uncertainty in ε_[LI]. Similar to Köhler et al. (2018), we also detect a state dependency with decreasing $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ towards colder climates for this dataset, more strongly so in the case of lower ε_[LI]. This state dependency is opposite to the one found in the CLIMBER-2 results (Sect. 3). The difference may be related either to the fact that fast climate feedbacks are too linear or that some slow feedbacks are underestimated in intermediate-complexity climate models like CLIMBER-2 (see Köhler et al., 2018, for a detailed discussion). At $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}=-\mathrm{2.04}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, the LGM value, $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ is only ${\mathrm{1.45}}_{-\mathrm{0.37}}^{+\mathrm{0.33}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$. The PALAEOSENS approach, which does not consider efficacy and is therefore equivalent to our approach using ε_[LI]=1, yields $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}=\mathrm{1.66}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ for PI and $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}=\mathrm{0.93}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ for the LGM (Fig. 4b). The specific paleoclimate sensitivities we find using the refined approach are hence generally larger than those obtained when neglecting efficacy differences. This is because, for the range of the impact of land ice changes on the LGM temperature anomaly implemented ($\mathit{\omega }=\mathrm{0.46}\pm \mathrm{0.14}$), the efficacy factor ε_[LI] is smaller than unity. In other words, these land ice changes contribute comparatively less per unit radiative forcing to the global temperature anomalies than the CO₂ changes.

Our inferred PI $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ is equivalent to an S^a of ${\mathrm{1.6}}_{-\mathrm{0.4}}^{+\mathrm{0.3}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ when only considering the uncertainty caused by the implemented range in ε_[LI] and to an S^a of ${\mathrm{1.6}}_{-\mathrm{0.2}}^{+\mathrm{0.1}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ when only considering the uncertainty in the conversion factor ϕ. The equivalent ECS is 5.8±1.3 K per CO₂ doubling when only considering the uncertainty caused by the implemented range in ε_[LI] and 5.8±0.6 K per CO₂ doubling when only considering the uncertainty in the conversion factor $\mathrm{\Delta }{R}_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$. The ECS we find is thus at the high end of the results of other approaches to obtain ECS (Knutti et al., 2017), e.g. the 2.0 to 4.3 K 95 % confidence range from a large model ensemble (Goodwin et al., 2018) and the 2.2 to 3.4 K 66 % confidence range from an emerging constraint from global temperature variability and CMIP5 (Cox et al., 2018). Hence, the low end of our ECS estimate is in the best agreement with these other estimates. This could mean that the relative influence of land ice changes on the LGM temperature anomaly is on the high side of or possibly higher than the 0.46±0.14 range we consider here. Alternatively, the conversion factor $\mathit{\varphi }=\mathrm{0.64}\pm \mathrm{0.07}$ we use to convert $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ to S^a is an overestimation, which could be caused by a larger-than-unity efficacy of long-term processes besides CO₂ and land ice changes. We have focused primarily on the effect of ε_[LI] on $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ in this analysis, and therefore we have, for simplicity, ignored uncertainties in the investigated proxy-inferred records themselves. A comprehensive description of these uncertainties and their influence on the calculated climate sensitivity can be found in Köhler et al. (2015a).

5 Conclusions

We have incorporated the concept of a constant efficacy factor (Hansen et al., 2005), which interrelates the global temperature responses to radiative forcing caused by land ice changes and CO₂ changes, into our framework of calculating specific paleoclimate sensitivity $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$. The aim of this effort has been to overcome the problem that land ice and CO₂ changes can lead to significantly different global temperature responses, even when they induce the same global-average radiative forcing. Firstly, we have assessed the usefulness of considering efficacy differences by applying our refined approach to results of 5 Myr CLIMBER-2 simulations (Stap et al., 2018a), where the separate effects of land ice changes and CO₂ changes can be isolated. In the results of these simulations, the error from assuming the efficacy factor to be constant in time is negligible. Thereafter, we have used our approach to reanalyse an 800 kyr proxy-inferred paleoclimate dataset (Köhler et al., 2015a). We have inferred a range in the land ice change efficacy factor ε_[LI] from the relative impact of land ice changes on the LGM temperature anomaly simulated by a 12-member climate model ensemble (Shakun, 2017). The thusly obtained efficacy factor ${\mathit{\epsilon }}_{\left[\mathrm{LI}\right]}={\mathrm{0.45}}_{-\mathrm{0.20}}^{+\mathrm{0.34}}$ is smaller than unity, implying that the impact on global temperature per unit of radiative forcing is less strong for land ice changes than for CO₂ changes. Consequently, our derived PI $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ of ${\mathrm{2.45}}_{-\mathrm{0.56}}^{+\mathrm{0.53}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ is ∼50 % larger than when efficacy differences are neglected. The equivalent S^a and ECS corresponding to this $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ are ${\mathrm{1.6}}_{-\mathrm{0.4}}^{+\mathrm{0.3}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ and 5.8±1.3 K per CO₂ doubling, respectively. The uncertainty in these estimates is only caused by the implemented range in ε_[LI].

Appendix A: Influence of the polar amplification factor

In the analysis performed in Sect. 4.2, we used a global temperature record that was obtained from northern high-latitude temperature anomalies using a polar amplification factor f_PA that varies from 2.7 in the coldest to 1.6 in the warmest conditions (Sect. 4.1). However, recent climate model simulations of the Pliocene using updated paleogeographic boundary conditions show that in warmer times polar amplification could have been nearly the same as in colder times (Kamae et al., 2016; Chandan and Peltier, 2017). We therefore repeat the analysis using the same range in ε_[LI] and the same dataset but with an applied constant f_PA=2.7 over the entire past 800 kyr to generate ΔT_g (ΔT_g2 in Köhler et al., 2015a).

Figure A1 The global temperature perturbations stripped of the inferred influence of land ice ($\mathrm{\Delta }{T}_{[-\mathrm{LI}]}^{\mathit{\epsilon }}$) calculated using Eq. (16) against $\mathrm{\Delta }{R}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ from the proxy-inferred paleoclimate dataset (Köhler et al., 2015a), using (a) ε_[LI]=0.79 (maroon dots), ε_[LI]=0.45 (cyan dots), and ε_[LI]=0.25 (green dots). Here, ε_[LI] is obtained from converting the multi-model assemblage of simulated relative influences of land ice changes on the LGM temperature anomaly (0.46±0.14) (Shakun, 2017). (b) Same as in panel (a) but using ε_[LI]=1 (grey dots), which is equivalent to the PALAEOSENS approach. The brown, blue, dark green (a), and black lines (b) represent second-order polynomial least-squares regressions through the data. Here, the global temperature anomalies are derived from the northern high-latitude temperature anomaly reconstruction assuming a constant polar amplification factor (f_PA) of 2.7, as opposed to the variable f_PA used in Fig. 4.

The constant polar amplification used here counteracts increasing state dependency towards low temperatures, as the temperature differences are no longer amplified by changing polar amplification. Hence, $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ is smaller at PI (${\mathrm{1.96}}_{-\mathrm{0.44}}^{+\mathrm{0.42}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ compared to ${\mathrm{2.45}}_{-\mathrm{0.56}}^{+\mathrm{0.53}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ using the variable f_PA) but diminishes less strongly towards colder conditions (Fig. A1a, cf. Fig. 4a). As before, the PALAEOSENS approach (equivalent to our approach using ε_[LI]=1) yields a lower PI $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}$ of 1.34 K W⁻¹ m² (Fig. A1b). The PI $S_{[{\mathrm{CO}}_{\mathrm{2}},\mathrm{LI}]}^{\mathit{\epsilon }}$ inferred here using our refined approach corresponds to an S^a of ${\mathrm{1.3}}_{-\mathrm{0.3}}^{+\mathrm{0.2}}\mathrm{K}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{m}}^{\mathrm{2}}$ and an ECS of ${\mathrm{4.6}}_{-\mathrm{1.3}}^{+\mathrm{1.0}}\mathrm{K}$ per CO₂ doubling.