Thermal evolution model

From energy conservation in a cooling convective system, the rate of change of the average internal temperature, T _i, is equal to the energy generated by the decay of radiogenic isotopes in the mantle, H, minus the surface heat loss, Q, as given by (Christensen Reference Christensen1985):

(1)$$ C_{\rm E} {\displaystyle{{\rm d}T_{i}} \over {{\rm d}t}} = H - Q $$

where C _E is the effective heat capacity of the whole Earth, taking into account the temperature gradient. This formulation does not discriminate between core and mantle contributions, which would require a parameterization of core–mantle interactions, but has the core contribution imbibed in it by the use of the whole Earth effective heat capacity (Labrosse and Jaupart Reference Labrosse and Jaupart2007). We will approximate the internal temperature by the effective mantle potential temperature, T _m.

Our formulation for H differs from the single term form commonly used in previous studies because we discriminate the four main decaying radiogenic heat sources, ²³⁵U, ²³⁸U, ²³²Th and ⁴⁰K (Labrosse and Jaupart Reference Labrosse and Jaupart2007; Turcotte and Schubert Reference Turcotte and Schubert2014):

(2)$$ H(t) = f_{\rm rm} \left[\sum_{i} C_{i} h_{i}\times \exp \left(-{\displaystyle{t\times \ln(2)} \over {t^{i}_{1/2}}}\right)\right] M_{\rm m} $$

where f _rm is the fraction of radiogenic material still present in the mantle, i denotes the isotope species (²³⁵U, ²³⁸U, ²³²Th and ⁴⁰K), C _i is the initial concentration of the radiogenic isotope i in the mantle, h _i is the energy release rate of the isotope i, t is the time, $t^{i}_{1/2}$ is the half-life of the isotope i and M _m is the mantle's mass. The difference between the two forms is small in the past, but becomes significant in the future due to the long half-life of ²³²Th. We used the estimates of Jaupart et al. (Reference Jaupart, Labrosse and Mareschal2007) and Arevalo Jr et al. (Reference Arevalo, McDonough and Luong2009) to set the past concentrations of the radiogenic isotopes in t = 0 in order to give the present total Earth heat budget of 46(2) TW, with 38.7 TW from the mantle, 21(4) TW from the radiogenic isotopes in the crust and the mantle and only 13.7 TW from the radiogenic isotopes in the mantle alone (Jaupart et al. Reference Jaupart, Labrosse and Mareschal2007; Arevalo Jr et al. Reference Arevalo, McDonough and Luong2009). See Table 1 for the initial concentrations and other parameters.

Table 1. Parameters for equation (2) using data from Arevalo Jr et al. (Reference Arevalo, McDonough and Luong2009)

The mantle would be richer in radiogenic elements if not for the continental crust, which concentrates radiogenic elements enough to produce 7.30(23) TW. We model this with the factor f _rm = 1 − (7.3/21)A _cc,N, where the 7.3/21 is the ratio of the radiogenic heat production of the continental crust to the total radiogenic heat production from the bulk silicate Earth today, and A _cc,N is the normalized continental crust area.

The rheology of our model uses an Arrhenius-type viscosity, with a dependence on the mean potential temperature of the mantle:

(3)$$ \eta = \eta_{0}\,{\rm exp}\left[{\displaystyle{E} \over {R}}\left({\displaystyle{1} \over {T_{\rm m}}}-{\displaystyle{1} \over {T_{0}}}\right)\right] $$

where η₀ is a reference viscosity in the reference temperature T ₀, E is the activation energy and R is the molar gas constant.

The plate-tectonics heat loss from the mantle is parameterized according to Conrad and Hager (Reference Conrad and Hager1999a,Reference Conrad and Hagerb) and Korenaga (Reference Korenaga2006):

(4)$$ Q = A \left[\displaystyle{C_{1} \alpha \rho_{\rm m} g (T_{\rm m} - T_{\rm s,0})^{3} (R_{\rm e} - R_{\rm i}) h} \over {C_{2} \eta + C_{3} \eta_{\rm L} (h/R_{\rm pb})^{3}} \right]^{{1}/{2}} $$

where A is a constant adjusted to fit the mantle convective heat loss of 38.7 TW, α is the coefficient of thermal expansion, ρ_m is the upper mantle mean density, g is the acceleration of gravity, T _s,0 is the reference surface temperature of 273 K, R _e and R _i are the external and internal radii of the mantle, respectively, h is the thickness of plate, η_L is the lithospheric viscosity and R _pb is the radius of curvature for plate bending. The remaining constants are given in Conrad and Hager (Reference Conrad and Hager1999a): C ₁ = π^−0.5, C ₂ ≈ 13.5 (Korenaga Reference Korenaga2006), and C ₃ ~ 2.5. This parameterization preserves the general shape of the formulation of McGovern and Schubert (Reference McGovern and Schubert1989) and Turcotte and Schubert (Reference Turcotte and Schubert2014), with a power law dependence on the mantle temperature, but with a different exponent, leading to a weaker dependence on the temperature, a behaviour more in line with the geophysical and geochemical constraints (Korenaga Reference Korenaga2006), as we will see below from the Urey ratio.

The mantle model is joined to the atmosphere model using the areal sea-floor spreading rate, S _r, as given in McGovern and Schubert (Reference McGovern and Schubert1989) and Franck (Reference Franck1998):

(5)$$ S_{\rm r} = {\displaystyle{(\gamma q) ^{2} \pi \kappa A_{\rm ob}} \over {4k^{2} (T_{\rm m} - T_{\rm s,0}) ^{2}}} $$

