On determining the point of no return in climate change

Earth’s global mean surface temperature has increased by about 1.0 C over the period 1880–2015. One of the main causes is thought to be the increase in atmospheric greenhouse gases. If greenhouse gas emissions are not substantially decreased, several studies indicate that there will be a dangerous anthropogenic interference with climate by the end of this century. However, there is no good quantitative measure to determine when it is “too late” to start reducing greenhouse gas emissions in order to avoid such dangerous interference. In this study, we develop a method for determining a so-called “point of no return” for several greenhouse gas emission scenarios. The method is based on a combination of aspects of stochastic viability theory and linear response theory; the latter is used to estimate the probability density function of the global mean surface temperature. The innovative element in this approach is the applicability to high-dimensional climate models as demonstrated by the results obtained with the PlaSim model.


Introduction
In the year 2100, which is as far away (or as close) as 1932 in the past, mankind will be living on an Earth with a different climate than today.At that time, we will know the 2100-mean Global Mean Surface Temperature (GMST) value and its increase, say ∆T , above the pre-industrial GMST value.
From the then available GMST records, it will also be known whether this change in GMST has been gradual or whether it was rather 'bumpy'.If the observational effort will continue as of today, there will also be an adequate observational record to determine whether the probability of extreme events (e.g., flooding, heat waves) has increased.
The outcomes of these future observations, to be made by future generations, will strongly depend on socio-economic and technological developments and political decisions which are made now and over the next decades.Fortunately, there is a set of tools available to inform decision makers: Earth System Models.These models come in different flavours, from global climate models (GCMs) providing details on the development of the ocean-atmosphere-ice-land system to integrated assessment models (AIMs) which also aim to describe the development of the broader socio-economic system.
During the preparation for the fifth assessment report (AR5) of the Intergovernmental Panel of Climate Change (IPCC), GCM studies have focussed on the climate system response to GHG changes as derived by AIMs from different socio-economic scenarios; the data from these simulations is gathered in the so-called CMIP5 archive (http://cmip-pcmdi.llnl.gov/cmip5/).
Depending on the representation of fast climate feedbacks in GCMs, determining their climate sensitivity, the CMIP5 models project a GMST increase ∆T of 2.5-4.5 • C over the period 2000-2100 (Pachauri et al., 2014).This does not mean that the actual measured value of ∆T in 2100 will be in this interval.For example, the GMST may be well outside this range because of current model errors which misrepresented the strength of a specific feedback.As a consequence, a transition might have occurred in the real climate system, which did not occur in any of the CMIP5 model simulations.Another possibility is that the GHG development eventually was far outside of the scenarios considered in CMIP5.
A crucial issue in 2100 will be whether a climate state has been reached where a dangerous anthropogenic interference (DAI) can be identified (Mann, 2009).In this case, present-day islands will have been swallowed by the ocean, extreme events have increased in frequency and magnitude (Smith and Schneider, 2009).These effects are then very inhomogeneously distributed over the Earth and have lead to enormous socio-economic consequences.If this is the case in 2100, then there is a point in time where we must have crossed the conditions for DAI.This time, marking the boundary of a 'safe' and 'unsafe' climate state, obviously depends on the metrics used to quantify the state of the complex climate system.
In very simplified views, this boundary is interpreted as a threshold on CO 2 concentration (Hansen et al., 2008) or on GMST.The latter, in particular the ∆T c = 2 • C threshold, has become an easy to communicate (and maybe therefore leading) idea to set mitigation targets for greenhouse gas reduction.Emission scenarios have been calculated (Rogelj et al., 2011) such that ∆T will remain below ∆T c .Although thresholds on GMST have been criticized for being very inadequate regarding impacts (Victor and Kennel, 2014), such a threshold forms the basis of policy making as is set forward in the Paris 2015 (COP21) agreement.
Suppose that measures are being taken to keep ∆T < ∆T c , does this mean that we are 'safe'?
The answer is a simple no, as regionally still DAI may have occurred such as the disappearance of island chains due to sea level rise (Victor and Kennel, 2014).Hence, attempts have been made to define what 'safe' means in a more general way, such the Tolerable Windows Approach (TWA) (Petschel-Held et al., 1999) and Viability Theory (VT) (Aubin, 2009).These approaches also deal with general control strategies to steer a system towards 'safety' when needed.On a more abstract level, both TWA and VT start by defining a desirable (or 'safe') subspace V of a state vector x in a general state space X.This subspace is characterized by constraints, such as thresholds on properties of x.For example, when x is a high-dimensional state vector of a GCM, such a threshold could be ∆T < ∆T c on GMST.When the time-development (or trajectory) of x is such that it moves outside the subspace V , a control is sought to steer the trajectory back into V .Note that this is an abstract formulation of the mitigation problem, when the amplitude of the emission of greenhouse gases is taken as control.Recently, Heitzig et al. (2016) have added more detail to regions in the space X which differ in their 'safety' properties and amount of flexibility in control to steer to 'safety'.
Given a certain desirable subspace of the climate system's state vector (e.g., to avoid DIA) and a suite of control options, (e.g., CO 2 emission reduction) it is important to know when it is too late to be able to steer the system to 'safe' conditions, say at the year 2100.In other words, when is the Point of No Return (PNR)?The TWA and VT approaches, and the theory in Heitzig et al. (2016), suffer from the 'curse of dimensionality' and cannot be used within CMIP5 climate models.For example, the optimization problems in VT and TWA lead to dynamic programming schemes which have up to now only been solved for model systems with low-dimensional state vectors.The approach in Heitzig et al. (2016) requires the computation of region boundaries in state space, which also becomes tedious in more than two dimensions.Hence, with these approaches it will be impossible to determine a PNR using reasonably detailed models of the climate system.
In this paper, we present an approach similar to TWA and VT, but one which can be applied to high-dimensional models of the climate system.Key in the approach is the estimation of the probability density function of the properties of the state vector x which determine the 'safe' subspace V .
The PNR problem is coupled to limitations in the control options (e.g. of emissions) and can be defined precisely using these options and stochastic viability theory.The methodology is presented in section 2 and just to illustrate the concepts, we apply the approach in section 3 to an idealized energy balance model with and without tipping behavior.In section 4, the application to a high-dimensional climate model follows, using data from the Planet Simulator (PlaSim, Fraedrich et al. (2005)).A summary and discussion in section 5 concludes the paper.

Methodology
Here we briefly describe the concepts we need from stochastic viability theory and then define the PNR problem, specifically in the climate change context.

Viable states
Viability theory studies the control of the evolution of dynamical systems to stay within certain constraints on the system's state vector (Aubin, 2009).Here we consider finite dimensional deterministic systems, with state vector x ∈ R d and vector field f : R d → R d , given by In the general formulation of viability theory a time-dependent input is also considered in the right hand side of Eq. ( 1) which can be used to control the path of the trajectory x(t) in state space.
For our purposes, we only need the concept of a viable state, which is related to constraints on the state vector defining a viable region V in state space, also called the viability constraint set.In the model (1) such a set can, for example, be defined by a threshold condition ||x|| < ||x * || An initial condition x 0 = x(t = 0) is called viable if x(t) ∈ V , for all 0 ≤ t ≤ t * , where t * is a certain end time.The set of all these initial conditions forms the viability kernel associated with V .
Stochastic extensions of viability theory consider finite dynamical systems defined by stochastic differential equations where X t ∈ R d is a multidimensional stochastic process, W t ∈ R n is a vector of n-independent standard Wiener processes and the matrix g ∈ R d×n describes the dependence of the noise on the state vector.The normalised probability density function (PDF) p(x, t) can be formally determined from the Fokker-Planck equation associated with Eq. (2).
A stochastic viability kernel V β consists of initial conditions X 0 for which the system has, for 0 ≤ t ≤ t * , a probability larger than a value β to stay in the viable region V (Doyen and De Lara, 2010).For example, in a one-dimensional version of Eq. (2) with state vector X t ∈ R and with a viable region V given by x ≤ x * a state X t is called viable, with tolerance probability β T , if and otherwise, X t is called non viable.

Linear response theory
In relatively idealized low-dimensional models (such as the energy balance model in section 3), the probability density functions can be easily computed by solving for the Fokker-Planck equation (see section 3.2).However, in order to find the temporal evolution of the PDF of the global mean surface temperature GMST under any CO 2 eq forcing in high-dimensional climate models, such as PlaSim in section 4, we will use linear response theory (LRT).With this theory, the effect of any small forcing perturbation on the system state can be calculated by running the climate model for only one forcing scenario (Ragone et al., 2014).
In LRT, the expectation value of an observable Φ, when forcing the system with a time-dependent function f (t), can be calculated by computing the convolution of a Green's function G Φ and the forcing f (t), according to To construct this Green's function, the property that the convolution in the time domain is the same as point-wise multiplication in the frequency domain is used.The Fourier transform of Eq. ( 4) is given by with χ Φ (ω), Φ f (ω) and f (ω) being the Fourier transforms of G Φ (t), Φ f (t) and f (t), respectively.Therefore, once the time evolution of the expectation value of an observable under a certain forcing is known, the Green's function of this observable can be constructed with Eq. ( 5) and consequently the linear response of the observable to any forcing can be calculated.

The Point of No Return problem
In the climate change context, scenarios of GHG increase and the associated radiative forcing have been formulated as Representative Concentration Pathways (RCPs).In Pachauri et al. (2014), there are four RCP scenarios (Fig. 1a) ranging from an increase in radiative forcing of 2.6 W m −2 (RCP2.6)at 2100 (with respect to 2000) to a forcing increase of 8.5 W m −2 (RCP8.5).
To define the PNR for each of these RCPs, a collection of mitigation scenarios on greenhouse gas emission has to be considered.These mitigation scenarios will lead to changes in GHG concentrations, described by functions F λ (t), where λ is a parameter.For instance, the collection F λ could result from mitigation measures that lead to an exponential decay to different stabilisation levels (measured in CO 2 equivalent, or CO 2 eq) within a certain time interval.An example of such a collection F λ is shown by the dashed and dotted red lines in Fig. 1b.The most extreme member of F λ is defined as the mitigation scenario (represented by a certain value of parameter λ) which has the steepest initial decrease at a certain time t (dashed curve in Fig. 1b).
Along the curve of a certain RCP scenario, there will be a point in time where action will be taken to reduce emissions of GHG; this is indicated by a time of action t b .Consider for example (Fig. 1b) that t b is chosen as the first year where the state vector X t is not viable anymore.A reduction in emissions is, however, not immediately followed by a decrease in CO 2 eq due to the long residence time of atmospheric CO 2 .In addition, there is a delay to take action in emission reduction due to technological, social, economic and institutional challenges.Hence, emission reduction will only start ∆t 1 years after t b .The CO 2 eq will, even after emissions have been reduced, also still increase over a time ∆t 2 .The time at which the CO 2 eq starts to reduce according to F λ is indicated by , where ∆t = ∆t 1 + ∆t 2 .Eventually, X t may become viable again and this point in time is indicated by t e (Fig. 1b).
For a given RCP scenario, tolerance probability β T , viable region V and collection F λ , we define the PNR (π t ) as the first year t c where, even when at that moment the most extreme CO 2 eq reduction scenario F λ applies, (a) either X t will be non viable for more than τ T years, where τ T is a set tolerance time, or (b) X t will be non viable in the year 2100.
The first PNR, which we will indicate below by π tol t , is based on limiting the amount of years that X t is non viable, since (during these years) society is exposed to risks from, for example, extreme weather events.The second PNR, which we will indicate below by π 2100 t imposes no restrictions on how long X t is non viable, but it is only based on that X t is non viable at the end on this century.
Hence, under the given set of mitigation options, it is guaranteed that the state will have left the viable region by the year 2100.We will use both PNR concepts in the results below.

Energy balance model
In this section, we illustrate the concepts and the computation of the PNR for an idealized energy balance model of Budyko-Seller type (Budyko, 1969;Sellers, 1969).We will also assume that the CO 2 eq can be directly controlled, and hence no carbon cycle model is needed to determine CO 2 eq from an emission reduction scenario.

Formulation
We use the stochastic extension of the model formulation as in Hogg (2008).The equation for the atmospheric temperature T t (in K) is given by The values and meaning of the parameters in Eq. ( 6) are given in Table 1.The first term in the right hand side of Eq. ( 6) represents the short-wave radiation received by the surface and α(T ) is the albedo function, given by This equation contains the effect of land ice on the albedo and continuous approximation of the Heaviside function.When the temperature T < T 0 , the albedo will be α 0 and when T > T 1 it will be α 1 and the albedo is linear in T for T ∈ [T 0 , T 1 ].The second term in the right hand side of Eq. ( 6) represents the effect of greenhouse gases on the temperature.It consists of a constant part (G), and a part (A ln C(t) C0 ) depending on the mean CO 2 eq concentration in the atmosphere (indicated by C(t)).The third term in the right hand side of Eq. ( 6) expresses the effect of long-wave radiation on the temperature and the last term represents noise with a constant standard deviation σ s .The standard value of σ s chosen as 3% of the value of G/c T , hence about 0.3K/year.The variance in CO 2 concentration originates mostly from seasonal variations, and the 3% is on the high side.Nevertheless, we still use this value, because if we take values smaller than 3% the PDF of the GMST will almost be a delta function and concepts can not be illustrated clearly.

Results: stochastic viability kernels
When using the global mean CO 2 eq concentration C in Eq. ( 6) as a time-independent control parameter, a bifurcation diagram can be easily (numerically) calculated for the deterministic case (σ s = 0).
In Fig. 2, such diagrams are plotted of C versus the equilibrium temperature T for two values of α 1 .To obtain realistic values for the temperature, the equilibrium temperature equilibria are shifted upwards by 30 K.This is done by substituting T with T − 30 and adapting the right hand side of Eq. ( 6) such that the new temperature is a steady state.This is obviously a bit artificial here, but we justify it by our aim to only wanting to illustrate the methodology; results from more realistic models will follow in section 4 below.The diagram corresponding to α 1 = 0.2 (Fig. 2a) has two saddle-node bifurcations which are absent for α 1 = 0.45 (Fig. 2b).From now on, the energy balance model with α 1 = 0.45 and α 1 = 0.2 will be called the monostable and bistable case, respectively.For σ s = 0, we explicitly determine the normalised PDF p(x, t).Rewriting Eq. ( 6) as with ), the Fokker-Planck equation of Eq. ( 8) is given by This differential equation is solved numerically for p(x, t) under any prescribed function C(t) with boundary conditions p(x u , t) = p(x l , t) = 0, where x l = 270 K and x u = 335 K, and an initial condition p(x, 0) (specified below) satisfying We first show stochastic viability kernels for each initial condition T 0 and C 0 , where C 0 is an initial CO 2 eq concentration and T 0 is the expectation value of the initial PDF of T t .As starting time, we take the year 2030 and suppose that the climate system will be forced by a certain RCP scenario from 2030 till 2200.For every C 0 , the original RCP scenario from Fig. 1a is adjusted such that its time development remains the same, but it has C 0 as CO 2 eq concentration in 2030.The PDF of the GMST p(x, t = 0) (t = 0 refers to the year 2030) has a prescribed variance (defined by σ 2 s ) and expectation value T 0 .
In Fig. 3, the stochastic viability kernels are plotted for the energy balance model forced by the RCP4.5 scenario and a viable region V defined by T ≤ 293 K.The results for the monostable and bistable cases are plotted in Fig. 3a and Fig. 3b, respectively.The colors indicate for each combination of T 0 and C 0 in which stochastic viability kernel the initial state (C 0 , T 0 ) is located.For example, consider the bistable case and an initial condition of T 0 = 288 K and C 0 = 400 ppmv, then this initial condition is in the kernel V β with β ≥ 0.9.This means that, with a probability larger than 0.9, a trajectory of the model starting at (C 0 , T 0 ) will remain viable up to the year 2200, where C follows the RCP4.5 scenario.The white areas contain initial conditions that are in a stochastic viability kernel V β with β < 0.5.
The sensitivity of the stochastic viability kernels with respect to RCP scenario, threshold defining the viable region V and amplitude of the noise σ s was also investigated (results not shown).The behaviour is as one can expect in that the area of the kernels becomes smaller (larger) when noise is larger (smaller), when the threshold temperature is smaller (larger) and when the radiative forcing associated with the RCP scenario is more (less) severe.For example for the RCP6.0 scenario, each combination of T 0 and C 0 (same range as in Fig. 3) is in a V β with β < 0.5 for both mono-and bistable cases.

