Extreme sensitivity and climate tipping points

A climate state close to a tipping point will have a degenerate linear response to perturbations, which can be associated with extreme values of the equilibrium climate sensitivity (ECS). In this paper we contrast linearized (`instantaneous') with fully nonlinear geometric (`two-point') notions of ECS, in both presence and absence of tipping points. For a stochastic energy balance model of the global mean surface temperature with two stable regimes, we confirm that tipping events cause the appearance of extremes in both notions of ECS. Moreover, multiple regimes with different mean sensitivities are visible in the two-point ECS. We confirm some of our findings in a physics-based multi-box model of the climate system.


Introduction
The equilibrium climate sensitivity (ECS) is widely used as a measure for expected future global warming. Following Charney's definition [7], the ECS is the increase in global mean surface temperature (GMST) per radiative forcing change after the fast-acting feedback processes in the Earth System reach equilibrium. Fast-acting means here that those processes are faster than the timehorizon for global mean temperature evolution that interests us, typically taken to be 100 years [30].
The value of ECS remains not very well constrained, as the expected warming per doubling of atmospheric CO 2 still contains a considerable uncertainty of 1.5 -4.5 • C [21]. In fact, this range has not changed much since first ECS estimates based on energy balance arguments, despite enormous developments in climate modelling, and improved observational methods [22]. In particular, large temperature changes and dangerous climate change as a consequence of increased atmospheric CO 2 cannot be excluded. For example, climate observations from the instrumental period have not narrowed down the range of expected climate change mainly because of uncertainties in quantification of the forcing [31]. Recently there have been concerted attempts to constrain ECS using emergent constraints in climate models [9]. Some sources of uncertainty for ECS lie in the classical measurement or model uncertainty, although in particular for observations the quantification of the applied forcing generally contains the largest uncertainties. Furthermore, by going back in time further than the instrumental period (e.g. last millennium, glacial cycles or even millions of years of palaeoclimate data) the uncertainty in both the forcing and the global mean temperature response becomes significant [22]. There has been a debate as to the cause of the uncertainty and extremes of the ECS distribution, in particular the long tail towards high sensitivity values. On the one hand, [29] suggest it is an inevitable consequence of nonlinear transformation of normally distributed feedbacks that appear in the denominator when calculating ECS. On the other hand, [39] suggest they are a sign of 'tipping points' owing to nonlinearities in the system -this has generated a lively debate. In this paper we highlight (a) the notion of ECS can usefully be generalised to a truly nonlinear geometric notion: the two-point sensitivity and (b) the distribution of ECS values is a valuable tool for characterising both state-dependent (feedback) dynamics and tipping of the climate system.
The climate system exhibits both internal and forced variability on many timescales. The consequence is that any 'equilibrium' is only relative to fixing part of the feedback processes that are internal to the climate system, in particular the 'slower' part. This requires an assumption that a time scale separation (into fast and slow processes) exists and the time scale of interest sits between fast and slow. For climate model simulations and observations of the last century there might be a time horizon where this is a reasonable assumption [31], but as we include palaeoclimate data and model simulations into the estimate of ECS this assumption needs to be carefully evaluated. In particular, methods to estimate ECS from palaeoclimate data or models differ from those of (short) climate model simulations; the latter generally derive ECS from the decay of the energy imbalance at the top of the atmosphere induced by a instantaneous doubling or quadrupling of CO 2 [18]; palaeoclimate reconstructions instead make the assumption that the reconstructed climate is in (energetic, short time scale) equilibrium and compare different of these 'equilibria' to each other for estimating ECS. Without compensating for slow feedback processes, palaeoclimate records give the so-called Earth System Sensitivity (ESS) that includes the effect of slow processes and boundary conditions (e.g. geography, vegetation and land ice) [27]. If estimates of these slow processes are available then ECS can be estimated from the ESS under an assumption of time scale separation [30,36].
Note that the ECS is usually thought of as a linearized response of the GMST to perturbations in the radiative balance of the earth. Next to incoming (short-wave) and outgoing (long-wave) radiation, feedback processes in the climate system play an important role in determining the ECS. In their sum these feedback processes tend to enhance ECS (net positive feedback) and associated time scales vary from fast to very slow. Examples of fast feedback processes include cloud feedbacks, water vapour feedback and sea-ice processes. The strength of each of these feedback processes depends on the background (long-term mean) climate state [38] and it is therefore not surprising that the sum of the fast feedback processes varies over time in particular when considering climate states far back in time and under very different boundary conditions. For example, from palaeoclimate records together with ice sheet modelling, it has been found that ECS varies considerably between glacial and interglacial states [23,15,24]. Both present-day and Palaeogene climate model simulations suggest state-dependence of ECS due to feedback processes [6,28]. For example, Transient Climate Response (TCR) considers the deep ocean warming as a slow process.
Abrupt climate shifts have occurred in the past climate system and therefore seem likely to occur in the future for a variety of reasons [26,33]. In the present climate system, potential tipping elements have been identified some of which may have a considerable impact on future values of GMST [26,14]. Even those tipping elements that have little affect on the GMST may cause significant regional damage and/or contribute to global mean climate change by triggering cascades of transitions involving other tipping elements [11]. Across such an abrupt transition there is a breakdown of the assumption of a linear response to perturbations, suggesting that the ECS does not adequately represent the temperature response to radiative perturbations [37,4]. In practice, when deriving ECS from palaeoclimate time series, which include abrupt transitions, these shifts may lead to extreme values of ECS.
In this paper we show that more general notions of ECS can be useful in understanding the response of a climate state to changes in radiative forcing -in addition to an 'instantaneous' linearised notion of ECS we explore 'incremental' and 'two-point' climate sensitivities that are distributions related to dynamic properties of the climate sytem: they characterise the geometry of the dynamics and are not simply estimators of a 'mean ECS'. In fact the distributions of sensitivities reflect the intrinsic uncertainty due to climate system dynamics. The paper is organized as follows: Section 2 introduces these notions of ECS and relates them to the underlying climate dynamics. We illustrate these concepts using a global energy balance model in Section 3. In particular, we relate properties of extremes of the ECS to the presence of tipping points and multistability. Section 4 finishes with a discussion of conclusions and some challenges for the future. Appendix A extends results of [36] and examines extremes of sensitivity associated with tipping points in a more realistic physics-based multi-box model of the glacial cycles by Gildor and Tziperman [17].

Sensitivities and the climate attractor
In order to understand variability, abrupt transitions and response to perturbations we consider the climate system as a high-dimensional multiscale complex dynamical system whose evolving trajectories form a climate attractor. The ECS can be defined on this attractor and regimes or states may be identified where a linear approximation of the response may be reasonable. Tipping points visible in the GMST will show up as large but occasional shifts between different 'climate regimes' of the attractor, or indeed different attractors. We visualise the attractor by projection onto climate variables relevant for determining ECS, i.e. the GMST T and the radiative forcing R per unit area [36]. Consider the energy balance model where the left hand side represents the rate of change of the global mean surface temperature T (with specific heat capacity c T ) and on the right hand side R forcing is the (external) radiative forcing (including changes in CO 2 ), R slow (R fast ) is the radiative perturbation due to all slow (fast) feedback processes within the climate system and R OLW is the outgoing longwave radiation, respectively. Following the formalism of [30], the specific climate sensitivity is which equals the Charney sensitivity S if ∆R slow is the sum of all slow feedback processes contributing to the ECS (and under the assumption of time scale separation). In practise, only some of the slow processes are accessible from palaeoclimate records (e.g. only land ice), in which case the specific climate sensitivity is only an approximation of the Charney sensitivity [30] (e.g. S [CO 2 ,LI] is the specific climate sensitivity considering only land ice changes as slow feedback). This ECS gives a linear prediction for change in temperature: For a specific energy balance model including regime shifts we can explicitly calculate ECS for the different regimes, see section 3. We note that several other authors have highlighted the need to improved notions of ECS: this includes [8] who propose to use a measure-based approach to understand climate sensitivity and [12] who consider conditional climate sensitivities constrained by temperature, coupled with resilience measures for switching to other regimes.

Observation of the climate attractor
We consider the climate system as a high dimensional dynamical system that evolves along trajectories x(t) according to a smooth flow where x ∈ X is in some high dimensional state space and ϕ t (x 0 ) evolves the initial state x 0 along by a time t. The global mean temperature T : X → R and radiative forcing R : X → R are considered to be observables of the underlying dynamical system on X. We assume that the dynamics of x are stationary, i.e. that there is a natural probability measure M on X such that typical trajectories x(t) of (4) satisfy for typical x 0 and any integrable observable p : X → R, i.e. the long-time average of p can be computed using an ergodic hypothesis, by averaging over the measure M in phase space. This implies that, for any open set A ⊂ X, the long-term average proportion of time a typical trajectory spends in A is M (A). In [36] it is supposed there is a stationary measure µ of points in the (∆R [CO 2 ,LI] , T )-plane according to how often they are visited over asymptotically long times, i.e. for any measurable subset A ⊂ R 2 we define for typical initial condition, where χ A is the indicator function, χ A (∆R, T ) = 1 if (∆R, T ) ∈ A and = 0 otherwise. Note that applying (5) with In other words, the measure µ is simply a projection of a natural measure M on the 'climate attractor' onto the two observables (R [CO 2 ,LI] , T ). In general we will consider throughout LI] i.e. the change in radiative forcing relative to some fixed reference levelR [CO 2 ,LI] usually chosen as the level during the pre-industrial climate.