where q is the average heat flux (Q divided by Earth's area), A _ob is the area of the ocean basins, κ is the thermal diffusivity and k is the thermal conductivity. The factor γ is meant to take into account the heterogeneity of the surface areal heat flow and the higher mantle heat flux of the ocean basins and below the mid-ocean ridge when compared with the average mantle heat flux. Its value was chosen so to give an areal-spreading rate closer to the observed.

The fixed reference surface temperature value of 0^°C is a good approximation for present day surface temperatures, when the deep oceans are cold. In the future, due to the expected warming of the planet, this approximation may introduce errors of ~ 12% for Q, ~ 8% for S _r and ~ 17% for $P_{{\rm CO}_{2}}$ when T _s = 373 K. However, at this surface temperature, $P_{{\rm CO}_{2}}$ already would be so low that it would not be significant.

In order to determine the area of the ocean basins, the surface of the planet has been divided into continental crust (A _cc) and oceanic crust (the ocean basins, A _ob). The time needed to build the actual continental volume is matter of considerable debate. Most of the models agree that the volume of the continental crust was smaller in the past, but its growing rate greatly varies from model to model. With the caveat that there are large uncertainties in the parameters and a variety of models, we will rely on the formulation of Rosing et al. (Reference Rosing, Bird, Sleep and Bjerrum2010) given by a simple sigmoid function, with zero growth in the future, and assuming a direct correlation between continental crust volume and surface continental crust area:

(6)$$ A_{\rm cc,N} = {\displaystyle{1} \over {1 + {\rm exp} [-2 (t-2)]}} $$

(7)$$ A_{\rm ob,N} = 1 - f_{\rm cc} A_{\rm cc,N} $$

(8)$$ A_{\rm ob} = A_{\rm E} A_{\rm ob,N} $$

where A _ob,N is the normalized ocean basins area, f _cc is the continental crust fraction of Earth's surface and A _E is the surface area of the Earth. With this modelling, the bulk of the continental crust was formed in the period from 3.0 to 1.5 Gyr ago and will be stationed in the actual value in the future.

Water exchange

The water reservoirs in the mantle and on the surface are expected to change with time due to gas exchanges. Estimates of the current abundance of mantle water are in a wide range: from 0.1 to 11 modern oceans of 1.4 × 10²¹ kg (Liu Reference Liu1988; Ahrens Reference Ahrens1989; Jambon and Zimmermann Reference Jambon and Zimmermann1990; Smyth Reference Smyth1994; Dai and Karato Reference Dai and Karato2009; Genda Reference Genda2016; Kurokawa et al. Reference Kurokawa, Foriel, Laneuville, Houser and Usui2018). Some parameterized models of mantle gas exchanges have used values around three modern oceans (McGovern and Schubert Reference McGovern and Schubert1989; Franck and Bounama Reference Franck and Bounama1995, Reference Franck and Bounama1997; Franck Reference Franck1998; Franck et al. Reference Franck, Block, von Bloh, Bounama, Schellnhuber and Svirezhev2000a, Reference Franck, Kossacki, von Bloth and Bounama2002). In our standard model, we assume an initial value of two ocean masses in the mantle and one on the surface.

We model the change rate of the mantle water reservoir with time, dM _mw/dt, using the rate of gain of water due to regassing of water from the surface reservoir, [dM _mw/dt] _r, and the rate of loss of water to the surface due to degassing, [dM _mw/dt] _d, as:

(9)$$ {\displaystyle{{\rm d} M_{\rm mw}} \over {{\rm d}t}} = \left[ {\displaystyle{{\rm d} M_{\rm mw}} \over {{\rm d}t}}\right]_{\rm r} - \left[{\displaystyle{{\rm d} M_{\rm mw}} \over {{\rm d}t}} \right]_{{\rm d}} $$

(10)$$ \left[ \displaystyle{{\rm d}M_{\rm mw}} \over {{\rm d}t}\right]_{\rm r} = f_{\rm bas}\rho_{\rm bas} d_{\rm h} f_{\rm r} S_{\rm r} $$

(11)$$ \left[ \displaystyle{{\rm d} M_{\rm mw}} \over {{\rm d}t}\right]_{\rm d} = \rho_{\rm mw} d_{\rm m} f_{{\rm d}}S_{\rm r} $$

where f _bas is the mean mass fraction of water in the basalt layer of the ocean crust, ρ_bas is the density of basalt, d _h is the depth of hydration in the oceanic crust, f _r is the fraction of water that actually enters the mantle and is not lost in back-arc volcanism (this value is adjusted in order to give the same amount of water we have today), ρ_mw is the density of water in the mantle (the mass of water in the mantle divided by the mantle volume), d _m is the melt generation depth, the depth bellow mid-ocean ridges where water can degass and f _d is analogous to f _r, but for the degassing, the fraction of water in the volume S _rd _m that actually degasses per unit of time.

Our description of d _h and f _d follows Cowan and Abbot (Reference Cowan and Abbot2014):

(12)$$ d_{\rm h} = {\rm min} \left[ d_{\rm h,E}\left( {\displaystyle{P} \over {P_{\rm E}}}\right) ^{\sigma}, d_{\rm b}\right] $$

(13)$$ f_{\rm d} = {\rm min}\left[ f_{\rm d,E}\left( {\displaystyle{P} \over {P_{\rm E}}}\right) ^{-\mu}, 1\right] $$

(14)$$ P = g \rho_{\rm w} d_{\rm oc} $$

(15)$$ d_{\rm oc} = {\displaystyle{M_{\rm ow}} \over {\rho_{\rm w}A_{\rm ob}}} $$

where P is the pressure at bottom of the ocean, P _E is the current value of the pressure, d _b is the thickness of basaltic oceanic crust, d _h,E is the nominal value of the hydration depth, f _d,E is the nominal degree of melt degassing, σ and μ describe the pressure dependence, ρ_w is the density of water, d _oc is the effective depth of the oceans and M _ow is the mass of water in the surface reservoir.

A possible extreme scenario for the gas exchanges is one in which a very large amount of water is degassed by the mantle, increasing sea levels to the point of covering the continents to the top of the highest possible mountain. This case could be very disruptive not only for terrestrial life on the planet but also for the silicate–carbon cycle, because without a considerable land area to promote rock weathering and CO₂ recycling, the planetary thermostat buffering could not work properly.

The depth of the oceans, d _oc, is calculated using the approximation of Cowan and Abbot (Reference Cowan and Abbot2014), in which the tops of continents are at the sea level, that is, the continental freeboard is much smaller than the depth of the ocean basins. For the present-day Earth, this value is d _oc = 3.94 km in our model. Cowan and Abbot (Reference Cowan and Abbot2014) also derived the maximum depth an ocean could have before covering the highest mountain of a planet, which, in the case of Earth, is $d_{\rm oc}^{\rm max} = 11.4\,{\rm km}$, and we will use this value as the upper limit for ocean depth habitability. Our choice of the lower limit follows Abe et al. (Reference Abe, Abe-Ouch, Sleep and Zahnle2011), which considers a depth in the order of tens of centimetres, $d_{\rm oc}^{\rm min} = 0.5\,{\rm m}$, appropriate for a very dry planet but with still enough water to be habitable.

Assuming that the system is closed against other forms of water loss (water loss due to photolysis and hydrogen escape to space will be discussed below), we can calculate the amount of water in the mantle and on the surface.

See Table 2 for a list of the parameters and constants used in the thermal evolution model and in the water exchange model.

Table 2. Parameters and constants for the thermal evolution and gaseous exchange models

B86 = Bevis (Reference Bevis1986), CH99a = Conrad and Hager (Reference Conrad and Hager1999a), CH99b = Conrad and Hager (Reference Conrad and Hager1999b), CA14 = Cowan and Abbot (Reference Cowan and Abbot2014), K06 = Korenaga (Reference Korenaga2006), K10 = Korenaga (Reference Korenaga2010), MS89 = McGovern and Schubert (Reference McGovern and Schubert1989), S80 = Schubert et al. (Reference Schubert, Stevenson and Cassen1980), S81 = Stacey (Reference Stacey1981).

Atmospheric model

The first component of the climate model is the evolution of the energy flux reaching the atmosphere from the Sun. We used the main sequence stellar luminosity parameterization of Rushby et al. (Reference Rushby, Claire, Osborn and Watson2013), considering data of Baraffe et al. (Reference Baraffe, Chabrier, Allard and Hauschildt1998), with dependence on stellar mass, M, and stellar age, t. We had to introduce a small correction to the fit in the form of a multiplicative constant, so we could have exactly one solar luminosity for an one solar mass star at t = 4.57 Gyr:

(16)$$ \eqalign{& \hskip-16.5pt L(M,t) \cr & = 1.00195 \times (-2.245 +0.7376 \times t + 16.03 \times M \cr & -0.02348 \times{ t}^{2}-4.596 \times tM -44.2 \times{M}^{2} \cr & +0.1212 \times{t}^{2}M+10.5 \times t{M}^{2}+59.23 \times{M}^{3}\cr & -0.2047 \times{t}^{2}{M}^{2} -10.43 \times t{M}^{3}-38.59 \times{M}^{4} \cr & +0.1132 \times{t}^{2}{M}^{3}+3.82 \times t{M}^{4}+10.46 \times{M}^{5} )} $$

This parameterization is more complicated than most of the forms for the evolution of solar luminosity in the literature, but closely reproduces the results of most of them over most of the time interval of interest. Significant deviations are found in Kasting et al. (Reference Kasting, Whitmire and Reynolds1993), with a larger increase in luminosity with time, and Gough (Reference Gough1981) and Franck et al. (Reference Franck, Kossacki and Bounama1999), with higher luminosities for times greater than ~8 Gyr. The flux, S, at the Earth position is obtained by dividing the luminosity, L, by the area of a sphere with a radius of one astronomical unit. The present-day value of S is S ₀ = 1368 W m⁻².

Following Lenton (Reference Lenton2000), Franck et al. (Reference Franck, Kossacki, von Bloth and Bounama2002) and Levenson (Reference Levenson2011), we consider a quasi-grey atmosphere (Chamberlain Reference Chamberlain1980; Chamberlain and Hunten Reference Chamberlain and Hunten1990) to obtain the mean surface temperature, T _s:

(17)$$ \sigma{T_{\rm s}}^{4} = {\displaystyle{ (1-A_{\rm p}) S} \over {4}}\left( 1 + {\displaystyle{3} \over {4}} \tau\right) $$

where σ is the Stefan–Boltzmann constant, A _p is the effective planetary albedo, S is the incoming flux at the top of the atmosphere and τ is the optical depth of the atmosphere, which we assume to be separately dependent on the concentrations of the two greenhouse gases considered, H₂O and CO₂:

(18)$$ \tau = \tau \;({\rm H}_{2}{\rm O}) + \tau\; ({\rm CO}_{2}) $$

We assumed the optical depths due to each gas to be a function of its partial pressure, following Lenton (Reference Lenton2000), which used data of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993). However, we decided to make our own fit, because the original fits of Lenton (Reference Lenton2000) were for the narrow temperature range from 0 to 40^°C:

(19)$$ \tau_{{\rm H}_{2}{\rm O}} = 0.0278 ( P_{{\rm H}_{2}{\rm O}})^{0.4079} $$

(20)$$ \tau_{{\rm CO}_{2}} = 1.276 (P_{{\rm CO}_{2}})^{0.602} $$

where $P_{{\rm H}_{2}{\rm O}}$ is the water saturation vapour pressure and $P_{{\rm CO}_{2}}$ is the carbon dioxide partial pressure. The $\tau _{{\rm H}_{2}{\rm O}}$ is more confidently applicable in the range 250 to ~400 K, and the $\tau _{ {\rm CO}_{2}}$ in the range 0.0–2.0 bar. The Clausius–Clapeyron equation gives the water pressure in function of the temperature and relative humidity:

(21)$$ P_{{\rm H}_{2}{\rm O}} = HP_{0} \exp \left[-{\displaystyle{L_{0}} \over {R}}\left( {\displaystyle{1} \over {T_{\rm s}}}-{\displaystyle{1} \over {T_{\rm s,0}}} \right) \right] $$

where H is the mean relative humidity, P ₀ is a reference value for the water vapour saturation curve in the reference temperature and L ₀ is the latent heat per mole of water. The partial pressure of CO₂ is dependent of the formulation of the weathering rate, as we will see below.

The planetary albedo is given by the parameterization of Williams and Kasting (Reference Williams and Kasting1997), a polynomial of the superficial albedo, a _s, the surface temperature, T _s, the cosine of the solar zenith angle, μ = cos(Z) and the CO₂ the partial pressure, $P_{{\rm CO}_{2}}$, which can be confidently applicable in the range 0 < a _s < 1, 190 K < T _s < 360 K, 0 < μ < 1, $10^{-5} \; {\rm bar} \lt P_{{\rm CO}_{2}} \lt 10 \; {\rm bar}$, but is well behaved and can be used beyond that temperature range. For temperatures higher than ~ 460 K we used the planetary albedo function for the Sun of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993). The original fit in Williams and Kasting (Reference Williams and Kasting1997) is divided into two parts at 280 K to reduce errors, but this creates a discontinuity between the two parts. Because 280 K is not so far from the current Earth's average temperature of 287 K, this discontinuity may cause jumps in the results of the planetary albedo and surface temperature. In order to prevent this, we smoothly join the two fits using two sigmoid functions:

(22)$$ A_{\rm p} = A_{\rm p,1}\times {\rm Sigmoid}_{1} + A_{\rm p,2}\times {\rm Sigmoid}_{2} $$