Results: Point of No Return
Again for illustration purposes, we assume that reduction of the emissions will have an immediate effect of the CO 2 eq, such that effectively the CO 2 eq is controlled.We choose the collection F λ to consist of mitigation scenarios that exponentially decay to the preindustrial CO 2 eq concentration, which is 280 ppmv.For this exponential decay, we consider different e-folding times τ d .The most extreme scenario has an exponential decay within 50 years, which corresponds to an e-folding time of τ d = 9 years.Hence, the collection F λ is given by (for τ d ≥ 9) In this equation, t c is the time at which the scenario is applied and C tc the associated CO 2 eq concentration at that moment.
Next, we determine PNR values π tol t for the energy balance model when it is forced by the four different RCP scenarios using a tolerance probability of β T = 0.9 and a tolerance time of τ T = 20  years.The π tol t values for a system forced with the RCP4.5, RCP6.0 and RCP8.5 scenarios are shown in Fig. 4 for both the monostable and bistable cases.As expected, the more extreme the RCP scenario, the earlier the PNR.This can be easily explained by the fact that when the CO 2 eq concentration is rising faster, the temperature will get non viable earlier.Consequently, the PNR will be earlier, since the GMST is only allowed to be non viable for at most τ T years.When the model is forced with RCP2.6, there is no PNR for both models.The reason for this is that the CO 2 eq concentration will remain low throughout the whole period and consequently the temperature will stay viable.The value of π tol t of the bistable case is for each scenario earlier than the value of the monostable case.This can be clarified by the fact that the PDF of the temperature in the bistable case will leave the viable region at a lower CO 2 eq concentration because of the existence of nearby equilibria.
The sensitivity of π tol t versus the tolerance time τ T and the tolerance probability β T was also investigated and the results are as expected (and therefore not shown).A longer tolerance time will shift π tol t to later times.For example, for the RCP4.5 scenario π tol t = 2071, 2088 and 2116 for τ T = 0, 20 and 50 years for the bistable case (for fixed β T = 0.9).With a fixed τ T = 20 years, the value of π tol t shifts to smaller values when the tolerance probability is increased.For example, for β T = 0.80 and 0.99, the values of π tol t are 2127 and 2058, respectively, for the bistable case (for β T = 0.9, π tol t = 2088, see Fig. 4).at each time interval, a χ 2 test was used to analyse the PDFs.For each time, the value of χ 2 > 0.05 and therefore the assumption that the PDF of the GMST is normally distributed appears justified.
The Green's functions for the expectation value and variance of GMST have been calculated with the instantaneously doubling CO 2 profile and the associated ensemble.From the ensemble, at each point in time the expectation value and variance are calculated to obtain the temporal evolution of these two variables.Subsequently, we have found the Green's functions using (5).To check whether these Green's functions perform well, we compared the temporal evolution of the expectation value and variance of the GMST under the 1 % forcing (calculated with (4)) with those directly generated with PlaSim (Fig. 5).The expectation value determined with LRT is close to the one directly generated by PlaSim.However, the variance of the ensemble generated by PlaSim is a lot noisier than the one calculated with LRT.Although the Green's function of the variance provides only a rough approximation, it has the right order of magnitude and we will use it to calculate the variance of the GMST for other forcing scenarios.

