Articles | Volume 27, issue 3
Research article
08 Jul 2020
Research article |  | 08 Jul 2020

Anthropocene climate bifurcation

Kolja Leon Kypke, William Finlay Langford, and Allan Richard Willms

This article presents the results of a bifurcation analysis of a simple energy balance model (EBM) for the future climate of the Earth. The main focus is on the following question: can the nonlinear processes intrinsic to atmospheric physics, including natural positive feedback mechanisms, cause a mathematical bifurcation of the climate state, as a consequence of continued anthropogenic forcing by rising greenhouse gas emissions? Our analysis shows that such a bifurcation could cause an abrupt change to a drastically different climate state in the EBM, which is warmer and more equable than any climate existing on Earth since the Pliocene epoch. In previous papers, with this EBM adapted to paleoclimate conditions, it was shown to exhibit saddle-node and cusp bifurcations, as well as hysteresis. The EBM was validated by the agreement of its predicted bifurcations with the abrupt climate changes that are known to have occurred in the paleoclimate record, in the Antarctic at the Eocene–Oligocene transition (EOT) and in the Arctic at the Pliocene–Paleocene transition (PPT). In this paper, the EBM is adapted to fit Anthropocene climate conditions, with emphasis on the Arctic and Antarctic climates. The four Representative Concentration Pathways (RCP) considered by the IPCC (Intergovernmental Panel on Climate Change) are used to model future CO2 concentrations, corresponding to different scenarios of anthropogenic activity. In addition, the EBM investigates four naturally occurring nonlinear feedback processes which magnify the warming that would be caused by anthropogenic CO2 emissions alone. These four feedback mechanisms are ice–albedo feedback, water vapour feedback, ocean heat transport feedback, and atmospheric heat transport feedback. The EBM predicts that a bifurcation resulting in a catastrophic climate change, to a pre-Pliocene-like climate state, will occur in coming centuries for an RCP with unabated anthropogenic forcing, amplified by these positive feedbacks. However, the EBM also predicts that appropriate reductions in carbon emissions may limit climate change to a more tolerable continuation of what is observed today. The globally averaged version of this EBM has an equilibrium climate sensitivity (ECS) of 4.34 K, near the high end of the likely range reported by the IPCC.

1 Introduction

Today, there is widespread agreement that the climate of the Earth is changing, but the precise trajectory of future climate change is still a matter of debate. Recently there has been much interest in the possibility of “tipping points” (or bifurcation points) at which abrupt changes in the Earth climate system occur (see Brovkin et al.1998; Ghil2001; Alley et al.2003; Seager and Battisti2007; Lenton et al.2008; Ditlevsen and Johnsen2010; Lenton2012; Ashwin et al.2012; Barnosky et al.2012; Drijfhout et al.2015; Bathiany et al.2016; North and Kim2017; Steffen et al.2018; Dijkstra2019; Wallace-Wells2019). Section 12.5.5 in IPCC (2013) gives an overview of such potential abrupt changes. At such points, a small change in the forcing parameters (whether anthropogenic or natural forcings) may cause a catastrophic change in the state of the system. In order to prepare for future climate change, it is of great importance to know if such abrupt changes can occur and, if so, when and how they will occur. Some authors have suggested that conventional general circulation models (GCMs) may be “too stable” to provide reliable warning of these sudden catastrophic events (Valdes2011), and that the study of paleoclimates may be a better guide to how abrupt changes may occur (Zeebe2011). Steffen et al. (2018) asked the fundamental question: “Is there a planetary threshold in the trajectory of the Earth System that, if crossed, could prevent stabilization in a range of intermediate temperature rises?”

In this paper, a simple energy balance model (EBM) is used to investigate the possible occurrence of such a threshold (tipping point or bifurcation point) in the climate of the Anthropocene. Energy balance models assume that the climate is in an equilibrium state, for which “energy in equals energy out” at each point of the Earth's surface and atmosphere. Thus, time dependence is eliminated from the model, greatly simplifying the analysis. In the literature, EBMs have played an important role in understanding climate and climate change (Budyko1968; Sellers1969; Sagan and Mullen1972; North et al.1981; Thorndike2012; Kaper and Engler2013; Dijkstra2013; Payne et al.2015; Hartmann2016; North and Kim2017; Ghil and Lucarini2020). Many of those EBMs have exhibited bifurcations. The present paper presents an EBM with more accurate diagnostic equations than the early EBMs, built upon the basic laws of geophysics and including nonlinear feedback processes that amplify anthropogenic CO2 forcing. A rigorous mathematical bifurcation analysis of this EBM has been presented in Kypke and Langford (2020). That analysis gave a mathematical proof of the existence of a cusp bifurcation in the EBM, complete with a determination of the centre manifold and of the universal unfolding parameters, as functions of the relevant physical parameters. The existence of the cusp bifurcation implies the coexistence of two distinct stable equilibrium climate states (bistability), as well as the existence of hysteresis; that is, two abrupt one-way transitions between these two states (via fold bifurcations) exist in the EBM. The present paper extends those results from the paleoclimate model in Kypke and Langford (2020) to the Anthropocene climate model studied here. This Anthropocene model gives predictions of climate changes in the 21st century and beyond.

One advantage of an EBM over a more complex GCM is that it facilitates the exploration of specific cause and effect relationships, as particular climate forcing factors are varied or ignored. Another advantage of an EBM is that rigorous mathematical analysis often can prove the existence of certain behaviours, such as bistability and bifurcations, that could only be surmised from numerical evidence, or missed, in more complicated models. Four versions of the EBM are considered here: a globally averaged temperature model and three regional models corresponding to Arctic, Antarctic, and tropical climates.

This EBM was validated in Dortmans et al. (2019), where it was applied to known paleoclimate changes. It successfully “predicted” the abrupt glaciation of Antarctica at the Eocene–Oligocene transition (EOT) and the abrupt glaciation of the Arctic at the Pliocene–Pleistocene transition (PPT), using bifurcation analysis. The transitions in the model were congruent with the abrupt cooling, from warm, equable “hothouse” climate conditions to a cooler state like the climate of today, with ice-capped poles, as indelibly recorded in the geological record at the EOT and PPT.

In adapting the previous paleoclimate EBM to the climate of the Anthropocene, this paper explores the possibility of a “reversal” of those two paleoclimate transitions; that is, a transition from today's climate with ice-capped poles to an equable hothouse climate state, such as existed in the Pliocene and earlier. It provides new mathematical evidence suggesting that catastrophic climate change in polar regions is inevitable in the coming decades and centuries if current anthropogenic forcing continues unabated. The EBM also suggests that if appropriate mitigation strategies are adopted (as recommended by the IPCC), such an outcome can be avoided.

Figure 1A visualization of the energy balance model (EBM). Here Q is the incident solar radiation. A fraction ξA of Q is absorbed by the atmosphere and another fraction ξR is reflected by clouds into space. The resulting solar forcing that strikes the surface is FS=(1-ξA-ξR)Q. The surface has albedo α, which means that αFS is reflected back to space and the remaining energy (1−α)FS is absorbed by the surface. The surface emits longwave radiation of intensity IS, of which a fraction ηIS is absorbed by greenhouse gases in the atmosphere and the remainder (1−η)IS escapes to space. The atmosphere emits longwave radiation IA, of which a fraction βIA goes downward to the surface and the remaining fraction (1−β)IA escapes to space. The three forcing terms (FA, FO, FC) represent atmospheric heat transport, ocean heat transport and vertical conduction/convection heat transport, respectively. Values of these and other parameters are given in Table 1.


Table 1Summary of variables and parameters used in the EBM. The values of the physical constants (ξa, ξR, kC, kW, GC, GW1, GW2) are as determined in Kypke (2019) and Dortmans et al. (2019).

Download Print Version | Download XLSX

The EBM of this paper has been kept as simple as possible, while incorporating the nonlinear physical processes that are essential to our exploration of bifurcation behaviour. In that sense, it follows in the tradition of simple energy balance models of Budyko (1968), Sellers (1969), North et al. (1981), and others. However, this EBM is but the first step in the authors' study of a hierarchy of nonlinear models of increasing complexity. That hierarchy is outlined in the concluding Sect. 4.

2 An energy balance model for climate change

The EBM is a simple two-layer model, with layers corresponding to the surface and atmosphere, respectively; see Fig. 1, which is based on Payne et al. (2015), Trenberth et al. (2009), and Wild et al. (2013). The symbols in Fig. 1 are defined in the caption of Fig. 1 and in Table 1. In the EBM of Fig. 1, the surface receives shortwave radiant energy FS=(1-ξA-ξR)Q from the sun, where Q is the incident solar radiation and ξA and ξR are the fractions of Q absorbed by the atmosphere and reflected back to space by clouds, respectively. The values of ξA and ξR are obtained from Trenberth et al. (2009) and Wild et al. (2013) (see the appendix and Table 1). At the surface, a fraction αFS is reflected back to space, where α is the surface albedo, and the remainder, (1−α)FS, is absorbed by the surface. The surface re-emits longwave radiant energy of intensity IS=σTS4 (Stefan–Boltzmann law) upward into the atmosphere. The atmosphere contains greenhouse gases that absorb a fraction η of the radiant energy IS from the surface, and the remainder, (1−η)IS, escapes to space. The atmosphere re-emits radiant energy of total intensity IA. Of this radiation IA, a fraction βIA is directed downward to the surface, and the remaining fraction (1−β)IA goes upward and escapes to space. Balancing the energy flows represented in Fig. 1 leads directly to the following dynamical system:


