Investigating climate tipping points under various emission reduction and carbon capture scenarios with a stochastic climate model^††thanks: Accepted for publication in Proc. Roy. Soc. A

It is highly likely that anthropogenic greenhouse gases are responsible for more than half of the increase in the global mean surface temperature between 1951 to 2010 [1]. Therefore, reducing the atmospheric concentration of greenhouse gases must be a central component of any climate change mitigation strategy. This reduction can be achieved in two ways. One involves reducing CO₂ emissions through alternative energy sources, increased fuel efficiency, and reduced consumption [2]. A second approach is expanding carbon sinks, e.g., by employing carbon capture and storage (CCS) technologies, which capture CO₂ from large emitting sources (such as factories) or directly from the atmosphere and prepares it for long-term storage [3, 4].

What complicates these matters is the possible existence of climate tipping points, i.e., climate regimes where small changes significantly alter the future state of the system [5, 6, 7]. This tipping behavior can involve the sudden disruption of climatological and ecological processes such as melting of ice sheets or large-scale death of rainforests [8]. Even if emission reduction and carbon capture strategies lead to long-term CO₂ reduction, its transient growth can trigger a climate tipping point with adversely irreversible impact. Therefore, mitigation strategies not only have to ensure long-term reduction of greenhouse gasses, but also ensure that their transient response remains below the critical tipping point levels.

The purpose of the present study is to quantify emission reduction and carbon capture scenarios that mitigate climate tipping points. To this end, we use a Budyko–Sellers-type model of the climate [9, 10]. This energy balance model assumes that the Earth’s radiation output is balanced by radiation input, and represents key aspects of glacial-interglacial climate transitions [11]. It relies on the hypothesis that glacial cycles are triggered by variations in the Earth’s orbit which alter the incoming solar radiation to the earth. These orbital cycles are amplified by limited feedback between greenhouse gases and temperature [12]. The resulting governing equations are a one-way coupling between the global mean surface temperature and CO₂ concentration, where the temperature is affected by CO₂ concentration through the greenhouse effect.

We model the evolution of the CO₂ concentration by a stochastic delay differential equation (SDDE), which takes into account the carbon emission and capture rates. This model allows us to parameterize the possible reduction in CO₂ emission rates as well as the capability of carbon sinks to absorb atmospheric CO₂. In particular, we show that, even when the emission rates are reduced, the CO₂ concentration exhibits a transient growth which may instigate a climate tipping point. We investigate the response of global mean surface temperature to various emission reduction and carbon capture scenarios, identifying the scenarios which will mitigate the climate tipping point transition.

We use two stochastic differential equations for the climate system, modeling the global mean surface temperature, $T$, and average concentration of CO₂, $C$, in the atmosphere. The equation for temperature is a Budyko-Sellers-type model, derived from the balance between incoming and outgoing radiations [9, 10]. The resulting temperature model reads,

where the temperature $T$ is in units of Kelvin, $C$ is the concentration of CO₂ in parts per million (ppm), and $c_{T}$ is the thermal inertia in units of $\mbox{Jm}^{-2}\mbox{K}^{-1}$. The term $Q_{0}(1-\alpha(T))$ represents short-wave radiation from the Sun, where $Q_{0}$ is the solar input in units of $\mbox{Wm}^{-2}$. The multiplier $\alpha(T)$ denotes temperature dependent albedo that accounts for the light reflecting off the Earth surface. The term $S+A\ln(C/C_{p})$ models the effect of greenhouse gases, where $A$ (in units of $\mbox{Wm}^{-2}$) is the direct forcing of CO₂ and determines the sensitivity of the climate equilibrium [29]. The parameter $S$ represents the trapping of outgoing radiations by greenhouse gases [11]. The constant $C_{p}$ denotes the preindustrial concentration of CO₂ in units of ppm. Finally, the last term $-\epsilon\sigma T^{4}$ represents long-wave radiation from the Sun, where $\sigma T^{4}$ represents the outgoing long wave radiation that is modified by the emissivity $\epsilon$.

