Anthropocene climate bifurcation

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 saddlenode 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-Pliocenelike 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.


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;Ghil, 2001;Alley et al., 2003;Seager and Battisti, 2007;Lenton et al., 2008;Ditlevsen and Johnsen, 2010;Lenton, 2012;Ashwin et al., 2012;Barnosky et al., 2012;Drijfhout et al., 2015;Bathiany et al., 2016;North and Kim, 2017;Steffen et al., 2018;Dijkstra, 2019;Wallace-Wells, 2019). 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 (Valdes, 2011), and that the study of paleoclimates may be a better guide to how abrupt changes may occur (Zeebe, 2011). 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 (Budyko, 1968;Sellers, 1969;Sagan and Mullen, 1972;North et al., 1981;Thorndike, 2012;Kaper and Engler, 2013;Dijkstra, 2013;Payne et al., 2015;Hartmann, 2016;North and Kim, 2017;Ghil and Lucarini, 2020). 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 CO 2 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.
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.

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  Table 1. In the EBM of Fig. 1, the surface receives shortwave radiant energy F S = (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 αF S is reflected back to space, where α is the surface albedo, and the remainder, (1 − α)F S , is absorbed by the surface. The surface re-emits longwave radiant energy of intensity I S = σ T 4 S (Stefan-Boltzmann law) upward into the atmosphere. The atmosphere contains greenhouse gases that absorb a fraction η of the radiant energy I S from the surface, and the remainder, (1 − η)I S , escapes to space. The atmosphere re-emits radiant energy of total intensity I A . Of this radiation I A , a fraction βI A is directed downward to the surface, and the remaining fraction (1 − β)I A goes upward and escapes to space. Balancing the energy flows represented in Fig. 1 leads directly to the following dynamical system: 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 The surface has albedo α, which means that αF S is reflected back to space and the remaining energy (1 − α)F S is absorbed by the surface. The surface emits longwave radiation of intensity I S , of which a fraction ηI S is absorbed by greenhouse gases in the atmosphere and the remainder (1 − η)I S escapes to space. The atmosphere emits longwave radiation I A , of which a fraction βI A goes downward to the surface and the remaining fraction (1 − β)I A escapes to space. The three forcing terms (F A , F O , F C ) 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. where Eq. (1) represents surface energy balance and Eq.
(2) represents atmosphere energy balance. Here c S and c A 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: F A is atmospheric heat transport, F O is ocean heat transport (both horizontally), and F C is conductive/convective heat transport, vertically from the surface to the atmosphere. In Eqs. (1) and (2) there is an asymmetry, in that temperature T S is used to represent the state of the surface in Eq. (1), but radiant energy intensity I A 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 (I S = σ T 4 S and I A = σ T 4 A , where = 0.9 is the emissivity of dry air). The use of T S 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 I A instead of T A 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 T A for the atmosphere. However, we can define a single value of radiant energy intensity I A , 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 I A instead of T A 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 (T R = 273.15 K). 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 CO 2 and H 2 O. 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 (ICAO, 1993), 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 CO 2 is independent of temperature, but the concentration of H 2 O decreases with altitude as the temperature decreases, according to the Clausius-Clapeyron relation (Pierrehumbert, where µ is the concentration of CO 2 in the atmosphere, measured in molar parts per million; δ is the relative humidity of water vapour (0 ≤ δ ≤ 1); γ = /T R is the nondimensionalized lapse rate (ICAO, 1993); 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 k C and k W are absorption coefficients for CO 2 and H 2 O, 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: 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 ω = /T R = 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).

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 landbased 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 F C has been modified to take into account both sensible and latent heat transport (Kypke, 2019). See the appendix, where the following formula is obtained (here in nondimensional form).
where H 1 and H 2 are nondimensional constants: In Eqs. (4) and (5), at equilibrium (i.e. d· dŝ = 0), the state variable i A is easily eliminated, leaving a single equation with a single state variable τ S ,

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 Langford, 2020), 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 : R 2 × R 4 → R 2 . Then an equilibrium point (τ S ,ī A ) of Eqs. (4) and (5), at which where ρ represents four physical parameters that may be varied in the model, 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 (F O , µ) is of primary interest. Equilibrium points (τ S ,ī A ) satisfying Eq. (12) have been computed in Kypke (2019). Having computed the equilibrium point (τ S ,ī 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, e 1 , e 2 , 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 e 1 , e 2 in Eq. (17), below. Then in (u, v) coordinates, the 2-D system Eq. (14) becomes where and Recall that the parameters µ and δ, representing CO 2 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.

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 F O , 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 CO 2 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.

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.5 W 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.0 W m −2 and RCP 4.5 corresponds to +4.5 W 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.6 W 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.
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 (IPCC, 2013). 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 CO 2 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 CO 2 pathway.   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 F A = 104 W 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.

EBM for the Anthropocene Arctic
In Fig. 5a 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 (F O , µ) 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 (F O , µ) 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.

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 F O . 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 F O = 10 W m −2 , F A = 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.90 • C (τ S = 0.8942). 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 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 CO 2 concentrations (µ) by year, up to year 2100. Then both use these CO 2 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 c S and c A 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.5 • C to the new value of +22.5 • C. 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.

Ocean and atmosphere heat transport feedback
In Fig. 6, the only forcing parameter that was changing was the CO 2 concentration (assumed due to anthropogenic forcing). Now we incorporate changes to ocean and atmosphere heat transport, which may amplify the effects of increasing CO 2 alone.
There is evidence that ocean heat transport, F O , 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 Table 2. Atmosphere and ocean heat fluxes into the Arctic as simulated in Koenigk and Brodeau (2014) 8.5). These authors found that, in a stable climate state that ensures global energy conservation, F O and F A 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 F O and F A may be in phase. We assume the latter situation in this paper; see Fig. 7b. In our model, climate forcings F O and F A 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 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 × 10 13 m 2 . This leads to atmospheric and ocean heat fluxes into the Arctic as summarized in Table 2. Because the change in F A is small relative to the changes in µ and F O , F A is kept constant at an intermediate value of 104 W m −2 in Figs. 5 to 7a. Figure 7a shows the change in Arctic surface temperature (relative to the year 2000 temperature of −28.90 • C) for the four RCPs with the ocean heat flux F O increasing linearly on each RCP until the year 2100, as projected in Koenigk and Brodeau (2014), using their data for F O in Table 2 but with constant F A = 104. Beyond the year 2100, until 2300, the ocean heat transport, F O , , 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 F O 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.3 • C. The other three RCPs remain below the freezing temperature. Figure 7b shows the same scenario as in 7a, with increasing µ and F O but with the atmospheric heat transport, F A , 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.  Koenigk and Brodeau (2014); see Table 2. (a) With constant F A = 104, the jump in temperature for RCP 8.5 now occurs nearly 40 years earlier than for the case of constant F O shown in Fig. 6. (b) The same as (a) except that in addition to increasing F O the atmosphere heat transport F A 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.