where Eq. (1) represents surface energy balance and Eq. (2) represents atmosphere energy balance. Here cS and cA are specific heat rate constants derived in Kypke (2019) and Kypke and Langford (2020) and listed in Table 1. There are three heat transport terms: FA is atmospheric heat transport, FO is ocean heat transport (both horizontally), and FC is conductive/convective heat transport, vertically from the surface to the atmosphere.

In Eqs. (1) and (2) there is an asymmetry, in that temperature TS is used to represent the state of the surface in Eq. (1), but radiant energy intensity IA is used instead of temperature to represent the state of the atmosphere in Eq. (2). Note that either temperature variables or energy intensity variables could have been used in either equation if we assume the Stefan–Boltzmann law (IS=σTS4 and IA=ϵσTA4, where ϵ=0.9 is the emissivity of dry air). The use of TS as state variable in surface Eq. (1) is the obvious choice, since the surface temperature is the most important variable in the EBM. However, the choice of IA instead of TA in Eq. (2) is less obvious. The atmosphere has thickness. In the actual atmosphere, temperature decreases with height above the surface at a rate called the lapse rate. Therefore, there is not just one value of temperature TA for the atmosphere. However, we can define a single value of radiant energy intensity IA, corresponding to the total energy emitted vertically by the slab of atmosphere, and use this instead of temperature in the energy balance Eq. (2). A second reason for the use of IA instead of TA in Eq. (2) is that this facilitates the use of the ICAO Standard Atmosphere lapse rate, as explained in the paragraph following Eqs. (4) and (5) below, and that allows for a more realistic representation of the greenhouse gas behaviour of water vapour.

The first step of the analysis of system Eqs. (1) and (2) is a rescaling of temperature T (in kelvin) to a new nondimensional temperature τ with τ=1 corresponding to the freezing temperature of water (TR=273.15K). Then all variables and parameters in the system can be made nondimensional by the scalings


where s is dimensionless time and χ is the dimensionless heat rate constant. The surface and atmosphere energy balance Eqs. (1) and (2) in nondimensional variables are then


In Fig. 1, the atmosphere is shown to be a uniform slab, even though the actual atmosphere is not a uniform slab. The essential nonlinear processes in the atmosphere, which the model must capture, are the heating effects of the greenhouse gases CO2 and H2O. According to the Beer–Lambert law, the absorptivity of these gases is determined by their optical depths. Therefore, in the model we substitute for optical depth in the slab, the values of optical depth that these gases would have in the ICAO Standard Atmosphere (ICAO1993), which is a good approximation to the actual atmosphere. In the ICAO model, the rate of change of temperature with altitude is assumed to be a negative constant −Γ, called the ICAO lapse rate; see Table 1. The concentration of CO2 is independent of temperature, but the concentration of H2O decreases with altitude as the temperature decreases, according to the Clausius–Clapeyron relation (Pierrehumbert2010). Then the optical depth of H2O is obtained by integration of the Clausius–Clapeyron relation from the surface to the tropopause, resulting in Eq. (6). In this way, the simple slab model has greenhouse effects close to those of these two gases in the actual atmosphere, where the temperature is not constant but decreases with altitude. This use of the ICAO lapse rate differentiates the present EBM from all previous EBMs in the literature. The calculation gives the absorptivity η due to greenhouse gases, in the atmosphere Eqs. (2) or (5), as


where μ is the concentration of CO2 in the atmosphere, measured in molar parts per million; δ is the relative humidity of water vapour (0δ1); γ=Γ/TR is the nondimensionalized lapse rate (ICAO1993); and Z is the tropopause height. (Since methane acts similarly to carbon dioxide as a greenhouse gas, it may be assumed that μ includes also the effects of methane.) Equation (6) is derived using fundamental laws of atmospheric physics: the Beer–Lambert law, the ideal gas law and the Clausius–Clapeyron equation; see Dortmans et al. (2019) for more details. In Eq. (6),


are dimensionless physical constants determined in Dortmans et al. (2019), and kC and kW are absorption coefficients for CO2 and H2O, respectively; see Table 1.

In the surface Eqs. (1) or (4), α is the albedo of the surface (0α1) and α depends strongly on temperature τS near the freezing point (τS=1). Typical values of surface albedo are 0.6–0.9 for snow, 0.4–0.7 for ice, 0.2 for crop land, and 0.1 or less for open ocean. Therefore, in the high Arctic, as the ice-cover melts, the albedo will transition from a high value such as αc=0.7 for snow/ice to a low value such as αw=0.08 for open ocean. Some authors have assumed this change in albedo to be a discontinuous step function (Dortmans et al.2018). However, all variables in this EBM have annually averaged values. As the Arctic thaws, the annual average albedo will transition gradually, over a number of years, from its high value for year-round ice-covered surface to its low value for year-round open ocean. Therefore, in this paper we use a smoothly varying albedo function, which better models this gradual transition from high to low albedo, as the mean temperature rises through the freezing point (τS=1). It is modelled by the hyperbolic tangent function:

(8) α ( τ S , ω ) = 1 2 ( [ α w + α c ] + [ α w - α c ] tanh ( τ S - 1 ω ) ) .

Here the parameter ω controls the steepness of this switch function. Analysis of polar ice data in recent years confirms that this function gives a good fit to the decline of ice cover and albedo in the Arctic with ω=Ω/TR=0.01 (Pistone et al.2014; Dortmans et al.2019). The dependence of Arctic Ocean sea-ice thickness on surface albedo parametrization in models has been investigated in Björk et al. (2013), where alternative albedo schemes are compared. The nonlinear dependence of albedo on temperature, as in Eq. (8), has been shown to lead to hysteresis behaviour (Stranne et al.2014; Dortmans et al.2019).

2.1 Refinement of the paleoclimate EBM to an Anthropocene EBM

Paleoclimate data are difficult to obtain and in general can only be inferred from proxy data. The situation is different for the Anthropocene. There is now an abundance of land-based and satellite climate data. Therefore, the EBM in this paper can be refined to take advantage of the additional data. The appendix details the changes made in this EBM from that which was presented in Dortmans et al. (2019) and Kypke and Langford (2020), to improve its accuracy for the Anthropocene. These changes do not change the fundamental behaviour of the EBM, including the existence of bifurcations, but they do make the numerical predictions more reliable. Table 1 of this paper may be compared with the corresponding Table 1 of Kypke and Langford (2020) to see how parameter values have been updated.

The total absorptivity η, given in Eq. (6) for paleoclimates, is now modified to reflect the fact that clouds absorb a fraction ηCl of the outgoing longwave radiation. In this paper


see the appendix.

The vertical heat transport term FC has been modified to take into account both sensible and latent heat transport (Kypke2019). See the appendix, where the following formula is obtained (here in nondimensional form).


where H1 and H2 are nondimensional constants:


In Eqs. (4) and (5), at equilibrium (i.e. dds^=0), the state variable iA is easily eliminated, leaving a single equation with a single state variable τS,


2.2 Cusp bifurcation in the EBM

In this subsection, we outline the proof that the cusp bifurcation, which was proven to exist in the Paleoclimate EBM (Kypke and Langford2020), in fact persists in this Anthropocene EBM Eqs. (4) and (5). Therefore, the conclusions of that paper carry over to this paper. Readers not interested in these mathematical details may skip this subsection.

The right-hand sides of Eqs. (4) and (5) can be represented by a single vector function F:R2×R4R2. Then an equilibrium point (τ¯S,i¯A) of Eqs. (4) and (5), at which dτSdt=diAdt=0, is a solution of

(12) F ( τ ¯ S , i ¯ A ; ρ ) = 0 ,

where ρ represents four physical parameters that may be varied in the model,

(13) ρ { μ , δ , F O , ω } .

See Table 1 for definitions of these parameters. Since the codimension of the cusp bifurcation is only two, there is some redundancy in the choice of these four parameters. For application to future climates, the parameter pair (FO,μ) is of primary interest. Equilibrium points (τ¯S,i¯A) satisfying Eq. (12) have been computed in Kypke (2019). Having computed the equilibrium point (τ¯S,i¯A) satisfying Eq. (12), the system may be translated to the origin (x,y)=(0,0), in new state variables defined by


and Eqs. (4) and (5) become