(23)$$ {\rm Sigmoid}_{1} = {\displaystyle{1} \over {1 + \exp [-0.15 (T_{\rm s} -280\,{\rm K} )]}} $$

(24)$$ {\rm Sigmoid}_{2} = 1 - {\rm Sigmoid}_{1} $$

(25)$$ \eqalign{ A_{\rm p,1} & = -6.8910 \times 10^{-1} + 1.0460a_{\rm s} + 7.8054 \times 10^{-3}T_{\rm s} \cr & \quad - 2.8373 \times 10^{-3}P_{{\rm CO}_{2}} -2.8899 \times 10^{-1}\mu \cr & \quad - 3.7412 \times 10^{-2}a_{\rm s}P_{{\rm CO}_{2}} - 6.3499 \times 10^{-3}\mu P_{{\rm CO}_{2}} \cr & \quad + 2.0122 \times 10^{-1}a_{\rm s}\mu -1.8508 \times 10^{-3}a_{\rm s} T_{\rm s} \cr & \quad+ 1.3649 \times 10^{-4}\mu T_{\rm s} +9.8581 \times 10^{-5} P_{{\rm CO}_{2}}T_{\rm s} \cr & \quad + 7.3239 \times 10^{-2}{a_{\rm s}}^{2} -1.6555 \times 10^{-5}{T_{s}}^{2} \cr & \quad + 6.5817 \times 10^{-4}{P_{{{\rm CO}}_{2}}}^{2} + 8.1218 \times 10^{-2}{\mu}^{2} } $$

(26)$$ \eqalign{ A_{\rm p,2} & = 1.1082 + 1.5172a_{\rm s} - 5.7993 \times 10^{-3}T_{\rm s} \cr & + 1.9705 \times 10^{-2}P_{{\rm CO}_{2}} - 1.8670 \times 10^{-1}\mu \cr & - 3.1355 \times 10^{-2}a_{s}P_{{\rm CO}_{2}} - 1.0214 \times 10^{-2}\mu P_{{\rm CO}_{2}} \cr & + 2.0986 \times 10^{-1}a_{\rm s}\mu - 3.7098 \times 10^{-3}a_{\rm s}T_{\rm s} \cr & - 1.1335 \times 10^{-4}\mu T_{\rm s} + 5.3714 \times 10^{-5}P_{{\rm CO}_{2}}T_{\rm s} \cr & + 7.5887 \times 10^{-2}{a_{\rm s}}^{2} + 9.2690 \times 10^{-6}{T_{\rm s}}^{2} \cr & -4.1327 \times 10^{-4}{P_{{\rm CO}_{2}}}^{2} + 6.3298 \times 10^{-2}{\mu}^{2} } $$

Because our model is zero-dimensional, with no latitudinal variations, and results more closely associated with global averages, we used the fits with Z = 60^°. Another problem is that humidity and cloud effects are hard to model. Most of the studies in the literature circumvent this difficulty using a saturated atmosphere (H = 1), a representative fixed surface albedo, a combination of the effects of the surface albedo and clouds adjusted to give the actual surface temperature of 287 K, and no feedback between surface temperature and surface albedo, like the appearance of high albedo materials (ice) in low temperatures. This results in a high surface albedo and a rather low planetary albedo, which is one of the limitations of the present model. Figure 1 compares some planetary albedo functions found in the literature. We can see that, although all of them present a parabolic-like shape, there are considerable shifts in the values of the albedo and this can greatly impact the surface temperature.

Fig. 1. Comparison between planetary albedo functions given in the literature and the fit adopted in this work (equation (22)). Here, the Williams and Kasting (Reference Williams and Kasting1997)'s albedo assumes a _s = 0.2, Z = 60^° and $P_{{\rm CO}_{2}}$ = 0.

Atmosphere escape

Atmospheric escape occurs when enough energy is imparted into the particles allowing them to achieve velocities greater than the escape velocity. The energy input can be done via the thermal energy of the atmosphere; direct absorption of stellar energy, primary by extreme ultraviolet (EUV), in the outer part of a planet's atmosphere; collision with energetic particles of the stellar wind; and kinetic energy from impact of large bodies (Pierrehumbert et al. Reference Pierrehumbert2010). Here, we will consider the loss of water by diffusion in the atmosphere reaching high altitudes where photolysis and escape can occur (Abe et al. Reference Abe, Abe-Ouch, Sleep and Zahnle2011). The maximum diffusion-limited escape flux is given by Abe et al. (Reference Abe, Abe-Ouch, Sleep and Zahnle2011):

(27)$$ F_{\rm dl} = f_{\rm t}(i){\displaystyle{b_{ i{\rm a}} (m_{\rm a} - m_{i}) g} \over {k T_{\rm strat}}} $$

(28)$$ b_{ i{\rm a}} = 1.9 \times 10^{19}{\left( {\displaystyle{T_{\rm strat}} \over {300\,{\rm K}}}\right)}^{0.75}\,{\rm cm}^{-1}\,{\rm s}^{-1} $$

(29)$$ T_{\rm strat} = 167\,{\rm K} { \left[{\displaystyle{S_{\rm eff} (1 - A_{\rm p})} \over {0.316}}\right]}^{0.25} $$

where f _t(i) represents the total mixing ratio of the element i in all forms at the homopause, b _ia is the binary diffusion coefficient of i in air, m _a and m _i are the mean mass of an ‘air molecule’ and the mass of the considered element i, respectively, k is the Boltzmann constant, T _strat is the stratosphere temperature (Kasting et al. Reference Kasting, Whitmire and Reynolds1993) and S _eff is the effective solar flux in dimensionless units (normalized to today flux). We assume that H₂O is the most important carrier of hydrogen and approximate f _t(i) by the mixing ratio of H₂O in the stratosphere given by equation (30).

Our prescription of the mixing ratio of the stratosphere follows figures 1c and 3 of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993), considering two temperature ranges for higher accuracy:

(30)$$f^{\rm strat}_{{\rm H}_{2}{\rm O}} = \left \{\matrix {0.43411 - {\displaystyle{0.43411} \over \displaystyle{1+\exp [ 0.123431 ( T_{\rm s} -375.476\,{\rm K})] }}, \hfill \cr {\rm if} \; {T_{\rm s}} \lt 400\,{\rm K} \hfill \cr 0.896558 - {\displaystyle{0.896558} \over \displaystyle{1 + \exp [ 0.042219 (T_{\rm s}-399.142\,{\rm K})] }}, \hfill \cr {\rm if} \; T_{\rm s} \geq 400\,{\rm K}\hfill }\right.$$

The energy-limited escape flux is given by (Watson et al. Reference Watson, Donahue and Walker1981; Erkaev et al. Reference Erkaev, Kulikov, Lammer, Selsis, Langmayr, Jaritz and Biernat2007; Abe et al. Reference Abe, Abe-Ouch, Sleep and Zahnle2011):

(31)$$ F_{\rm el} = {\displaystyle{\pi {R_{\rm E}}^{3} \varepsilon S_{\rm EUV}} \over {GM_{\rm E} m_{i}}} $$

where R _E is the Earth's radius, ε is the efficiency of the heating of gas by the stellar high energy flux, S _EUV, G is the gravitational constant and M _E is Earth's mass. The value of the efficiency coefficient depends on the spectral region and time, but usually is in the range from 0.1 to 0.6. Values around 0.2 are frequently found (Watson et al. Reference Watson, Donahue and Walker1981; Penz et al. Reference Penz, Micela and Lammer2008; Lammer et al. Reference Lammer, Odert, Leitzinger, Khodachenko, Panchenko, Kulikov, Zhang, Lichtenegger, Erkaev, Wuchterl, Micela, Penz, Biernat, Weingrill, Steller, Ottacher, Hasiba and Hanslmeier2009; Kuramoto et al. Reference Kuramoto, Umemoto and Ishiwatari2013; Zahnle and Catling Reference Zahnle and Catling2017).

The high energy photons in the EUV and X-ray part of the spectrum dissociate H₂O molecules and are absorbed by H and H₂ promoting their thermal escape. However, the exact spectral range of the stellar high energy flux varies from work to work. Here, we used the EUV in the range of 10–120 nm plus the Lyman-alpha line at 121.6 nm. With this definition, we have (Ribas et al. Reference Ribas, Guinan, Gudel and Audard2005):

(32)$$ S_{\rm EUV}(t) = 29.7 \times {t}^{-1.23} + 19.2 \times {t}^{-0.72} $$

where t is the time in Gyr, $S_{\rm XEUV}$ is given in erg s⁻¹ cm⁻² and stabilizes in a plateau for t < 0.1 Gyr. The first power law corresponds to the 10–120 nm range, and the second one to the Lyman-alpha line.

A comparison between F _dl and F _el shows that the limiting factor will be the diffusion for low mixing ratios and high EUV fluxes, but the energy limit will become more important for a high mixing ratio, because is unlikely that all the water will be broken and the hydrogen heated up. We will always choose the lower of the two values: min[F _dl, F _el].

Weathering model

We followed the approach of Walker et al. (Reference Walker, Hays and Kasting1981) for the global carbon cycle and the formulation of the continental weathering flux, W, based on the kinetics of weathering reactions between CO₂ and silicate rocks and runoff estimates with some additional parameters:

(33)$$ {\displaystyle{W} \over {W_{0}}} = B_{\rm f}{\left( {\displaystyle{P_{{\rm H}_{2}{\rm O}}} \over {P^{0}_{{\rm H}_{2}{\rm O}}}}\right) }^{ a}{\left( {\displaystyle{P_{{\rm CO}_{2}}} \over {P^{0}_{{\rm CO}_{2}}}}\right) }^{ b} \exp\left[ - {\displaystyle{E_{\rm a}} \over {R}}\left( {\displaystyle{1} \over {T_{\rm s}}} - {\displaystyle{1} \over {T_{\rm sat}}}\right) \right] $$

Here, B _f is a factor for abiotic dampening of weathering with formulation given in the Biosphere Feedback section, a and b are constants, E _a is the activation energy of the weathering reaction, T _sat is the saturation temperature and the indices 0 indicate present values of the variables, except for CO₂ partial pressure value, which is the pre-industrial level of 0.280 mbar. We did not introduced a dependency on the fraction of exposed dry land because Abbot et al. (Reference Abbot, Cowan and Ciesla2012) indicate that weathering behaviour does not depend strongly on land fraction for partially ocean-covered planets and that even a small land area may be enough to maintain the weathering thermostat. The B _f factor refers to increase of continental weathering due to vascular plants, so the weathering rates before ~400 million years ago should be lower by the inverse of that factor. Weathering rates increase with increasing $P_{{\rm CO}_{2}}$ or T _s because both factors speed up the weathering reactions. The $P_{{\rm H}_{2}{\rm O}}$ term models the change in runoff due to change in the saturation water vapour pressure: the hotter, the wetter; more rain, more weathering.

The values for a and b vary considerably in the literature: 0.3 ≥ a ≥ 1.0 and 0.25 ≥ b ≥ 1.0. We choose conservatively the more frequently used values a = 0.65 and b = 0.5 (Pierrehumbert et al. Reference Pierrehumbert2010).

Following Walker et al. (Reference Walker, Hays and Kasting1981), the average rate of silicate weathering must be nearly proportional to the rate release of CO₂ from the mantle, which is proportional to S _r, to maintain a balanced flow of carbon in the atmosphere–ocean–mantle system. Equalizing W/W ₀ to $S_{\rm r}/S_{\rm r}^{0}$ and isolating $P_{{\rm CO}_{2}}$ give the equation for CO₂ partial pressure:

(34)$$\eqalign {P_{{\rm CO}_{2}} & = P^{0}_{{\rm CO}_{2}} \left\{\left({\displaystyle{1} \over {B_{\rm f}}}\right) \left( {\displaystyle{S_{\rm r}} \over {S_{\rm r}^{0}}}\right) \left( {\displaystyle{P_{{\rm H}_{2}{\rm O}}} \over {P^{0}_{{\rm H}_{2}{\rm O}}}}\right) ^{-a} \right. \cr & \left. \times \exp \left[{\displaystyle{E_{\rm a}} \over {R}} \left({\displaystyle{1} \over {T_{\rm s}}} - {\displaystyle{1} \over {T_{\rm sat}}}\right) \right] \right\}^{{1}\over{b}} } $$

This derivation for $P_{{\rm CO}_{2}}$ is an approximation much simpler than to follow fluxes of carbon through reservoirs, but depends on the condition of equilibrium between the flux of carbon to and out of the atmosphere. We tried to ensure this condition using in our calculations temporal steps larger than the time required to stabilize the carbon–silicate cycle, around 10⁵ years (Volk Reference Volk1987; Colbourn et al. Reference Colbourn, Ridgwell and Lenton2015). This condition may be broken if the fraction of exposed dry land vanishes or equals one, because there would be no land to be weathered and maintain the cycle or there would be no surface water to lubricate the plates and maintain intense plate tectonics.

Biosphere feedback

The influence of life on Earth is evident in the atmosphere through the out of thermodynamic equilibrium high concentrations of O₂ and CH₄, and not only that, but also through more weathering caused by vascular plants pumping more CO₂ into the ground, where chemical reactions take place, changes in humidity and changes in the surface albedo, as in the Daisyworld of Watson and Lovelock (Reference Watson and Lovelock1983). According to the Gaia hypothesis of Lovelock and Margulis (Reference Lovelock and Margulis1974), life may be necessary for a planet to be fully habitable, changing the environmental conditions to make it more stable and habitable in the long run. Here, we adopt a very simple biosphere model (Volk Reference Volk1987) of a primary biological productivity, with Π, the biomass generation rate, depending on T _s, $P_{{\rm CO}_{2}}$ and t, having feedback only on $P_{{\rm CO}_{2}}$ (and indirectly on T _s):

(35)$$ {\displaystyle{\Pi} \over {\Pi _{0}}} = \Pi _{\rm T} \times \Pi _{\rm P}\times \Pi _{\rm t} $$

(36)$$ \Pi _{\rm T} = 1-{\left[ {\displaystyle{ (T_{\rm s}-273\,{\rm K})-T_{\rm prod}} \over {T_{\rm prod}}}\right] }^{2} $$

(37)$$ \Pi _{\rm P} = {\displaystyle{P_{{\rm CO}_{2}}-P_{\rm min}} \over {P_{1/2}+\left( P_{{\rm CO}_{2}}-P_{\rm min}\right) }} $$

(38)$$ \Pi _{\rm t} = {\displaystyle{1} \over {1 + \exp\left[ -50 ( t-t_{\rm vasc}) \right] }} $$

(39)$$ B_{\rm f} = {\displaystyle{\Pi} \over {\Pi _{0}}} \left( 1 - {\displaystyle{1} \over {B}}\right) + {\displaystyle{1} \over {B}} $$

where Π₀ is the recent value of the biological productivity, T _prod is the temperature for the maximum productivity, P _min is the minimum partial pressure of CO₂ needed to maintain the biosphere, P _1/2 is the partial pressure of CO₂ at which the productivity is half the maximum value only considering the pressure term, t _vasc is the epoch when vascular plants appeared and B is a constant describing the abiotic weathering dampening. For clearness, we divided the productivity in three terms; the first one is the temperature dependence, a parabolic-like curve with roots in 0 and 50^°C and maximum in T _prod; the second term is the CO₂ partial pressure dependence, ranging from 0 at P _min to 1 for higher pressures; the third term is a sigmoid function in order that the productivity is zero before the appearance of vascular plants. Of course, the real productivity would be non-zero before t _vasc, but, here, we intend to distinguish the effect of vascular plants on the weathering, and not to model the entire biosphere. By the same reason, we consider 50^°C as the maximum temperature for the biological productivity, even when considering a temperature upper limit for the biosphere as a whole of 100^°C. In the literature, the CO₂ partial pressure lower limit P _min varies from 0.15 mbar for C3 plants (Lovelock and Whitfield Reference Lovelock and Whitfield1982; Caldeira and Kasting Reference Caldeira and Kasting1992; Lenton and von Bloh Reference Lenton and von Bloh2001) to 0.01 mbar in the case of C4 plants (Caldeira and Kasting Reference Caldeira and Kasting1992; Franck et al. Reference Franck, Block, von Bloh, Bounama, Schellnhuber and Svirezhev2000a,Reference Franck, Block, von Bloh, Bounama, Schellnhuber and Svirezhevb, Reference Franck, Bounama and von Bloh2006). We decided to investigate both limits. The parameters used in the atmospheric, weathering and biosphere models are found in Table 3.

Table 3. Constants and parameter of the atmospheric, weathering and biosphere models.

Considering all the criteria above, our planet will be considered habitable if it has an effective ocean depth between 0.5 m < d _oc < 11.5 km, a CO₂ partial pressure higher than 0.15 mbar for C3 plants and 0.01 mbar for C4 plants, and a mean superficial temperature in the range 273 K < T _s < 373 K. As seen below, the last criterion will be the most important one.

Deep past Earth's climate

The general trend is of a hotter and more active geosphere in the past and a colder, more viscous and less active mantle in the future. The parameterization we used implies a more stable heat flux through time than the classical models (Fig. 2(d)), being only ~40% higher 4 Gyr ago than today, compared to more than 2.4 times higher of other models (Franck Reference Franck1998; Sandu et al. Reference Sandu, Lenardic and McGovern2011). This impacts the spreading rate S _r and, indirectly, the $P_{{\rm CO}_{2}}$ and the water exchange.

Fig. 2. Most important variables of the model. Data points are geochemical estimates. Dashed lines are from other geophysical models.

The lower solar luminosity in the Archean would not be enough to maintain warm temperatures with the same atmosphere as today and Earth should be fairly colder at that time. However, some geochemical (Karhu and Epstein Reference Karhu and Epstein1986; Knauth et al. Reference Knauth and Lowe2003; Robert and Chaussidon Reference Robert and Chaussidon2006) and geophysical (Franck et al. Reference Franck, Block, von Bloh, Bounama, Schellnhuber and Svirezhev2000a, Reference Franck, Bounama and von Bloh2006; Schwartzman Reference Schwartzman2002) studies indicate very high temperatures (> 320 K) in the Archean. Since these temperatures are for the sea water, air temperatures (the surface temperature in our model) may have been even higher. As we can see in Fig. 3, the surface temperatures are above 273 K during almost all the Earth's history, but always below the current temperature of 287 K.

Fig. 3. Comparison of surface and water temperature given by models and geochemical estimates. The thick solid red curve is the prediction of our standard model, the shaded red region, the model uncertainty. Dashed, dotted and dot-dashed thin lines are results for geophysical and atmospheric models. Boxes indicate uncertainties of geochemical estimates, arrows, geochemical upper limits. The black dot with a white border marks the recent value of 287 K at 4.5 Gyr. A colour version of this figure is available in the online journal.

Walker et al. (Reference Walker, Hays and Kasting1981) comment that the runoff estimates they used to model the weathering feedback mechanism may be ‘contaminated’ by the effect of biota, promoting more weathering than expected in an abiotic Earth. Schwartzman and Volk (Reference Schwartzman and Volk1989) indicated that the weathering today may be 10, 100 or even 1000 times greater than in the past, before the appearance of vascular plants and specialized microbes. We modelled this possible effect of life on the weathering with a very conservative factor of B = 2. A larger factor would increase CO₂ partial pressures and temperatures in the past, but this would be incompatible with the geochemical data for CO₂ concentrations in the atmosphere (Fig. 4). If CO₂ alone could not account for this effect without contradicting the geochemical estimates, other factors may be at play.

Fig. 4. Comparison of CO₂ partial pressure given by models and geochemical estimates. The thick solid red curve is the prediction of our standard model, the shaded red region, the model uncertainty. Dashed thin lines are results for geophysical and atmospheric models. Data points are geochemical estimates. Boxes indicate uncertainties of geochemical estimates, arrows, geochemical upper or lower limits. The black dot with a white border marks the recent value of 0.280 mbar at 4.5 Gyr. The two lines of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993) define minimum and maximum concentrations to attain, respectively, 273 and 373 K. A colour version of this figure is available in the online journal.