Our choice of the temperature dependent albedo function $\alpha(T)$ shown in Fig. 1 is similar to [30], although we use a different set of parameters as reported in Table 1. In particular, we define

which transitions smoothly between albedo parameters $\alpha_{1}$ and $\alpha_{2}$, where $\alpha_{1}$ represents the current global albedo, while $\alpha_{2}$ represents the global albedo of the Earth. The case $\alpha_{1}>\alpha_{2}$ corresponds to the melting of ice on the Earth’s surface, whereas the case $\alpha_{1}<\alpha_{2}$ represents the formation of ice. Here, we only consider the case $\alpha_{1}>\alpha_{2}$ which conforms to the present trends.

The melting of ice depends on the temperature levels. The function $\Sigma(T)$ is chosen to reflect this temperature dependence, so that $\alpha(T)$ transitions smoothly between the albedo $\alpha_{1}$ at temperature $T_{1}$, to the threshold of $\alpha_{2}$ at temperature $T_{2}$. More precisely, we define

where

is a smooth approximation of the unit Heaviside function. The parameter $T_{\alpha}$ controls the transition rate between temperatures $T_{1}$ and $T_{2}$.

Following [30], we model the unresolved subgrid processes by a stochastic term $\eta_{T}\mathrm{d}W_{T}$, where $\eta_{T}>0$ is a constant amplitude and $W_{T}(t)$ denotes the standard Wiener process. The subscript $T$ is used to distinguish the stochastic processes affecting the temperature from those affecting CO₂, to be described shortly. Adding this stochastic term, the temperature model reads

where $F(T,C)$ is defined in equation (1). Note that, although the driving stochastic force is a Gaussian white noise, the response $T$ itself is non-Gaussian and non-white owing to the nonlinear nature of $F(T,C)$.

To close the equations, we also need to model the evolution of CO₂. Dijkstra and Viebahn [29] prescribe $C(t)$ as an explicit function of time. Ashwin and von der Heydt [30] model the CO₂ evolution as a random walk confined to an interval $[C_{1},C_{2}]$. To mimic the current trends and to allow for possible emission reduction and carbon capture, we model the CO₂ evolution as a stochastic delay differential equation,

where $\beta(t)$ is a time-dependent emission rate, $G$ is a time-independent carbon capture rate, and $\tau$ is a time delay. The uncertainties in CO₂ concentration are modeled with the stochastic term $\eta_{C}\mathrm{d}W_{C}$, where $\eta_{C}$ is a constant noise intensity and $W_{C}(t)$ is a standard Wiener process.

The source term in equation (6) models the CO₂ emission into the atmosphere. We allow for a time-dependent rate $\beta(t)$ to model reduction or enhancement of the CO₂ emissions. The sink term, on the other hand, models the CO₂ captured from the atmosphere either through natural phenomena (e.g., photosynthesis [31]) or through artificial technologies (e.g., chemical looping [32]). For the CO₂ sinks, we allow for a time delay $\tau$ to model the time between carbon capture and its effect being felt in the atmospheric concentration. This type of delay is typical in control theory, where there is often a lag between a modifying action and its actualization [33, 34]. If this delay is non-existent or negligible, one can set $\tau=0$. Here, we set $\tau=1$ year, representing the time it takes for CO2 data to be updated and the carbon capture strategies adjusted correspondingly. Nonetheless, we have varied the time delay in the interval $1-20$ years (not shown here) and only observed variations of approximately $10$ ppm in the CO₂ evolution, which do not significantly alter the results.

Equations (5) and (6) form a closed set of equations modeling the global mean temperature $T$ and CO₂ concentration $C$. Our main goal here is to investigate whether a combination of emission reduction and carbon capture may avert a possible upcoming climate tipping point. To this end, we consider various scenarios in terms of emission reduction and carbon capture as discussed in section 3.

The parameter values used in the model are listed in Table 1 and their choice is discussed in Appendix B. For a prescribed CO₂ concentration, the temperature model (5) exhibits three distinct regimes as shown in Fig. 2. Regime 1: If $C<C_{1}=378$ ppm, the temperature has a single stable equilibrium satisfying $F(T,C)=0$.

