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.


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 2 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 icecapped. The major climate forcing factors (solar constant, orbital parameters, CO 2 concentration and locations of the continents) were all very similar to today (Fedorov et al., 2006(Fedorov et al., , 2010Haywood 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., 2006Fedorov et al., , 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); 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 2 concentration worldwide (DeConto et al., 2008;Goldner et al., 2014;Pagani et al., 2011). The EBM of this paper includes both CO 2 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., 1981Barron et al., , 1995Sloan and Barron, 1992). In a similar vein, Huber and Caballero (2011) 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 Hutchinson et al., 2018;Lunt et al., 2016Lunt et al., , 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 the-ory 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 2 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 2 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.
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 2 increases with the density of the atmosphere but is independent of temperature, while the greenhouse warming of H 2 O increases with temperature but is independent of the density (or partial pressure) of the other gases present. The greenhouse warming of methane (CH 4 ) acts in a similar fashion to that of CO 2 ; therefore, CH 4 can be incorporated into the CO 2 concentration. As an increase in CO 2 concentration causes climate warming, this warming causes an increase in evaporation of H 2 O into the atmosphere, which further increases the climate warming beyond that due to CO 2 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 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 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 2 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).

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  Table 1 and Sect. 2.1.
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.
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.

Energy balance
The model consists of two energy balance equations, one for the atmosphere and one for the surface. In Fig. 1, the socalled 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 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 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 where A 1 and A 2 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 = 173.2 Wm −2 and at the Equator is 418.8 Wm −2 (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 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: 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. (5).
The emission of long-wave radiation, I , from a body is governed by the Stefan-Boltzmann law, I = σ T 4 , 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 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 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 (2).
Here instead, we use (T S , I A ) as the two independent variables, and then we formally let T A be defined by Eq. (7). 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 After normalization, the freezing temperature of water is represented by τ = 1 and the atmosphere and surface energy balance Eqs. (1)-(6) simplify to where, from Appendix A, 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 0.8 ≤ τ ≤ 1.2, which corresponds approximately to a range in more familiar degrees Celsius of −54 • C ≤ T ≤ +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.

Optical depth and the Beer-Lambert law
The goal of this section is to define the absorptivity parameter η in the EBM Eq. (9) (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 µ, δ, k C , k 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 where k (m 2 kg −1 ) is the absorption coefficient of the material, ρ (kg m −3 ) is the density of the absorbing substance such as CO 2 and z (m) is distance along the path. The differential Eq. (13) may be integrated from z = 0 (the surface) to z = Z (the tropopause) to give 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 (9) 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 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 so that (17) 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. (14), and then combine them using Eq. (16). This calculation is presented for the case of CO 2 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 .

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

The grey gas approximation
Although it is well-known that gases like CO 2 and H 2 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 2 and H 2 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.

Carbon dioxide
The concentration of CO 2 in the atmosphere is usually expressed as a ratio, in molar parts per million (ppm) of CO 2 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 2 is approximately mm C = 44×10 −3 kg mol −1 . 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 −1 , respectively. Neglecting other trace gases in the atmosphere, a weighted average gives the molar mass of the dry atmosphere as mm A = 29 × 10 −3 kg mol −1 . Therefore, the CO 2 concentration µ measured in molar ppm is converted to mass concentration in kg ppm by multiplication by the ratio mm C /mm A ≈ 1.52. If ρ A (z) is the density of the atmosphere at altitude z in kg m −3 , then the mass density of CO 2 at the same altitude, with molar concentration µ ppm, is It is known that CO 2 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 2 decreases at exactly the same rate, according to Eq. (18). Substituting Eq. (18) into Eq. (14) determines the optical depth λ C of CO 2 λ C = 1.52 µ Now consider a vertical column of air, of unit cross section, from surface to tropopause. The integral in Eq. (19) 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 Z 0 ρ(A) dz = P A /g, where g is acceleration due to gravity. Therefore, the optical depth λ C of CO 2 in the actual atmosphere in Eq. (19) is With λ C so determined, it follows from Eq. (15) that 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 2 is G C = 1.166 × 10 −3 . The dry atmosphere EBM is obtained by assuming there is no water vapour and no clouds in the atmosphere. Hence η = η C and δ = 0. We assume also that ξ R = ξ A = 0. Making these substitutions in the EBM Eqs. (9)-(10) yields Due to the nonlinearity of the ice-albedo function in 5, there can be in fact up to three points of intersection as shown in Fig. 2a. If the CO 2 level µ decreases sufficiently, the warm state (τ S > 1) disappears, while as CO 2 increases, the frozen state (τ S < 1) may disappear. The CO 2 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. (22)-(23) 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.

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 2 , the concentration of water vapour (H 2 O) in the atmosphere varies widely with location and altitude. This is because the partial pressure of H 2 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 sat W (T ). The actual partial pressure of water vapour is then a fraction, δ, of this saturated value, where δ is called relative humidity. While P sat W (T ) varies greatly with T in the atmosphere, the relative humidity δ is comparatively constant. If the actual P W (T ) exceeds the saturated value P sat W (T ) (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), 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. (24) and (25). We may use the ideal-gas law in the form P W = ρ W R W T 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. (24) and (25) gives Transforming to the dimensionless temperature τ = T /T R as in Sect. 2.1, this becomes 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. (15) and (14), Here k W is the absorption coefficient of water vapour; see Appendix A5. In order to evaluate the integral in Eq. (28), we need to know how ρ W varies with height z. We have shown that ρ W is a function of temperature, given by Eqs. (26) or (27). 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 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 = 6.49×10 −3 K m −1 (ICAO, 1993). Using this assumption, the variation in temperature with height is given as where the normalized lapse rate is γ = /T R = 2.38 × 10 −5 m −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 (30) may be used to change the variable of integration in Eq. (28) for the optical depth of water vapour, from z to τ . The result is Now substitute Eq. (27) into the integral Eq. (31) and simplify to where the greenhouse gas parameters G W1 and G W2 for water vapour are defined as Finally, using the Beer-Lambert law, the absorptivity η W of water vapour in Eq. (28) is determined by its optical depth, and is now a function of the surface temperature, 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 2 kg −1 and the greenhouse gas parameters for H 2 O are G W1 = 17.89 and G W2 = 12.05. Figure 3a shows that the function η W (τ S ) in Eq. (34) 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. (10) increases past τ S = 1, then the absorptivity of water vapour η W (τ S ), acting as a greenhouse gas in Eq. (9), increases rapidly, thus amplifying the heating of the atmosphere by the radiation I S = σ T 4 S 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 green- .08, and Z = 9 km). (a) Absorptivity η W (τ S ) as given by Eq. (34) with relative humidity δ = 0.4, 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 .

Combined CO 2 and H 2 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 2 as in Eq. (21) and η W is the absorptivity of H 2 O as in Eq. (28), 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. (17), The full nondimensional two-layer EBM is therefore specified by Eqs. (9)- (12) and (35), and the parameter values in Table 1.

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 2 itself in the atmosphere. However, if the concentration of CO 2 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 2 alone.

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.

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 2 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., , 2013, as has the Pliocene Model Intercomparison Project (PlioMIP; Haywood et al., 2011Haywood et al., , 2016Zhang et al., 2013). Advances in the extraction and interpretation of proxy data have given a clearer picture of the warm Pliocene climate 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  Peltier (2017, 2018), and other references therein. Using a fully coupled atmosphereocean 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 2 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 2 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 2 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., 2006Fedorov et al., , 2010Zhang 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 2 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 30million-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 2 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 Wolfe et al., 2017;Huber and Caballero, 2011); also, CO 2 concentration µ was higher than today. Estimates of Eocene CO 2 concentration µ vary from 1000-1500 ppm (Pagani et al., 2005(Pagani et al., , 2006 to 490 ppm (Wolfe et al., 2017). For this EBM, we set mid-Eocene CO 2 at µ = 1000 (Pagani et al., 2005). Both temperature and CO 2 concentration have decreased steadily but not monotonically, with many fluctuations from their Eocene values to pre-industrial modern values. The overall decrease in CO 2 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 O = 50 W m −2 , because the change in F O is relatively small. The blue atmosphere equilibrium curves represent values of CO 2 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 saddlenode 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 2 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.
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 t = 50(1 − ν) Ma.
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 2 concentration µ decreased more significantly. Therefore, we express both µ and F O as decreasing functions of bifurcation parameter ν µ = 1000 − 730 · ν ppm 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 preindustrial values (IPCC, 2013). Equation (37) 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 saddlenode bifurcation. This bifurcation occurs at approximately ν = 0.91, which corresponds to F O = 50.9 Wm −2 and µ = 336, which is in good agreement with the determination (Seki et al., 2010) that the warm Pliocene pCO 2 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. (36), 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 2 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(Fedorov et al., , 2010Steph 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).

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 Lear et al., 2008;Miller et al., 2008;Scher et al., 2011Scher et al., , 2015Ladant et al., 2014;Ruddiman, 2014). In the mid-Cretaceous, the continent of Antarctica was in the South Pa-cific 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 −2 (Barron et al., 1981), to a much lower value, estimated here to be 30 Wm −2 . 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 2 level in the Antarctic is the same as elsewhere. For this EBM, we estimate the early Eocene global CO 2 level as µ = 1100 ppm, (Pagani et al., 2005(Pagani et al., , 2006Cronin, 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 2 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 2 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 2 concentration in the atmosphere. Since both CO 2 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 shows the solutions of the EBM equations for Antarctic parameter values, including δ = 0.67, F A = 45 Wm −2 , Q = 173.2 Wm −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 2 concentration µ is decreased from the top to bottom blue lines. In Fig. 5b the reverse is true; the CO 2 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 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 2 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: Here, ν = 0 corresponds to the high early Eocene values of µ = 1100 ppm and F O = 100 Wm −2 , while ν = 1 corresponds to the low Oligocene-Miocene boundary values of µ = 400 ppm and F O = 30 Wm −2 . Thus, in Eq. (39), 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 A = 45 Wm −2 . The solar radiation is Q = 173.2 Wm −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 6a shows a saddle-node bifurcation at ν = 0.779, corresponding to forcing parameter values µ = 555 ppm and F O = 45.5 W m −2 in the model, and to geological time about 30 Ma, assuming the linear time relation in Eq. (38). 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  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. (39), and, on the magenta curve, µ is held fixed at the Eocene value while F O decreases as in Eq. (39). No saddle-node bifurcation destroying the warm equilibrium occurs in either scenario. 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 2 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;Pa-gani et al., 2011;Scher et al., 2011Scher et al., , 2015. The green curve represents a scenario in which µ decreases as in Eq. (39) 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. (39). In neither case does a glaciation event occur. The analysis of this paper implies that significant decreases in both CO 2 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 Lear et al., 2008; 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 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 2 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.

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.
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 ν. µ = 1130 − 860 · ν, The mid-Cretaceous CO 2 concentration µ = 1130 ppm is as determined by Fletcher et al. (2008), see Table 3. The ocean heat transport F O = −55 Wm −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.

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., 1981Barron et al., , 1995Sloan 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 2 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 2 concentrations between 600 and 2400 ppmv, and  Table 2. (a) The blue curves represent the atmosphere equilibrium equation for µ = 300, 600, 900 and 1200 from bottom to top. (b) Bifurcation diagram for the tropics EBM. Here, CO 2 concentration µ and ocean heat transport F O are both decreasing, with increasing ν, (0 ≤ ν ≤ 1), as given by Eq. (40). As ν varies, both the warm equilibrium state (τ S > 1) and the frozen equilibrium state (τ S < 1) persist; there are no saddle-node bifurcations.
tropical Atlantic upper ocean temperatures between 33 and 42 • C. Based on fossil samples, Fletcher et al. (2008) estimate mid-Cretaceous atmospheric CO 2 concentrations of 1130 ppmv and we choose their estimate for use in this EBM, see Table 3.
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. (9)-(10), 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 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 2 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 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).  Barron et al. (1981) 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., 1981Barron et al., , 1995Sloan and Barron, 1990) and the early Eocene 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.

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 2 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,  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.
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 2 concentrations, over the geological time period studied here (Cronin, 2010;Pagani et al., 1999Pagani et al., , 2005. The model of this paper assumes a smooth decline in CO 2 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 2 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 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 onedimensional 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 .
Code availability. The original MATLAB code for the computations that support this research is available in the appendix of Dortmans (2017).

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 B.

A1 Solar radiation
The value of the annually averaged solar radiation at either pole is Q = 173.2 Wm −2 and at the Equator is 418.8 Wm −2 (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 = 340 W m −2 . Of this, 79 W m −2 is directly absorbed by the atmosphere, and 76 W m −2 is reflected by the atmosphere back into space (both primarily due to clouds). Hence we define the two solar radiation fractions 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, 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. (5), and reproduced here: The parameter determines the steepness of the transition between α C and α W . We have chosen ω ≡ /T R = 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 shows the nondimensional energy balance Eq. (9) (in blue) and Eq. (10) (in red) for a range of values of the parameter ω = /T 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) andTrenberth et al. (2009) indicate that the globally averaged amount of infrared radiation being emitted by the surface of the Earth is 398 W m −2 . Since the EBM assumes the Earth is a black-body radiator, the Stefan-Boltzmann law dictates that the corresponding temperature is 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 4 while the latter is obtained from averaging actual temperatures. We use the above value for T avg S 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 com- Figure A1. Graphs of the nondimensionalized EBM Eqs. (9)-(10). The red curves represent the surface EBM Eq. (10). The smooth red curves represent α(T S ) as in Eq. (A3), with ω ≡ T R = 0.001, 0.01, 0.02, 0.5, on successively flatter red curves. The blue curve represents the atmosphere Eq. (9). There are up to three intersections of the blue curve Eq. (9) with a given red curve Eq. (10) in this figure. Other parameters in the EBM are α C = 0.7, α W = 0.08, δ = 0.67, µ = 400, F A = 45 Wm −2 , F O = 50 Wm −2 and Z = 9 km. bines 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) andTrenberth et al. (2009), the magnitude of this forcing is F C = 104 W m −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 avg S . In the past several decades, ET has increased by 0.6 W m −2 per decade (Wang et al., 2010). Extrapolating back 100 years, this would be 6 W m −2 . 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 S = T avg S of about 6. If the de-pendence of F C on T S was simply linear, then, given that F C (T avg S ) = 104, this slope would predict that F C was negative for T S < T avg S − 17.3 ≈ −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 or equivalently, where A 1 and A 2 are constants. In terms of the nondimensional variables and parameters Eq. (8), this equation becomes The constants A 1 and A 2 are chosen so that the forcing at the global average surface temperature T avg S is 104, and so that the forcing at the freezing point is 20 % of this value, that is, F C (T R ) = 0.2F C (T avg S ). Thus (A8) Figure A2 shows a plot of Eq. (A5) for these parameter values. With these choices of the constants, the slope at T avg S is 5.9, agreeing well with the above argument that it should be about 6.

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. (13) 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 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. 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 −2 of which 40.1 W m −2 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 : water vapour 50 %, clouds 25 % and CO 2 (plus other gases) 25 %. Furthermore, they determined that these ratios remain unchanged, even after a doubling of CO 2 . Therefore, setting x = η C = η Cl and η W = 2x in Eq. (17) we get This cubic has one real root,x, x ≈ 0.3729.
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 η Cl =x.
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 η C =x in Eq. (21) yieldinĝ x = 1 − e −µ·1.52×10 −6 k C P A /g ⇒k C = g ln(1 −x) µ(1.52 × 10 −6 )P A .
With this value for k C , the coefficient G C , defined in Eq. (20) has the value G C = 1.52 × 10 −6 k C P A g = 1.166 × 10 −3 .
To calibrate k W we proceed as follows. The latent heat of vaporization of water is L v = 2.2558 × 10 6 m 2 s −2 and the ideal-gas constant specific to water vapour is R W = 461.5 m 2 s −2 K −1 . Thus from Eq. (33) we have 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. (34) we set η W = 2x withx given by Eq. (A9), τ S to the present normalized global average temperature T avg S /T 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 γ = /T R = 2.38 × 10 −5 m −1 .
Using these values and inverting Eq. (34) to isolate k W we get With this value for k W , the second greenhouse gas parameter for water vapour is

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, β = 1 2 in Eq.
(2). 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 169.9 + 29.9 W m −2 upward to space and emit 340.3 W m −2 downward to the Earth. Therefore, in order to have the model in this paper represent the atmosphere more realistically, we set β = 340.3 340.3 + 169.9 + 29.9 ≈ 0.63.
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 2 concentration µ. It is usually defined as the change T in the global mean temperature 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 = 340 Wm −2 (Wild et al., 2013;Trenberth et al., 2009). They indicate that 161 W m −2 of sunlight is absorbed by the Earth while 24 W m −2 is reflected back on average. Thus we take α W = 24 161 + 24 = 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 δ = δ avg = 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).