An energy balance model for paleoclimate transitions

A new energy balance model (EBM) is presented and is used to study paleoclimate transitions. While most previous EBMs only dealt with the globally averaged climate, this new EBM has three variants: Arctic, Antarctic and tropical climates. The EBM incorporates the greenhouse warming effects of both carbon dioxide and water vapour, and also includes ice–albedo feedback and evapotranspiration. The main conclusion to be inferred from this EBM is that the climate system may possess multiple equilibrium states, both warm and frozen, which coexist mathematically. While the actual climate can exist in only one of these states at any given time, the EBM suggests that climate can undergo transitions between the states via mathematical saddle-node bifurcations. This paper proposes that such bifurcations have actually occurred in Paleoclimate transitions. The EBM is applied to the study of the Pliocene paradox, the glaciation of Antarctica and the so-called warm, equable climate problem of both the mid-Cretaceous Period and the Eocene Epoch. In all cases, the EBM is in qualitative agreement with the geological record.

1 Introduction

For approximately 75 % of the last 540 million years of the paleoclimate history of the Earth, the climate of both polar regions was mild and free of permanent ice-caps (Cronin, 2010; Crowley, 2000; Hubert et al., 2000). Today, both the North and South poles are ice-capped; however, there is overwhelming evidence that these polar ice-caps are melting. The Arctic is warming faster than any other region on Earth. The formation of the present-day Arctic and Antarctic ice-caps occurred abruptly, at widely separated times in the geological history of the Earth. This paper explores some of the underlying mechanisms and forcing factors that have caused climate transitions in the past, with a focus on transitions that were abrupt. The understanding gained here will be applied in a subsequent paper to the important problem of determining whether anthropogenic climate change may lead to an abrupt transition in the future.

We present a new two-layer energy balance model (EBM) for the climate of the Earth. General knowledge of climate and of climate change has been advanced by many studies employing simple EBMs (Budyko, 1968; Kaper and Engler, 2013; McGehee and Lehman, 2012; North et al., 1981; Payne et al., 2015; Sagan and Mullen, 1972; Sellers, 1969; Stap et al., 2017). In general, these EBMs facilitate exploration of the relationship between specific climate forcing mechanisms and the resulting climate changes. The EBM presented here includes a more accurate representation of the role of greenhouse gases in climate change than has been the case for previous EBMs. The model is based on fundamental principles of atmospheric physics, such as the Beer–Lambert law, the Stefan–Boltzmann law, the Clausius–Clapeyron equation and the ideal-gas equation. In particular, the modelling of water vapour acting as a greenhouse gas in the atmosphere, presented in Sect. 2.3.3, is more physically accurate than in previous EBMs and it shows that water vapour feedback is important in climate change. Also, ice–albedo feedback plays a central role in this EBM. The nonlinearity of this EBM leads to bistability (existence of multiple stable equilibrium states), to hysteresis (the climate state realized in the model depends on the past history) and to bifurcations (abrupt transitions from one state to another). Forcing factors that are included explicitly in the model include insolation, CO₂ concentration, relative humidity, evapotranspiration, ocean heat transport and atmospheric heat transport. Other factors that may affect climate change, such as geography, precipitation, vegetation and continental drift, are not included explicitly in the EBM but are present only implicitly in so far as they affect those included factors.

During the Pliocene Epoch, 2.6–5.3 Ma, the climate of the Arctic region of Earth changed abruptly from ice-free to ice-capped. The major climate forcing factors (solar constant, orbital parameters, CO₂ concentration and locations of the continents) were all very similar to today (Fedorov et al., 2006, 2010; Haywood et al., 2009; Lawrence et al., 2009; De Schepper et al., 2015). Therefore, it is difficult to explain why the early Pliocene climate was so different from that of today. That problem has been called the Pliocene paradox (Cronin, 2010; Fedorov et al., 2006, 2010). Currently, there is great interest in mid-Pliocene climate as a natural analogue of the future warmer climate expected for Earth, due to anthropogenic forcing. In recent years, significant progress has been achieved in understanding Pliocene climate, e.g. by Haywood et al. (2009); Salzmann et al. (2009); Ballantyne et al. (2010); Steph et al. (2010); von der Heydt and Dijkstra (2011); Zhang and Yan (2012); Zhang et al. (2013); Sun et al. (2013); De Schepper et al. (2013); Tan et al. (2017); Chandan and Peltier (2017, 2018) and Fletcher et al. (2018); see also references therein. This paper proposes that bistability and bifurcation may have played a fundamental role in determining the Pliocene climate.

The waning of the warm ice-free paleoclimate at the South Pole, leading eventually to abrupt glaciation of Antarctica at the Eocene–Oligocene transition (EOT) about 34 Ma, is believed to have been caused primarily by two major geological changes (although other factors played a role). One is the movement of the continent of Antarctica from the Southern Pacific Ocean to its present position over the South Pole, followed by the development of the Antarctic Circumpolar Current (ACC), thus drastically reducing ocean heat transport to the South Pole (Cronin, 2010; Scher et al., 2015). The second is the gradual draw down in CO₂ concentration worldwide (DeConto et al., 2008; Goldner et al., 2014; Pagani et al., 2011). The EBM of this paper includes both CO₂ concentration and ocean heat transport as explicit forcing factors. The results of this paper suggest that the abrupt glaciation of Antarctica at the EOT was the result of a bifurcation that occurred as both of these factors changed incrementally; and furthermore that both are required to explain the abrupt onset of Antarctic glaciation.