Incremental and two-point sensitivities
In order to predict the temperature at some fixed future time-horizon ∆t in response to a change in radiative forcing, we consider the quotient (2) for fixed changes in time. In this case we can view the distribution of what we call incremental sensitivities as the spread of trajectories from the currently estimated values of (∆R, T ), where from now on we write the radiative forcing corrected for slow feedback as R. On the other hand we can consider a time-independent choice of pairs of points on the climate attractor to obtain two-point sensitivities.
Let us assume the current state at time t = 0 of the climate system is given by a measure σ 0 on phase space X (denoted by point A in Fig. 1). Note that this will always be a measure rather than a point because of lack of knowledge of sub-grid parametrized processes (e.g. [34]) but it will project onto the current values (∆R 0 , T 0 ) = (∆R(x), T (x)) for all x in the support of σ 0 . As time progresses, this state will spread to give a measure at time t that is for any A ⊂ X ( Fig. 1 shows a trajectory in black and others from the ensemble starting at A in grey). The incremental sensitivity for a time interval ∆t is then . Over long time, if there is decay of correlations and mixing of trajectories on the climate attractor [35,34] then σ ∆t → M in the weak sense as ∆t → ∞, and so we expect the distribution of long-term incremental sensitivities for ∆t → ∞ become time-independent for typical trajectories within the attractor: where Note that (7) means that the distribution of long-term sensitivities starting at (∆R 0 , T 0 ) can be written in terms of the geometry of the projected measure µ where we define the two-point sensitivity as The distribution of long-term incremental sensitivities (10) for a generic choice of the initial climate state suggests [36] a time-independent notion of climate sensitivity that can be found by picking pairs of points (R 0,1 , T 0,1 ) independently distributed according to µ and evaluating (2). This means that for any A ⊂ R we can use µ to assign a probability to the sensitivity being in A: with S ∞ 0,1 defined as in (11). This gives, in some sense, a maximal set of possibilities for the sensitivities in that it compares the observables T and ∆R over all possible time points and possible trajectories of the system. This is comparable to the conditional climate sensitivity of [12] except rather than dividing into regimes, they restrict to deviations of temperature at most δ T from T 0 . In the case that the sensitivity is fixed at S 0 , note that S ∞ 0,1 is a Dirac δ-distribution centred at S 0 . Figure 1: Schematic diagram to demonstrate the instantaneous, incremental and two-point sensitivities for an ensemble of trajectories starting at point A that evolves on the climate attractor. The equilibrium of a mean energy balance model is shown as a red curve. For a specific trajectory (shown in black) in the ensemble the slopes of AC and AD correspond to incremental sensitivities (for fixed ∆t) or two-point sensitivities (for varying ∆t). The instantaneous sensitivity is the slope of the tangent to the closest solution in the equilibrium model at B. Note that the cold regime T < T thr and the warm regime T > T thr have different asymptotic (instantaneous) sensitivities corresponding to slopes EB and DF respectively.