For Eq. (14) to have a steady-state bifurcation at the equilibrium point (0,0), the Jacobian J of F in Eq. (12) must have a zero eigenvalue, λ1=0, at that point. (A Hopf bifurcation is not possible in this system.) For stability, the second eigenvalue satisfies λ2<0. The corresponding eigenvectors, e1,e2, form an eigenbasis. A linear transformation takes (x,y) coordinates to eigenbasis coordinates (u,v), where (x,y)=T(u,v), and the columns of T are the normalized eigenvectors e1,e2 in Eq. (17), below. Then in (u,v) coordinates, the 2-D system Eq. (14) becomes



(16) = f C 0 + 4 η 0 τ ¯ S 3 + η 0 τ ¯ S 4 , k = β χ , ϕ = 1 + k


(17) e 1 = 1 , e 2 = - k 1 .

Recall that the parameters μ and δ, representing CO2 concentration and water vapour relative humidity, respectively, enter into Eq. (15) through the function η defined in Eq. (9). For more details, see Kypke (2019), Kypke and Langford (2020), and Kuznetsov (2004).

If the Centre Manifold Theorem applies to Eq. (15), then there exists a flow-invariant centre manifold, which is tangent to the u axis. The applicability of this theorem has been verified, and the centre manifold has been computed for the present Anthropocene EBM as was done for the paleoclimate EBM in Kypke (2019) and Kypke and Langford (2020). Details are omitted here. A phase portrait for Eq. (15) in a neighbourhood of the cusp equilibrium point, together with a portion of this centre manifold (in red), is shown in Fig. 2 in (u,v) coordinates. In this figure, trajectories quickly collapse to the centre manifold around the equilibrium point (0,0), as predicted by the Centre Manifold Theorem. The cusp equilibrium manifold for Eq. (15) in normal form is shown in Fig. 3. Here ζ1 and ζ2 are the normal form unfolding parameters for the cusp bifurcation. Note the coexistence of three equilibrium points with different values of u (two stable and the middle one unstable) inside the cusp-shaped region, but only one equilibrium point (stable) outside of that region.

Figure 2Phase portrait of system Eq. (15), together with a portion of the centre manifold (in red), in the (u,v) eigenbasis coordinates. Parameter values are those at the computed cusp point. The yellow dot marks the cusp equilibrium point. Note the rapid approach to the centre manifold from initial points not on the centre manifold in contrast to the slow evolution along the centre manifold.


Figure 3Cusp bifurcation diagram for the EBM in normal form coordinates. (a) Graph of the equilibrium surface with normal form unfolding parameters (ζ1,ζ2). (b) Projection of this surface in 3-D onto the (ζ1,ζ2) plane. The blue semi-cubic parabola represents fold bifurcations, and it separates the (ζ1,ζ2) plane into two open 2-D regions. Inside the cusp region, that is between the two branches of the semi-cubic parabola, there exist three equilibrium solutions u: two stable and the middle one unstable. Outside of the semi-cubic parabola there exists only one unique equilibrium solution u and it is stable.


3 Anthropocene climate forecasts

In this section, the EBM of Sect. 2 is applied to the present and future climates of the Earth to investigate the possibility of climate bifurcations (or tipping points) in the Anthropocene. The principal parameters chosen to be explored are carbon dioxide concentration μ, ocean heat transport FO, and relative humidity δ. The EBM is adapted locally to three separate regions, namely the Arctic, Antarctic, and the tropics, in Sect. 3.2, 3.3, and 3.4, respectively.

Carbon dioxide production due to human activities has been well documented as a driver of climate change in the Anthropocene. Projections of future atmospheric CO2 levels are available under various future scenarios; we follow the Representative Concentration Pathways of IPCC (2013), which is described in Sect. 3.1. Ocean heat transport is a difficult quantity to predict, as many different factors influence the transport of heat to various regions of the world via the oceans. Changes in temperature can change ocean heat transport which in turn affects local temperatures. This is ocean heat transport feedback, which is explored in Sect. 3.2.2. Similarly, the role of atmospheric heat transport feedback is investigated briefly in Sect. 3.2.2. In addition, water vapour is a powerful greenhouse gas, with a positive feedback effect that is investigated in Sect. 3.2.3.

In Sect. 3.5, a globally averaged model is considered, mainly for the purpose of determining the global equilibrium climate sensitivity (ECS) of this EBM, for comparison with the ECS of other models as reported in IPCC (2013) and Priostosescu and Huybers (2017); see Sect. 3.6.

3.1 Representative Concentration Pathways (RCPs)

The IPCC has standardized four Representative Concentration Pathways (RCPs), which are used for projections of future carbon dioxide concentrations; see van Vuuren et al. (2011) and Box TS.6 in IPCC (2013). These RCPs are scenarios for differing levels of anthropogenic forcings on the climate of the Earth and represent differing global societal and political “storylines”. The scenarios are named RCP 8.5, RCP 6.0, RCP 4.5, and RCP 2.6, after their respective peak radiative forcing increases in the 21st century. That is, in the year 2100, RCP 8.5 will reach its maximum radiative forcing due to anthropogenic emissions of +8.5W m−2 relative to the year 1750. This scenario is understood as one where emissions continue to rise and are not mitigated in any way. RCP 6.0 corresponds to +6.0W m−2 and RCP 4.5 corresponds to +4.5W m−2 relative to 1750. These are stabilization scenarios, where greenhouse gas emissions level off to a target amount by the end of the century. Finally, RCP 2.6 corresponds to +2.6W m−2 in 2100, relative to 1750. This is a mitigation scenario, where strong steps are taken to eliminate the increase and eventually reduce anthropogenic greenhouse gas emissions. Figure 4 shows the carbon dioxide concentrations projected to the year 2500 in the IPCC scenarios for the four RCPs. The carbon dioxide increase is relatively moderate for RCPs 2.6 and 4.5 and even decreasing eventually for RCP 2.6. The increase for RCP 6.0 is larger, and it is drastic for RCP 8.5.

Figure 4Carbon dioxide concentration μ, as projected by the IPCC for each of the four RCP hypothetical scenarios. This figure is generated from data at IIASA (2019).


These RCPs represent hypothetical forcings due to human activity up to the year 2100. Beyond 2100, they assume a “constant composition commitment”, where emissions are kept constant, which serves to stabilize the scenarios beyond 2100 (IPCC2013). Of course, emissions could continue to increase or be greatly reduced (“zero emissions commitment”) at some point in the future. However, the constant emission commitment dataset provided in IPCC (2013), and shown in Fig. 4, is what is utilized in this work. In the following sections we enter the CO2 concentrations μ shown in Fig. 4, one at a time along each of the four IPCC RCPs, into the versions of our EBM for the Arctic. Antarctic, etc., and we then let the EBM climate evolve quasi-statically along each CO2 pathway.

3.2 EBM for the Anthropocene Arctic

Figure 5a is the manifold of equilibrium states for the EBM parameterized by (FO,μ) for the case of an Arctic climate under Anthropocene conditions. Figure 5b is the projection of this manifold onto the parameter plane, showing the fold bifurcation lines as boundaries between coloured regions. Parameter values are as in Table 1, with FA=104W m−2 and δ=0.6 taken as the nominal values for the modern Arctic throughout this section, except in Fig. 8. These figures were computed as in Kypke and Langford (2020). The cusp point, seen in Fig. 3, still exists but is not visible in Fig. 5, because it is outside of the range of parameters included in the figure.

Figure 5Energy balance model of the modern Arctic. Parameter values are as in Table 1. The red dots locate today's Arctic climate conditions. Ocean heat transport FO increases from 0 to 50 W m−2 and carbon dioxide concentration μ increases from 0 to 2000 ppm. (a) 3-D plot of equilibrium manifold. (b) Projection on the (FO,μ) parameter plane.


In Fig. 5a, today's climate state lies on the lower (blue) portion of the equilibrium manifold, as shown by the red dot. The upper (yellow) portion represents a coexisting warm equable climate state, similar to the climate of Earth in the Pliocene and earlier. The intermediate (green) portion represents an unstable (and unobservable) climate state.

Similarly, in Fig. 5b, the blue area represents unique cold stable states, yellow represents unique warm stable states, and the green region is the overlap region between the two fold bifurcations, where both warm and cold states coexist. Hence, on moving in the (FO,μ) parameter space starting from the blue region, crossing the green region, and into the yellow region, there would be no observable change in climate on crossing from blue to green; however, crossing the boundary from green to yellow would cause a catastrophic jump from cold to warm stable climate states. Alternatively, if the (FO,μ) parameter values moved from the yellow, through the green, into the blue region in Fig. 5b, there would be no abrupt change in climate state on crossing the yellow–green boundary, but a sudden transition to a cold state would occur at the green–blue boundary. This behaviour is called hysteresis.

3.2.1 Arctic climate for the four RCPs

The paramount question considered in this paper can now be phrased as follows. Can a bifurcation leading to a warmer and more equable climate state be expected in the EBM if it is allowed to evolve along one of the four RCPs in Fig. 4? In Fig. 5b, this bifurcation would correspond to crossing the line of fold bifurcations separating the green and yellow regions on increasing μ and possibly FO.