In the mid-Cretaceous Period (100 Ma) the climate of the entire Earth was much more equable than it is today. This means that, compared to today, the pole-to-Equator temperature gradient was much smaller and also the summer/winter variation in temperature at mid- to high-latitudes was much less (Barron, 1983). The differences in forcing factors between the Cretaceous Period and modern times appeared to be insufficient to explain this difference in the climates. Eric Barron called this the warm, equable Cretaceous climate problem and he explored this problem in a series of pioneering papers (Barron, 1983; Barron et al., 1981, 1995; Sloan and Barron, 1992). In a similar vein, Huber and Caballero (2011) and Sloan and Barron (1990, 1992) observed that the early Eocene (56–48 Ma) encompasses the warmest climates of the past 65 million years, yet climate modelling studies had difficulty explaining such warm and equable temperatures. Therefore, this situation was called the early Eocene warm, equable climate problem. The EBM of the present paper suggests that the mathematical mechanism of bistability provides plausible answers to both the mid-Cretaceous and the early Eocene equable climate problems. In fact, from the perspective of this simple model, these two are mathematically the same problem. Therefore, we call them collectively the warm, equable Paleoclimate problem. Recent progress in paleoclimate science has succeeded in narrowing the gap between proxy and general circulation model (GCM) estimates for the Cretaceous climate (Donnnadieu et al., 2006; Ladant and Donnadieu, 2016; O'Brien et al., 2017) and for the Eocene climate (Baatsen et al., 2018; Hutchinson et al., 2018; Lunt et al., 2016, 2017).

The principal contribution of this paper is a simple climate EBM, based on fundamental physical laws, that exhibits bistability, hysteresis and bifurcations. We propose that these three phenomena have occurred in the paleoclimate record of the Earth and they help to explain certain paleoclimate transitions and puzzles as outlined above. A key property of this EBM is that its underlying physical principles are highly nonlinear. As is well known, nonlinear equations can have multiple solutions, unlike linear equations, which can have only one unique solution (if well-posed). In our EBM, the same set of equations can have two or more coexisting stable solutions (bistability); for example, an ice-capped solution (like today's climate) and an ice-free solution (like the Cretaceous climate), even with the same values of the forcing parameters. The determination of which solution is actually realized by the planet at a given time is dependent on past history (hysteresis). Changes in forcing parameters may drive the system abruptly from one stable state to another, at so-called “tipping points”. In this paper, these tipping points are investigated mathematically, and are shown to be bifurcation points, which are investigated using mathematical bifurcation theory. Bifurcation theory tells us that the existence of bifurcation points is preserved (but the numerical values may change) under small deformations of the model equations. Thus, even though this conceptual model may not give us precise quantitative information about climate changes, qualitatively there is good reason to believe that the existence of the bifurcation points in the model will be preserved in similar more refined models and in the real world.

Previous climate models exhibiting multiple equilibrium states have indicated that bifurcations may cause abrupt climate transitions (Budyko, 1968; Sellers, 1969; North et al., 1981; Paillard, 1998; Alley et al., 2003; Rial et al., 2004; Lindsay and Zhang, 2005; Ferreira et al., 2010; Thorndike, 2012; Payne et al., 2015; Stap et al., 2017). None of these authors employ the recently developed mathematical theory of bifurcations to the extent used in this paper. Abrupt transitions in Quaternary glaciations have been studied by Ganopolski and Rahmstorf (2001); Paillard (2001a, b); Calov and Ganopolsky (2005) and Robinson et al. (2012). These glacial cycles are strongly influenced by orbital forcings, which are not within the scope of this paper. Scientists have long considered that abrupt transitions (or bifurcations) in the Atlantic Meridional Overturning Circulation (AMOC) were possible, and that this could contribute to abrupt climate change (Stommel, 1961; Rahmstorf, 1995; Ganopolski and Rahmstorf, 2001; Lenton et al., 2012). This type of change in the AMOC is sometimes called a Stommel bifurcation. However, that phenomenon is also outside the scope of the EBM of this paper, which does not include ocean geography or meridional dependence. (An extension of this EBM, to a partial differential equation (PDE) spherical shell model, is planned.)

With regard to the EOT, a variety of indicators, including the analysis of fossil plant stomata (Steinthorsdottir et al., 2016), imply that decreasing atmospheric pCO₂ in the Eocene preceded the large shift in oxygen isotope records that characterize the EOT. It was hypothesized in Steinthorsdottir et al. (2016) that at the EOT, a certain threshold of pCO₂ was crossed, resulting in an abrupt change of climate mode. Possible mechanisms for the threshold that they hypothesized are studied in Sect. 3.2.

Even for changes as recent as the mid-Holocene, there is a debate, for example, over whether the abrupt desertification of the Sahara is due to a bifurcation, as suggested in Claussen et al. (1999) using an Earth model of intermediate complexity (EMIC), or is a transient response of the AMOC to a sudden termination of freshwater discharge to the North Atlantic, as proposed in Liu et al. (2009), using a coupled atmosphere–ocean GCM (AOGCM). Similarly, for the glaciation of Greenland at the Pliocene–Pleistocene transition, recent work in Tan et al. (2018), using a coupled GCM–ice-sheet model, shows good agreement with proxy records without the need for bifurcations. The model of this paper does not adapt to these two situations, which are localized and away from the poles where the axis of symmetry restricts the dynamics and facilitates the analysis presented here.

Further investigation of climate changes, using a range of climate models from EBM and EMIC to GCM and AOGCM, is warranted to clarify the underlying mechanisms of abrupt climate change. Very sophisticated GCMs, which include many 3-D processes, are only able to run a few climate trajectories; while EBMs and EMICs may explore more possibilities and investigate climate transitions (tipping points) with major simplifications and with less effort. More rigorous mathematical analysis is possible on small models, which may then suggest lines of inquiry on large models. The Earth climate system is extremely complex. For best results, a hierarchy of climate models is necessary.

The energy balance model presented here is conceptual and qualitative. It contains many simplifying assumptions and is not intended to give a detailed description of the climate of the Earth with quantitative precision. It complements but does not replace more detailed GCMs. Geographically, this model is as simple as possible. It follows in a long tradition of slab models of the atmosphere. Previous slab models represented the atmosphere of the Earth as a globally averaged uniform slab at a single temperature T. The temperature T is determined in those models by a global energy balance equation of the form energy in = energy out. Such models are unable to differentiate between different climates at different latitudes; for example, if the polar climate is changing more rapidly than the tropical climate. In this new model, the forcing parameters of the slab atmosphere are chosen to represent each one of three particular latitudes: Arctic, Antarctic or tropics. Each of these regions is represented by its own slab model, with its own forcing parameters and its own surface temperature T_S. In addition, each region has its own variable I_A representing the intensity of the radiation re-emitted by the atmosphere. The two independent variables I_A and T_S are determined in each model by two energy balance equations, expressing energy balance in the atmosphere and energy balance at the surface, respectively. In this way, the different climate responses of these three regions to their respective forcings can be explored.

Table 1 Summary of variables and parameters used in the model. Many of these parameters have standard textbook values (e.g. $T_{\mathrm{R}},\mathit{\sigma },L_{\mathrm{v}},R_{\mathrm{W}}$). The standard lapse rate Γ is from ICAO (1993). Other parameters such as k_C and k_W are derived in Appendix Appendix A.

The role played by greenhouse gases in climate change is a particular focus of this model. Greenhouse gases trap heat emitted by the surface and are major contributors to global warming. The very different roles of the two principal greenhouse gases in the atmosphere, carbon dioxide and water vapour, are analyzed here in Sects.  2.3.2 and 2.3.3, respectively. The greenhouse warming effect of CO₂ increases with the density of the atmosphere but is independent of temperature, while the greenhouse warming of H₂O increases with temperature but is independent of the density (or partial pressure) of the other gases present. The greenhouse warming of methane (CH₄) acts in a similar fashion to that of CO₂; therefore, CH₄ can be incorporated into the CO₂ concentration. As an increase in CO₂ concentration causes climate warming, this warming causes an increase in evaporation of H₂O into the atmosphere, which further increases the climate warming beyond that due to CO₂ alone (this is true both in the model and in the real atmosphere). This effect is known as water vapour feedback. The energy balance model presented here is the first EBM to incorporate these important roles of the greenhouse gases in such detail.

The paper concludes with two appendices. In Appendix Appendix A, model parameters that are difficult or impossible to determine for paleoclimates are calibrated using the abundant satellite and surface data available for today's climate. In addition, justification is given for parameter values chosen for the model. In Appendix Appendix B the paleoclimate model of this paper is adapted to modern-day conditions and its equilibrium climate sensitivity (ECS) is determined. Here, ECS is the change in global mean temperature produced by a doubling of CO₂ in the model, starting from the pre-industrial value of 270 ppm. For this EBM, the ECS is determined to be ΔT=3.3 ^∘C, which is at the high end of the range accepted by the IPCC (IPCC, 2013).

2 The energy balance climate model

In this EBM, the atmosphere and surface are each assumed to be in energy balance. Short-wave radiant energy from the sun is partly reflected by the atmosphere back into space, a small portion is absorbed directly by the atmosphere and the remainder passes through to the surface. The surface reflects some of this short-wave energy (which is assumed to escape to space) and absorbs the rest, re-emitting long-wave radiant energy of intensity I_S upward into the atmosphere. The atmosphere is modelled as a slab, with greenhouse gases, that absorbs a fraction η of the radiant energy I_S from the surface. The atmosphere re-emits radiant energy of total intensity I_A. Of this radiation I_A, a fraction β is directed downward to the surface, and the remaining fraction (1−β) goes upward and escapes to space.

This model is based on the uniform slab EBM used in Payne et al. (2015), modified as shown in Fig. 1. In our case, the “slab” is a uniform column of air of unit cross section, extending vertically above the surface to the tropopause, and located either at a pole or at the Equator. The symbols in Fig. 1 are defined in Table 1.

Figure 1 A visualization of the energy balance model. Symbols are defined in Table 1 and Sect. 2.1.

This section presents the mathematical derivation of the EBM. Readers interested only in the climate applications of the model may skip this section and go directly to Sect. 3. A preliminary version of this EBM was presented in a conference proceedings paper (Dortmans et al., 2018). The present model incorporates several important improvements over that model. The differences between that model and the one presented here are indicated where appropriate in the text. Furthermore, the previous paper (Dortmans et al., 2018) considered only the application of the model to the Arctic climate and the Pliocene paradox; it did not study Antarctic or tropical climate or the Cretaceous warm, equable climate problem, as does the present paper.

2.1 Energy balance

The model consists of two energy balance equations, one for the atmosphere and one for the surface. In Fig. 1, the so-called forcings are shown as arrows, pointing in the direction of energy transfer. From Fig. 1, the energy balance equations for the atmosphere and surface are given, respectively, by

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{0}& {\displaystyle }=F_{\mathrm{A}}+F_{\mathrm{C}}+{\mathit{\xi }}_{\mathrm{A}}Q+\mathit{\eta }{I}_{\mathrm{S}}-I_{\mathrm{A}},\mathrm{and}\text{(2)}& {\displaystyle }\mathrm{0}& {\displaystyle }=F_{\mathrm{O}}-F_{\mathrm{C}}+(\mathrm{1}-\mathit{\alpha })F_{\mathrm{S}}-I_{\mathrm{S}}+\mathit{\beta }{I}_{\mathrm{A}}.\end{array}$$

Symbols and parameter values for the model are defined in Table 1. Appendix A4 provides derivations and justification for the values of the empirical parameters. The forcings F_O and F_A represent ocean and atmosphere heat transport, respectively, and are specified as constants for each region of interest. Heat transport by conduction/convection from the surface to the atmosphere is denoted F_C. This quantity will be largely dependent on surface temperature, T_S. As described in Appendix Appendix A we have modelled it as a hyperbola that is mostly flat for temperatures below freezing, and grows roughly linearly for temperatures above freezing so that

$$\begin{array}{}\text{(3)}& F_{\mathrm{C}}=A_{\mathrm{1}}\left(T_{\mathrm{S}}-T_{\mathrm{R}}\right)+\sqrt{A_{\mathrm{1}}^{\mathrm{2}}(T_{\mathrm{S}}-T_{\mathrm{R}}{)}^{\mathrm{2}}+A_{\mathrm{2}}^{\mathrm{2}}},\end{array}$$

where A₁ and A₂ are constants. Since the model is concerned with temperatures around the freezing point of water, we set this as a reference temperature, T_R=273.15 K.

The annually averaged intensity of solar radiation striking a surface parallel to the Earth's surface but at the top of the atmosphere is Q. The value of Q at either Pole is $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$ and at the Equator is 418.8 Wm⁻² (McGehee and Lehman, 2012; Kaper and Engler, 2013). A fraction ξ_R of this short-wave radiation is reflected by the atmosphere back into space and a further fraction ξ_A is directly absorbed by the atmosphere; the remainder penetrates to the surface. See Appendix A1 for the derivation of values for ξ_R and ξ_A. The solar radiation striking the surface of the Earth is

$$\begin{array}{}\text{(4)}& F_{\mathrm{S}}=(\mathrm{1}-{\mathit{\xi }}_{\mathrm{R}}-{\mathit{\xi }}_{\mathrm{A}})Q.\end{array}$$

The surface albedo is the fraction, α, of this solar radiation that reflects off the surface back into space. Thus the solar forcing absorbed by the surface is (1−α)F_S, and the solar radiation reflected back to space is αF_S. Typical values of the 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. In this paper we introduce a smoothly varying albedo given by the hyperbolic tangent function:

$$\begin{array}{}\text{(5)}& \mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left([{\mathit{\alpha }}_{\mathrm{W}}+{\mathit{\alpha }}_{\mathrm{C}}]+[{\mathit{\alpha }}_{\mathrm{W}}-{\mathit{\alpha }}_{\mathrm{C}}]\tanh \left({\displaystyle \frac{T_{\mathrm{S}}-T_{\mathrm{R}}}{\mathrm{\Omega }}}\right)\right),\end{array}$$

where α_C and α_W are the albedo values for cold and warm temperatures, respectively, and the parameter Ω determines the steepness of the transition between α_C and α_W. See Appendix A2 for a full explanation of Eq. (3).

The emission of long-wave radiation, I, from a body is governed by the Stefan–Boltzmann law, I=ϵσT⁴, where ϵ is the emissivity, σ is the Stefan–Boltzmann constant, and T is temperature. The surface of the Earth acts as a black-body radiator; so for the radiation emitted from the surface ϵ=1, thus

$$\begin{array}{}\text{(6)}& I_{\mathrm{S}}=\mathit{\sigma }{T}_{\mathrm{S}}^{\mathrm{4}}.\end{array}$$

Previous authors have postulated an idealized uniform atmospheric temperature T_A for the slab model, so that the intensity of radiation emitted by the atmosphere, I_A, is

$$\begin{array}{}\text{(7)}& I_{\mathrm{A}}=\mathit{ϵ}\mathit{\sigma }{T}_{\mathrm{A}}^{\mathrm{4}}.\end{array}$$

The emissivity is ϵ=0.9 since the atmosphere is an imperfect black-body radiator. A uniform temperature for the atmosphere does not exist in the real world, where T_A varies strongly with height, unlike T_S, which has a single value. Previous two-layer EBMs have used (T_S,T_A) as the two independent variables in the two energy balance Eqs. (1) and (1). Here instead, we use (T_S,I_A) as the two independent variables, and then we formally let T_A be defined by Eq. (5). The fraction of I_A that reaches the surface is β=0.63; see Appendix A6.

The parameter η represents the fraction of the infrared radiation I_S from the surface that is absorbed by the atmosphere and is called absorptivity. The major constituents of the atmosphere are nitrogen and oxygen and these gases do not absorb any infrared radiation. The gases that do contribute to the absorptivity η are called greenhouse gases. Chief among these are carbon dioxide and water vapour. The contribution of these two greenhouse gases to η are analyzed in Sects. 2.3.2 and 2.3.3, respectively. Although both contribute to warming of the climate, the underlying physical mechanisms of the two are very different. In general, the contribution to η from water vapour is a function of temperature. Another major contributor to absorption is the liquid and solid water in clouds. We model this portion of the absorption as constant, since we do not include any data on cloud cover variation. However, we experimented with making this portion vary with surface temperature, and the results were not qualitatively different than those presented here.

Nondimensional temperatures

We rescale temperature by the reference temperature T_R=273.15 K and define new nondimensional temperatures and new parameters

$$\begin{array}{ll}{\displaystyle }{\mathit{\tau }}_{\mathrm{S}}& {\displaystyle }={\displaystyle \frac{T_{\mathrm{S}}}{T_{\mathrm{R}}}},i_A={\displaystyle \frac{I_{\mathrm{A}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}},q={\displaystyle \frac{Q}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}},f_O={\displaystyle \frac{F_{\mathrm{O}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}},f_A={\displaystyle \frac{F_{\mathrm{A}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}},\\ \text{(8)}& {\displaystyle }& {\displaystyle }{f}_C={\displaystyle \frac{F_{\mathrm{C}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}},a_{\mathrm{1}}={\displaystyle \frac{A_{\mathrm{1}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{3}}}},a_{\mathrm{2}}={\displaystyle \frac{A_{\mathrm{2}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}},\mathit{\omega }={\displaystyle \frac{\mathrm{\Omega }}{T_{\mathrm{R}}}}.\end{array}$$

After normalization, the freezing temperature of water is represented by τ=1 and the atmosphere and surface energy balance Eqs. (1)–(4) simplify to

$$\begin{array}{}\text{(9)}& {\displaystyle }{i}_A& {\displaystyle }=f_A+f_C\left({\mathit{\tau }}_{\mathrm{S}}\right)+{\mathit{\xi }}_{\mathrm{A}}q+\mathit{\eta }\left({\mathit{\tau }}_{\mathrm{S}}\right){\mathit{\tau }}_{\mathrm{S}}^{\mathrm{4}},\text{(10)}& {\displaystyle }{i}_A& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}\left[f_C\left({\mathit{\tau }}_{\mathrm{S}}\right)-f_O-(\mathrm{1}-\mathit{\alpha }({\mathit{\tau }}_{\mathrm{S}}\left)\right)(\mathrm{1}-{\mathit{\xi }}_{\mathrm{R}}-{\mathit{\xi }}_{\mathrm{A}})q+{\mathit{\tau }}_{\mathrm{S}}^{\mathrm{4}}\right],\end{array}$$

where, from Appendix Appendix A,

$$\begin{array}{}\text{(11)}& {\displaystyle }{f}_C\left({\mathit{\tau }}_{\mathrm{S}}\right)& {\displaystyle }=a_{\mathrm{1}}\left({\mathit{\tau }}_{\mathrm{S}}-\mathrm{1}\right)+\sqrt{a_{\mathrm{1}}^{\mathrm{2}}({\mathit{\tau }}_{\mathrm{S}}-\mathrm{1}{)}^{\mathrm{2}}+a_{\mathrm{2}}^{\mathrm{2}}},\text{(12)}& {\displaystyle }\mathit{\alpha }\left({\mathit{\tau }}_{\mathrm{S}}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left([{\mathit{\alpha }}_{\mathrm{W}}+{\mathit{\alpha }}_{\mathrm{C}}]+[{\mathit{\alpha }}_{\mathrm{W}}-{\mathit{\alpha }}_{\mathrm{C}}]\tanh \left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{S}}-\mathrm{1}}{\mathit{\omega }}}\right)\right).\end{array}$$

The range of surface temperatures T_S observed on Earth is restricted to an interval around the freezing point of water, 273.15 K, and therefore the nondimensional temperature τ lies in an interval around τ=1. In this paper, we assume $\mathrm{0.8}\le \mathit{\tau }\le \mathrm{1.2}$, which corresponds approximately to a range in more familiar degrees Celsius of $-\mathrm{54}{}^{\circ }\mathrm{C}\le T\le +\mathrm{54}$ ^∘C. Another reason for the upper limit on temperature is that the Clausius–Clapeyron law used in Sect. 2.3.3 fails to apply at temperatures above the boiling point of water.

2.2 Optical depth and the Beer–Lambert law

The goal of this section is to define the absorptivity parameter η in the EBM Eq. (6) (or Eq. 1) in such a way that the atmosphere in the uniform slab model will absorb the same fraction η of the long-wave radiation I_S from the surface as does the real nonuniform atmosphere of the Earth. This absorption is due primarily to water vapour, carbon dioxide and clouds. Previous energy balance models have assigned a constant value to η, often determined by climate data. In the present EBM, η is not constant but is a function of other more fundamental physical quantities, such as $\mathit{\mu },\mathit{\delta },k_{\mathrm{C}},k_{\mathrm{W}}$ and T. This function is determined by classical physical laws. In this way, the present EBM adjusts automatically to changes in these physical quantities, and represents a major advance over previous EBMs.

The Beer–Lambert law states that when a beam of radiation (or light) enters a sample of absorbing material, the absorption of radiation at any point z is proportional to the intensity of the radiation I(z) and also to the concentration or density of the absorber ρ(z). This bilinearity fails to hold at very high intensity of radiation or high density of absorber, neither of which is the case in the Earth's atmosphere. Whether this law is applied to the uniform slab model or to the nonuniform real atmosphere, it yields the same differential equation

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}z}}=-k\mathit{\rho }\left(z\right)I\left(z\right),\end{array}$$

where k (m² kg⁻¹) is the absorption coefficient of the material, ρ (kg m⁻³) is the density of the absorbing substance such as CO₂ and z (m) is distance along the path. The differential Eq. (6) may be integrated from z=0 (the surface) to z=Z (the tropopause) to give

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{I_{\mathrm{T}}}{I_{\mathrm{S}}}}=e^{-{\int }_{\mathrm{0}}^{Z}k\mathit{\rho }\left(z\right)\mathrm{d}z}\equiv e^{-\mathit{\lambda }},\mathrm{where}\mathit{\lambda }\equiv \underset{\mathrm{0}}{\overset{Z}{\int }}k\mathit{\rho }\left(z\right)\mathrm{d}z.\end{array}$$

Here I_T≡I(Z) is the intensity of radiation escaping to space at the tropopause z=Z, and λ is the so-called optical depth of the material. Note that λ is dimensionless. The absorptivity parameter η in Eqs. (1) and (6) represents the fraction of the outgoing radiation I_S from the surface that is absorbed by the atmosphere (not to be confused with the absorption coefficient k). It follows from the Beer–Lambert law that η is completely determined by the corresponding optical depth parameter λ; that is

$$\begin{array}{}\text{(15)}& \mathit{\eta }={\displaystyle \frac{I_{\mathrm{S}}-I_{\mathrm{T}}}{I_{\mathrm{S}}}}=\mathrm{1}-e^{-\mathit{\lambda }}.\end{array}$$

For a mixture of n attenuating materials, with densities ρ_i, absorption coefficients k_i and corresponding optical depths λ_i, the Beer–Lambert law extends to

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{I_{\mathrm{T}}}{I_{\mathrm{S}}}}=e^{-{\mathrm{\Sigma }}_{i=\mathrm{1}}^n{\int }_{\mathrm{0}}^Z{k}_i{\mathit{\rho }}_i\left(z\right)\mathrm{d}z}\equiv e^{-{\mathrm{\Sigma }}_{i=\mathrm{1}}^n{\mathit{\lambda }}_i},\end{array}$$

so that

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathit{\eta }=\mathrm{1}-e^{-{\mathrm{\Sigma }}_{i=\mathrm{1}}^n{\mathit{\lambda }}_i}=\mathrm{1}-\prod _{i=\mathrm{1}}^n{e}^{-{\mathit{\lambda }}_i}=\mathrm{1}-\prod _{i=\mathrm{1}}^n(\mathrm{1}-{\mathit{\eta }}_i),\end{array}$$