Some recent studies defy the idea of a super-hot Archean Earth, indicating colder maximum temperatures (< 315 K) or even lower ones (Sleep and Hessler Reference Sleep and Hessler2006; Hren et al. Reference Hren, Tice and Chamberlain2009; Blake et al. Reference Blake, Chang and Lepland2010; de Wit et al. Reference de Wit and Furnes2016). Combining this with other possible phenomena that could alter Earth's temperature, for example, surface albedo (Rosing et al. Reference Rosing, Bird, Sleep and Bjerrum2010), cloud properties (Goldblatt and Zahnle Reference Goldblatt and Zahnle2010; Rondanelli and Lindzen Reference Rondanelli and Lindzen2012), total atmospheric pressure (Li et al. Reference Li, Pahlevan, Kirschvink and Yung2009) and, of course, other greenhouse gases that could be present at that time in different concentrations (Pavlov et al. Reference Pavlov, Kasting, Brown, Rages and Freedman2000; Haqq-Misra et al. Reference Haqq-Misra, Domagal-Goldman, Kasting and Kasting2008; Wordsworth and Pierrehumbert Reference Wordsworth and Pierrehumbert2013; Byrne and Goldblatt Reference Byrne and Goldblatt2014) and with not only one mechanism, but a combination of several processes we could relieve this problem (Wolf and Toon Reference Wolf and Toon2014). Solving the problem of Earth's paleoclimate will have also astrobiological implications, because it could help set the outer limit of the habitable zone.