Figure 6 shows the increase in surface temperature in the Arctic region, using historical data from the year 1900 to the present, and then the EBM forecasts up to the years 2100 and 2300, holding μ to each of the four RCPs in Fig. 4, and with constants FO=10Wm-2, FA = 104 W m−2, and δ=0.6. The temperature change shown is relative to the Arctic temperature of the EBM in the year 2000, which was −28.90C (τS=0.8942).

Figure 6Arctic surface temperature change, relative to the year 2000 average temperature of -28.90C, for each of the four RCPs with constant FO=10Wm-2 and with δ=0.6. In (a), the EBM temperature change is projected to year 2100, and in (b) it is projected to year 2300, following the assumptions of the RCP pathways in Fig. 4. Note the dramatic jump in temperature on RCP 8.5, resulting from a saddle-node bifurcation, predicted near the year 2170. The vertical arrow is not an actual trajectory of the dynamical system but rather represents the transition to a new equilibrium state that occurs after the disappearance of the saddle node.


Figure 6a may be compared to the results in Figs. AI.8 and AI.9 in IPCC (2013). Those figures use the RCPs of Fig. 4 and an ensemble of climate models to forecast surface temperature changes for the Arctic up to the year 2100 for land and sea separately. IPCC Fig. AI.8 displays the winter months of December–February and AI.9 is for June–August. Figure 6 of this paper does not distinguish surface covering, so it is given as an annual average value. Figure 6a shows predictions of the EBM to the year 2100, which is the same time frame as for the IPCC GCM projections. Both use the RCP scenarios to determine hypothetical CO2 concentrations (μ) by year, up to year 2100. Then both use these CO2 concentrations as input to the respective models (EBM or GCM) to determine predicted climate changes for the same period. They are in good agreement if a weighted average of the sea and land temperature changes are considered and if the winter months are more representative of an annual value for the Arctic climate.

Supported by the agreement until year 2100 between Fig. 6a and the IPCC reported values of temperature change on the RCP paths, the EBM forecast was then extended to the year 2300 (see Fig. 6b), which uses the same parameter values as Fig. 6a. It exhibits a saddle-node bifurcation for RCP 8.5 near the year 2170. Following the disappearance of the cooler equilibrium state in this saddle-node bifurcation, bifurcation theory tells us that there exists a neighbourhood in the state space, where the saddle node once existed, inside of which trajectories move slowly through a so-called “ghost equilibrium” that is a remnant of the saddle node. (The transit time has inverse square-root dependence on the unfolding parameter in that neighbourhood.) Outside of that neighbourhood, trajectories evolve with velocity determined by cS and cA in Eqs. (1) and (2) to the upper stable equilibrium point. (The vertical arrow in Fig. 6b represents that transition but is not an actual solution of the dynamical system, and similarly for the vertical arrows in Figs. 7 and 10.) This bifurcation on RCP 8.5 results in a drastic increase in temperature: the mean Arctic temperature jumps by +26.5C to the new value of +22.5C. This is warmer than the mean annual Arctic temperature in the Pliocene but is consistent with what existed in the Eocene (Willard et al.2019). Because of simplifications made in this EBM, these numbers should not be taken literally as quantitatively accurate forecasts; however, the qualitative existence of a dramatic increase in temperature due to a bifurcation must be taken seriously. The topological methods employed in this work ensure that bifurcation in this model is a mathematically persistent feature that will be manifest in all “nearby” models; see Kypke and Langford (2020).

The other three RCPs show no such jump in Fig. 6, and indeed all three stay well below freezing. However, it must be borne in mind that the IPCC assumptions (used here) have all four RCPs levelling off to a target value by the end of this century; see Fig. 4. That may be overly optimistic.

3.2.2 Ocean and atmosphere heat transport feedback

In Fig. 6, the only forcing parameter that was changing was the CO2 concentration (assumed due to anthropogenic forcing). Now we incorporate changes to ocean and atmosphere heat transport, which may amplify the effects of increasing CO2 alone.

There is evidence that ocean heat transport, FO, into the Arctic is increasing. For example, Koenigk and Brodeau (2014) project ocean heat transport above 70 N to increase from 0.15 PW in 1860 to 0.2 and 0.3 PW in 2100, for RCP 2.6 and 8.5, respectively. At the same time, they find in their model that atmospheric heat transport decreases slightly, from 1.65 PW in 1850 to 1.6 PW (1.5 PW) for RCP 2.6 (RCP 8.5). These authors found that, in a stable climate state that ensures global energy conservation, FO and FA tend to be out of phase; see, for example, the coupled climate model in Koenigk and Brodeau (2014). However, Yang et al. (2016) show that in a more realistic situation when the climate is perturbed by both heat and freshwater fluxes, the changes in FO and FA may be in phase. We assume the latter situation in this paper; see Fig. 7b.

In our model, climate forcings FO and FA are expressed as power per unit area (W m−2) determined as follows (see Table 2). First, the surface area of the Arctic region is estimated. The Arctic is taken to be the surface of the Earth above the 70th parallel; as such the surface area is

(18) Arctic Surface Area = 2 π R 2 ( 1 - cos θ ) ,

where R is the radius of the Earth, 6371 km, and θ is 90 minus the latitude. Hence, the surface area of the Earth above 70 N is approximately 1.538×1013m2. This leads to atmospheric and ocean heat fluxes into the Arctic as summarized in Table 2. Because the change in FA is small relative to the changes in μ and FO, FA is kept constant at an intermediate value of 104 W m−2 in Figs. 5 to 7a.

Table 2Atmosphere and ocean heat fluxes into the Arctic as simulated in Koenigk and Brodeau (2014), using the global coupled climate model EC-Earth.

Download Print Version | Download XLSX

Figure 7Arctic surface temperature change projected to year 2300 (relative to year 2000 temperature of -28.90C), with linearly increasing ocean heat transport FO, interpolating the data in Koenigk and Brodeau (2014); see Table 2. (a) With constant FA=104, the jump in temperature for RCP 8.5 now occurs nearly 40 years earlier than for the case of constant FO shown in Fig. 6. (b) The same as (a) except that in addition to increasing FO the atmosphere heat transport FA also increases, as in Yang et al. (2016), linearly from 104 to 129 W m−2. Now upward transitions occur on both RCP 8.5 and 6.0, as indicated by the arrows.


Figure 7a shows the change in Arctic surface temperature (relative to the year 2000 temperature of -28.90C) for the four RCPs with the ocean heat flux FO increasing linearly on each RCP until the year 2100, as projected in Koenigk and Brodeau (2014), using their data for FO in Table 2 but with constant FA=104. Beyond the year 2100, until 2300, the ocean heat transport, FO, , is held constant at its 2100 value. In this scenario, the onset of the jump (via a fold bifurcation) to a warm equable climate is advanced dramatically. The bifurcation for RCP 8.5 occurs in the year 2118, more than 40 years earlier than was the case with a constant FO in Fig. 6. The temperatures before and after the jump in 2118 between two stable states are 4.6 and +24.7 C, respectively, which is a sudden increase of 29.3C. The other three RCPs remain below the freezing temperature.

Figure 7b shows the same scenario as in 7a, with increasing μ and FO but with the atmospheric heat transport, FA, also increasing linearly from 104 W m−2 in the year 2000 to 129 W m−2 in 2100, being constant thereafter. In this case, the saddle-node bifurcation occurs even earlier for RCP 8.5, and a new bifurcation appears for RCP 6.0. Both of these changes make mitigation more challenging.

3.2.3 Water vapour feedback

Overall, water vapour is known to be the most powerful greenhouse gas in the atmosphere (Dai2006; Pierrehumbert2010; IPCC2013). Warming of the surface causes evaporation of more water vapour, which causes further greenhouse warming and a further rise in surface temperature. Thus, water vapour amplifies the warming due to other causes. This is called water vapour feedback. Empirical studies such as Dai (2006) show that the increase in surface specific humidity H with surface temperature T is close to that predicted by the Clausius–Clapeyron equation as in Eq. (6) (outside of desert regions). The relative humidity, RH or δ, changes little with surface temperature, even as the specific humidity H increases (Serreze et al.2012). For paleoclimates, Jahren and Sternberg (2003) have described an Eocene Arctic rain forest with RH estimated to be δ=0.67. Modern data, from a variety of sources, suggest similar values of Arctic RH. For example, Shupe et al. (2011) report Arctic RH at the surface over 0.7 and atmospheric RH at 2.5 km altitude near 0.6, with relatively small seasonal and spatial variations.

Therefore, in the EBM Eqs. (1) and (2), it is assumed that the greenhouse warming effect of water vapour is determined by the Clausius–Clapeyron relation as in Eq. (6) and the RH δ remains constant. We assumed a value of δ=0.60 for the Arctic atmosphere in the previous section.