Sensitivities and climate regimes
By partitioning the climate attractor into a number of regimes, we can condition the sensitivities on staying within a regime, or undergoing a transition between regimes. By making an optimal partition of the attractor projected into (∆R, T ) space we can hope to find localised distributions of sensitivities for pairs in the same regime. As in [36] we consider these sensitivities conditional on climate regime by partitioning µ into two distributions corresponding to being in a cold (C) or warm (W) state. In our case we set for any measurable A ⊂ R 2 and some threshold temperature T thr . As in [36] we define distributions of conditional sensitivities by Conditional sensitivities for changes of regime are for example The distribution of sensitivities (12) is then a sum of the four conditional sensitivities Moreover, (14) means that we have a symmetry The distributions of S CW and S W C correspond to choices of pairs across the two regimes: these distributions are associated with 'tipping between regimes'. Even though the two-point sensitivities may measure states very far apart in time, we find that extreme values of the sensitivity are usually associated with choice of points from two different regimes.

Sensitivity and tipping in climate models
To illustrate the notions of instantaneous and two-point sensitivities, we consider a conceptual energy balance model: a more complex model is briefly discussed in Appendix A. We consider a variant of the Budyko-Ghil-Sellers energy balance model [5,16,32] for GMST. This model builds on [12,39] and has multiple regimes with state-dependent sensitivity in each. It is a special case of (1) for global mean surface temperature T (t) with atmospheric CO 2 concentration C(t) as a parameter: For this equation, Q 0 represents the solar input modulated by the temperature-dependent albedo α(T ). The change in radiative forcing due to atmospheric greenhouse gases is where A = 5.35 Wm −2 is the direct forcing effect of CO 2 and C 0 represents pre-industrial CO 2 levels. Finally, the outgoing long wave radiation σT 4 is modified by a temperature-dependent emissivity 0 < (T ) < 1.
We consider a temperature-dependent emissivity decreasing from one plateau to a lower one because of changes in water vapour and cloud feedbacks. There are other choices [39], but for simplicity we assume here where 1 and 2 are the limit emissivities for low and high temperatures respectively, T 0 is the threshold and E > 0 corresponds to the (wide) range of temperatures over which there is variation (see Figure 2b). Note that water vapour feedback is sometimes included in the CO 2 term, resulting in an additional constant and modified A [19,12]. Here we separate radiative forcing due to CO 2 and temperature-dependent water feedbacks in emissivity. As in [12], we assume that the albedo varies with temperature due to changes in land-ice feedback processes: we assume there are threshold temperatures T 1 < T 2 associated with changes of albedo α(T ) and define a function that switches from 0 for T < T 1 to 1 for T > T 2 : H(T ) is approximately a Heaviside unit step function and we use a smooth approximation H(T ) = (1 + tanh(T / H ))/2 as in [12]. As in that paper, we write the albedo so that it changes smoothly from a higher albedo α 1 in the presence of more ice surface (T < T 1 ) to a lower α 2 in the presence of more ocean surface (T > T 2 ), see Figure 2a. Note that [12] consider a global transition from ice-covered to ocean-covered earth -here we model a large but regional change in ice cover and therefore have a smaller contrast in global albedo between the two states; our choice of parameters might be more realistic for albedo variations between glacial and interglacial states.
We add a stochastic term to (15) that represents unresolved subgrid processes with amplitude η T : The parameters listed in Table 1 are used, except where specified. Note that the deterministic equilibria of (15) are at F (T, C) = 0, which gives From (18), this means we have equilibria at Figure 2 illustrates temperature dependence of albedo and emissivity as well as the resulting equilibrium forcing ∆R = A ln(Γ(T )/C 0 ) needed to give this temperature. Note there is a unique equilibrium for each T , but not necessarily for each C: as discussed in [12,39] there are three branches of equilibria for a range of C: for the parameters used there is bistability in the region denoted using the red lines in Figs. 3 and 5. We can define the instantaneous sensitivity 1 as is the total feedback factor in this model: S corresponds to the slope of the tangent of the equilibrium (non-stochastic) model (see Fig. 1, point B). The sensitivities on the stable branches differ due to both nonlinearity of black body radiation and change in emissivity. Due to the choice of albedo and emissivity changes in this model, α is nonzero only in the bistable regime and is nonzero only in the temperature range where it varies  Table 1.
(a) temperature-dependence of albedo α(T ) and (b) emissivity (T ); (c) CO 2 and (d) radiative forcing levels necessary to give temperature equilibria, corresponding to (18) and (19), respectively. Note the region of multistability, and temperature-dependence of the sensitivity corresponding to slopes in the bottom right figure. Figure 2. State-dependence between glacial and interglacial states has been detected in estimates of specific climate sensitivities from different palaeo-data, suggesting lower sensitivity during cold periods than during warm periods (e.g. [38,15] who estimate a close approximation of the Charney sensitivity and find warm (interglacial) climate states to be about 60% more sensitive than cold (glacial) states). We can compute the curvature of (19) as Note that this is small except near the folds at T ≈ T 1 and T ≈ T 2 with maximum absolute value d 2 dT 2 ∆R of order −1 H . This confirms that the saddle-node becomes non-smooth in the limit H → 0. The second derivative gives the size of the quadratic correction a in Zaliapin et al. [39]; very large values of the slope near the saddle nodes correspond to the run-away climate observed in [4].
For the model (17), the atmospheric CO 2 concentration is a parameter for the energy balance dynamics. We explore this by considering a 'wandering' CO 2 profile such that γ(t) := ln(C(t)) undertakes a Brownian motion with growth in variance η γ per unit time between reflecting limit values. More precisely, we consider CO 2 dynamics governed by soft reflecting boundary conditions at ln C min and ln C max : where and we use parameters K = 10 −7 s −1 , C min = 10 2 ppm, C max = 10 3 ppm, η C = 2 × 10 −6 s −1/2 .
Clearly there are common causes of variability of temperature and CO 2 and so in general there will be strong correlations between the noise terms W T and W γ ; for convenience we assume here that they are uncorrelated. In most studies of climate sensitivity, carbon cycle processes are not treated as feedbacks but rather as forcing (external to the system); this assumption corresponds in our model to CO 2 driving temperature changes with no direct impact of temperature on CO 2 .
The parameter η C determines the timescale of wandering of the CO 2 : we consider cases where this is slower than, or comparable to the the timescale of evolution of T . Figure 3(b) shows a time series for a typical simulation of (17) with wandering CO 2 (22) and parameters as in Table 1, while the corresponding time series of T is shown in Figure 3(a); we see as C crosses thresholds (for the bistable regime as calculated in (20)) the state of the system tips between warm and cold states. Global mean temperature T vs C and the relative radiative forcing ∆R [CO 2 ] are shown in Figure 2 (c,d) for the non-stochastic model and in Figure 4 for the stochastic model, respectively. Observe the region of bistability around ∆R [CO 2 ] = 0 (Figure 4b) that switches rapidly between cool high albedo and warm low albedo states via saddle-node bifurcations. There is approximate linearity away from these tipping points, but with different mean slopes.