Results: Point of No Return under CO 2 eq control
We first consider the case without a carbon cycle model, again assuming that the CO 2 eq concentration can be controlled directly.The scenarios F λ chosen for use in PlaSim are exponentially decaying to different stabilisation levels (varying between 400 and 550 ppmv, see Edenhofer et al. (2010)).
This stabilization level is taken as the parameter λ.We assume that stabilisation happens within 100 years, which corresponds to an e-folding time τ d of about 25 years; the mitigation scenarios F λ are then given by  where t c is again the time at which the mitigation scenario is applied and C tc the associated CO 2 eq concentration.The most extreme mitigation scenario in F λ in terms of CO 2 eq decrease is the one that stabilises at a CO 2 eq concentration of 400 ppmv.
We next determine the PNR π 2100 t by requiring that the GMST must be viable in 2100 using a tolerance probability of β T = 0.90.Furthermore, the viable region is set at T ≤ 16.15 • C, which corresponds to temperatures less than 2 • C above the preindustrial GMST.The values of π 2100 t for all the RCP scenarios are plotted in Fig. 6a.Solid curves show the RCP scenarios while dashed curves present the most extreme scenario F λ .For RCP8.5, π 2100 t is 10 years earlier than for RCP6.0, since the CO 2 eq concentration increases much faster for the RCP8.5 scenario.The mitigation scenario after the point of no return, represented by the dashed line, is the same for all RCP scenarios.This is related to our definition of π 2100 t , where it is required that the GMST is viable in 2100.The mitigation scenario that is plotted is the ultimate scenario that guarantees this.It indicates that for each CO 2 scenario the associated π 2100 t is given by the intersection of that CO 2 eq scenario and the mitigation scenario.This is because it is considered that an exponential decay to 400 ppmv within 100 years is always possible, no matter the CO 2 eq concentration at t c .However, when this concentration becomes too high, this mitigation scenario is not very realistic anymore.
The influence of the tolerance probability on π 2100 t for the RCP4.5 scenario is plotted in Fig. 6b, where we only consider a tolerance probability of 0.8, 0.9 and 0.99.When the tolerance probability is higher, it takes longer before the GMST will be viable again and thus the PNR π 2100 t will be earlier.
However, the differences are very small, since the mitigation scenarios that guarantee viability in 2100 for the different tolerance probabilities are very close.