where ${\mathit{\eta }}_i=\mathrm{1}-e^{-{\mathit{\lambda }}_i}$. Equation (10) is the key to solving the problem posed in the first sentence of this subsection. For the ith absorbing material in the slab model, we set its optical depth λ_i to be equal to the value of the optical depth that this gas has in the Earth's atmosphere as given by Eq. (7), and then combine them using Eq. (9). This calculation is presented for the case of CO₂ in Sect. 2.3.2 and for water vapour in Sect. (2.3.3). For the third absorbing material in our model, clouds, we assume a constant value η_Cl.

2.3 Greenhouse gases

The two principal greenhouse gases are carbon dioxide (CO₂) and water vapour (H₂O). Because they act in different ways, we determine the absorptivities η_C and η_W, and optical depths λ_C and λ_W of CO₂ and H₂O separately, and then combine their effects, along with the absorption due to clouds, η_Cl, using the Beer–Lambert law for mixtures, Eq. (10). Methane acts similarly to CO₂ and can be included in the optical depth for CO₂. Other greenhouse gases have only minor influence and are ignored in this paper.

2.3.1 The grey gas approximation

Although it is well-known that gases like CO₂ and H₂O absorb infrared radiation I_S only at specific wavelengths (spectral lines), in this paper the grey gas approximation is used; that is, the absorption coefficient k_C or k_W is given as a single number averaged over the infrared spectrum (Pierrehumbert, 2010). The thesis by Dortmans (2017) presents a survey of values in the literature for the absorption coefficients k_C and k_W of CO₂ and H₂O, respectively, in the grey gas approximation. The values used in this paper are given in Table 1; and are derived as described in Appendix A5.

Figure 2 Dry atmosphere EBM Eqs. (15)–(15) (with $F_{\mathrm{A}}=\mathrm{65}{\mathrm{Wm}}^{-\mathrm{2}}$, $F_{\mathrm{O}}=\mathrm{50}{\mathrm{Wm}}^{-\mathrm{2}}$, $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, α_W=0.08 and Z=9 km). (a) CO₂ values μ=400, 800, 1200, 1600, 2000; from bottom to top blue curves. A unique cold equilibrium exists for μ=400 (bottom blue curve), multiple equilibria exist for μ=800, 1200 and 1600, a unique warm equilibrium exists for μ=2000 (top blue curve). (b) Bifurcation diagram for the same parameters as in (a) showing the result of solving the EBM for τ_S as a function of the forcing parameter μ, holding all other forcing parameters constant. Two saddle-node (or fold) bifurcations occur at approximately μ=681 and μ=1881. Three distinct solutions exist between these two values of μ. As μ decreases from right to left in this figure, if the surface temperature starts in the warm state (τ_S>1), then it abruptly “falls” from the warm state to the frozen state when the left bifurcation point is reached. Similarly, starting on the frozen state (τ_S<1), if μ increases sufficiently, then the state will jump upward to the warm branch when the right bifurcation point is reached. Both of these transitions are irreversible (one way).

2.3.2 Carbon dioxide

The concentration of CO₂ in the atmosphere is usually expressed as a ratio, in molar parts per million (ppm) of CO₂ to dry air, and written as μ. There is convincing evidence that μ has varied greatly in the geological history of the Earth, and has decreased slowly over the past 100 million years; however, today μ is increasing due to human activity. The value before the industrial revolution was μ=270 ppm, but today μ is slightly above 400 ppm.

Although traditionally μ is measured as a ratio of molar concentrations, in practice both the density ρ and the absorption coefficient k are expressed in mass units of kilograms. Therefore, before proceeding, μ must be converted from a molar ratio to a ratio of masses in units of kilograms. The mass of one mole of CO₂ is approximately $m{m}_C=\mathrm{44}\times {\mathrm{10}}^{-\mathrm{3}}$ kg mol⁻¹. The dry atmosphere is a mixture of 78 % nitrogen, 21 % oxygen and 0.9 % argon, with molar masses of 28, 32 and 40 g mol⁻¹, respectively. Neglecting other trace gases in the atmosphere, a weighted average gives the molar mass of the dry atmosphere as mm${}_{\mathrm{A}}=\mathrm{29}\times {\mathrm{10}}^{-\mathrm{3}}$ kg mol⁻¹. Therefore, the CO₂ concentration μ measured in molar ppm is converted to mass concentration in kg ppm by multiplication by the ratio mm${}_C/{\mathrm{mm}}_{\mathrm{A}}\approx \mathrm{1.52}$. If ρ_A(z) is the density of the atmosphere at altitude z in kg m⁻³, then the mass density of CO₂ at the same altitude, with molar concentration μ ppm, is

$$\begin{array}{}\text{(18)}& {\mathit{\rho }}_C\left(z\right)=\mathrm{1.52}{\displaystyle \frac{\mathit{\mu }}{{\mathrm{10}}^{\mathrm{6}}}}{\mathit{\rho }}_{\mathrm{A}}\left(z\right)\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}.\end{array}$$

It is known that CO₂ disperses rapidly throughout the Earth's atmosphere, so that its concentration μ may be assumed independent of location and altitude (IPCC, 2013). As the density of the atmosphere decreases with altitude, the density of CO₂ decreases at exactly the same rate, according to Eq. (11). Substituting Eq. (11) into Eq. (7) determines the optical depth λ_C of CO₂

$$\begin{array}{}\text{(19)}& {\mathit{\lambda }}_C=\mathrm{1.52}{\displaystyle \frac{\mathit{\mu }}{{\mathrm{10}}^{\mathrm{6}}}}{k}_{\mathrm{C}}\underset{\mathrm{0}}{\overset{Z}{\int }}{\mathit{\rho }}_{\mathrm{A}}\left(z\right)\mathrm{d}z.\end{array}$$

Now consider a vertical column of air, of unit cross section, from surface to tropopause. The integral in Eq. (12) is precisely the total mass of this column. The atmospheric pressure at the surface is P_A, and this is the total weight of the column. Therefore ${\int }_{\mathrm{0}}^Z\mathit{\rho }\left(A\right)\mathrm{d}z=P_{\mathrm{A}}/g$, where g is acceleration due to gravity. Therefore, the optical depth λ_C of CO₂ in the actual atmosphere in Eq. (12) is

$$\begin{array}{}\text{(20)}& {\mathit{\lambda }}_C=\mathit{\mu }{G}_C,\text{where}{G}_C\equiv \mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{6}}{k}_{\mathrm{C}}{\displaystyle \frac{P_{\mathrm{A}}}{g}}.\end{array}$$

With λ_C so determined, it follows from Eq. (8) that

$$\begin{array}{}\text{(21)}& {\mathit{\eta }}_C=\mathrm{1}-e^{-{\mathit{\lambda }}_C}=\mathrm{1}-\exp (-\mathit{\mu }{G}_C).\end{array}$$

As listed in Table 1 and derived in Appendix A5, the calibrated value for k_C is 0.07424, and therefore the value for the greenhouse gas parameter for CO₂ is $G_C=\mathrm{1.166}\times {\mathrm{10}}^{-\mathrm{3}}$.

The dry atmosphere EBM is obtained by assuming there is no water vapour and no clouds in the atmosphere. Hence η=η_C and $\mathit{\delta }=\mathrm{0}.$ We assume also that ${\mathit{\xi }}_{\mathrm{R}}={\mathit{\xi }}_{\mathrm{A}}=\mathrm{0}$. Making these substitutions in the EBM Eqs. (6)–(6) yields

$$\begin{array}{}\text{(22)}& {\displaystyle }{i}_A& {\displaystyle }=f_A+f_C\left({\mathit{\tau }}_{\mathrm{S}}\right)+[\mathrm{1}-\exp (-\mathit{\mu }{G}_C\left)\right]\cdot {\mathit{\tau }}_{\mathrm{S}}^{\mathrm{4}},\text{(23)}& {\displaystyle }{i}_A& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}\left(f_C\left({\mathit{\tau }}_{\mathrm{S}}\right)-f_O-[\mathrm{1}-\mathit{\alpha }({\mathit{\tau }}_{\mathrm{S}}\left)\right]q+{\mathit{\tau }}_{\mathrm{S}}^{\mathrm{4}}\right).\end{array}$$

Due to the nonlinearity of the ice–albedo function in 3, there can be in fact up to three points of intersection as shown in Fig. 2a. If the CO₂ level μ decreases sufficiently, the warm state (τ_S>1) disappears, while as CO₂ increases, the frozen state (τ_S<1) may disappear. The CO₂ level μ is quite elevated in Fig. 2, but the effects of water vapour as a greenhouse gas and of clouds have been ignored.

Figure 2b introduces a bifurcation diagram (Kuznetsov, 2004), in which the two EBM Eqs. (15)–(15) have been solved for the surface temperature τ_S, which is then plotted as a function of the parameter μ. Figure 2b shows the three distinct solutions τ_S of the EBM: a warm solution (τ_S>1), a frozen solution (τ_S<1), and a third solution that crosses through τ_S=1 and connects the other two solutions in saddle-node bifurcations (Kuznetsov, 2004). Stability analysis (Dortmans, 2017) shows that the warm and frozen solutions are stable (in a dynamical systems sense), while the third solution (denoted by a dashed line) is unstable.

2.3.3 Water vapour and the Clausius–Clapeyron equation

In this section, we determine the absorptivity of water vapour as a function of temperature, η_W(τ_S), using fundamental physical laws including the Clausius–Clapeyron equation, the ideal-gas law and the Beer–Lambert law (Pierrehumbert, 2010). We also assume the idealized lapse rate of the International Standard Atmosphere (ISA) as defined by ICAO (ICAO, 1993). Unlike CO₂, the concentration of water vapour (H₂O) in the atmosphere varies widely with location and altitude. This is because the partial pressure of H₂O varies strongly with the local temperature. In fact, it is bounded by a maximum saturated value, which itself is a nonlinear function of temperature, $P_{\mathrm{W}}^{\mathrm{sat}}\left(T\right)$. The actual partial pressure of water vapour is then a fraction, δ, of this saturated value,

$$\begin{array}{}\text{(24)}& P_{\mathrm{W}}\left(T\right)=\mathit{\delta }{P}_{\mathrm{W}}^{\mathrm{sat}}\left(T\right),\mathrm{0}\le \mathit{\delta }\le \mathrm{1},\end{array}$$

where δ is called relative humidity. While $P_{\mathrm{W}}^{\mathrm{sat}}\left(T\right)$ varies greatly with T in the atmosphere, the relative humidity δ is comparatively constant. If the actual P_W(T) exceeds the saturated value $P_{\mathrm{W}}^{\mathrm{sat}}\left(T\right)$ (i.e. δ>1), the excess water vapour condenses out of the atmosphere and falls as rain or snow.

The saturated partial pressure at temperature T is determined by the Clausius–Clapeyron equation (Pierrehumbert, 2010),

$$\begin{array}{}\text{(25)}& P_{\mathrm{W}}^{\mathrm{sat}}\left(T\right)=P_{\mathrm{W}}^{\mathrm{sat}}\left(T_{\mathrm{R}}\right)\exp \left({\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{W}}}}\left[{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{R}}}}-{\displaystyle \frac{\mathrm{1}}{T}}\right]\right),\end{array}$$

where T_R is the reference temperature, here chosen to be the freezing point of water (273.15 K); L_v is the latent heat of vaporization of water and R_W is the ideal-gas constant for water, see Table 1. The actual partial pressure of water vapour at relative humidity δ and temperature T is then given by combining Eqs. (15) and (16).

We may use the ideal-gas law in the form P_W=ρ_WR_WT to convert the partial pressure P_W of water vapour at temperature T to mass density ρ_W of water vapour at that temperature. Substituting into Eqs. (15) and (16) gives

$$\begin{array}{}\text{(26)}& {\mathit{\rho }}_{\mathrm{W}}\left(T\right)=\mathit{\delta }{\displaystyle \frac{P_{\mathrm{W}}^{\mathrm{sat}}\left(T_{\mathrm{R}}\right)}{R_{\mathrm{W}}T}}\exp \left({\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{W}}}}\left[{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{R}}}}-{\displaystyle \frac{\mathrm{1}}{T}}\right]\right).\end{array}$$

Transforming to the dimensionless temperature $\mathit{\tau }=T/T_{\mathrm{R}}$ as in Sect. 2.1.1, this becomes

$$\begin{array}{}\text{(27)}& {\mathit{\rho }}_{\mathrm{W}}\left(\mathit{\tau }\right)=\mathit{\delta }\left({\displaystyle \frac{P_{\mathrm{W}}^{\mathrm{sat}}\left(T_{\mathrm{R}}\right)}{R_{\mathrm{W}}{T}_{\mathrm{R}}}}\right){\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\exp \left({\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{W}}{T}_{\mathrm{R}}}}\left[{\displaystyle \frac{\mathit{\tau }-\mathrm{1}}{\mathit{\tau }}}\right]\right).\end{array}$$

The Beer–Lambert law in Sect. 2.2 implies that the absorptivity of a greenhouse gas η_i is completely determined by its optical depth λ_i. For water vapour, from Eqs. (8) and (7),

$$\begin{array}{}\text{(28)}& {\mathit{\eta }}_{\mathrm{W}}=\mathrm{1}-e^{-{\mathit{\lambda }}_W}\mathrm{where}{\mathit{\lambda }}_W\equiv \underset{\mathrm{0}}{\overset{Z}{\int }}{k}_{\mathrm{W}}{\mathit{\rho }}_{\mathrm{W}}\left(z\right)\mathrm{d}z.\end{array}$$

Here k_W is the absorption coefficient of water vapour; see Appendix A5. In order to evaluate the integral in Eq. (19), we need to know how ρ_W varies with height z. We have shown that ρ_W is a function of temperature, given by Eqs. (17) or (18). Therefore, we need an expression for the variation in temperature T with height z. Under normal conditions, the temperature T decreases with height in the troposphere. This rate of decrease is called the lapse rate Γ, and is defined as

$$\begin{array}{}\text{(29)}& \mathrm{\Gamma }\equiv -{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}.\end{array}$$

Normally, Γ is positive and is close to constant in value from the surface to the tropopause. The International Civil Aviation Organization has defined, for reference purposes, the ISA, in which Γ is assigned the constant value $\mathrm{\Gamma }=\mathrm{6.49}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{K}{\mathrm{m}}^{-\mathrm{1}}$ (ICAO, 1993). Using this assumption, the variation in temperature with height is given as