Among the different gases that may contribute to a larger greenhouse effect in the Archean, methane, CH₄, a trace gas today, but a powerful greenhouse gas, could be very important. It is produced by methanogenic microorganisms, and may have reached concentrations of ~500 ppmv (Pope et al. Reference Pope, Bird and Rosing2012), and, as long as the ratio CH₄/CO₂ remained less than about 0.1, it could contribute to a considerable warming of the atmosphere (Haqq-Misra et al. Reference Haqq-Misra, Domagal-Goldman, Kasting and Kasting2008).

Higher surface temperatures in the past due to other factors would cause, in our model, an increase in weathering and a decrease in CO₂ concentrations, which would be more in line with the geochemical estimates around t = 2.0 Gyr, and would allow a B factor larger than 2, as indicated by Schwartzman and Volk (Reference Schwartzman and Volk1989). Mean temperatures low as 273 K may not have transformed Earth in a snowball. Ice-free regions around the equator would be possible at a mean temperature of ~260 K. Only when the ice advances over latitudes lower than ~ 30^°, the feedback between albedo and temperature would cause a global glaciation (Tajika Reference Tajika2003).

However, since our model aims to be a first approximation with only two greenhouse gases, CO₂ and H₂O, a fixed surface albedo and no cloud or albedo feedback, the predicted temperatures above 273 K are satisfactory. In the future of Earth, when CO₂ concentrations are expected to decrease and temperatures to rise due to the increasing solar luminosity, our treatment will be more appropriated.

We did not predict a fast and intense degasification of the mantle in the past as shown by other models (Fig. 5(a)). However, a modest net exchange of water between mantle and surface is present with a maximum of 23% gain of surface water at t = 1.5 Gyr. This is related to the mantle temperature, the mean heat flux and the spreading rate in the sense that a less variable heat flux implies a more stable spreading rate and a mild volatile (water) flux in the past. As a result, the surface water mass is almost constant in the last 0.5 Gyr.

Fig. 5. Evolution of mass and effective depth of the ocean. (a) Comparison of oceans masses from models and geochemical estimates. Our standard model results for the surface reservoir are the thick solid blue line, for the mantle, the thick solid brown line. Dashed blue and red lines refer to hydrogen loss efficiencies ε = 0.2 and 1.0. Dot-dashed thin lines are for geophysical models. (b) Our standard model results for the oceanic effective depth is the thick solid blue line, for hydrogen loss efficiencies ε = 0.2 and 1.0, the dashed blue lines. The labels on the curves stand for: A11 = Abe et al. (Reference Abe, Abe-Ouch, Sleep and Zahnle2011), B01 = Bounama et al. (Reference Bounama, Franck and von Bloh2001), F98 = Franck (Reference Franck1998), F02 = Franck et al. (Reference Franck, Kossacki, von Bloth and Bounama2002), FB97 = Franck and Bounama (Reference Franck and Bounama1997), H09 = Hren et al. (Reference Hren, Tice and Chamberlain2009), P12 = Pope et al. (Reference Pope, Bird and Rosing2012).

Estimates of ocean masses in the past are scarce and based on oxygen isotopic composition, which can be interpreted as variations in water temperature or changes in ocean volume (Pope et al. Reference Pope, Bird and Rosing2012). The later could be caused by water loss due to hydrogen escape or by net flux of water into the mantle from the surface. According to Pope et al. (Reference Pope, Bird and Rosing2012) both causes may be possible in the history of Earth. Because our model does not predict relevant hydrogen escape in the past, all changes in ocean volume would be due to net exchanges of water between mantle and surface.