Results: Point of No Return under emission control
Finally, we consider the more realistic case where emissions are controlled and a carbon model converts emissions to CO 2 eq.A simple carbon model relating emissions E to concentrations C is given by where C CO2,0 is the initial concentration.The Green's function for CO 2 is taken directly from Joos et al. ( 2013): where the parameters are shown in Table 2.The quantity E CO2 is the CO 2 emission in ppm yr −1 that has been converted from GtC yr −1 using the Carbon molecular weight as E CO2 [ppm yr −1 ] = γ E CO2 [GtC/yr], with γ = 0.469 69 ppm GtC −1 .The emissions for the RCP scenarios are taken from Meinshausen et al. (2011) 1 .The carbon model underestimates CO 2 levels for very high emission scenarios as it does not include saturation of natural CO 2 sinks.
Following Table 8.SM.1 of Myhre et al. (2013) we obtain the changes in radiative forcing compared to preindustrial (in W m −2 ) due to changes in CO 2 as where C 0 is the pre-industrial (1750) CO 2 concentration.
We use the same PlaSim ensemble of instantaneous CO 2 doubling runs again to determine a Green's function that relates radiative forcing changes to temperature changes as where G T is the data-based function determined from LRT.The total radiative forcing is taken as where we introduce a scaling constant A to correct for the high climate sensitivity of the PlaSim model compared to typical CMIP5 models.Based on trial runs attempting to reconstruct mean CMIP5 RCP temperature trajectories with RCP CO 2 emissions we choose A = 0.6.
For PlaSim, the Green's function G T , as determined through LRT, is well approximated by a one-time scale exponential: with b 1 = 0.25 K W −1 m 2 yr −1 and τ b1 = 4.69 yr.We determine a Green's function for the temperature variance in the same way.
where τ e = 25 yr is the e-folding timescale of the emission reduction that we keep constant.Using The counter-intuitive lowering of the curve for RCP2.6 (also slightly for RCP4.5) is due to very fast emission reductions in these RCP scenarios.So starting emission reduction at later times may lead to lower total emissions (and hence, temperatures).The CO 2 concentration for the four RCP scenarios as computed by the model (solid) and following the exponential mitigation starting from the Point of no Return (dashed) are shown in Fig. 7b.Note how emissions 'still in the pipeline' lead to CO 2 increases even after the reduction is initiated.Note that this approach does not factor in the uncertainty in the carbon model as we do not have a Green's function propagating the carbon uncertainty through the temperature model.Including this would very likely increase the variance in the pdf and move the Point of No Return to an earlier year.On the other hand, the PLASIM variance is quite small, so the 90% threshold is not vastly different from the mean.2014) stated with high confidence that: "Without additional mitigation efforts beyond those in place today, and even with adaptation, warming by the end of the 21st century will lead to high to very high risk of severe, widespread and irreversible impacts globally".If no measures are taken to reduce GHG emissions during this century and neither will there be any new technological developments that can reduce GHGs in the atmosphere, it is likely that the GMST will be 4 • C higher than the preindustrial GMST at the end of the 21st century (Pachauri et al., 2014).Consequently, it is important that anthropogenic emissions are regulated and significantly reduced before widespread and irreversible impacts occur.It would help motivate mitigation to know when it is 'too late'.
In this study we have defined the concept of the Point of No Return (PNR) in climate change more precisely, using stochastic viability theory and a collection of mitigation scenarios.For an energy balance model, as in section 3, the probability density function could be explicitly computed and hence stochastic viability kernels could be determined.The additional advantage of this model is that one can easily construct a bistable regime, so one can investigate the effects of tipping behavior on the PNR.We used this model (where the assumption was made that CO 2 could be controlled directly instead of through emissions), to illustrate of concept of PNR which is based on a tolerance time for which the climate state is non viable.For the RCP scenarios considered, one finds that the PNR is smaller in the bistable than in the monostable regime of this model.The occurrence of possible transitions to warm states in this model indeed cause the PNR to be 'too late' earlier.
The determination of the PNR in the high-dimensional PlaSim climate model, however, shows the key innovation in our approach, i.e. the use of linear response theory (LRT) to estimate the probability density function of the GMST.PlaSim was used to compute another variant of a PNR based on only requiring that the climate state is viable in the year 2100.Hence, the PNR here is the time such that no allowed mitigation scenario can be chosen to keep GMST below a certain threshold at year 2100 with a specified probability.In the PlaSim results, we used a viability region that was defined as GMSTs lower than 2 • C above the pre-industrial value, but with our methodology, the PNR can be easily determined for any threshold defining the viable region.The more academic case where we assume that GHG levels can be controlled directly provides PNR (for RCP4.5, RCP6.0 and RCP8.5) values around 2050 (section 4.2).However, the more realistic case where the emissions are controlled (section 4.3) and a carbon model is used, reduces the PNR for these three RCP scenarios by about 30 years.The reason is that there is a delay between the decrease in GHG gas emissions and concentrations.
Although our approach provides new insights into the PNR in climate change, we recognize there is potential to substantial further improvement.First of all, the PlaSim model has a too high climate sensitivity compared to CMIP5 models.Although in the most realistic case (section 4.3), we somehow compensate for this effect, it would be much better to apply the LRT approach to CMIP5 simulations.Second, in the LRT approach, we assume the GMST distributions to be Gaussian.This is well justified in PlaSim, as can be verified from the PlaSim simulations, but it may not be the case for a typical CMIP5 model.Third, for the more realistic case in section 4.3, we do not capture the uncertainties in the carbon model and hence in the radiative forcing.
A large ensemble such as that available for PlaSim is not available (yet) for any CMIP5 model.However, we have recently applied the same methodology to two CMIP5 ensembles of models, i.e. a 34 member ensemble of abrupt CO 2 quadrupling and a 35 member ensemble of smooth 1% CO 2 increase per year.The CO 2 quadrupling ensemble was used to derive the Green's function and then the 1% CO 2 increase ensemble was used as a check on the resulting response.The probability density function of GMST increase is close to Gaussian for the 1% CO 2 increase ensemble but clearly deviates from a Gaussian distribution for the 4x CO 2 forcing ensemble, particularly at later times.Although the ensemble is relatively small and the models within the ensemble are different (but many are related), the results for the LRT determined GMST response (Aengenheyster, 2017) are surprisingly good.This indicates that the methodology has a high potential to be successfully applied to results of simulations of CMIP5 (and in the future CMIP6) models.The applicability of LRT to other observables than GMST can in principle be performed but the results may be less useful (e.g.due to non-Gaussian distributions).
Because PlaSim is highly idealized compared to a typical CMIP5 model one cannot attribute much importance to the precise PNR values obtained for the PlaSim model as in Fig. 7.However, we think that our approach is general enough for handling many different political and socio-economical scenarios combined with state-of-the-art climate models when adequate response functions of CMIP5 models have been determined (e.g. using LRT).Hence, it will be possible to make better estimates of the PNR for the real climate system.We therefore hope that eventually these ideas on the PNR in climate change will become part of the decision-making process during future discussions about climate change.