Regime 2: For CO₂ concentrations $C\in(C_{1},C_{2})$ a bifurcation takes place, whereby two additional equilibria are born. One of the new equilibria is unstable (dashed line in Fig. 2) whereas the other one is stable. As a result, in this regime the temperature model is bistable. Regime 3: If $C>C_{2}=478.6$, the model switches back to a single stable equilibrium. However, in this regime, the global mean temperature is significantly higher than the equilibrium temperature in regime 1.

In our model, the CO₂ concentration is not constant, instead it evolve according to the SDDE (6). We choose the initial conditions $T(0)=15\degree$C and $C(0)=410\,$ppm which reflect the current global mean temperature and CO₂ concentration, respectively [35, 36, 37]. The model parameters, listed in Table 1, are chosen such that the current climate lies in bistable regime 2 near the lower branch of equilibria. If the current trends continue, i.e., the CO₂ concentration keeps increasing, the temperature $T$ also continues to increase gradually. Eventually, one reaches the tipping point $C=C_{2}$ where a slight increase in CO₂ concentration leads to a dramatic increase in the temperature by about six degrees in Celsius. Our goal is to determine the emission reduction and carbon capture scenarios that would avert this catastrophic climate tipping point.

To this end, we allow for a time-dependent emission rate $\beta(t)$ as shown in figure 3. It contains three adjustable periods. First is a period of inaction, up to time $t=t_{1}$, where the emission rate remains constant at its current level $a$. It is followed by a reduction period, $t_{1}<t<t_{2}$ where the emission rate decreases linearly toward its terminal value $b=a/n$. This reduction continues for a period of $\Delta t=t_{2}-t_{1}$ years until it plateaus at time $t=t_{2}$. The period $t>t_{2}$ constitutes the terminal stage where the emission rate remains at the constant level $b$.

Therefore, the emission rate has the general form

where $b=a/n$ for an integer $n\geq 2$. The current emission rate $a$ is estimated from the data in Ref. [38]. They estimate that there are $11.8\pm 0.9$ gigatons of carbon emissions each year from industrial processes and land usage. Converting this to concentration yields $5.556\pm 0.4237$ ppm added to the atmosphere each year. The current emission rate $a$ is then estimated from $aC(0)=5.556/(365\cdot 24\cdot 60\cdot 60)$ ppm/s, where $C(0)\simeq 410$ ppm is the current CO₂ concentration. This yields $a=4.2971\cdot 10^{-10}\;\mbox{s}^{-1}$. In the following sections, we study the evolution of the global mean temperature $T$ under various choices of the remaining free parameters, $t_{1}$, $t_{2}$ and $n$.

To conclude the choice of parameters, we need to specify the parameters $G$ and $\tau$ in the CO₂ model (6). Recall that the parameter $G$ denotes the capacity of carbon sinks to remove CO₂ from the atmosphere. Assuming that carbon capture takes place only due to natural phenomena, the value of $G$ can be estimated from the data in Ref. [38]. They estimate that approximately $6.5$ gigatonnes of carbon are taken from the atmosphere each year from natural processes such as photosynthesis or absorption by the oceans. Converting this to concentration yields $3.0603$ ppm removed from the atmosphere each year. The parameter $G$, then can be estimated by $GC(0)=3.0603/(365\cdot 24\cdot 60\cdot 60)$ ppm/s, which yields $G=2.3669\cdot 10^{-10}\;\mbox{s}^{-1}$.

Finally, we allow for a delay $\tau$ in the CO₂ sinks, which can be set to zero if no such delay exists. We have not been able to determine this parameter from the available literature. The results in Section 5 are reported for the delay time $\tau$ equal to one year. However, we examined delay times up to three years and did not observe a significant effect on the results.

The system of equations (5) and (6) is a one-way coupling between the temperature $T$ and the CO₂ concentration $C$, whereby the CO₂ concentration evolves independently of the temperature. As discussed in section 3, if the CO₂ concentration increases beyond $C_{2}=478.6\,$ppm, the climate system undergoes a tipping point transition leading to an abrupt increase of the global mean surface temperature (see Fig. 2).