Water vapour feedback
Overall, water vapour is known to be the most powerful greenhouse gas in the atmosphere (Dai, 2006;Pierrehumbert, 2010;IPCC, 2013). 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 CO 2 follows RCP 8.5. Here, the reference temperatures in year 2000 are as follows: red curve −28.899 • C, green curve −29.418 • C, and blue curve −30.320 • C. 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.

Anthropocene Arctic EBM summary
This EBM for the Anthropocene Arctic predicts that a bifurcation will occur, leading to catastrophic warming of the Arctic if CO 2 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 CO 2 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 CO 2 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. On all three curves, CO 2 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.

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 (Hartmann, 2016). Hence, F O 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 (Pirazzini, 2004). Atmospheric heat transport is F A = 97 W 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.
For the Antarctic, Fig. 9a is the equilibrium manifold for the energy balance model parameterized by (F O , µ), 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 F O = 12 W m −2 and a carbon dioxide concentration less than µ = 3000 ppm. 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.78 • C for the year 2000, and the value of ocean heat transport into the Antarctic is assumed to be F O = 14 W m −2 , an annual mean for the sea-ice zone (approximately 70 • S) 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.  increases linearly up to the year 2100, where it has a value of F O = 20 W 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 Brodeau, 2014). The first thing to notice is that a hysteresis loop now exists. With increasing CO 2 on RCP 8.5, there is a jump between stable states, from T = +29.1 • C to T = +56.5 • C, occurring in year 2225, where T is the change in surface temperature relative to the year 2000 temperature of −32.78 • C. 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-towarm" 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.

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.8 W m −2 and relative humidity δ = 0.8. Heat transports F A = −38 W m −2 and F O = −39 W m −2 are both negative, because heat is transported away from the Equator towards the poles (Hartmann, 2016). 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 F O , F A values and very low carbon dioxide levels, the climate state known as "snowball Earth" (Kaper and Engler, 2013;Pierrehumbert, 2010) 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 CO 2 . Water vapour is a much more effective greenhouse gas than carbon dioxide (Pierrehumbert, 2010). 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 CO 2 concentration has little additional greenhouse warming effect. Figure 11 reveals relatively low increases in temperature for the tropics compared to the poles, as CO 2 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 2 • C, in all four RCP scenarios.

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 F O = 0 and F A = 0, are simplified as follows.
Here α is as defined in Eq. (8), η is as in Eq. (9), and f C 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 CO 2 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 24 185 = 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−24 340 = 0.2235 = 22.35 % is reflected by clouds and 79 340 = 0.2324 = 23.24 % is absorbed by the atmosphere. The planetary boundary layer (PBL) altitude is 600 m (Ganeshan and Wu, 2016). 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 C D = 1.5 × 10 −3 .

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 CO 2 levels (IPCC, 2013;Knutti et al., 2017;Priostosescu and Huybers, 2017). 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 CO 2 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). 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.

ECS for the globally averaged EBM
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  , 2017). 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.

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, F O values are kept constant at their minimal values: 10 W m −2 for the Arctic, 14 W m −2 for the Antarctic, and −39 W 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 CO 2 concentrations since the industrial revolution, these ECS values are consistent with observations to date.

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 CO 2 continue on their current pathway. However, if the increase in atmospheric CO 2 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 (Pierrehumbert, 2010), 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.
where C DL 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, F C , which replaces the F C term that was introduced in Dortmans et al. (2019). This term is defined here as a function of surface temperature T S .
Here, P 0 is the atmospheric pressure at the surface (Z = 0) and R A is the ideal gas constant specific to dry air. This equation is scaled by 1 σ T 4 R to nondimensionalize it, creating in Eq. (10), where T R = 273.15 K is the reference temperature. As this f C 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 F C , used by the authors in Dortmans et al. (2019), was a simple functional form calibrated to empirical data. The result was a relationship between F C and T S that is quantitatively very similar to that given by Eq. (A3). In Kypke and Langford (2020), F C 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 https://doi.org/10.5281/zenodo.3834344 (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 (Kypke, 2019). 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.