Figure 1 .
Figure 1.(a) CO2eq trajectories of the RCP scenarios used by the IPCC in CMIP5.(b) The solid red curve represents a typical RCP scenario.At the time t b the climate state becomes non-viable, while at t = tc a CO2eq reduction F λ applies; at time te, the climate state is viable again.

Figure 2 .
Figure 2. Bifurcation diagram of the deterministic energy balance model for α1 = 0.45 ((a), monostable model) and α1 = 0.2 ((b), bistable model).The solid curve represents a stable equilibrium while a dashed curve represents an unstable equilibrium.

Figure 3 .
Figure 3.The stochastic viability kernels for the monostable and bistable cases forced by the RCP4.5 scenario.The viable region is defined as T ≤ 293K and is indicated by the red dashed line.This plots show, for each combination of T0 and C0, in which stochastic viability kernel these initial values are located.The numbers in the colourbar stand for the β in V β .For convenience, the bifurcation diagram of the deterministic model is also shown.

Figure 4 .
Figure 4.The PNR π tol t for a system forced with different RCP scenarios, tolerance probability βT = 0.9 and tolerance time τT = te − t b = 20 years.The triangles indicate the point of no return for the bistable case andthe squares for the monostable case.The dotted line is the most extreme scenario of F λ with an exponential decay to 280 ppmv and an e-folding time of 9 years.Note that for both cases there is no PNR when the model is forced with the RCP2.6 scenario.