Extreme sensitivities, early warning signals and tipping between regimes
The energy balance model (17) with wandering CO 2 (22) gives a framework in which one can test correlation between extreme values of sensitivity and tipping between climate regimes, as  well as testing other possible early warning signals for a tipping event. Indeed, we find a rise in instantaneous sensitivity seems to act as a good precursor in cases of slow variation of CO 2 . Figure 5(left) shows the variation of instantaneous sensitivity and two early warning signals for tipping between regimes for the wandering variation of CO 2 concentration ln C(t) relative to T  Figure 3. The red line corresponds to the steady solution of (15) with dependence on CO 2 . timescale, with η C = 2 × 10 −6 s −1/2 ; this means that CO 2 variations are comparable in timescale to T fluctuations. By considering the nearest equilibrium point on C = Γ(T ) for fixed C (18), we evaluate the instantaneous sensitivity S using (21) and plot this in middle panel. Observe there is a qualitative change in S before and after tipping events. There is a clear precursor and then singularity as the tipping point is crossed. Note that in this case the instantaneous sensitivity depends only on the current CO 2 level and the nearest branch -fluctuations in T around the branch do not affect S, while fluctuations in C do. There are also apparent 'false alarm': for example, the fluctuations of S around 3 kyr and 18 kyr. Figure 5(left, panels d,e) show detrended estimates of sd(t) and AR1(t) [40] using moving averages with length τ = 500yrs. Both do not show any precursors before tipping events. By contrast, Figure 5(right) shows early warning signals for tipping between regimes for slower variation of CO 2 concentration with η C = 5 × 10 −7 s −1/2 , where T evolves faster than C. Unlike the case on the left, this slower switching gives a precursor of increasing AR1. In both cases there are increasing fluctuations of instantaneous sensitivity S. Note however, that the instantaneous S we consider here is based on our model equations, and will be more difficult to access from complex model realisations or observations.