$$\begin{array}{}\text{(30)}& T\left(z\right)=T_{\mathrm{S}}-\mathrm{\Gamma }z\mathrm{or}\mathit{\tau }\left(z\right)={\mathit{\tau }}_{\mathrm{S}}-\mathit{\gamma }z,\end{array}$$

where the normalized lapse rate is $\mathit{\gamma }=\mathrm{\Gamma }/T_{\mathrm{R}}=\mathrm{2.38}\times {\mathrm{10}}^{-\mathrm{5}}{\mathrm{m}}^{-\mathrm{1}}$. The tropopause height Z ranges from 8 to 11 km at the poles and 16 to 18 km at the Equator, based on satellite measurements (Kishore et al., 2006). Therefore, both T(Z) and τ(Z) are positive. For this paper we will take the height to be Z=9 km at the poles and Z=17 km at the Equator.

Equation (21) may be used to change the variable of integration in Eq. (19) for the optical depth of water vapour, from z to τ. The result is

$$\begin{array}{}\text{(31)}& {\mathit{\lambda }}_W=k_{\mathrm{W}}\underset{\mathrm{0}}{\overset{Z}{\int }}{\mathit{\rho }}_{\mathrm{W}}\left(\mathit{\tau }\right(z\left)\right)\mathrm{d}z={\displaystyle \frac{k_{\mathrm{W}}}{\mathit{\gamma }}}\underset{{\mathit{\tau }}_{\mathrm{S}}-\mathit{\gamma }Z}{\overset{{\mathit{\tau }}_{\mathrm{S}}}{\int }}{\mathit{\rho }}_{\mathrm{W}}\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }.\end{array}$$

Now substitute Eq. (18) into the integral Eq. (22) and simplify to

$$\begin{array}{}\text{(32)}& {\mathit{\lambda }}_W\left({\mathit{\tau }}_{\mathrm{S}}\right)=\mathit{\delta }{G}_{\mathrm{W}\mathrm{2}}\underset{{\mathit{\tau }}_{\mathrm{S}}-\mathit{\gamma }Z}{\overset{{\mathit{\tau }}_{\mathrm{S}}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\exp \left(G_{\mathrm{W}\mathrm{1}}\left[{\displaystyle \frac{\mathit{\tau }-\mathrm{1}}{\mathit{\tau }}}\right]\right)\mathrm{d}\mathit{\tau },\end{array}$$

where the greenhouse gas parameters G_W1 and G_W2 for water vapour are defined as

$$\begin{array}{}\text{(33)}& G_{\mathrm{W}\mathrm{1}}\equiv {\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{W}}{T}_{\mathrm{R}}}}\mathrm{and}{G}_{\mathrm{W}\mathrm{2}}\equiv {\displaystyle \frac{k_{\mathrm{W}}{P}_{\mathrm{W}}^{\mathrm{sat}}\left(T_{\mathrm{R}}\right)}{R_{\mathrm{W}}{T}_{\mathrm{R}}\mathit{\gamma }}}.\end{array}$$