The Clausius–Clapeyron equation tells us that below the freezing temperature (τS=1) the concentration of water vapour is very low and therefore it has very little greenhouse effect. However, above freezing, if a source of water is available (e.g. oceans), then the concentration of water vapour and its greenhouse warming effect increase rapidly. This is shown clearly in Fig. 8, where the three curves show different levels of relative humidity, δ, but all assume that CO2 follows RCP 8.5. Here, the reference temperatures in year 2000 are as follows: red curve -28.899C, green curve -29.418C, and blue curve -30.320C. On each of the curves of Fig. 8, the dashed portions with negative slope are unstable, while portions with positive slope are asymptotically stable (in the sense of Liapunov). Bistability (the coexistence of two stable solutions) occurs sooner when water vapour is present. The lower continuous blue curve with δ=0 shows no thawing (τS<1) in this range.

Figure 8Arctic surface temperature change with increasing relative humidity of water vapour, δ=0.0, 0.4, 0.6, and with fixed FO and FA, projected to year 2300. (Temperature change is relative to the year 2000 temperature; see text for details.) On all three curves, CO2 is increasing in time according to RCP 8.5. The red curve is essentially the same as that shown in Fig. 6 with δ=0.6. The blue and green curves have water vapour fixed at δ=0.0 and δ=0.4, respectively. For temperatures significantly below freezing, water vapour makes little contribution to temperature change. However, above freezing (τS>1), the greenhouse warming effect of water vapour is dramatic.


3.2.4 Anthropocene Arctic EBM summary

This EBM for the Anthropocene Arctic predicts that a bifurcation will occur, leading to catastrophic warming of the Arctic if CO2 emissions continue to increase along RCP 8.5 without mitigation. This is true in the model even if ocean and atmosphere heat transport remain unchanged. The amplifying effects of ocean and atmosphere heat transport can make this abrupt climate change become even more dramatic and occur even earlier. Water vapour feedback further intensifies global warming. However, the EBM predicts that RCPs with reduced CO2 emissions (due, for example, to effective mitigation strategies if introduced soon enough) may avert the drastic consequences of such a bifurcation.

Further work on Anthropocene Arctic climate modelling will include the effects of other positive feedback mechanisms, e.g. the greenhouse effects of methane and CO2 that will be released as the permafrost thaws in the Arctic and the Hadley cell feedback that may occur as global circulation patterns change. With such additional amplification in the Arctic taken into account, and no mitigation strategies in place, a saddle-node bifurcation resulting in a transition to a warmer Arctic climate state may occur even earlier than predicted by the present model.

3.3 EBM for the Anthropocene Antarctic

It has long been recognized that the climate of the Southern Hemisphere is generally colder than that of the Northern Hemisphere for a number of reasons (Feulner et al.2013). In particular the Antarctic is colder than the Arctic. The Antarctic climate is affected by the Antarctic Circumpolar Current (ACC), which flows freely west to east around Antarctica in the Southern Ocean, unimpeded by continental barriers. The ACC blocks the poleward heat transport from the warm oceans of the Southern Hemisphere (Hartmann2016). Hence, FO is restricted to be below 20 W m−2 in the Antarctic EBM. Additionally, the cold surface albedo αC=0.8 is greater for the Antarctic than the Arctic, because the snow and ice that covers the continent is more pure than that in the Arctic. Cloud albedo is reduced, with a value of 7 % as opposed to the value of 12.12 % for the Arctic (Pirazzini2004). Atmospheric heat transport is FA=97W m−2, as determined in Zhang and Rossow (1997). Finally, the Antarctic region is much drier than the Arctic; hence, a relative humidity of δ=0.4 is used.

Figure 9Energy balance model for the Antarctic. (a) 3-D equilibrium manifold, (b) Projection of the fold bifurcations. The red dot locates today's climate conditions on the cold surface.


For the Antarctic, Fig. 9a is the equilibrium manifold for the energy balance model parameterized by (FO,μ), and Fig. 9b is the projection of the fold bifurcations onto the parameter plane. In Fig. 9b, the yellow area represents a warm stable climate, the blue area a cold stable climate, and the green area represents the overlap region (between the two fold bifurcations) where bistability exists. It can be seen in Fig. 9 that a bifurcation from a cold state to a warm state in the EBM cannot occur for an ocean heat transport value of less than FO=12W m−2 and a carbon dioxide concentration less than μ=3000ppm.

Figure 10a shows the temperature change, following each of the RCP curves in the Antarctic, extended to the year 2300. The reference temperature is -32.78C for the year 2000, and the value of ocean heat transport into the Antarctic is assumed to be FO=14W m−2, an annual mean for the sea-ice zone (approximately 70S) of the Antarctic, as determined in Wu et al. (2001). This scenario does not exhibit a bifurcation from a cold climate state to a warm state. This suggests that for the Antarctic, a change in μ alone is not sufficient to cause a hysteresis loop to exist, between two coexisting stable states, in the context of modern and near-future carbon dioxide concentrations.

Figure 10Antarctic surface temperature change projected to year 2300, relative to year 2000 temperature of -32.78C; (a) on each of the four RCPs with constant FO=14W m−2 and (b) with increasing ocean heat transport FO; see text for details. The upward arrow indicates the transition to a warmer equilibrium climate state after the saddle-node bifurcation.


Figure 10b presents the Antarctic model for values of FO that increase with time as μ increases. The value of FO is kept constant at 5 W m−2 until the year 2000, after which time it increases linearly up to the year 2100, where it has a value of FO=20W m−2, after which it is held constant again. This increase might represent an increase in sea levels, caused by thawing of the Arctic, loss of the Antarctic ice shelves, and subsequent increase in ocean heat transport (Koenigk and Brodeau2014). The first thing to notice is that a hysteresis loop now exists. With increasing CO2 on RCP 8.5, there is a jump between stable states, from ΔT=+29.1C to ΔT=+56.5C, occurring in year 2225, where ΔT is the change in surface temperature relative to the year 2000 temperature of -32.78C. The above-freezing average temperatures on the upper warm branch of equilibrium states are consistent with estimates of Antarctic temperatures in the Eocene (Passchier et al.2013). Such a transition would imply melting of the Antarctic ice cap and a drastic rise in sea levels around the world. The return bifurcation from warm to cold (with time reversed) is visible in Fig. 10b. The “cold-to-warm” transition occurs at a later time in the Antarctic than in the Arctic and at a higher temperature. The differences between the Antarctic and Arctic bifurcations in the EBM are due to the differences in ocean heat transport and ice albedos. The difference could be larger if other factors are taken into account, e.g. the Antarctic has thicker ice and hence more heat is required to melt enough ice to cause a change in albedo. This could be represented with a larger value of ω in the tanh switch function; then, a greater temperature change is required for the function to switch from a cold albedo value to a warm one.

3.4 EBM for the Anthropocene tropics

Next, the EBM is adapted to model the climate of the tropics by choosing parameter values that are annual mean, zonally averaged values at the Equator. This gives insolation Q=418.8W m−2 and relative humidity δ=0.8. Heat transports FA=-38W m−2 and FO=-39W m−2 are both negative, because heat is transported away from the Equator towards the poles (Hartmann2016). The shortwave cloud cooling ξR (the albedo of the clouds) is also greater in the tropics, and the surface has a lower albedo. The value ξR=22.35 % is determined in the appendix from the global energy budgets of Trenberth et al. (2009) and Wild et al. (2013).

As the tropics have annual average temperatures well above the freezing point of water, ice–albedo feedback is absent in the tropics and a bifurcation from a cold stable state to a warm state can not occur under Anthropocene conditions. However, if forced to low FO, FA values and very low carbon dioxide levels, the climate state known as “snowball Earth” (Kaper and Engler2013; Pierrehumbert2010) is a possibility. That scenario is not explored in this paper, as it is not relevant to the Anthropocene.

The large relative humidity of δ=0.8 in the tropics serves to mitigate the radiative forcing of increasing CO2. Water vapour is a much more effective greenhouse gas than carbon dioxide (Pierrehumbert2010). The atmosphere of the tropics contains more water vapour (the product of greater relative humidity and warmer temperatures). The total atmospheric longwave absorption η is given in Eq. (9). In the EBM of the tropics, the water vapour content is so high that ηW (and thus total absorptivity η) is almost “maxed out” at η≈1. Hence the total absorptivity is dominated by the contribution due to water vapour, and an increase in CO2 concentration has little additional greenhouse warming effect.

Figure 11 reveals relatively low increases in temperature for the tropics compared to the poles, as CO2 concentration increases along the RCPs. Both the absence of a bifurcation and the mitigation due to a large existing water vapour greenhouse forcing taken together cause the temperature increase relative to the year 2000 to be less than 2C, in all four RCP scenarios.

Figure 11Tropical surface temperature changes for the four RCP scenarios, forecast to year 2300 relative to year 2000 (25.19 C), with constant FO and FA. The temperature change is relatively small in the tropics.