Figure 5 .
Figure 5. (a) The expectation value and (b) variance of GMST generated by PlaSim (orange) and determined through LRT (blue) for the 1% CO2 concentration increase.

Figure 6 .
Figure 6.(a) The PNR π 2100 t for the RCP2.6,RCP4.5, RCP6.0 and RCP8.5 scenarios for a tolerance probability of βT = 0.9 and ∆t = 0.The solid lines represent the RCP scenarios and the dashed line the most extreme scenario from F λ .Note that these dashed lines coincide.(b) The point of no return for RCP4.5 for different tolerance probabilities.
the carbon model we compute the instantaneous CO 2 concentrations for each such scenario and use the Green's functions for GMST mean and variance to determine the PDF at year 2100 for each starting year t b .Assuming Gaussian distributions (as mentioned, this is well satisfied for the original PlaSim ensemble) we can then easily determine the temperature threshold below which 90% of the values fall.The first year for which this threshold is above 2 K gives π 2100 t .Note that the value of t c (in Fig.1b) where the CO 2 starts to decrease is determined by the coupled carbon-climate model.The warming in 2100 predicted by our simple climate model when starting exponential CO 2 emission reduction in a given year is shown in Fig.7a.The intersections between the RCP curves (solid color) and the dashed line (giving 2 K warming) provide values of π 2100 t .Values do not differ much for the different RCP (4.5, 6.0 and 8.5) scenarios and are before 2030.RCP2.6 does not have a Point of No Return as its emission scenario is sufficient to keep the warming safely below 2 K.

Figure 7 .
Figure 7. (a) Warming in 2100 when starting exponential CO2 emission reduction in a given year.(b) CO2 concentration for the four RCP scenarios as computed by the model (solid) and following exponential mitigation starting from the Point of no Return (dashed).

Table 1 .
Value and meaning of the parameters in the energy balance model given by Eq. (6).

Table 2 .
Model Parameters.No units are given for dimensionless parameters At a year t b > 2005 we start reduction of emissions at an exponential rate, i.e. for t > t b the emissions follow