Finally, using the Beer–Lambert law, the absorptivity η_W of water vapour in Eq. (19) is determined by its optical depth, and is now a function of the surface temperature,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathrm{W}}\left({\mathit{\tau }}_{\mathrm{S}}\right)=\mathrm{1}-\exp \left[-{\mathit{\lambda }}_W\left({\mathit{\tau }}_{\mathrm{S}}\right)\right]\\ \text{(34)}& {\displaystyle }& {\displaystyle }=\mathrm{1}-\exp \left[-\mathit{\delta }{G}_{\mathrm{W}\mathrm{2}}\underset{{\mathit{\tau }}_{\mathrm{S}}-\mathit{\gamma }Z}{\overset{{\mathit{\tau }}_{\mathrm{S}}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\exp \left(G_{\mathrm{W}\mathrm{1}}\left[{\displaystyle \frac{\mathit{\tau }-\mathrm{1}}{\mathit{\tau }}}\right]\right)\mathrm{d}\mathit{\tau }\right].\end{array}$$

The definite integral in this expression is easily evaluated numerically in the process of solving the atmosphere and surface EBM equations. As given in Table 1 and derived in Appendix A5, the calibrated value of k_W is 0.05905 m² kg⁻¹ and the greenhouse gas parameters for H₂O are

$$G_{\mathrm{W}\mathrm{1}}=\mathrm{17.89}\text{and}{G}_{\mathrm{W}\mathrm{2}}=\mathrm{12.05}.$$

Figure 3 Climate EBM (Eqs. 6–6) including greenhouse effect of water vapour only, μ=0 (with F_A=45, F_O=50, $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, α_W=0.08, and Z=9 km). (a) Absorptivity η_W(τ_S) as given by Eq. (25) with relative humidity $\mathit{\delta }=\mathrm{0.4},\mathrm{0.6}$ and 0.8 from bottom to top blue lines. Below freezing (τ_s<1) water vapour has little influence as a greenhouse gas, but absorptivity increases rapidly to approach η=1 for τ>1. (b) Bifurcation diagram for EBM with μ=0, and relative humidity increasing from 0 to 1. Note the accelerated increase in surface temperature τ_S above τ_S=1, compared to that in Fig. 2b, due to the rapid increase in η_W.

Figure 3a shows that the function η_W(τ_S) in Eq. (25) increases rapidly from near 0 to near 1 as τ_S increases past 1, and it steepens towards a step function as δ increases. This implies that if the surface temperature τ_S in the surface energy balance Eq. (6) increases past τ_S=1, then the absorptivity of water vapour η_W(τ_S), acting as a greenhouse gas in Eq. (6), increases rapidly, thus amplifying the heating of the atmosphere by the radiation $I_{\mathrm{S}}=\mathit{\sigma }{T}_{\mathrm{S}}^{\mathrm{4}}$ from the surface. Energy balance requires a corresponding increase in the radiation I_A transmitted from the atmosphere back to the surface, further increasing the surface temperature τ_S. This is a positive feedback loop called water vapour feedback. This positive feedback is manifested in Fig. 3b in the additional increase in T_S on the warm branch, beyond that due to the ice–albedo bifurcation. Compare this with Fig. 2, where no water vapour is present. This rapid nonlinear change is due to the relatively large size of the greenhouse constant $G_{\mathrm{W}\mathrm{1}}=\frac{L_{\mathrm{v}}}{R_{\mathrm{W}}{T}_{\mathrm{R}}}=\mathrm{17.89}$ in the exponent of the Clausius–Clapeyron Eq. (18), as it reappears in Eq. (25).

2.3.4 Combined CO₂ and H₂O greenhouse gases

The combined effect of two greenhouse gases is determined by the Beer–Lambert law as shown in Sect. 2.2. If η_C is the absorptivity of CO₂ as in Eq. (14) and η_W is the absorptivity of H₂O as in Eq. (19), then the combined absorptivity η of these two is obtained by adding the two corresponding optical depths λ_C and λ_W. The overall absorptivity of the atmosphere including these two greenhouse gases and the (assumed) constant effect of clouds is, by Eq. (10),

$$\begin{array}{ll}{\displaystyle }\mathit{\eta }\left({\mathit{\tau }}_{\mathrm{S}}\right)& {\displaystyle }=\mathrm{1}-(\mathrm{1}-{\mathit{\eta }}_C)(\mathrm{1}-{\mathit{\eta }}_{\mathrm{W}})(\mathrm{1}-{\mathit{\eta }}_{\mathrm{Cl}})\\ {\displaystyle }& {\displaystyle }=\mathrm{1}-\exp \left[-{\mathit{\lambda }}_C-{\mathit{\lambda }}_W\left({\mathit{\tau }}_{\mathrm{S}}\right)\right](\mathrm{1}-{\mathit{\eta }}_{\mathrm{Cl}})\\ {\displaystyle }& {\displaystyle }=\mathrm{1}-\exp [-\mathit{\mu }\cdot G_C-\mathit{\delta }{G}_{\mathrm{W}\mathrm{2}}\\ \text{(35)}& {\displaystyle }& {\displaystyle }\underset{{\mathit{\tau }}_{\mathrm{S}}-\mathit{\gamma }Z}{\overset{{\mathit{\tau }}_{\mathrm{S}}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\exp \left(G_{\mathrm{W}\mathrm{1}}\left[{\displaystyle \frac{\mathit{\tau }-\mathrm{1}}{\mathit{\tau }}}\right]\right)\mathrm{d}\mathit{\tau }\left]\right(\mathrm{1}-{\mathit{\eta }}_{\mathrm{Cl}}).\end{array}$$

The full nondimensional two-layer EBM is therefore specified by Eqs. (6)–(6) and (25), and the parameter values in Table 1.

2.4 Positive feedback mechanisms

The above analysis shows that there are two highly nonlinear positive feedback mechanisms in this EBM. Both serve to amplify an increase (or decrease) in surface temperature T_S near the freezing point as follows. Consider the case of rising temperature. The first feedback is ice–albedo feedback, due to a change in albedo in the surface equilibrium equation. As illustrated in Fig. A1, if the surface temperature increases slowly through the freezing point, causing a drop in albedo, there is a large increase in solar energy absorbed, leading to an abrupt increase in temperature. The second is water vapour feedback in the atmosphere equilibrium equation. As shown in Fig. 3, if the surface temperature continues to increase above freezing, then the absorptivity of water vapour increases dramatically, strengthening the greenhouse effect for water vapour and further increasing the temperature. Both of these mechanisms act independently of the concentration of CO₂ itself in the atmosphere. However, if the concentration of CO₂ goes up, causing a rise in temperature, then each of these two positive feedback mechanisms can amplify the increase in temperature that would occur due to CO₂ alone.

3 Applications of the energy balance model

The goal of this section is an exploration of the underlying causes of abrupt climate changes that have occurred on Earth in the past 100 million years, using, as a tool, the EBM developed in Sect. 2. The most dramatic climate changes occurred in the two polar regions of the Earth. The climate of the tropical region of the Earth has changed relatively little in the past 100 million years. In this section, the EBMs for the two polar regions are applied to the Arctic Pliocene paradox, the glaciation of Antarctica, the warm, equable Cretaceous problem and the warm, equable Eocene problem. For comparison, a control EBM with parameter values set to those of the tropics shows no abrupt climate changes in the tropical climate for the past 100 million years.

3.1 EBM for the Pliocene paradox

The Arctic region of the Earth's surface has had ice cover year-round for only the past few million years. For at least 100 million years prior to about 3 Ma, the Arctic had no permanent ice cover, although there could have been seasonal snow in the winter. Recently, investigators have found plant and animal remains, in particular on the farthest northern islands of the Canadian Arctic Archipelago, which demonstrate that there was a wet temperate rainforest there for millions of years, similar to that now present on the Pacific Northwest coast of North America (Basinger et al., 1994; Greenwood et al., 2010; Struzik, 2015; West et al., 2015; Wolfe et al., 2017). The relative humidity has been estimated at 67 % (Jahren and Sternberg, 2003), and this value has been chosen for δ in the Arctic EBM of this section. The change from an ice-free to an ice-covered condition in the Arctic occurred abruptly, during the late Pliocene and early Pleistocene. This is sometimes called the Pliocene–Pleistocene transition (PPT 3.0–2.5 Ma; Willeit et al., 2015; Tan et al., 2018). It has been a longstanding challenge to explain this dramatic change in the climate; however, significant progress has been made recently for the case of the Greenland ice sheet (GIS). The authors of Willeit et al. (2015) have solved the “inverse problem” for the GIS by finding a schematic pCO₂ concentration scenario from an ensemble of transient simulations using an EMIC that gives the best fit to data from 3.2 to 2.4 Ma, taking into account the obliquity cycle. Meanwhile, Tan et al. (2018) have used an AOGCM asynchronously coupled with sophisticated ice sheet models to reproduce the waxing and waning of the GIS across the PPT, obtaining good qualitative agreement with ice rafted debris data reconstructions.

Currently, there is great interest in the mid-Pliocene climate, because it is the most recent paleoclimate that resembles the future warmer climate now predicted for the Earth. Significant progress in understanding the Pliocene climate has been achieved in recent years. The Pliocene Research, Interpretation and Synoptic Mapping (PRISM) project of the US Geological Survey has contributed to this goal (Dowsett et al., 2011, 2013, 2016), as has the Pliocene Model Intercomparison Project (PlioMIP; Haywood et al., 2011, 2016; Zhang et al., 2013). Advances in the extraction and interpretation of proxy data have given a clearer picture of the warm Pliocene climate (Haywood et al., 2009; Salzmann et al., 2009; Ballantyne et al., 2010; Steph et al., 2010; Seki et al., 2010; Bartoli et al., 2011; De Schepper et al., 2013; Knies et al., 2014; O'Brien et al., 2014; Brierley et al., 2015; Fletcher et al., 2018). At the same time, computer models have achieved closer agreement with proxy data, see for example Haywood et al. (2009); Dowsett et al. (2011); Zhang and Yan (2012); Sun et al. (2013); Willeit et al. (2015); Tan et al. (2017, 2018); Chandan and Peltier (2017, 2018), and other references therein. Using a fully coupled atmosphere–ocean GCM, Lunt et al. (2008) considered the following forcing factors contributing to late Pliocene glaciation: decreasing carbon dioxide concentration, closure of the Panama seaway, end of a permanent El Niño state, tectonic uplift and changing orbital parameters. They concluded that falling CO₂ levels were primarily responsible for the formation of the Greenland ice sheet in the late Pliocene. Recently, Tan et al. (2018) have strengthened this conclusion.

During the Pliocene Epoch, important forcing factors that determine climate were very similar to those of today. The Earth orbital parameters, the CO₂ concentration, solar radiation intensity, position of the continents, ocean currents and atmospheric circulation all had values close to the values they have today. Yet, in the early/mid-Pliocene, 3.5–5 million years ago, the Arctic climate was much milder than that of today. Arctic surface temperatures were 8–19 ^∘C warmer than today and global sea levels were 15–20 m higher than today, and yet CO₂ levels are estimated to have been 340–400 ppm, about the same as 20th century values (Ballantyne et al., 2010; Csank et al., 2011; Tedford and Harington, 2003). As mentioned in the Introduction section, the problem of explaining how such dramatically different climates could exist with such similar forcing parameter values has been called the Pliocene paradox (Cronin, 2010; Fedorov et al., 2006, 2010; Zhang and Yan, 2012).

Another interesting fact concerning polar glaciation is that although both poles have transitioned abruptly from ice-free to ice-covered conditions, they did so at very different geological times. The climate forcing conditions of Earth are highly symmetric between the two hemispheres, and for most of the past 200 million years (or more) the climates of the two poles have been similar. However, there was an anomalous interval of about 30 million years, from the EOT, 34 Ma, to the early Pliocene, 4 Ma, when the Antarctic was largely ice covered but the Arctic was largely land ice free. Because CO₂ disperses rapidly in the Atmosphere, its concentration μ must be the same everywhere at any given time. Therefore we seek a forcing factor other than μ to account for this 30-million-year period of broken symmetry. One obvious difference is geography. Since the Eocene, the South Pole has been land-locked in Antarctica, while the North Pole has been in the Arctic Ocean. Therefore, our two EBMs for the North Pole and South Pole have very different values for ocean heat transport F_O. We will show that this difference is sufficient to account for the gap of 30 million years between the Antarctic and Arctic glaciations.

Our Pliocene Arctic EBM brackets the Pliocene Epoch between the mid-Eocene (50 Ma) and pre-industrial modern times, and it models the effects of the slow decrease in both CO₂ and ocean heat transport F_O in the Arctic over this long time interval. In this Arctic model, abrupt glaciation of the Arctic is inevitable, due to the existence of a bifurcation point.

During the Eocene (56–34 Ma), temperatures were much higher than today, especially in the Arctic (Greenwood et al., 2010; Wolfe et al., 2017; Huber and Caballero, 2011); also, CO₂ concentration μ was higher than today. Estimates of Eocene CO₂ concentration μ vary from 1000–1500 ppm (Pagani et al., 2005, 2006) to 490 ppm (Wolfe et al., 2017). For this EBM, we set mid-Eocene CO₂ at μ=1000 (Pagani et al., 2005). Both temperature and CO₂ concentration have decreased steadily but not monotonically, with many fluctuations from their Eocene values to pre-industrial modern values. The overall decrease in CO₂ concentration observed since the Eocene may be attributed to decreased volcanic activity, increased absorption and sequestration by vegetation and the oceans, continental erosion, and other sinks.

The changes in ocean heat transport to the Arctic are more complicated and derive from many factors only summarized here (De Schepper et al., 2015; Haug et al., 2004; Knies et al., 2014; Lunt et al., 2008; Zhang et al., 2013). There was a slow drop in global sea level, in large part due to the gradual accumulation of vast amounts of water in the form of ice and snow on Antarctica. It has been estimated that the total amount of ice today in the Antarctic is equivalent to a change in sea level of about 58 m (Fretwell et al., 2013; IPCC, 2013). This drop in sea level likely reduced the flow of warm tropical ocean water into the Arctic. Against this background, several other factors came into play due to changing geography. In the Eocene, the North Atlantic Ocean was not always connected to the Arctic, but the Turgai Sea existed between Europe and Asia and connected the warm Indian Ocean to the Arctic, until about 29 Ma (Briggs, 1987). By the Oligocene, the Turgai Sea had closed and the North Atlantic had opened between Greenland and Norway, forming a deep-water connection to the Arctic Ocean. The Bering Strait opened and closed. During the Pliocene the formation of the Isthmus of Panama about 3.5 Ma cut off a warm equatorial current that had existed between the Atlantic and Pacific, at least since Cretaceous times. On the Atlantic side of the isthmus, the sea water became warmer, and became more saline due to evaporation. The Gulf Stream carried this warm salty water to western Europe. One might expect the Gulf Stream to transport more heat into the Arctic. However, some believe that the Gulf Stream actually contributed to glaciation in the Arctic (Haug et al., 2004; Bartoli et al., 2005) as follows. Evaporation from the Gulf Stream waters contributed to rainfall across northern Europe and Siberia, increasing the flow of fresh water in rivers emptying into the Arctic Ocean. This reduced the salinity, and hence the density, of the Arctic Ocean waters. The Gulf Stream waters in the North Atlantic now cooler, and denser due to high salinity, were forced downward by the less dense Arctic waters, which began to flow into the North Atlantic on the surface. The denser Gulf Stream waters returned southward as a deep ocean current, without having conveyed much heat to the Arctic. Meanwhile the low salinity Arctic surface water, with a higher freezing temperature, began to freeze, resulting in higher albedo and accelerating Arctic glaciation. In large measure, these changing geographical factors partially cancelled each other in their contributions to ocean heat transport to the Arctic. In the EBM, we summarize all of the above heat transport mechanisms by specifying a slow overall decrease in ocean heat transport to the Arctic, represented by the single forcing parameter F_O.

Figure 4a shows graphs of the surface and atmosphere equilibrium equations for the Arctic Pliocene EBM, for varying values of μ. The figure shows only one surface equilibrium curve (red), with $F_{\mathrm{O}}=\mathrm{50}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, because the change in F_O is relatively small. The blue atmosphere equilibrium curves represent values of CO₂ concentration μ falling from 1100 to 200 ppm, from the top to bottom blue curves. It is clear that there may exist up to three points of intersection of a given atmosphere equilibrium curve (blue) with the surface equilibrium curve (red); namely, a warm equilibrium state τ_s>1, a frozen equilibrium state τ_S<1 and a third (intermediate) solution, which is always unstable (when it exists). As μ decreases and the warm equilibrium state τ_s>1 approaches the local minimum on the red S-shaped curve, the unstable intermediate equilibrium state moves down the middle branch of the red S-shaped curve. When they meet, these two equilibria coalesce, then disappear, via a saddle-node bifurcation. Beyond this saddle-node bifurcation, only one equilibrium state remains, which is the stable frozen state on the left in the Figure. Dynamical systems theory tells us that following this bifurcation, the system will transition rapidly to that frozen equilibrium state. The paleoclimate record shows that CO₂ concentration was trending downward for millions of years before and during the Pliocene. Therefore, Fig. 4a predicts that an abrupt drop in temperature to a frozen state would be inevitable if this trend continued far enough.

Figure 4 Pliocene Arctic EBM Eqs. (6)–(6). Parameter values δ=0.67, $F_{\mathrm{A}}=\mathrm{45}{\mathrm{Wm}}^{-\mathrm{2}}$, $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, α_W=0.08 and Z=9 km. (a) CO₂ takes values μ=200, 500, 800, 1100 and 1400 ppm from bottom to top on the blue curves with fixed $F_{\mathrm{O}}=\mathrm{50}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ on the red curve. The warm equilibrium state τ_S>1 disappears as μ decreases on successively lower blue curves. (b) Bifurcation diagram for the Pliocene paradox. Here, CO₂ concentration μ and ocean heat transport F_O decrease simultaneously with increasing ν, ($\mathrm{0}\le \mathit{\nu }\le \mathrm{1}$), as given by Eq. (26). As ν increases, the warm equilibrium solution (τ_S>1) disappears in a saddle-node bifurcation at approximately ν=0.91, corresponding to forcing parameters μ=336 ppm and F_O=50.9. The equilibrium temperatures at this value of ν are 3.9 and −26 ^∘C, and geological time about 4.5 Ma. To the right of this point, only the frozen equilibrium state exists. Between about ν=0.12 and ν=0.91, the frozen and warm equilibrium states coexist, separated by the unstable intermediate state.

In order to explore this downward trend further, we bracket the Pliocene Epoch between the mid-Eocene Epoch (50 Ma) and the pre-industrial modern era (300 years ago), and define a surrogate time variable ν by

$$\begin{array}{}\text{(36)}& t=\mathrm{50}(\mathrm{1}-\mathit{\nu })\mathrm{Ma}.\end{array}$$

As reviewed above, it is believed that ocean heat transport F_O decreased modestly over this time period, mainly due to the drop in global sea level, while the CO₂ concentration μ decreased more significantly. Therefore, we express both μ and F_O as decreasing functions of bifurcation parameter ν

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\mu }=\mathrm{1000}-\mathrm{730}\cdot \mathit{\nu }\mathrm{ppm}\\ \text{(37)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{O}}=\mathrm{60}-\mathrm{10}\cdot \mathit{\nu }{\mathrm{Wm}}^{-\mathrm{2}}.\end{array}$$

Here, ν=0 corresponds to estimated mid-Eocene values of μ and F_O (Pagani et al., 2005; Wolfe et al., 2017; Barron et al., 1981), while ν=1 corresponds to modern pre-industrial values (IPCC, 2013). Equation (26) defines μ and F_O as linear functions of ν. In the real world, neither μ nor F_O decreased linearly. This is not an obstacle for our bifurcation analysis. What is important is that, somewhere between ν=0 and ν=1, a bifurcation point is crossed.

Figure 4b is a bifurcation diagram, which shows the dependence of surface temperature τ_S on this bifurcation parameter ν. Note that for ν=0 (mid-Eocene values, 50 Ma) only the warm equilibrium state exists. At about ν=0.116 (44 Ma) both warm and frozen equilibrium states exist. However, as ν increases toward ν=1, the warm equilibrium state disappears in a saddle-node (or fold) bifurcation, leaving only the frozen equilibrium state to the right of the saddle-node bifurcation. This bifurcation occurs at approximately ν=0.91, which corresponds to $F_{\mathrm{O}}=\mathrm{50.9}{\mathrm{Wm}}^{-\mathrm{2}}$ and μ=336, which is in good agreement with the determination (Seki et al., 2010) that the warm Pliocene pCO₂ was between 330 and 400 ppm, similar to today. The temperatures of the warm state and frozen state, at the bifurcation value of ν=0.91 in the EBM, are +3.9 and −26 ^∘C, respectively. The surrogate time of the bifurcation point, ν=0.91, corresponds to a geological time of t=4.5 Ma, from Eq. (25), which is close to the time of glaciation of the Arctic in the geological record.

The EBM plotted in Fig. 4 provides a plausible explanation for the Pliocene paradox. For millions of years up to the mid-Pliocene, while the Arctic temperature remained above freezing on the warm solution branch in Fig. 4b, the climate change was incremental. Then the slowly acting physical forcings of decreasing CO₂ concentration and decreasing ocean heat transport, F_O, were amplified by the mechanisms of ice–albedo feedback and water vapour feedback, both of which act very strongly when the temperature crosses the freezing point of water. The EBM suggests that when the freezing temperature was approached, the Arctic climate changed abruptly via a saddle-node bifurcation as in Fig. 4b, from a warm state to a frozen state.

Permanent El Niño

Another explanation has been proposed for the Pliocene paradox. There is convincing evidence that at the beginning of the Pliocene, there was a permanent El Niño condition in the tropical Pacific Ocean (Fedorov et al., 2006, 2010; Steph et al., 2010; Cronin, 2010; von der Heydt and Dijkstra, 2011). However, some have disputed this finding (Watanabe et al., 2011). It has been suggested that a permanent El Niño condition could explain the warm early Pliocene, and that the onset of the El Niño–Southern Oscillation (ENSO) was the cause of sudden cooling of the Arctic during the Pliocene. Today, it is known that ENSO can influence weather patterns as far away as the Arctic. However, the present authors propose that bistability and bifurcation provide a more satisfactory explanation for the Pliocene paradox, and suggest that the concurrent change in ENSO may have been a consequence, not the cause, of the changing Pliocene Arctic climate (work in progress).

3.2 EBM for the glaciation of Antarctica

Antarctica in the mid-Cretaceous Period was ice free, and it remained so for the remainder of the Cretaceous Period and into the Paleocene and early Eocene. However, recent investigations in paleoclimate science have shown that there was an abrupt drop in temperature and an onset of glaciation in Antarctica at the EOT about 34 Ma (Katz et al., 2008; Lear et al., 2008; Miller et al., 2008; Scher et al., 2011, 2015; Ladant et al., 2014; Ruddiman, 2014). In the mid-Cretaceous, the continent of Antarctica was in the South Pacific Ocean and the South Pole was located in open ocean, warmed by South Pacific Ocean currents (Cronin, 2010). Then, the continent of Antarctica moved poleward and began to encroach upon the South Pole towards the end of the Cretaceous, fully covering the South Pole before the end of the Eocene (Briggs, 1987). The diminishing marine influence on the South Pole coincided with the onset of cooling of Antarctica about 45 Ma (mid-Eocene). The opening of both the Drake Passage (between South America and Antarctica) and of the Tasmanian Gateway (between Antarctica and Australia) near the end of the Eocene (34 Ma) led to the development of the ACC, which further isolated the South Pole from warm ocean heat transport and accelerated the cooling and glaciation of Antarctica (Cronin, 2010). Therefore, from the early Eocene to the Oligocene, the ocean heat transport to the South Polar region decreased significantly, from near the Cretaceous value, about 100 Wm⁻² (Barron et al., 1981), to a much lower value, estimated here to be 30 Wm⁻². The fact that the glaciation of the Antarctic took place 30 million years before the glaciation of the Arctic is very likely due primarily to the much larger decrease in ocean heat transport that took place at the South Pole.

Because the atmosphere is well mixed, at any given time the CO₂ level in the Antarctic is the same as elsewhere. For this EBM, we estimate the early Eocene global CO₂ level as μ=1100 ppm, (Pagani et al., 2005, 2006; Cronin, 2010); decreasing to approximately modern levels by the end of the Oligocene, 23 Ma, that is to μ=400 ppm.

In the recent literature, there has been a discussion of whether the primary cause of the onset of Antarctic glaciation is the slow decline in CO₂ concentration or the decrease in poleward ocean heat transport due to the opening of ocean gateways and the development of the ACC (DeConto et al., 2008; Goldner et al., 2014; Pagani et al., 2011; Scher et al., 2015). In fact, Ladant et al. (2014) state “The reasons for this greenhouse–icehouse transition have long been debated, mainly between the 'tectonic-oceanic' hypothesis and the CO₂ hypothesis”. In the work by Scher et al. (2015), it is observed that the onset of the ACC coincided with major changes in global ocean circulation, which probably contributed to the drawdown in CO₂ concentration in the atmosphere. Since both CO₂ concentration and ocean heat transport F_O are explicit parameters in the EBM presented here, they can be varied independently in the model to investigate which one played a primary role.

Figure 5 EBM for the glaciation of Antarctica (with δ=0.67, $F_{\mathrm{A}}=\mathrm{45}{\mathrm{Wm}}^{-\mathrm{2}}$, $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, α_W=0.15, Z=9 km). (a) Graphs of EBM equations, with F_O fixed at 90 W m⁻² in the surface equation (red line), and with CO₂ concentration μ in the atmosphere equation increasing (from bottom to top blue lines) $\mathit{\mu }=\mathrm{200},\mathrm{400},\mathrm{600},\mathrm{800},\mathrm{1000}$ and 1200 ppm. (b) Graphs of the same EBM equations but with CO₂ concentration μ fixed at 1100 ppm (blue line) and with ocean heat transport F_O decreasing (from bottom to top red lines) 80, 60, 40, 20 and 0 W m⁻².

Figure 5 shows the solutions of the EBM equations for Antarctic parameter values, including δ=0.67, $F_{\mathrm{A}}=\mathrm{45}{\mathrm{Wm}}^{-\mathrm{2}}$, $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, α_W=0.15, and Z=9 km. Points of intersection of one red line (surface equation) and one blue line (atmosphere equation) are the equilibrium solutions of the system. As noted earlier, there can be up to three equilibrium solutions. In Fig. 5a the ocean heat transport F_O is held fixed, while CO₂ concentration μ is decreased from the top to bottom blue lines. In Fig. 5b the reverse is true; the CO₂ concentration μ is held fixed while the ocean heat transport F_O is decreased from bottom to top red lines. For a sufficiently small value of F_O, the warm solution disappears in Fig. 5b. This occurs when the warm solution meets the unstable intermediate solution in a saddle-node bifurcation. Beyond this bifurcation, only the frozen solution (τ_S<1) exists. In contrast, reducing μ in Fig. 5a does not affect the existence of the warm equilibrium, only the existence of the cold equilibrium.

Figure 6 displays bifurcation diagrams in which the surface temperature τ_S, which has been determined as a solution of the EBM, is plotted while parameters in the EBM are varied smoothly, from the early Eocene values, 55 million years ago, to the late Oligocene values, 23 million years ago. These dates are chosen to bracket the time of the glaciation of Antarctica. First we introduce a bifurcation parameter ν, which acts as a surrogate time variable. The bifurcation parameter ν is related to geological time by

$$\begin{array}{}\text{(38)}& t=\mathrm{55}-\mathrm{32}\cdot \mathit{\nu }\mathrm{Ma}.\end{array}$$

Thus, ν=0 corresponds to the early Eocene, 55 Ma, and ν=1 corresponds to the late Oligocene, 23 Ma. During this time, it is known that both CO₂ concentration μ and ocean heat transport F_O were decreasing. To study the effects of the simultaneous reduction of μ and F_O, we express both as simple linear functions of the bifurcation parameter ν, as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\mu }=\mathrm{1100}-\mathrm{700}\cdot \mathit{\nu },\\ \text{(39)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{O}}=\mathrm{100}-\mathrm{70}\cdot \mathit{\nu },\mathrm{0}\le \mathit{\nu }\le \mathrm{1}.\end{array}$$

Here, ν=0 corresponds to the high early Eocene values of μ=1100 ppm and $F_{\mathrm{O}}=\mathrm{100}{\mathrm{Wm}}^{-\mathrm{2}}$, while ν=1 corresponds to the low Oligocene–Miocene boundary values of μ=400 ppm and $F_{\mathrm{O}}=\mathrm{30}{\mathrm{Wm}}^{-\mathrm{2}}$. Thus, in Eq. (27), as ν increases from 0 to 1, the climate forcing factors μ and F_O decrease linearly from their Eocene values to late Oligocene values. Strictly speaking, these forcings did not change linearly in time; however, the important fact is that overall they were decreasing, For our bifurcation results the decrease does not need to be linear in the EBM. The atmospheric heat transport parameter F_A is held constant at $F_{\mathrm{A}}=\mathrm{45}{\mathrm{Wm}}^{-\mathrm{2}}$. The solar radiation is $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, the warm albedo value is α_W=0.15 and the tropopause height is Z=9 km. Other parameters are as in Table 1.

Figure 6 Bifurcation diagrams for the glaciation of Antarctica (with δ=0.67, F_A=45, $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$, α_W=0.15, Z=9 km). (a) This shows a saddle-node bifurcation at ν=0.779, corresponding to forcing parameter values μ=555 ppm and $F_{\mathrm{O}}=\mathrm{45.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in Eq. (27). As ν increases through this saddle-node bifurcation point, the climate will transition abruptly from a warm state to a frozen state. The warm state and frozen state temperatures coexisting at the bifurcation point are +3.5 and −22.1 ^∘C. (b) Two superimposed bifurcation diagrams as in (a), except that, on the green curve, F_O is held fixed at its Eocene values while μ deceases as in Eq. (27), and, on the magenta curve, μ is held fixed at the Eocene value while F_O decreases as in Eq. (27). No saddle-node bifurcation destroying the warm equilibrium occurs in either scenario.

Figure 6a shows a saddle-node bifurcation at ν=0.779, corresponding to forcing parameter values μ=555 ppm and $F_{\mathrm{O}}=\mathrm{45.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in the model, and to geological time about 30 Ma, assuming the linear time relation in Eq. (26). This corresponds closely to best estimates of the timing of glaciation in the geological record of about 34 Ma at the EOT. This quantitative agreement should not be taken very seriously, as it is largely due to assumptions made in the modelling. However, the existence of the bifurcation, implying an abrupt glaciation of Antarctica, can be taken seriously, as this is a robust (structurally stable) property of the model. The warm state and frozen state temperatures coexisting at the bifurcation point are +3.5 and −22.1 ^∘C. As ν increases past the bifurcation point ν=0.779, the warm climate state ceases to exist, and the climate system transitions (“falls”) rapidly from the saddle-node point to the frozen state. Thus the EBM reproduces the abrupt transition to the glaciation of Antarctica, which is seen in the geological record at the EOT.

Figure 6b explores the relative importance of decreasing CO₂ concentration μ and decreasing ocean heat transport F_O in the glaciation of Antarctica, a subject that has been much debated in the literature (DeConto et al., 2008; Goldner et al., 2014; Ladant et al., 2014; Miller et al., 2008; Pagani et al., 2011; Scher et al., 2011, 2015). The green curve represents a scenario in which μ decreases as in Eq. (27) but F_O is held fixed at its Eocene value, and the magenta curve represents a case in which μ is held fixed at its Eocene value while F_O decreases according to Eq. (27). In neither case does a glaciation event occur. The analysis of this paper implies that significant decreases in both CO₂ concentration μ and ocean heat transport F_O are required to achieve a saddle-node bifurcation, reproducing the observed transition to a frozen Antarctic state at the EOT.

While the glaciation of Antarctica is an accepted fact in paleoclimate science (Miller et al., 2008; Lear et al., 2008; Pagani et al., 2011; Scher et al., 2011; Goldner et al., 2014; Scher et al., 2015), the suddenness of the climate change that occurred in Antarctica near the Eocene–Oligocene transition (34 Ma), a time when the forcing parameters were changing slowly, has been difficult to explain. The bifurcation analysis presented here gives a simple but plausible explanation for the abruptness of this event. Furthermore, this EBM supports the hypothesis that both falling CO₂ concentration μ and decreasing ocean heat transport F_O (due to gateway openings and development of the ACC) are essential to an explanation of the sudden glaciation of Antarctica at the EOT.

3.3  EBM for the tropics

In the tropics, many of the values of the forcing parameters are different from their values in the Arctic and Antarctic, see Table 2. The geological record shows little change in the tropical climate over the past 100 million years, other than minor cooling. Even when the Arctic climate changed dramatically in the late Pliocene, the tropical climate changed very little.

Figure 7a shows solutions of the EBM equations for tropical parameter values. Note that at no value of the forcings used here does the warm equilibrium point approach a saddle-node bifurcation point. Thus, our EBM is in agreement with the geological record.

Figure 7 Tropics EBM, using Eqs. (6)–(6) and parameter values as in Table 2. (a) The blue curves represent the atmosphere equilibrium equation for $\mathit{\mu }=\mathrm{300},\mathrm{600},\mathrm{900}$ and 1200 from bottom to top. (b) Bifurcation diagram for the tropics EBM. Here, CO₂ concentration μ and ocean heat transport F_O are both decreasing, with increasing ν,  ($\mathrm{0}\le \mathit{\nu }\le \mathrm{1}$), as given by Eq. (27). As ν varies, both the warm equilibrium state (τ_S>1) and the frozen equilibrium state (τ_S<1) persist; there are no saddle-node bifurcations.

Table 2 Summary of parameters used in the tropics EBM. Relative humidity δ is higher in the tropics than the Arctic, and the forcings F_O and F_A are negative instead of positive. The insolation Q is as determined by McGehee and Lehman (2012). Tropopause height Z is from Kishore et al. (2006). Parameters not listed here are as in Table 1.

A bifurcation diagram is constructed for the tropics EBM, spanning mid-Cretaceous to modern pre-industrial times, similar to those for the Antarctic and Arctic glaciations. Here ν=0 corresponds to mid-Cretaceous values and ν=1 corresponds to modern pre-industrial values. We let μ and F_O decrease linearly with ν.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\mu }=\mathrm{1130}-\mathrm{860}\cdot \mathit{\nu },\\ \text{(40)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{O}}=-\mathrm{55}+\mathrm{29}\cdot \mathit{\nu },\mathrm{0}\le \mathit{\nu }\le \mathrm{1}.\end{array}$$

The mid-Cretaceous CO₂ concentration μ=1130 ppm is as determined by Fletcher et al. (2008), see Table 3. The ocean heat transport $F_{\mathrm{O}}=-\mathrm{55}{\mathrm{Wm}}^{-\mathrm{2}}$, from the tropics in the mid-Cretaceous, is from Barron et al. (1981). Astrophysicists have determined that solar luminosity is slowly increasing with time (Sagan and Mullen, 1972). For the mid-Cretaceous, insolation was approximately 1 % less than it is today (Barron, 1983). This difference is considered too small to be significant in our model. The bifurcation diagram for the tropics EBM is shown in Fig. 7b. Note that no bifurcations occur for parameter values relevant to the tropics. This is in agreement with paleoclimate records that show little change in tropical climate, even when polar climates changed dramatically.

3.4 EBM for the Cretaceous and Eocene warm, equable climate problems

One hundred million years ago, in the mid-Cretaceous Period, the climate of Earth was much more equable than today. “More equable” means that the pole-to-Equator temperature gradient was much smaller, and also the seasonal summer/winter temperature variations were much smaller. The climate in the tropics was only slightly warmer than today, but the climate at both poles was much warmer than today. An abundance of plant and animal life thrived under these conditions, from the Equator to both poles, including of course dinosaurs. The question of how this globally ice-free climate could have been maintained was explored in pioneering work of Barron and coworkers (Barron, 1983; Barron et al., 1981, 1995; Sloan and Barron, 1992). Barron called this the warm, equable Cretaceous climate problem. Early GCMs of the 1980s adjusted to mid-Cretaceous forcing parameter values, had difficulty giving good agreement with climate proxies. In order to obtain polar temperatures in agreement with mid-Cretaceous values, these simulations typically assumed increased CO₂ levels more than 4 times modern levels; but then the tropical temperatures predicted by the models were too high. Today, the issue of reproducing an equable Equator-to-pole temperature gradient remains a challenge for modellers.

There have been many investigations of the correlation between early climate models and proxy data for the Cretaceous climate (Pagani et al., 2005; Bice et al., 2006; Donnnadieu et al., 2006; Cronin, 2010; Ruddiman, 2014). Recent studies have succeeded in narrowing the gap (Craggs et al., 2012; Bowman et al., 2014; O'Brien et al., 2017; Lunt et al., 2016; Ladant and Donnadieu, 2016). Based on ocean drilling samples, Bice et al. (2006) estimate Cretaceous CO₂ concentrations between 600 and 2400 ppmv, and tropical Atlantic upper ocean temperatures between 33 and 42 ^∘C. Based on fossil samples, Fletcher et al. (2008) estimate mid-Cretaceous atmospheric CO₂ concentrations of 1130 ppmv and we choose their estimate for use in this EBM, see Table 3.

Fletcher et al. (2008)Jahren and Sternberg (2003)Barron et al. (1981)Barron et al. (1981)

Table 3 Summary of parameters used in the polar mid-Cretaceous model. The main difference in Cretaceous climate between the two poles was that ocean heat transport was much higher in the Antarctic. Parameters not listed here are the same as in Table 1. The values of F_A and F_O shown here are estimated from Barron et al. (1981).

Paleoclimate data from the mid-Cretaceous show little difference in climate between the two warm poles at that time (Barron, 1983; Barron et al., 1981). The major difference in forcings between the two poles at that time was that ocean heat transport F_O was much higher in the Antarctic than the Arctic, as shown in Table 3. This is due to the location of the South Pole in open ocean during the Cretaceous, as the continent of Antarctica had not yet drifted to its present position over the South Pole. Estimates for atmospheric transport are also different (Barron et al., 1981). Therefore, we model the mid-Cretaceous climate at the two poles using Eqs. (6)–(6), with the values from Table 3. The EBM equilibrium curves are shown in Fig. 8, where it is clear that a warm (τ_S>0) and a frozen (τ_S<0) equilibrium state exists at both poles for the given parameter values. The Antarctic warm equilibrium state is warmer than in the Arctic because of the higher ocean heat transport F_O to the Antarctic. From Sect. 3.3, the tropical region of Earth also had a warm equilibrium state under mid-Cretaceous conditions. Therefore, the EBM of this paper implies the existence of a pole-to-pole warm, equable Cretaceous climate, as is seen in the geological record.

Figure 8 also implies that a frozen equilibrium state is mathematically possible at each pole during the Cretaceous. In that case, the tropics could remain in the warm state, thus giving the mathematical possibility of a Cretaceous climate that is not equable, but has warm tropics and ice-covered poles, like today's climate. This may help explain why some computer simulations, originally designed for today's climate conditions, found a mathematically existing Cretaceous climate that resembled today's climate more than the warm, equable climate that physically existed on Earth in the Cretaceous.

The climate of the Eocene Epoch (≈55 million years ago) was the warmest of the past 65 million years, but a little cooler than the Cretaceous (Sloan and Barron, 1992; Pagani et al., 2005; Cronin, 2010; Huber and Caballero, 2011; Hutchinson et al., 2018). Both poles were ice-free and the pole-to-Equator temperature gradient was much smaller than today. As for the mid-Cretaceous, early computer simulations, based on modern climate conditions, had difficulty in reproducing the early Eocene warm, equable climate (Sloan and Barron, 1990, 1992; Sloan and Rea, 1995; Jahren and Sternberg, 2003; Huber and Caballero, 2011). This discrepancy has been called the early Eocene warm, equable climate problem. Here we combine these two “problems” under one name as the warm, equable Paleoclimate problem.

The main difference in forcings between early Eocene conditions and those of the mid-Cretaceous is that the global CO₂ concentration μ may have been a little lower and the ocean heat transport F_O to the Antarctic was less. Referring to Fig. 8, this means that the blue and green atmosphere equilibrium curves move downward slightly, and the magenta Antarctic surface equilibrium curve moves upward slightly. The orange Arctic surface equilibrium curve does not change significantly. With these small changes, all of the equilibrium climate states (intersection points) persist in Fig. 8. The figure for the Eocene EBM (not shown here) is topologically the same as Fig. 8 for the mid-Cretaceous EBM,

Thus, the EBM predicts that a warm ice-free equilibrium state exists (mathematically) at both poles in the mid-Cretaceous and early Eocene. From Sect. 3.3, the equatorial region also has a warm equilibrium state for Eocene conditions. Therefore, this EBM study supports the existence of a pole-to-pole warm, equable climate in both the mid-Cretaceous and the early Eocene. Additionally, in the EBM there is the coexistence of a nonequable climate, with cold poles, under exactly the same forcing conditions, which suggests the following plausible solution to the “warm, equable paleoclimate problem” in both the mid-Cretaceous (Barron, 1983; Barron et al., 1981, 1995; Sloan and Barron, 1990) and the early Eocene (Sloan and Barron, 1990, 1992; Sloan and Rea, 1995; Jahren and Sternberg, 2003; Huber and Caballero, 2011). While the Earth's climate existed in the warm, equable climate state of the EBM, computer simulations of Barron and others may have correctly computed the coexisting nonequable solution.

Figure 8 Solutions of the mid-Cretaceous EBM for both poles, using parameter values in Table 3. The blue and orange curves represent the EBM for the Arctic while the green and magenta curves represent the Antarctic. The two Z-shaped curves represent the surface equilibrium equations for the Arctic (upper orange curve) and the Antarctic (lower magenta curve). The only differences in these two are that ocean heat transport F_O is much greater in the Antarctic than the Arctic, and the atmospheric transport, F_A, is smaller. Both poles support both warm and frozen equilibrium paleoclimate states. The warm equilibrium states at the two poles are marked by stars.

4 Conclusions and future work

This paper presents a new energy balance model (EBM) for the climate of Earth, one that elucidates the distinctive roles of carbon dioxide and water vapour as greenhouse gases, and also the role of ice–albedo feedback, in climate change. Nonlinearity of the EBM leads to multiple solutions of the mathematical equations and to bifurcations that represent transitions between coexisting stable equilibrium states. This EBM sheds new light on several important problems of paleoclimate science; namely, the Pliocene paradox, the abrupt glaciation of Antarctica, and the warm, equable mid-Cretaceous and early Eocene climate problems. Predictions of the EBM are in qualitative agreement with the paleoclimate record.

There has been a range of values of paleoclimate forcing parameters in the literature, in particular for mid-Cretaceous and Eocene CO₂ concentrations and temperatures. The specific choices made here, while informed by proxies, are somewhat arbitrary. However, this fact does not affect the validity of our main conclusions. The coexistence of multiple solutions and the bifurcations demonstrated in the EBM are robust phenomena. That is, the existence of these multiple solutions and bifurcations will persist over a range of values of the forcing parameters. The main conclusion of this paper that sudden and significant changes in climate have occurred, even while forcing parameters were changing very gradually, follows from the existence of mathematical bifurcations in the EBM, not from particular choices of the forcing parameters.

As the paleoclimate record becomes clearer, there is growing evidence for small but rapid fluctuations in some parameter values, including CO₂ concentrations, over the geological time period studied here (Cronin, 2010; Pagani et al., 1999, 2005). The model of this paper assumes a smooth decline in CO₂ concentration. However, this does not invalidate our main conclusions. The theory of “stochastic bifurcation” (Namachchivaya, 1990; Arnold et al., 1996) tells us that if stochastic noise is added parametrically to a deterministic bifurcation problem, then typically the location of the bifurcation (in terms of the bifurcation parameter) may change; but the existence of the bifurcation is preserved.

Further work on this EBM is in progress. Having demonstrated the applicability of the EBM to known paleoclimate transitions, this EBM is now being applied to anthropogenic climate change with the goal of predicting the future climate effects of continued increases in CO₂ concentration.

The most significant question is whether a bifurcation is to be expected in our future. This bifurcation would lead to climate warming even stronger than that currently underway, leading possibly to a warm, equable climate, at least in the Northern Hemisphere, similar to the climate of Earth before the Pliocene.

The equilibrium climate sensitivity (ECS) of the EBM, adapted to present-day satellite data, has been calculated and is presented in Appendix Appendix B of this paper.

This scalar EBM will be generalized to a two-point boundary value problem (BVP) in the altitude variable z, using the Schwarzschild equations to replace the ICAO International Standard Atmosphere (ISA) approximation for lapse rate. That change will more accurately model the behaviour of the greenhouse gases in the atmosphere. Next, that one-dimensional EBM BVP will be incorporated into a generalization of the spherical shell PDE model of Lewis and Langford (2008) and Langford and Lewis (2009). We anticipate that this 3-D zonally symmetric Navier–Stokes Boussinesq model will confirm the fundamental predictions relating to bistability and bifurcation of the present simple EBM. Also, it will enable the study of a third positive feedback mechanism (in addition to the two studied in this paper); namely, Hadley cell convection feedback, which influences atmospheric heat transport F_A.

Appendix A: Empirical calibration of EBM parameters

Some of the parameters used in the EBM have standard values, available in textbooks and reported in Table 1. Others, however, are determined empirically in this appendix. Modern values of these parameters may be determined from today's abundant satellite and land-based data. It is assumed that the values determined here are also valid for paleoclimates.

The primary source for these empirical calibrations are the data presented by Wild et al. (2013, Fig. 4) and Trenberth et al. (2009, Fig. 3). The data from these authors are globally averaged values. We are applying our model to specific regions (the Arctic, Antarctic or tropics) and therefore adjust some of these values as discussed below. We also employ the model with globally averaged values to give an equilibrium climate sensitivity (ECS) calculation in Appendix Appendix B.

A1 Solar radiation

The value of the annually averaged solar radiation at either pole is $Q=\mathrm{173.2}{\mathrm{Wm}}^{-\mathrm{2}}$ and at the Equator is 418.8 Wm⁻² (McGehee and Lehman, 2012; Kaper and Engler, 2013). Values from Fig. 4 of Wild et al. (2013) indicate the globally averaged solar radiation is $Q=\mathrm{340}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. Of this, 79 W m⁻² is directly absorbed by the atmosphere, and 76 W m⁻² is reflected by the atmosphere back into space (both primarily due to clouds). Hence we define the two solar radiation fractions

$$\begin{array}{}\text{(A1)}& {\mathit{\xi }}_{\mathrm{A}}={\displaystyle \frac{\mathrm{79}}{\mathrm{340}}}=\mathrm{0.2324}\text{and}{\mathit{\xi }}_{\mathrm{R}}={\displaystyle \frac{\mathrm{76}}{\mathrm{340}}}=\mathrm{0.2235}.\end{array}$$

Since we incorporate no information on varying cloud cover in the model, and since there is very little such information for paleoclimates, we assume these values are constant around the globe. Clearly cloud cover is correlated with surface temperature. We also tried making ξ_A and ξ_R vary with T_S in a manner similar to how we model heat conduction/convection, F_C below, but found no qualitative differences in the results. For this reason, and to keep the model simple, we have kept these values of ξ_A and ξ_R as global constants in our model.

A2 Surface albedo

The Earth's albedo, α, varies considerably depending on the surface features. Typical values 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. In previous papers, including Dortmans et al. (2018), the polar albedo α is assumed to jump between two discrete values, a cold albedo α_C for the ice/snow-covered surface, when below the freezing temperature, and a warm albedo α_W corresponding to land or open ocean above the freezing temperature; that is,

$$\begin{array}{}\text{(A2)}& \mathit{\alpha }=\left\{\begin{array}{ll}{\mathit{\alpha }}_{\mathrm{C}}& \text{if}{T}_{\mathrm{S}}\le T_{\mathrm{R}},\\ {\mathit{\alpha }}_{\mathrm{W}}& \text{if}{T}_{\mathrm{S}}>T_{\mathrm{R}}.\end{array}\right.\end{array}$$

This discontinuous albedo function is conceptually simple but it is not an accurate representation of what would actually happen if the polar region cooled from ice-free to ice-covered conditions. Recall that this model represents the annually averaged climate. As the polar region cools, there will be a transition period in which warm ice-free summers get shorter and cold ice-covered winters get longer. The annually averaged albedo, therefore, would not jump abruptly from low to high constant values as in Eq. (A2); it would transition more smoothly between the summer and winter extreme values. Therefore, in this paper we introduce a more realistic (and smooth) sigmoidal albedo given by the hyperbolic tangent function, Eq. (3), and reproduced here:

$$\begin{array}{}\text{(A3)}& \mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left([{\mathit{\alpha }}_{\mathrm{W}}+{\mathit{\alpha }}_{\mathrm{C}}]+[{\mathit{\alpha }}_{\mathrm{W}}-{\mathit{\alpha }}_{\mathrm{C}}]\tanh \left({\displaystyle \frac{T_{\mathrm{S}}-T_{\mathrm{R}}}{\mathrm{\Omega }}}\right)\right).\end{array}$$

The parameter Ω determines the steepness of the transition between α_C and α_W. We have chosen $\mathit{\omega }\equiv \mathrm{\Omega }/T_{\mathrm{R}}=\mathrm{0.01}$. With this choice, the albedo changes by 80 % (that is, it goes from 0.1 to 0.9 of the difference between α_C and α_W) over a span of 6 C in annual mean temperature, which is realistic; see Fig. A1.

Figure A1 Graphs of the nondimensionalized EBM Eqs. (6)–(6). The red curves represent the surface EBM Eq. (6). The smooth red curves represent α(T_S) as in Eq. (A3), with $\mathit{\omega }\equiv \frac{\mathrm{\Omega }}{T_{\mathrm{R}}}=\mathrm{0.001},\mathrm{0.01},\mathrm{0.02},\mathrm{0.5}$, on successively flatter red curves. The blue curve represents the atmosphere Eq. (6). There are up to three intersections of the blue curve Eq. (6) with a given red curve Eq. (6) in this figure. Other parameters in the EBM are α_C=0.7, α_W=0.08, δ=0.67, μ=400, $F_{\mathrm{A}}=\mathrm{45}{\mathrm{Wm}}^{-\mathrm{2}}$, $F_{\mathrm{O}}=\mathrm{50}{\mathrm{Wm}}^{-\mathrm{2}}$ and Z=9 km.

Figure A1 shows the nondimensional energy balance Eq. (6) (in blue) and Eq. (6) (in red) for a range of values of the parameter $\mathit{\omega }=\mathrm{\Omega }/T_{\mathrm{R}}$ in Eq. (A3). The blue (atmosphere) curve may have up to three intersections with a red (surface) curve, implying the existence of multiple equilibrium states.

A3 Average global surface temperature

Both Wild et al. (2013) and Trenberth et al. (2009) indicate that the globally averaged amount of infrared radiation being emitted by the surface of the Earth is 398 W m⁻². Since the EBM assumes the Earth is a black-body radiator, the Stefan–Boltzmann law dictates that the corresponding temperature is

$$\begin{array}{}\text{(A4)}& T_{\mathrm{S}}^{\mathrm{avg}}={\left({\displaystyle \frac{\mathrm{398}}{\mathit{\sigma }}}\right)}^{\mathrm{1}/\mathrm{4}}=\mathrm{289.45}\mathrm{K}.\end{array}$$

This temperature is equivalent to 16.3 ^∘C, which is a little higher than the accepted average surface temperature of about 14 ^∘C. This is not surprising since the former is obtained by essentially averaging T⁴ while the latter is obtained from averaging actual temperatures. We use the above value for $T_{\mathrm{S}}^{\mathrm{avg}}$ in the calibrations described below, since it is consistent with the EBM.

A4 Heat convection/conduction and evapotranspiration

Evapotranspiration (ET) is the transport of water from the surface to the atmosphere in the form of water vapour. It combines the effects of evaporation from the surface and transpiration by plants. Globally, the largest contributor is evaporation from the surface of the oceans and the main determining factor there is the ocean surface temperature. The ET process also transports heat from the surface to the atmosphere, in the form of latent heat. Recently Wang et al. (2010) have shown that global ET has been increasing over the past several decades.

The forcing term F_C in the EBM represents transport of heat away from the surface to the atmosphere by conduction plus convection plus change of state of water. The most important of these is the upward transport of latent heat. Surface water evaporates, taking heat from the surface. As warm moist air rises and cools, the water vapour condenses, releasing its latent heat into the surrounding atmosphere. According to both Wild et al. (2013) and Trenberth et al. (2009), the magnitude of this forcing is $F_{\mathrm{C}}=\mathrm{104}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. However, this is a globally averaged value and, due to strong dependence on temperature, is not likely valid at temperatures other than the globally averaged surface temperature $T_{\mathrm{S}}^{\mathrm{avg}}$. In the past several decades, ET has increased by 0.6 W m⁻² per decade (Wang et al., 2010). Extrapolating back 100 years, this would be 6 W m⁻². In the past century the average global surface temperature has risen about 1 ^∘C. Thus it is reasonable to assume the dependence of F_C on the surface temperature T_S should have a slope near $T_{\mathrm{S}}=T_{\mathrm{S}}^{\mathrm{avg}}$ of about 6. If the dependence of F_C on T_S was simply linear, then, given that $F_{\mathrm{C}}\left(T_{\mathrm{S}}^{\mathrm{avg}}\right)=\mathrm{104}$, this slope would predict that F_C was negative for $T_{\mathrm{S}}<T_{\mathrm{S}}^{\mathrm{avg}}-\mathrm{17.3}\approx -\mathrm{3}$ ^∘C. Since negative values for F_C are unreasonable, and since there likely is still some heat transport even near freezing, instead of a linear function, we have modelled the dependence of F_C on temperature as a hyperbolic function that is nearly zero for temperatures below freezing, and increases roughly linearly for temperatures above freezing. Specifically, we use the function

$$(T_{\mathrm{S}}-T_{\mathrm{R}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{A}_{\mathrm{1}}}}\left(F_{\mathrm{C}}-{\displaystyle \frac{A_{\mathrm{2}}^{\mathrm{2}}}{F_{\mathrm{C}}}}\right),$$

or equivalently,

$$\begin{array}{}\text{(A5)}& F_{\mathrm{C}}\left(T_{\mathrm{S}}\right)=A_{\mathrm{1}}(T_{\mathrm{S}}-T_{\mathrm{R}})+\sqrt{A_{\mathrm{1}}^{\mathrm{2}}(T_{\mathrm{S}}-T_{\mathrm{R}}{)}^{\mathrm{2}}+A_{\mathrm{2}}^{\mathrm{2}}},\end{array}$$

where A₁ and A₂ are constants. In terms of the nondimensional variables and parameters Eq. (6), this equation becomes

$$\begin{array}{}\text{(A6)}& f_C\left({\mathit{\tau }}_{\mathrm{S}}\right)=a_{\mathrm{1}}({\mathit{\tau }}_{\mathrm{S}}-\mathrm{1})+\sqrt{a_{\mathrm{1}}^{\mathrm{2}}({\mathit{\tau }}_{\mathrm{S}}-\mathrm{1}{)}^{\mathrm{2}}+a_{\mathrm{2}}^{\mathrm{2}}}.\end{array}$$

The constants A₁ and A₂ are chosen so that the forcing at the global average surface temperature $T_{\mathrm{S}}^{\mathrm{avg}}$ is 104, and so that the forcing at the freezing point is 20 % of this value, that is, $F_{\mathrm{C}}\left(T_{\mathrm{R}}\right)=\mathrm{0.2}{F}_{\mathrm{C}}\left(T_{\mathrm{S}}^{\mathrm{avg}}\right)$. Thus

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{A}_{\mathrm{2}}=\mathrm{0.2}\cdot \mathrm{104}=\mathrm{20.8},\\ \text{(A7)}& {\displaystyle }& {\displaystyle }\text{and}{A}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}(T_{\mathrm{S}}^{\mathrm{avg}}-T_{\mathrm{R}})}}\left(\mathrm{104}-{\displaystyle \frac{A_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{104}}}\right)=\mathrm{3.063},\end{array}$$

or in terms of the nondimensional parameters,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{a}_{\mathrm{1}}={\displaystyle \frac{A_{\mathrm{1}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{3}}}}=\mathrm{2.650},\\ \text{(A8)}& {\displaystyle }& {\displaystyle }\text{and}{a}_{\mathrm{2}}={\displaystyle \frac{A_{\mathrm{2}}}{\mathit{\sigma }{T}_{\mathrm{R}}^{\mathrm{4}}}}=\mathrm{6.590}\times {\mathrm{10}}^{-\mathrm{2}}.\end{array}$$

Figure A2 shows a plot of Eq. (A5) for these parameter values. With these choices of the constants, the slope at $T_{\mathrm{S}}^{\mathrm{avg}}$ is 5.9, agreeing well with the above argument that it should be about 6.

Figure A2 Variation of evapotranspiration/convection with surface temperature. The hyperbolic function is tangent to the horizontal as the temperature decreases, and grows almost linearly with a slope of about 6 as temperature increases.

A5 Atmospheric absorption of infrared radiation

Two of the empirical parameters in the EBM are k_W and k_C, defined in Sect. 2.2 as the grey gas absorption coefficients for water vapour and carbon dioxide, respectively. More precisely, the absorption coefficient k of a gas is its proportionality coefficient in the Beer–Lambert law, as defined in Eq. (6) of Sect. 2.2. In the grey gas approximation, it is assumed that the absorption coefficient k is constant for all frequencies of incident radiation, and depends only on the total intensity of that radiation. In reality though, k is not an absolute physical constant because gases absorb radiation at discrete spectral lines. The total amount of energy absorbed depends on the distribution of energy in the radiation at each of these spectral wavelengths. Values of k_W and k_C for the particular spectrum of incident radiation experienced by the atmosphere may be obtained empirically from modern era satellite and land-based measurements. We assume that k_W and k_C for the atmosphere do not change with time, even over the geological timescales considered here.

Atmospheric absorption, η, of infrared radiation emitted by the Earth is, in this EBM, attributed to the two primary greenhouse gases: water vapour and carbon dioxide, and to the liquid and solid water making up clouds. According to Trenberth et al. (2009), at present as a global average, the Earth emits 398 W m⁻² of which 40.1 W m⁻² passes through the atmosphere to space. Thus approximately 90 % of the radiation is absorbed in the atmosphere, that is η=0.9. The present-day absorption factor η can be attributed as follows (Schmidt et al., 2010): water vapour 50 %, clouds 25 % and CO₂ (plus other gases) 25 %. Furthermore, they determined that these ratios remain unchanged, even after a doubling of CO₂. Therefore, setting $x={\mathit{\eta }}_C={\mathit{\eta }}_{\mathrm{Cl}}$ and η_W=2x in Eq. (10) we get

$$\mathrm{0.9}=\mathrm{1}-(\mathrm{1}-x)(\mathrm{1}-\mathrm{2}x)(\mathrm{1}-x)=\mathrm{2}{x}^{\mathrm{3}}-\mathrm{5}{x}^{\mathrm{2}}+\mathrm{4}x.$$

This cubic has one real root, $\widehat{x}$,

$$\begin{array}{}\text{(A9)}& \widehat{x}\approx \mathrm{0.3729}.\end{array}$$

In this model, since we have no data on cloud cover for paleoclimates, we make the assumption that absorption due to clouds is a constant globally and temporally, and set

$$\begin{array}{}\text{(A10)}& {\mathit{\eta }}_{\mathrm{Cl}}=\widehat{x}.\end{array}$$

We also tried varying η_Cl with temperature using a hyperbolic function like that used for F_C in Eq. (A5); however, we found no qualitative change in the results and therefore, for simplicity, have left η_Cl as a constant. To calibrate the absorption coefficient k_C for the greenhouse gas carbon dioxide, we set ${\mathit{\eta }}_C=\widehat{x}$ in Eq. (14) yielding

$$\begin{array}{ll}{\displaystyle }\widehat{x}& {\displaystyle }=\mathrm{1}-e^{-\mathit{\mu }\cdot \mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{6}}{k}_{\mathrm{C}}{P}_{\mathrm{A}}/g}\\ {\displaystyle }⟹& {\displaystyle }{k}_{\mathrm{C}}={\displaystyle \frac{g\ln (\mathrm{1}-\widehat{x})}{\mathit{\mu }(\mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{6}})P_{\mathrm{A}}}}.\end{array}$$

Using the present-day value of μ=400 ppm, the value of $\widehat{x}$ from Eq. (A7), the standard atmospheric pressure value of $P_{\mathrm{A}}=\mathrm{101.3}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{Pa}$, and the acceleration due to gravity $g=\mathrm{9.8}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$ we get

$$\begin{array}{ll}{\displaystyle }{k}_{\mathrm{C}}& {\displaystyle }={\displaystyle \frac{\left(\mathrm{9.8}\right)\ln (\mathrm{1}-\widehat{x})}{\left(\mathrm{400}\right)(\mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{6}})(\mathrm{101.3}\times {\mathrm{10}}^{\mathrm{3}})}}\\ \text{(A11)}& {\displaystyle }& {\displaystyle }=\mathrm{0.07424}{\mathrm{m}}^{\mathrm{2}}{\mathrm{kg}}^{-\mathrm{1}}.\end{array}$$