The main objective of this paper is to determine emission reduction and carbon capture scenarios that ensure the mitigation of this tipping point. Extreme events, such as tipping point transitions, have been the subject of much research, with an emphasis on their causal mechanisms [39, 40, 41], probabilistic quantification [42, 43, 44, 45], and data-driven prediction [46, 47, 48, 18]. Only recently, control strategies for mitigating extreme events have been proposed [46, 34, 49]. In particular, Farazmand [34] proposes a time-delay feedback control for mitigating noise-induced transitions in multistable systems. This control strategy relies on the stationary equilibrium density of the system. Given the time-dependent emission rate $\beta(t)$ and the linearity of the CO₂ model (6), it does not possess such an equilibrium density. Consequently, the framework of Ref. [34] is not applicable here.

Therefore, we take a different approach here and derive a quantitative upper bound for the CO₂ concentration which determines whether the climate tipping point transitions can be averted. To this end, we consider the CO₂ model (6) without the stochastic term,

Note that (8) is a delay differential equation which requires the initial condition $C(s)$ to be specified as a function over the interval $[-\tau,0]$. Since we consider a relatively short delay $\tau$ of one year, the CO₂ concentration does not change significantly over this time interval. As a result, we assume the constant initial condition $C(s)=C_{0}$ for all $s\in[-\tau,0]$, where $C_{0}=410$ ppm is the current level of CO₂ concentration. The following theorem provides an upper bound for the solutions $C(t)$. Adding the stochastic term $\eta_{C}\mathrm{d}W_{C}$ only leads to small fluctuations around this upper bound, without fundamentally altering the results.

Recall that the climate tipping point occurs if the CO₂ levels exceed $C_{2}=478.6$ ppm. Therefore, ensuring that the upper bound (9) is uniformly below $C_{2}$ is a sufficient condition for avoiding the tipping point. This upper bound encapsulates the competition between CO₂ emission rate $\beta(t)$ and the carbon capture rate $G$. The negative-valued function $\gamma(t)$, appearing in the exponent of the upper bound, is proportional to the carbon capture rate $G$. But it also depends on the emission rate $\beta(t)$ due to the delay $\tau$. This leads to a non-trivial dependence of the upper bound on the carbon emission and capture rates. In section 5, we investigate the shape of this upper bound for various parameter values.

The upper bound (9) also provides a sufficient condition for the asymptotic decay of CO₂ concentration.

We emphasize that, to avert the climate tipping point, it is not sufficient for the CO₂ concentration to decay asymptotically. As we show in Section 5, even a transient growth of CO₂ that exceeds the critical threshold $C_{2}$ will lead to a tipping point transition.

In this section, we present the numerical results obtained from the model (5)-(6). First, we focus on the CO₂ model (6) and investigate the behavior of its upper bound (9). Recall that to avoid the climate tipping point it is sufficient for this upper bound to remain below the critical level $C_{2}$.

Figure 4 shows the upper bound (9) with the delay time $\tau=1$ year, emission reduction time span $\Delta t=50$ years, and the carbon capture rate $G=G_{0}=2.37\cdot 10^{-10}\;\mbox{s}^{-1}$. The initial CO₂ level is set at $C(0)=C_{0}=410\,$ppm, which is the estimated CO₂ concentration in the year 2019 [37]. Figure 4 represents the optimistic scenario where emission reductions begin immediately ($t_{1}=0$) and continue to decrease linearly for $\Delta t=50$ years. The figure shows three terminal emission rates $b=a/n$ with $n=2,5,10$. In every case, we observe a transient growth of CO₂ levels. But none of them reach the tipping point $C_{2}$ and they decay asymptotically.

On the other hand, figure 4 shows the case where the transient growth surpasses the tipping point $C=C_{2}$. In this figure, we fix the terminal CO₂ emission rate at $b=a/3$ with the time window for reaching this rate being $\Delta t=50$ years. We vary the number of years $t_{1}$ before the linear reduction in emission rate begins. If reductions begin in $t_{1}=5$ years, the tipping point can still be avoided. However, if $t_{1}\geq 10$ years, the tipping point may be reached and consequently the global mean surface temperature may sharply increase by about six degrees. These observations accentuate the need for immediate action on carbon emission reduction.