The increase in ocean depth (Fig. 5(b)) may be more related to the continental crust growth than with the larger ocean mass. As the continental crust swells up, oceans basins are left with a smaller area of the planet. Because the oceans appear do not have covered too much of the continental surface, the ocean basins become deeper to maintain approximately the same volume of water in a region with less surface area, increasing the ocean depth.

Earth's climate hereafter

Urey ratio.

The Urey ratio, Ur, a key parameter for the planetary thermal budget and state, is defined as the ratio of the heat production (due to radiogenic elements) to the surface heat loss: Ur = H/Q in our case. Ur > 1 implies a net heat gain and an increase in temperature in the mantle, and Ur < 1, a net heat loss and a decrease in temperature in the mantle. Our model gives a present Ur value of 0.35 (Table 4). This is considerably lower than the ~0.8 of classical parameterized convective models (Davies Reference Davies1980; Schubert et al. Reference Schubert, Stevenson and Cassen1980; Richter Reference Richter1984; McGovern and Schubert Reference McGovern and Schubert1989; Franck and Bounama Reference Franck and Bounama1995; Sandu et al. Reference Sandu, Lenardic and McGovern2011), but consistent with the ~0.34 of the geochemical estimates (Arevalo Jr et al. Reference Arevalo, McDonough and Luong2009) and more recent parameterized models from the literature (Korenaga Reference Korenaga2006; Lenardic et al. Reference Lenardic, Cooper and Moresi2011; Korenaga Reference Korenaga2012).

Table 4. Value of the most relevant variables for our model at t = 4.5 Gyr, today.

For comparison, the present Urey ratio of Mars is estimated to be ~0.59 (Plesa et al. Reference Plesa, Tosi, Grott and Breuer2015), but uncertainties regarding basic aspects of the planet are so large that it is not possible a much deeper analysis. For most of the other terrestrial bodies of the inner Solar System, there are only inferences of the surface heat flux. Besides the Earth, only the Moon had in situ surface heat flux measurements. With the current InSight mission (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) on Mars, we will have high quality data about the interior of the planet. This may allow us to use our model to study the thermal evolution of Mars and its habitability.

Temperature.

Our model exhibits the trend shown by other models of increasing T _s in the next 1–3 Gyr, predicting T _s = 373 K in 1.63( + 0.14, − 0.05) Gyr from now. The temperature rise could be faster (O'Malley-James et al. Reference O'Malley-James, Greaves, Raven and Cockell2013) or slower (Abe et al. Reference Abe, Abe-Ouch, Sleep and Zahnle2011), but most of the models predict a similar behaviour of the future increase of temperature (Fig. 3). Part of this convergence may be due to the fact that some of them (ours included) are based direct or indirectly on the atmospheric models of Kasting (Reference Kasting1988) and Kasting et al. (Reference Kasting, Whitmire and Reynolds1993).

Our approximation of a mean temperature may mask latitudinal and altitudinal differences in temperature. Some forms of life could survive in cold refuge (poles and high mountains) for millions years more than the rest of the planet, before also succumb due to the high temperatures (O'Malley-James et al. Reference O'Malley-James, Greaves, Raven and Cockell2013). In addition, a ‘deep biosphere’, a subsurface region of a planet able to sustain mild temperatures, may extend the planetary habitability lifetime for some more time until the surface runaway greenhouse leads to temperatures so high that the entire crust will become uninhabitable (Gold Reference Gold1992; McMahon et al. Reference McMahon, O'Malley-James and Parnell2013). Inlight of these other possible habitable regions, our results can be viewed as conservative estimates, with some life still enduring after the epoch when T _s = 373 K.

Since Kasting et al. (Reference Kasting, Whitmire and Reynolds1993) and Kopparapu et al. (Reference Kopparapu, Ramirez, Kasting, Eymet, Robinson, Mahadevan, Terrien, Domagal-Goldman, Meadows and Deshpande2013) used the water loss and runaway greenhouse criteria to estimate the inner limit of the CHZ, we cannot compare our results directly with theirs. However, they present results for the variation of surface temperature in function of solar flux, allowing us to find the solar flux that gives T _s = 373 K. These values are shown in Fig. 6(c). Even if this procedure is only an approximation, the T _s = 373 epoch should occur shortly after the water loss limit and is close to our estimate.

Fig. 6. Comparison of the results from different models about the CO₂ partial pressure limit for C3 plants (a), C4 plants (b) and superficial temperature limit (c). Our results are the upper points in every panel and the uncertainty is the green shaded region. Horizontal bars indicate variable results from varying parameters of the models.

We can see that there are at least two clusters of data. In the first one, our results cluster with the prediction of the collapse of the biosphere due to high temperatures at about 1.5 Gyr. The second one (Schwartzman Reference Schwartzman2002; Abe et al. Reference Abe, Abe-Ouch, Sleep and Zahnle2011) predicts a much more remote collapse, at around 2.6 Gyr. Kopparapu et al. (Reference Kopparapu, Ramirez, Kasting, Eymet, Robinson, Mahadevan, Terrien, Domagal-Goldman, Meadows and Deshpande2013) present an updated version of the model of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993), but it is not clear to us why this implies an extreme scenario of an inner limit of the CHZ located very close to Earth and a very early end of the biosphere when compared to the results of Kasting et al. (Reference Kasting, Whitmire and Reynolds1993).

Most of the models do not explicitly consider clouds and their effects on the albedo and greenhouse effect, which probably push the results to earlier times. In the same way, most of them consider a water saturated atmosphere, which, in view of the difficulty to accurately model humidity in a zero- or one-dimensional model, may reasonably describe the moist and runaway greenhouse, but it still overpredicts the rate of the thermal runaway. Wolf and Toon (Reference Wolf and Toon2014), using a three-dimensional climate model, showed that the surface temperature would increase less rapidly when compared with the results of one-dimensional models. For a 15.5% increase in solar flux, temperature would reach a global mean of 312.9 K and relative humidity would be around only 50%. They did not extend their simulations to higher temperatures, but, as we can see from Fig. 3, the global trend of their results is more consistent with the curve of Abe et al. (Reference Abe, Abe-Ouch, Sleep and Zahnle2011), indicating a late collapse.

Considering our results and those from the literature, it seems difficult to derive a common and convergent estimative for the end of the biosphere due to the rise of temperature. Nonetheless, the end hardly will happen before ~1.5 Gyr. In this way, we may set the thermal end of the biosphere at around 1.6 Gyr, but, as mentioned above, this could be a conservative estimate, and could be pushed to later times.