With this value for k_C, the coefficient G_C, defined in Eq. (13) has the value

$$\begin{array}{}\text{(A12)}& G_C=\mathrm{1.52}\times {\mathrm{10}}^{-\mathrm{6}}{k}_{\mathrm{C}}{\displaystyle \frac{P_{\mathrm{A}}}{g}}=\mathrm{1.166}\times {\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

To calibrate k_W we proceed as follows. The latent heat of vaporization of water is $L_{\mathrm{v}}=\mathrm{2.2558}\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}$ and the ideal-gas constant specific to water vapour is $R_{\mathrm{W}}=\mathrm{461.5}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$. Thus from Eq. (24) we have

$$\begin{array}{}\text{(A13)}& G_{\mathrm{W}\mathrm{1}}={\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{W}}{T}_{\mathrm{R}}}}=\mathrm{17.89}.\end{array}$$

Dai (2006) indicates that an average value for relative humidity, δ, from 1976 to 2004 for the region of the Earth between 60^∘ S and 75^∘ N is 0.74. Over the same latitudes, Dai reports that averages over water and over land are 0.79 and 0.65, respectively, except that over deserts the humidity drops to 30 %–50 %. The polar regions are typically drier; however, there is a significant amount of water north of 75^∘ N and south of 60^∘ S, hence we have chosen δ^avg=0.74 as the average global relative humidity for purposes of calibrating k_W. In Eq. (25) we set ${\mathit{\eta }}_{\mathrm{W}}=\mathrm{2}\widehat{x}$ with $\widehat{x}$ given by Eq. (A7), τ_S to the present normalized global average temperature $T_{\mathrm{S}}^{\mathrm{avg}}/T_{\mathrm{R}}$, δ to the present average global relative humidity δ^avg, Z to an average global tropopause height of Z^avg=14000 m, and the normalized lapse rate to $\mathit{\gamma }=\mathrm{\Gamma }/T_{\mathrm{R}}=\mathrm{2.38}\times {\mathrm{10}}^{-\mathrm{5}}{\mathrm{m}}^{-\mathrm{1}}$. Using these values and inverting Eq. (25) to isolate k_W we get