3.5 EBM for globally averaged temperature

Changes in globally averaged temperature can be modelled more easily than changes in regional temperatures due to the fact that, in a globally averaged equilibrium model, overall net horizontal transport of energy, by the oceans and the atmosphere, are both zero. Thus, the two-layer EBM Eqs. (4) and (5), globally averaged with FO=0 and FA=0, are simplified as follows.


Here α is as defined in Eq. (8), η is as in Eq. (9), and fC is as in Eq. (10). Parameters ξA and ξR are as in Table 1 and Sect. 2.1.

Figure 12 shows the change in globally averaged equilibrium surface temperature relative to the year 2000 global average (τS=1.064, 17.59 C), as determined by the EBM Eq. (19) to the year 2300. It is assumed that CO2 evolves with time along each of the four RCPs defined in IPCC (2013) and displayed in Fig. 4. The other parameters, assumed constant, are as follows. The global relative humidity δ is fixed at a value of 0.74. This is determined from Dai (2006), where it lies at the lower end of a range of averages. Surface albedo is highly variable regionally, so a global average was calculated from Wild et al. (2013), much like the atmospheric shortwave absorption and the cloud albedo. From Fig. 1 of Wild et al. (2013), of the global average solar radiation of 185 W m−2 that reaches the surface, a portion of 24 W m−2 is reflected. Thus, the global average surface albedo is 24185=0.13=13%. The values for cloud albedo and atmospheric shortwave radiation are calculated as follows. The global average incident solar radiation Q at the top of the atmosphere is 340 W m−2, of which 100-24340=0.2235=22.35% is reflected by clouds and 79340=0.2324=23.24% is absorbed by the atmosphere. The planetary boundary layer (PBL) altitude is 600 m (Ganeshan and Wu2016). Finally, for the purposes of sensible and latent heat transport (see appendix), the wind speed U is 5 m s−1 (Nugent et al.2014) and the drag coefficient is CD=1.5×10-3.

Figure 12Change in globally averaged surface temperature, relative to year 2000 global average temperature of 17.59 C, calculated to the year 2300 for the EBM on each of the four RCPs.


3.6 Equilibrium climate sensitivity

Equilibrium climate sensitivity (ECS) is a useful and widely adopted tool used to estimate the effects of anthropogenic forcing in a given climate model. The ECS of a model is defined as the change in the globally averaged surface temperature, after equilibrium is obtained, in response to a doubling of atmospheric CO2 levels (IPCC2013; Knutti et al.2017; Priostosescu and Huybers2017). The starting carbon dioxide concentration is that of the preindustrial climate, taken to be μ=270 ppm. The doubled value is then μ=540 ppm. Since the Earth has not yet experienced a doubling of CO2 concentration since the industrial revolution, these numbers cannot be verified. Calculation of the global ECS for the EBM of this paper facilitates comparisons with other climate models as reported in IPCC (2013).

3.6.1 ECS for the globally averaged EBM

Table 3 gives both the nondimensional τS and the degree Celsius temperature values for the μ=270 climate, the μ=540 climate, and the temperature difference. This difference is the ECS of the global EBM of this paper.

Table 3ECS for the globally averaged EBM.

Download Print Version | Download XLSX

For the models used in IPCC (2013), ECS values lie within a likely range of 1.5 to 4.5 C. Values of less than 1 C are deemed to be extremely unlikely, and values greater than 6 C are very unlikely (IPCC2013). The value of 4.343 C calculated for this global EBM lies just inside the likely range, at the high end. Recent work gives evidence that statistical climate models based on historical data tend to lie on the lower end of likely ECS values, with a range of 1.5 to 3 C, whereas nonlinear GCMs tend to have larger ECS values (Priostosescu and Huybers2017). Therefore, as the global EBM presented in this paper is nonlinear and is based on physical rather than statistical modelling, it may be expected to fall on the side of larger ECS values.

3.6.2 Regional ECS values

Local ECS values may be determined for each of the three regional models, the Arctic, Antarctic, and the tropics, as defined in Sects. 3.2 to 3.4. These values are given in Table 4. In all cases, FO values are kept constant at their minimal values: 10 W m−2 for the Arctic, 14 W m−2 for the Antarctic, and −39W m−2, for the tropics. The regional ECS values are high, 7.95 and 7.54 C, respectively, for the Arctic and Antarctic, and low, 1.27 C for the tropics. Although the Earth has not yet experienced a doubling of CO2 concentrations since the industrial revolution, these ECS values are consistent with observations to date.

Table 4ECS values for each of the three regional EBMs. The ECS is much greater for the poles than for the tropics, in agreement with observations.

Download Print Version | Download XLSX

4 Conclusions and future work

The analysis of this paper shows that a cusp bifurcation can occur in an energy balance model (EBM), which could lead to hysteresis behaviour and an abrupt warming of the Anthropocene climate, to a climate state that is like nothing that has existed on Earth since the Pliocene. The model has been constructed from fundamental nonlinear processes of atmospheric physics. This bifurcation is most likely to occur in the Arctic climate. It would lead to catastrophic warming, if increases in atmospheric CO2 continue on their current pathway. However, if the increase in atmospheric CO2 is mitigated sufficiently, this bifurcation can still be avoided. Climate changes in the Arctic, Antarctic, and tropics are compared. The globally averaged equilibrium climate sensitivity (ECS) of the EBM is 4.34 C, which is at the high end of the likely range reported in IPCC (2013).

Future work will strengthen the conclusions of this paper by extending this simple EBM to more comprehensive models, which still allow rigorous bifurcation analysis to be performed. In the first generalization, the two-layer model will be replaced with a column model, with the atmosphere extending continuously from the surface to the tropopause. The ICAO Standard Atmosphere assumption will be replaced with a Schwarzschild radiation model of the atmosphere (Pierrehumbert2010), which will determine the lapse rate from the solution of a two-point boundary value problem. This Schwarzschild column model will be used to study, in addition to the positive feedback processes of this paper, the amplifying effects of methane from permafrost feedback in the Arctic.

The next generalization will combine the ideas of these energy balance models with a 3-D Navier–Stokes–Boussinesq partial differential equation (PDE) model, representing the convectively driven atmosphere as a fluid in a rotating spherical shell, presented in Lewis and Langford (2008) and Langford and Lewis (2009). In that model, surface temperatures were prescribed as boundary conditions; however, with the assumption of an energy balance constraint, the meridional surface temperature gradient will be determined implicitly. Meridional heat transport also will be determined in this model and may be compared with the poleward heat transports in the EBM of this paper. The code to solve this PDE model for a cusp bifurcation has been written. Later, guided by these results, the climate bifurcations found in these analytical models will be sought in an open-source general circulation model. This hierarchy of models is expected to add credibility to the prediction, presented here, that the Earth's climate system is capable of exhibiting dramatic topological changes (bifurcations) in the Anthropocene, leading to a climate state that resembles the pre-Pliocene climate of Earth.

Appendix A: Determination of parameters in the Anthropocene EBM

The determination of the physical parameters appearing in the EBM (Eqs. 4 and 5) is summarized here. In most cases, these parameters were determined in earlier papers of the authors dealing with paleoclimates (Dortmans et al.2019; Kypke and Langford2020). The focal point is the scaled, two-dimensional Eqs. (4) and (5) in Sect. 2.1, with α given by Eq. (8), η given by Eq. (9), and fC given by Eq. (10).

First consider the incoming solar radiation Q. A fraction ξA is absorbed by the atmosphere and a fraction ξR is reflected by the clouds, as seen in Fig. 1 and Eqs. (4) and (5). These fractions were determined in Kypke (2019) and appendix A1 of Dortmans et al. (2019), using data for the global average energy budget described in Trenberth et al. (2009), Wild et al. (2013); and Kim and Ramanathan (2012). The globally averaged values were ξA=0.2324 and ξR=0.2235 as listed in Table 1. For the polar regions, the albedo of clouds is less than elsewhere. Using data collected in the Surface Heat Budget of the Arctic Ocean (SHEBA) program (Intrieri et al,2002; Shupe and Intrieri2003), the revised polar value of ξR=0.1212 was determined in Kypke (2019).

Clouds in the atmosphere have a second main effect, the absorption of a fraction, ηCl, of the longwave radiation outgoing from the surface of the Earth (Hartmann2016). This effect warms the Earth's atmosphere. Through data in the SHEBA program, the longwave cloud forcing in the Arctic was found to be 51 W m−2 in the paper of Shupe and Intrieri (2003). Using this SHEBA data, Kypke (2019) determined the fraction ηCl=0.255, as used in Eq. (9); see Table 1.