We emphasize that in every case shown in figure 4 the CO₂ levels decay asymptotically. However, the more germane factor is whether the transient growth of CO₂ would exceed the climate tipping point $C=C_{2}$.

Next we consider the full model, i.e., the couple system of equations (5)-(6). Recall that this is a one-way coupling, where the temperature $T$ is affected by the CO₂ concentration through the green house effect. In contrast, the CO₂ concentration evolves independently of the temperature.

This system of stochastic delay differential equations is integrated numerically using the predictor-corrector scheme developed by Cao et al. [50]. The numerical integration is carried out after nondimensionalizing the equations by defining $\hat{C}=C/C_{p}$, $\hat{T}=T/T_{0}$, and $\hat{t}=t/t_{0}$, where $C_{p}=280\,$ppm denotes the preindustrial CO₂ level, $T_{0}=288\,$K denotes the emissivity threshold given by Ashwin and von der Heydt [30], and $t_{0}=10^{7}\,$s is an arbitrary time scale (approximately one-third of a year). The time step of the numerical integrator in the non-dimensional time $\hat{t}$ is $0.0086$, which is equivalent to one day in dimensional time.

To investigate the behavior of the climate system with regards to various emission reduction scenarios, we simulate the model over a 200 year time span for a range of parameters $n$ and $\Delta t$. Recall that $n$ determines the terminal emission rate $b=a/n$ and $\Delta t$ determines the time it takes to reach this plateau (see figure 3). For each set of parameters $(n,\Delta t)$, we simulate $10^{3}$ realizations of the stochastic climate model and compute the ensemble average of the global mean surface temperature $T$ and the CO₂ concentration. In all simulations, the initial conditions are set to $T(0)=15^{\circ}\,$Celsius and $C(0)=410\,$ppm, which reflect the estimated values in the year 2019 [35, 36].

Figure 5 shows the ensemble means of CO₂ concentration $C$ and temperature $T$ as a function of time and the parameters $(n,\Delta t)$. This figure has three panels as described below:

1. 1.

   Figure 5(a): $t_{1}=0$ and $G=G_{0}=2.37\cdot 10^{-10}\,\mbox{s}^{-1}$.

2. 2.

   Figure 5(b): $t_{1}=0$ and $G=G_{0}/2=1.18\cdot 10^{-10}\,\mbox{s}^{-1}$.

3. 3.

   Figure 5(c): $t_{1}=25\,$years and $G=G_{0}=2.37\cdot 10^{-10}\,\mbox{s}^{-1}$.

Figures 5(a)-(b) correspond to the scenario where emission reductions begin immediately ($t_{1}=0$). In panel (a), the carbon capture rate is equal to the empirical value $G_{0}=2.37\cdot 10^{-10}\,\mbox{s}^{-1}$. In this case, the CO₂ concentration increases transiently, but does not reach its critical value $C_{2}$. As a result, no tipping point phenomena is observed. In fact, the temperature does not exceed $16$ degrees Celsius and decays towards $12$ degrees asymptotically.

Panel (b) is identical to (a) except that we use a lower carbon capture rate $G=G_{0}/2=1.18\cdot 10^{-10}\,\mbox{s}^{-1}$. The carbon capture rate may decrease over time due to various factors, but most prominently due to deforestation [51, 52]. In this case, the CO₂ concentration still undergoes a transient growth before eventually decreasing. Because of the lower carbon capture rate, however, the transient growth surpasses the critical threshold $C_{2}$, and leads to a drastic increase in the global mean surface temperature. This tipping point does not occur for all carbon emission scenarios. If the transition period $\Delta t$ is short enough and the terminal emission rate $b=a/n$ is small enough, the tipping point can still be mitigated.