$$\begin{array}{ll}{\displaystyle }{k}_{\mathrm{W}}& {\displaystyle }=-\ln (\mathrm{1}-\mathrm{2}\widehat{x})\mathit{\gamma }{\left[\mathit{\delta }{\displaystyle \frac{P_{\mathrm{W}}^{\mathrm{sat}}\left(T_{\mathrm{R}}\right)}{R_{\mathrm{W}}{T}_{\mathrm{R}}}}\underset{T_{\mathrm{S}}^{\mathrm{avg}}-\mathit{\gamma }{Z}^{\mathrm{avg}}}{\overset{T_{\mathrm{S}}^{\mathrm{avg}}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\exp \left(G_{\mathrm{W}\mathrm{1}}\left[{\displaystyle \frac{\mathit{\tau }-\mathrm{1}}{\mathit{\tau }}}\right]\right)\mathrm{d}\mathit{\tau }\right]}^{-\mathrm{1}}\\ \text{(A14)}& {\displaystyle }& {\displaystyle }=\mathrm{0.05905}{\mathrm{m}}^{\mathrm{2}}{\mathrm{kg}}^{-\mathrm{1}}.\end{array}$$

With this value for k_W, the second greenhouse gas parameter for water vapour is

$$\begin{array}{}\text{(A15)}& G_{\mathrm{W}\mathrm{2}}\equiv {\displaystyle \frac{k_{\mathrm{W}}{P}_{\mathrm{W}}^{\mathrm{sat}}\left(T_{\mathrm{R}}\right)}{\mathit{\gamma }{R}_{\mathrm{W}}{T}_{\mathrm{R}}}}=\mathrm{12.05}.\end{array}$$

A6 Atmospheric emission of radiation

In previous slab models, including Dortmans et al. (2018) and Payne et al. (2015), it was assumed that the upward and downward radiation intensities from the atmosphere are equal; that is, $\mathit{\beta }=\frac{\mathrm{1}}{\mathrm{2}}$ in Eq. (1). This would be the case for an actual uniform slab; however, the real atmosphere is not uniform in temperature nor density, and in fact both the temperature and density of the atmosphere are much higher in value near the surface than near the tropopause. The net effect of this nonuniformity is that, of the total radiation intensity I_A emitted by the atmosphere, almost two-thirds goes back to the surface and only a little more than one-third escapes to space, according to satellite and surface data (Trenberth et al., 2009; Wild et al., 2013). From Trenberth et al. (2009), the atmosphere and clouds emit $\mathrm{169.9}+\mathrm{29.9}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ upward to space and emit 340.3 W m⁻² downward to the Earth. Therefore, in order to have the model in this paper represent the atmosphere more realistically, we set

$$\begin{array}{}\text{(A16)}& \mathit{\beta }={\displaystyle \frac{\mathrm{340.3}}{\mathrm{340.3}+\mathrm{169.9}+\mathrm{29.9}}}\approx \mathrm{0.63}.\end{array}$$

instead of 0.50.

Appendix B: Equilibrium climate sensitivity

Equilibrium climate sensitivity (ECS) is a useful measure of the sensitivity of a given climate model to an increase in CO₂ concentration μ. It is usually defined as the change ΔT in the global mean temperature $\stackrel{\mathrm{‾}}{T}$, resulting from a doubling of μ, starting from the accepted pre-industrial value μ=270 ppm (IPCC, 2013; Forster, 2016; Knutti et al., 2017; Proistosescu and Huybers, 2017). The ECS provides a first-order estimate of the amount of present-day global warming predicted by a given model as μ increases. Because ECS is a single number, it facilitates comparisons between models. It has been used extensively in the assessment reports of the IPCC (IPCC, 2013), where ΔT values in the range 1.5 to 4.5 ^∘C have been documented. Some have cited this wide range of ECS values as a sign of weakness of the IPCC methodology. Recently these results have been reconciled (Proistosescu and Huybers, 2017). They showed that a linear statistical analysis of historical data gives low estimates in the 1.5 to 3 ^∘C range, while nonlinear models give higher estimates.

In terms of our EBM, to apply it to a global average we make the following settings. The global mean insolation at the top of the atmosphere is $Q=\mathrm{340}{\mathrm{Wm}}^{-\mathrm{2}}$ (Wild et al., 2013; Trenberth et al., 2009). They indicate that 161 W m⁻² of sunlight is absorbed by the Earth while 24 W m⁻² is reflected back on average. Thus we take

$${\mathit{\alpha }}_{\mathrm{W}}={\displaystyle \frac{\mathrm{24}}{\mathrm{161}+\mathrm{24}}}=\mathrm{0.13}.$$

For the tropopause height Z we take an average value of Z^avg=14 km. For the global model, the ocean and atmospheric transport terms, F_O and F_A, will be zero.

The ECS value we calculate will be the difference between the warm equilibrium temperatures when μ=270 and μ=540. For these calculations we took the relative humidity to be $\mathit{\delta }={\mathit{\delta }}^{\mathrm{avg}}=\mathrm{0.74}$, in accord with the value obtained in Dai (2006). The graphs of the surface and atmosphere equilibrium equations are shown in Fig. B1. The equilibrium solutions of this global mean version of the EBM yield T_S(270)=14.3, T_S(540)=17.6 and ΔT=3.3 ^∘C. This is in excellent agreement with accepted values (IPCC, 2013).

Figure B1 Equilibrium climate sensitivity (ECS) for millennial data. The equilibrium points marked with stars show the change in T_S for a doubling of μ from 270 to 540 ppm; T_S(270)=14.3, T_S(540)=17.6 and ΔT=3.3 ^∘C.