The absorption coefficients kC for carbon dioxide and kW for water vapour in the atmosphere (see Table 1) were calculated using an empirical approach, based on the modern-day global energy budget (Trenberth et al.2009; Wild et al.2013). Figure 1 in Trenberth et al. (2009) provides the global mean surface radiation as 396 W m−2, along with an atmospheric window of 40 W m−2. This atmospheric window, 403960.1, is then equal to 1−η. Schmidt et al. (2010) provide percentage contributions of carbon dioxide and water vapour and clouds in an all-sky scenario, based on simulations using modern climate conditions from the year 1980. The calculated values for ηC and ηW are then used to calculate the corresponding optical depths, λC, and λW for the case of the modern atmosphere, and these are used to solve for kC and kW, which then appear in the GC and GW2 terms, respectively, in Table 1 and Eq. (7).

The vertical transport of sensible and latent heat is a difficult and complicated process to model, so many approximations are made to keep it within the scope of this work. For more details, see Kypke (2019). The heat transports are modelled via bulk aerodynamic exchange formulae describing fluxes between the surface and the atmosphere (Pierrehumbert2010; Hartmann2016). For the sensible heat flux, cp is the specific heat of the air being heated, and T is its temperature. The bulk aerodynamic formula for sensible heat (SH) is

(A1) SH = c p ρ C DS U ( T S - T A ) ,

where CDS is the drag coefficient for temperature and U is the mean horizontal wind velocity. The density of the atmosphere ρ is determined as a function of both surface temperature TS and altitude Z using the barometric formula, and a constant lapse rate Γ (ICAO1993) is used to determine the temperature gradient.

In the case of latent heat (LH), the moisture content is represented by Lvr, where Lv is the latent heat of vaporization of water and r is the mass mixing ratio of condensable air to dry air (Pierrehumbert2010). Then

(A2) LH = L v ρ C DL U ( r S - r A ) ,

where CDL is the drag coefficient for moisture. In the following, for simplicity, we set


The mass mixing ratio is equal to the saturation mixing ratio times the relative humidity. The saturation mixing ratio depends on the saturation vapour pressure, which is a function of temperature as given by the Clausius–Clapeyron Eq. (6). The sensible and latent heat transports are combined into a single term, FC, which replaces the FC term that was introduced in Dortmans et al. (2019). This term is defined here as a function of surface temperature TS.


Here, P0 is the atmospheric pressure at the surface (Z=0) and RA is the ideal gas constant specific to dry air. This equation is scaled by 1σTR4 to nondimensionalize it, creating fC=FCσTR4 in Eq. (10), where TR=273.15K is the reference temperature. As this fC represents energy moving from the surface to the atmosphere, it is subtracted from the surface Eq. (4) and added to the atmosphere Eq. (5). A different model of FC, used by the authors in Dortmans et al. (2019), was a simple functional form calibrated to empirical data. The result was a relationship between FC and TS that is quantitatively very similar to that given by Eq. (A3). In Kypke and Langford (2020), FC was set equal to zero for the paleoclimate Arctic and Antarctic models for simplicity, since both SH and LH are very small for temperatures that are below freezing.

Code and data availability

The MATLAB files, used to solve the model equations and to produce the figures in this article, have been posted online at (Kypke et al.2020).

Author contributions

The original conception of this research project was due to WFL. Most of the mathematical analysis and computations were performed by KLK, starting with his MSc Thesis (Kypke2019). ARW co-supervised the thesis and provided valuable contributions and insights throughout the course of the work. All authors contributed to the writing of the article.

Competing interests

The authors declare that they have no conflict of interest.


The authors thank Marcus Garvie for his support of this work and thank the two NPG reviewers for suggestions and constructive criticisms that have improved the quality of this paper.

Financial support

This research has been supported by the Natural Sciences and Engineering Research Council of Canada (grant no. RGPIN-2015-05090).

Review statement

This paper was edited by Stefano Pierini and reviewed by Michel Crucifix and Peter Ditlevsen.


Alley, R. B., Marotzke, J., Nordhaus, W. D., Overpeck, J. T., Peteet, D. M., Pielke Jr., A., Pierrehumbert, R. T., Rhines, P. B., Stocker, T. F., Talley, L. D., and Wallace, J. M.: Abrupt climate change, Science, 299, 2005–2010, 2003. a

Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system, Philos. T. Roy. Soc. A, 370, 1166–1184, 2012. a

Barnosky, A.D., Hadly, E. A., Bascompte, J., Berlow, E. L., Brown, J. H., Fortelius, M., Getz, W. M., Harte, J., Hastings, A., Marquet, P. A., Martinez, N. D., Mooers, A., Roopnarine, P., Vermeij, G., Williams, J. W., Gillespie, R., Kitzes, J., Marshall, C., Matzke, N., Mindell, D. P., Revilla, E., and Smith, A. B.: Approaching a state shift in Earth's biosphere, Nature, 486, 52–58, 2012. a

Barron, E. J., Thompson, S. L., and Schneider, S. H.: An ice-free Cretaceous? Results from climate model simulations, Science, 212, 10–13, 1981. 

Bathiany, S., Dijkstra, H., Crucifix, M., Dakos, V., Brovkin, V., Williamson, M. S., Lenton, T. M., and Scheffer, M.: Beyond bifurcation: using complex models to understand and predict abrupt climate change, Dyn. Stat. Clim. Syst., 1, 1–31,, 2016. a

Björk, G., Stranne, C., and Borenäs, K.: The sensitivity of the Arctic Ocean sea ice thickness and its dependence on the surface albedo parameterization, J. Climate, 26, 1355–1370, 2013. a

Brovkin, V., Claussen, M., Petoukhov, V., and Ganopolski, A.: On the stability of the atmosphere-vegetation system in the Sahara/Sahel region, J. Geophys. Res., 103, 31613–31624, 1998. a

Budyko, M. I.: The effects of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619, 1968. a, b

Dai, A.: Recent climatology, variability and trends in global surface humidity, J. Climate, 19, 3589–3606, 2006. a, b, c

Dijkstra, H. A.: Nonlinear Climate Dynamics, Cambridge University Press, New York, 2013. a

Dijkstra, H. A.: Numerical bifurcation methods applied to climate models: analysis beyond simulation, Nonlin. Process. Geophys., 26, 359–369, 2019. a

Ditlevsen, P. D. and Johnsen, S. J.: Tipping points: Early warning and wishful thinking, Geophys. Res. Lett., 37, L19703,, 2010. a

Dortmans, B., Langford, W. F., and Willms, A. R.: A conceptual model for the Pliocene paradox, in: Recent Advances in Mathematical and Statistical Methods, Springer Proceedings in Mathematics and Statistics vol. 259, edited by: Kilgour, D. M., Kunze, H., Makarov, R., Melnik, R., Wang, X., IV AMMCS International Conference, Waterloo, Canada, 20–25 August 2017, Springer Nature Switzerland AG, Gewerbestrasse 11, 6330 Cham, Switzerland, pp. 339–349, 2019. a

Dortmans, B., Langford, W. F., and Willms, A. R.: An energy balance model for paleoclimate transitions, Clim. Past, 15, 493–520,, 2019. a, b, c, d, e, f, g, h, i, j, k

Drijfhout, S., Bathiany, S., Beaulieu, C., Brovkin, V., Claussend, M., Huntingford, C., Scheffer, M., Giovanni Sgubin, G., and Swingedouw, D.: Catalogue of abrupt shifts in Intergovernmental Panel on Climate Change climate models, P. Natl. Acad. Sci. USA, 112, E5777–E5786,, 2015. a

Feulner, G., Rahmstorf, S., Levermann, A., and Volkwardt, S.: On the origin of the surface air temperature difference between the hemispheres in Earth's present-day climate, J. Climate, 26, 7136–7150, 2013. a

Ganeshan, M. and Wu, D. L.: The open-ocean sensible heat flux and its significance for Arctic boundary layer mixing during early fall, Atmos. Chem. Phys., 16, 13173–13184,, 2016. a

Ghil, M.: Hilbert problems for the geosciences in the 21st century, Nonlin. Process. Geophys., 8, 211–222, 2001. a

Ghil, M. and Lucarini, V.: The Physics of Climate Variability and Climate Change, arXiv 1910.00583v2, 26 February 2020. a

Hartmann, D. L.: Global Physical Climatology, 2nd edition, Elsevier, Amsterdam, 2016. a, b, c, d, e

International Institute for Applied Systems Analysis (IIASA): The Representative Concentration Pathways data for Fig. 4, available at:, last access: 2019. a

International Civil Aviation Organization: Manual of the ICAO Standard Atmosphere: extended to 80 kilometres, Third Edition, ICAO, 1000 Sherbrooke Street West Suite 400, Montreal, Quebec, Canada H3A 2R2, 1993. a, b, c

Intrieri, J. M., Fairall, C. W., Shupe, M. D., Persson, P. O. G., Andreas, E. L., Guest, P. S., and R. E. Moritz, R. E.: An annual cycle of Arctic surface cloud forcing at SHEBA, J. Geophys. Res., 107, 13–14, 2002. a