Two-point sensitivities and tipping
An approximation of the stationary density of the global attractor for the system (17,22) is shown in Figure 6. We classify the system regime as one of: where we choose a threshold T thr = 10 C between the two stable branches: see Figure 1. We simulate a single very long trajectory (5 × 10 5 years) of the energy balance model (15,17) with wandering CO 2 and use this to create a density plot of the climate attractor projected onto T vs  (20). (c) shows estimated instantaneous sensitivity from the nearest equilibrium of (17). Note the gradual rise and fluctuations in S on approach to the tipping point, and the two levels of S corresponding to the differing sensitivities of the two stable branches. (d,e) show standard deviation and AR1 coefficient: note that the AR1 coefficient seems to have predictive power only for the right column.
∆R, as shown in Figure 6. This is used to consider the two-point sensitivities and probabilities of tipping as in Figure 7. The left column is computed by sampling incremental sensitivities (8) from points for increments up to 20 kyr. The right column is computed by sampling 10 7 pairs of points from the distribution in Figure 6 and using the two-point sensitivity (11). We observe: • There is good qualitative agreement between the incremental sensitivities averaged over long delays and the two-point sensitivities sampled independently from the attractor. Indeed, the autocorrelation of the timeseries for T (not shown) has substantially decayed and has its first zeros around 20kyr.
• High probability of tipping (see Figure 7(b,d)) corresponds mostly to extremes of S that may be positive or negative S.
• Within the W and C regimes, the sensitivities are closely clustered but have different means for the W and C state. We can estimate these using average temperatures and (21) as S ≈ 0.79 K[W m 2 ] −1 ) for the W and S ≈ 0.55 for the C state, respectively.
• Note that there are relatively low-probability 'shoulders' of the distributions within-regime. These are due to the classification of regimes also including states that are in transitions: although the system is in transition, both points are still classified as the same regime.   Figure 7: (a,c) Conditional two-point sensitivities and (b,d) probabilities of tipping from warm (a,b) and cold (c,d) states, for the energy balance model (17) with wandering CO 2 (22). The left column is computed using a range of delays up to 20 kyr while the right column is the distribution of two-point sensitivity from 10 7 independently sampled pairs of points from the distribution in Figure 6. Note that sensitivities outside the horizontal range are rounded into the last bin.