In fact, as shown in figure 6(a), there is a nonlinear boundary in the parameter space $(n,\Delta t)$ which separates the tipping point transitions from its successful mitigation. For $n=2$, the tipping point occurs regardless of the transition time $\Delta t$. Therefore, the terminal carbon emission rate has to reduce to at least one-third of its current value to mitigate the climate tipping point. Larger $n$ allows for longer transition period $\Delta t$. In figure 6(b), we show the ensemble average of the temperature $T$ for three parameter values $(n,\Delta t)$. The shaded areas mark one standard deviation from the mean. The point on the boundary $(n=4,\Delta t=29)$ has a large standard deviation because the stochastic process $\eta_{c}\mathrm{d}W_{c}$ can kick the trajectory towards the tipping regime or away from it. The parameters $n=3$ and $\Delta t=29$ lead to the tipping behavior with high probability, where the temperature rises to about $24^{\circ}$ Celsius. In contrast, $n=6$ successfully averts the climate tipping point, after a slight transient increase in the global mean surface temperature.

Finally, figure 5(c) shows the case where emission reduction is delayed by 25 years ($t_{1}=25$), and the carbon capture rate stays at its current level ($G=G_{0}$). In this case, CO₂ concentration again exhibits a transient growth past the critical value $C_{2}$, and as a result, a drastic increase in global temperature ensues. Although the CO₂ concentration reduces asymptotically to as low as $150\,$ppm, the transient tipping point will have dire consequences, such as sea level rise and droughts. This observation reaffirms the need for immediate action towards significant emission reduction. Otherwise, the only alternative would be to significantly increase the CO₂ sinks $G$, through artificial carbon capture technologies.

We considered a stochastic climate model as a one-way coupling between CO₂ concentration $C$ and the global mean surface temperature $T$. The CO₂ concentration is governed by a stochastic delay differential equation, allowing for the modeling of various emission reduction and carbon capture scenarios. The temperature $T$ satisfies a stochastic differential equation derived from energy balance (Budyko–Sellers model). The temperature is coupled with the CO₂ equation to model the effect of greenhouse gasses.

The model has a tipping point behavior: if the CO₂ concentration exceeds the critical value of $C=478.6$ ppm, the global mean surface temperature will increase abruptly by about six degrees Celsius. The CO₂ model exhibits transient growth which allows for reaching this tipping point in finite time, even when CO₂ decreases asymptotically. We derived a tight upper bound for the CO₂ concentration, which provides sufficient conditions for mitigating the climate tipping point (see Theorem 4.1). This upper bound depends on several parameters such as the CO₂ emission rate $\beta(t)$, the carbon capture rate $G$, and the delay time $\tau$. However, since the upper bound is explicitly known, it can be analyzed numerically with low computational cost in order to determine the emission reduction and carbon capture scenarios which would mitigate the climate tipping point.

We examined various emission reduction scenarios by combining our analytic results with Monte Carlo simulations of the climate model. In particular, we find that the climate tipping point in our model can be averted if CO₂ emission reductions begin immediately and the emission rate decreases to one-third of its current level within 50 years. However, if these reductions are delayed by as much as ten years, the transient growth of the CO₂ concentration will exceed the tipping point value, leading to a drastic and abrupt increase in the global mean surface temperature. In the latter case, increasing the carbon capture rate could still avert the tipping point.

It remains to be seen whether our conclusions carry over to more sophisticated climate models, such as box models or general circulation models. We will investigate such models in future work. Although the temperature model would be more complex, our analysis of the CO₂ concentration can be used to derive sufficient conditions that guarantee tipping point evasion. Furthermore, deriving an invariant probability distribution for the CO₂ model is of great interest as it would quantify the probability of transitions which in turn leads to necessary and sufficient conditions for mitigating the climate tipping point. Finally, the CO₂ model itself can be improved, e.g., by using data-driven discovery methods to obtain simple models that agree with observational data [53].

The data and code for reproducing the results of this manuscript are publicly available at https://github.com/mfarazmand/ClimateMitigation

We are grateful to Prof. Tamás Tél for recommending reference [30], and Prof. Pedram Hassanzadeh for fruitful conversations. This work was partially supported by the National Science Foundation grant DMS-1745654.