Intergovernmental Panel on Climate Change: Climate Change 2013: The Physical Science Basis, in: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G. K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. M., Cambridge Univ. Press, Cambridge, UK and New York, NY, USA, available at: (last access: 11 March 2019), 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Jahren, A. H. and Sternberg, L. S. L.: Humidity estimate for the middle Eocene Arctic rain forest, Geology, 31, 463–466, 2003. a

Kaper, H. and Engler, H.: Mathematics and Climate, Society for Industrial and Applied Mathematics, Philadelphia, 2013. a, b

Kim, D. and Ramanathan, V.: Improved estimates and understanding of global albedo and atmospheric solar absorption, Geophys. Res. Lett., 39, L24704,, 2012. a

Knutti, R., Rugenstein, M. A. A., and Hegerl, G. C.: Beyond equilibrium climate sensitivity, Nat. Geosci., 10, 727–744, 2017. a

Koenigk, T. and Brodeau, L.: Ocean heat transport into the Arctic in the twentieth and twenty-first century in EC-Earth, Clim. Dynam., 42, 3101–3120, 2014. a, b, c, d, e, f

Kuznetsov, Y. A.: Elements of Applied Bifurcation Theory, 3rd edition, Applied Mathematical Sciences vol. 112, Springer-Verlag, New York, 2004. a

Kypke, K. L.: Topological Climate Change and Anthropogenic Forcing, M.Sc. Thesis, Department of Mathematics and Statistics, University of Guelph, Guelph, Ontario, Canada, 137 pp., 2019. a, b, c, d, e, f, g, h, i, j, k

Kypke, K. L. and Langford, W. F.: Topological climate change, Int. J. Bifurc. Chaos, 30, 2030005,, 2020. a, b, c, d, e, f, g, h, i, j, k, l

Kypke, K. L., Langford, W. F., and Willms, A. R.: MATLAB files for figures for “Anthropocene Climate Bifurcation” in Nonlinear Processes in Geophysics,, 2020. a

Langford, W. F. and Lewis, G. M.: Poleward expansion of Hadley cells, Can. Appl. Math. Quart., 17, 105–119, 2009. a

Lenton, T. M.: Arctic climate tipping points, AMBIO, 41, 10–22, 2012. a

Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping elements in the Earth's climate system, P. Natl. Acad. Sci. USA, 105, 1786–1793, 2008. a

Lewis, G. M. and Langford, W. F.: Hysteresis in a rotating differentially heated spherical shell of Boussinesq fluid, SIAM J. Appl. Dyn. Syst., 7, 1421–1444, 2008. a

North, G. R. and Kim, K.-Y.: Energy Balance Climate Models, Wiley–VCH, Weinheim, Germany, 2017. a, b

North, G. R., Cahalan, R. F., and Coakley, J. A.: Energy balance climate models, Rev. Geophys. Space Phys., 19, 91–121, 1981. a, b

Nugent, A. D., Smith, R. B., and Minder, J. R.: Wind speed control of tropical orographic convection, J. Atmos. Sci., 71, 2695–2712, 2014. a

Passchier, S., Bohaty, S. M., Jiménez-Espejo, F., Pross, J., Röhl, U., van de Flierdt, T., Escutia, C., and Brinkhuis, H.: Early Eocene to middle Miocene cooling and aridification of East Antarctica, Geochem. Geophy. Geosy., 14, 1399–1410, 2013. a

Payne, A. E., Jansen, M. F., and Cronin, T. W.: Conceptual model analysis of the influence of temperature feedbacks on polar amplification, Geophys. Res. Lett., 42, 9561–9570, 2015. a, b

Pierrehumbert, R. T.: Principles of Planetary Climate, Cambridge University Press, Cambridge, UK, 2010. a, b, c, d, e, f, g

Pirazzini, R.: Surface albedo measurements over antarctic sites in summer, J. Geophys. Res., 109, D20118,, 2004. a

Pistone, K., Eisenman, I., and Ramanathan, V.: Observational determination of albedo decrease caused by vanishing Arctic sea ice, P. Natl. Acad. Sci. USA, 111, 3322–3326, 2014. a

Proistosescu, C. and Huybers, P. J.: Slow climate mode reconciles historical and model-based estimates of climate sensitivity, Sci. Adv., 3, e1602821,, 2017. a, b, c

Sagan, C. and Mullen, G.: Earth and Mars: Evolution of atmospheres and surface temperatures, Science, 177, 52–56, 1972. a

Schmidt, G. A., Ruedy, R. A., Miller, R. L., and Lacis, A. A.: Attribution of the present day total greenhouse effect, J. Geophys. Res., 115, D20106,, 2010. a

Seager, R. and Battisti, D. S.: Challenges to our understanding of the general circulation: abrupt climate change, in: The Global Circulation of the Atmosphere: Phenomena, Theory, Challenges, edited by: Schneider, T. and Sobel, A. S., Princeton University Press, Princeton, NJ, pp. 331–371, 2007. a

Sellers, W. D.: A global climate model based on the energy balance of the Earth–atmosphere system, J. Appl. Meteorol., 8, 392–400, 1969. a, b

Serreze, M. C., Barrett, A. P., and Stroeve, J.: Recent changes in tropospheric water vapor over the Arctic as assessed from radiosondes and atmospheric reanalyses, J. Geophys. Res., 117, D10104,, 2012. a

Shupe, M. D. and Intrieri, J. M.: Cloud radiative forcing of the Arctic surface: The influence of cloud properties, surface albedo, and solar zenith angle, J. Climate, 17, 616–628, 2003. a, b

Shupe, M. D., Walden, Von P., Eloranta, E., Uttal, T., Campbell, J. R., Starkweather, S. M., and Shiobara, M.: Clouds at Arctic atmospheric observatories. Part I: occurrence and macrophysical properties, J. Appl. Meteorol. Clim., 50, 626–644, 2011. a

Steffen, W., Rockström, J., Richardson, K., Lenton, T. M., Folke, C., Liverman, D., Summerhayes, C. P., Barnosky, A. D., Cornell, S. E., Crucifix, M., Donges, J. F., Fetzer, I., Lade, S. J., Scheffer, M., Winkelmann, R., and Schellnhuber, H. J.: Trajectories of the Earth system in the Anthropocene, P. Natl. Acad. Sci. USA, 115, 8252–8259,, 2018. a, b

Stranne, C., Jakobsson, M., and Björk, G.: Arctic Ocean perennial sea ice breakdown during the Early Holocene Insolation Maximum, Quart. Sci. Rev., 92, 123–132, 2014. a

Thorndike, A.: Multiple equilibria in a minimal climate model, Cold Reg. Sci. Technol., 76, 3–7,, 2012. a

Trenberth, K. E., Fasullo, J. T., and Kiehl, J.: Earth's global energy budget, B. Am. Meteorol. Soc., 90, 311–323, 2009. a, b, c, d, e, f

Valdes, P.: Built for stability, Nat. Geosci., 4, 414–416, 2011. a

van Vuuren, D. P., Edmonds, J., Kainuma,M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J.-F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J. and Rose, S. K.: The representative concentration pathways: an overview, Clim. Change, 109, 5–31, 2011. a

Wallace-Wells, D.: The Uninhabitable Earth: Life After Warming, Penguin Random House LLC, New York, 2019. a

Wild, M., Folini, D., Schär, C., Loeb, N., Dutton, E. G. and König-Langlo, G.: The global energy balance from a surface perspective, Clim. Dynam., 40, 3107–3134, 2013. a, b, c, d, e, f, g

Willard, D. A., Donders, T. H., Reichgelt, T., Greenwood, D. R., Sangiorgi, F., Peterse, F., Nierop, K. G. J., Frieling, J., Schouten, S., and Sluijs, A.: Arctic vegetation, temperature, and hydrology during Early Eocene transient global warming events, Global Planet. Change, 178, 139–152, 2019.  a

Wu, X., Budd, W. F., Worby, A. P., and Allison, I.: Sensitivity of the Antarctic sea-ice distribution to oceanic heat flux in a coupled atmosphere-sea-ice model, Ann. Glaciol., 33, 577–584, 2001. a

Yang, H., Zhao, Y.-Y., and Liu, Z.-Y.: Understanding Bjerknes compensation in atmosphere and ocean transports using a coupled box model, J. Climate, 29, 2145–2160, 2016. a, b

Zeebe, R. E.: Where are you heading Earth?, Nat. Geosci. 4, 416–417, 2011. a

Zhang, Y.-C. and Rossow, W. B.: Estimating meridional energy transports by the atmospheric and oceanic circulations using boundary fluxes, J. Climate, 10, 2358–2373, 1997. a

Short summary
The climate of Earth is governed by nonlinear processes of geophysics. This paper presents energy balance models (EBMs) embracing these nonlinear processes which lead to positive feedback, amplifying the effects of anthropogenic forcing and leading to bifurcations. We define bifurcation as a change in the topological equivalence class of the system. We initiate a bifurcation analysis of EBMs of Anthropocene climate, which shows that a catastrophic climate change may occur in the next century.