Sensitivities in the absence of tipping
Considering the same model (15,17) with wandering CO 2 (22), we use different parameters to contrast the results in the previous sections with cases where there is no bistability. In particular we consider parameters as in Table 1 except for: • Default albedo contrast: α 1 = 0.52, α 2 = 0.47 (i.e. also as in Table 1).
• No albedo contrast: α 1 = α 2 = 0.495. Figure 8 shows the low and no albedo contrast cases, comparing to the default case Figure 5(left). For the low albedo contrast case there is no longer a region of bistability, but there is nontrivial variation of the instantaneous sensitivity along the attractor. For the no albedo contrast case, the instantaneous sensitivity is close to constant. The projection of the climate attractor into the (∆R, T ) plane is shown in Figure 9(a-c), while the corresponding distribution of two-point sensitivities in Figure 9(d-f). Observe the presence of non-unimodal distributions for (d,e) associated with regions with different two-point sensitivity, and clear skewness and tails again associated with the geometry of the measure in (a,b). Note the higher average in (e) corresponds to there being only a single regime in (b) that runs over a wide range of temperatures.
In physical terms, the skewness (and long tails) in (d,e) originate from the state-dependence and nonlinearity of feedbacks (i.e. non-constant feedback factors). The bistability of the two regimes with different feedbacks gives the two peaks in the distribution of Figure 9(d). However, Figure 9(e) still has two peaks: these originate from state dependence on the same attractor (Figure 9(b). For this low-albedo-contrast case, there is no 'tipping' but we still find very non-Gaussian distribution of S that comes from nonlinearities in the system that, in this case, do not produce tipping. Note that only in the no albedo contrast case (c) is there a plausible fit to Gaussian.

Discussion
We demonstrate that state-dependence and the presence of tipping points produces signatures in the distribution of instantaneous and two-point notions of ECS. We explore this using a global energy balance model where state-dependence and multistability originate from the dependence of both albedo and emissivity on temperature.
For the deterministic version of our model (15) with fixed CO 2 the changes in albedo mean there can be bistability between regimes, while the changes in emissivity contribute to different sensitivities within these regimes. The distribution of ECS comes from several sources -nonlinearities that result in tipping points and/or state-dependence of the feedbacks and sub-grid variability that we model here as stochastic perturbations. Such regime-dependent sensitivity and extremes associated with tipping points are also visible in the more complex Gildor and Tziperman model [17] also investigated in [36], as outlined in Appendix A.
For the stochastic model (17) with wandering CO 2 , regime-dependent sensitivity is visible as differences in slope of the stable regimes for the T vs ∆R plots (see Figure 4(b)). The densities of the stable regimes for the T vs ∆R plots (see Figure 6) show varying slopes and so conditional two-point sensitivities for the two regimes (see Figure 7(a) and (c)) can have peaks at different sensitivities. We compare several notions of sensitivity. These are the instantaneous sensitivity  associated with the slope of the equilibrium branch, incremental sensitivity associated with a fixed delay, and two-point sensitivity that compares arbitrary points on the climate attractor. The presence of tipping points gives extremes of sensitivities in that (i) there are large fluctuations of instantaneous sensitivities for the nearest equilibrium just before a tipping point (see Figure 5); and (ii) in the distributions of conditional two-point sensitivities that cross regimes (see Figures 7(d) and (f)). It is remarkable that the two-point sensitivities are so informative, given that they compare points that may be very far from each other in time.
There remains much work to be done to understand the relation between and limitations of these notions of ECS, and indeed ECS calculated in other ways, for example from instantaneous doubling of CO 2 : this will involve transient non-equilibrium processes due to ocean thermal inertia. Note that determining ECS from palaeoclimate time-series [37], the two-point notion clearly has an advantage that we are not limited by the time-resolution of the time series.
The energy balance model can be criticised as being very simple and hard to parametrize in terms of the various physical processes that contribute to albedo and emissivity. Moreover, we consider CO 2 in (17) purely as a forcing term which ignores known land surface and ocean processes where temperature is known to affect CO 2 balance. However, the model is complex enough to confirm that extremes in ensembles of computed climate sensitivities can indicate nearby tipping points. Computations presented in Appendix A confirm this picture in a more complex box-model for the glacial cycles, where the CO 2 is modelled dynamically.

Future perspectives
When (and how) extremes of sensitivity can be effective precursors of a tipping event will depend on a number of factors. In particular, the timescale of dynamics of the climate response needs to be faster than the timescale of changes in forcing. Figure 5(left) shows that as CO 2 variability is rapid, this results in tipping points with little precursor visible in changes to AR1, though it is visible in the instantaneous sensitivity. There may be 'rate-induced' tipping points [1,2] that appear when the timescale of the forcing interacts with that of the system. The size of the region of effective nonlinearity can also vary -for (17) this is affected by H which 'smooths' the fold bifurcations. Note also that although tipping points do give rise to extremes in the distribution of ECS, extremes do not necessarily indicate a tipping point.
Translating these results to more complex models will be difficult: the 'climate attractor' is harder to define in the presence of large recurrence times and a variety of nonlinear multiscale processes. In such cases, interpreting transitions as tipping points is a challenge; nonetheless, the palaeoclimate record does show a variety of large and sudden transitions [25]. For example, ice core/ocean core records indicate repeated sudden changes in (regional) surface temperature associated with glacial cycles [19] or Dansgaard-Oeschger [13] events as well as global transitions, for example the greenhouse-icehouse transition at the Eocene-Oligocene transition [10]. Although glacial cycles can be found in models such as [17] as relaxation oscillation with clear regimes, for climate reconstruction data these regimes are not so clear (e.g. [23,15]).
Recent work [33] suggests we are at a crossroads in terms of the future earth system state. On the one hand, looking at the palaeoclimate record for the last 1 million years suggests that we are overdue descent into an ice age. On the other hand, comparison of anthropogenic CO 2 emissions with the palaeo record suggest the next tipping point may be to much warmer 'hothouse' earth. A better understanding of improved indicators such as two-point ECS and what they say about the climate response to changes in greenhouse gases, together with a better understanding of hothouse earth climate states that may have existed in the past (e.g. the Palaeocene climate [20,3]) should help our understanding and guide future generations in their need to avoid dangerous climate change.

A Appendix: Tipping in a multi-box climate model with ocean biogeochemistry
To investigate the notion of two-point sensitivity in a more complex and physics based earth system model, we explore here the Gildor-Tziperman (2001) model [17] with Milaknovich forcing and biogeochemistry in the ocean. In this model the ocean and atmosphere consists of 4 latitudinal (zonally averaged) boxes each, with two layers in the ocean and one in the atmosphere. Within the boxes a variety of thermodynamic quantities (eg temperature T ) and species (eg ocean salinity S and atmospheric CO 2 ) are modelled. In addition there is dynamic land ice (slowly evolving) and sea ice (rapidly evolving) as well as fluxes that join these boxes. The system is sufficiently complex to allow modelling of the Pleistocene ice-age oscillations of land-ice in response to Milankovich forcing. The glacial-interglacial cycles appear in this model as internally generated self-sustained oscillations, which are then modified by the Milankovich forcing. More details are given in [17,36]. Figure 10 shows projections of a long trajectory (500 kyr) on the climate attractor of this more complex climate model onto the plane of global mean temperature T against (a) ∆R [CO 2 ] (b) ∆R [LI] and (c) ∆R [CO 2 ,LI] [36]. Observe that the all three projections of clearly show two climate regimes, a lower 'cold' state (corresponding to large amounts of sea ice and T < 12.28C) and an upper 'warm' state (corresponding to T > 12.28C). The projection on the combined radiative forcing of CO 2 and the slow feedback in land ice changes ∆R [LI] shows a clear slope and hence 'mean' sensitivities in both regimes (Figure 10c). Following the formalism of [30], the slopes in (a) should reflect the Earth System sensitivity, while (c) should give a good approximation of the Charney or equilibrium climate sensitivity. Note that in this model there is only one slow feedback to correct for, namely the land-ice albedo feedback. When projecting onto the ∆R [CO 2 ] plane ( Figure 10a) the cold regime appears very diffuse and with very high (or sometimes negative) Earth System sensitivities, suggesting that in the cold branch local climate dynamics are not entirely determined by CO 2 . Similarly, when projecting only on land ice changes ∆R LI (Figure 10b) the warm branch appears more diffuse, i.e. CO 2 dynamics seem to be more important than land ice dynamics on this branch.
The left column of Figure 11 (a-d) shows the distribution of two-point climate sensitivities for R [CO 2 ,LI] conditional on regime: this is comparable to Figure 7(left) in that one can observe (i) clearly localised distributions of ECS in (a,c) conditional on remaining within the W or C regime, (ii) a broader distribution in (c): this seems to be associated with the curvature of the C regime branch in Figure 10c), (iii) a clear association of tipping from W to C (b) or from C to W (d) being associated with extreme sensitivities. Note that for this model there is no energy balance model available and so it is not possible to compute the instantaneous sensitivity. For comparison we show in the right column of Figure 11 (e-h) the same distributions for the Earth System Sensitivities (ESS), which are not compensated for the slow feedback (in this model the land-ice albedo feedback). Observe that for both regimes there is a much broader distribution for the ESS (e and g) than for the ECS in (a) and (c). In particular in the cold regime, earth system sensitivities (g) are much higher than equilibrium sensitivities (c) because the land-ice albedo feedback is very strong in the cold (land-ice covered) states.  Figure 10: Plots of T vs ∆R for the Gildor-Tziperman model [17]; see [36]. (a) projection onto the (∆R [CO 2 ] , T )-plane, which should reflect the Earth System Sensitivity only considering CO 2 as forcing; (b) projection onto the (∆R [LI] , T )-plane, considering only the slow land-ice albedo feedback as forcing; (c) projection onto the (∆R [CO 2 ,LI] , T )-plane, considering both CO 2 and the slow feedbacks (in this model there is only the land-ice albedo feedback slow) as forcing. Following the formalism of [30] the slopes in this graph should reflect a good approximation of the Charney ECS.  Figure 11: Conditional sensitivities and probabilities of tipping for glacial cycles simulations with the Gildor-Tziperman model [17] (as in Figure 7 for the energy balance model). (a-d) Conditional sensitivities compensating for slow feedbacks (i.e. reflecting ECS) and probabilities of tipping (i.e. from the distribution shown in Figure 10c). (e-h) Conditional earth system sensitivities ESS (not compensated for slow feedbacks) and probabilities of tipping (i.e. from the distribution shown in Figure 10a).