The end of life on Earth is not the end of the world: converging to an estimate of life span of the biosphere?

Abstract

Environmental conditions have changed in the past of our planet but were not hostile enough to extinguish life. In the future, an aged Earth and a more luminous Sun may lead to harsh or even uninhabitable conditions for life. In order to estimate the life span of the biosphere we built a minimal model of the co-evolution of the geosphere, atmosphere and biosphere of our planet, taking into account temperature boundaries, CO₂ partial pressure lower limits for C3 and C4 plants, and the presence of enough surface water. Our results indicate that the end of the biosphere will happen long before the Sun becomes a red giant, as the biosphere faces increasingly more difficult conditions in the future until its collapse due to high temperatures. The lower limit for CO₂ partial pressure for C3 plants will be reached in 170(+ 320, − 110) Myr, followed by the C4 plants limit in 840(+ 270, − 100) Myr. The mean surface temperature will reach 373 K in 1.63(+ 0.14, − 0.05) Gyr, a point that would mark the extinction of the biosphere. Water loss due to internal geophysical processes will not be dramatic, implying almost no variation in the surface ocean mass and ocean depth for the next 1.5 billion years. Our predictions show qualitative convergence and some quantitative agreement with results found in the literature, but there is considerable scattering in the scale of hundreds of millions of years for all the criteria devised. Even considering these uncertainties, the end of the biosphere will hardly happen sooner than 1.5 Gyr.

Introduction

Life has a long history on Earth, existing for at least 3.7 Gyr (Ohmoto et al. 2004; Noffke et al. 2013; Bell et al. 2015). But, for how much longer is it going to last? In the far future, 7.6 billion years from now, the Sun will enter its red giant phase, increasing in radius and luminosity. The higher luminosity will calcinate the surface of our planet raising the temperatures to make it uninhabitable. Moreover, the Sun will grow larger engulfing the orbits of Mercury, Venus and Earth (Schröder and Smith 2008). Even considering the existence of a ‘deep biosphere’ (Gold 1992; McMahon et al. 2013), the crust may be too hot to be habitable much before that point, certainly marking the end of life on Earth. However, the end of life may come sooner, long before the Sun's exit of the main sequence, due to the increase of solar luminosity with time.

There are quantitative estimates of the biosphere life span since the 1980s. The work of Lovelock and Whitfield (1982) introduced a number of criteria for habitability, for example, a minimum concentration of carbon dioxide in the atmosphere in order to maintain plant photosynthesis, predicting the collapse of the photosynthetic megaflora in about 100 Myr. Later studies investigated the problem with the help of more detailed models, increasing the expected life span of the biosphere to 800 Myr (Caldeira and Kasting 1992), 1.2 Gyr (Lenton and von Bloh 2001) and 1.6 Gyr (Franck et al. 2006).

A more systematic and general approach to the problem uses the concept of the Circumstellar Habitable Zone (CHZ), formalized by Kasting et al. (1993). The location of the CHZ boundaries depends on the atmospheric model and the adopted criteria, but, conservatively, the inner limits are defined by water loss or the runaway greenhouse, and the outer limits are defined by the first carbon dioxide condensation or the maximum greenhouse effect. Basically, this is a ‘follow the liquid water’ approach, focusing on maintaining liquid water on the planetary surface in a sustainable way.

In the inner limit of the CHZ, the increasing solar luminosity introduces more energy into the Earth system and warm up the oceans, increasing vaporization and water vapour content in the atmosphere. Water vapour, being a greenhouse gas, increases the trapping of energy in the atmosphere and the temperatures, promoting more vaporization. Emission increases with increasing temperatures so equilibrium is possible until certain point. When the surface temperature reaches ~340 K, the water vapour content of the lower atmosphere will be ~20%, and loss of water due to photolysis and hydrogen escape may become relevant, and Earth will undergo a moist greenhouse. From there on, a runaway greenhouse is established with a rapid increase in temperature and the complete loss of the oceans (Kasting et al. 1993; Kopparapu et al. 2013).

In our model, the main factor threatening life is the increase of temperature, which is related to the decrease in CO₂ partial pressure (due to weathering), the loss of surface liquid water (due to vaporization and hydrogen escape into space), and natural upper temperature limits for many processes in living cells. Different temperature limits appear in the literature: 340 K (Godolt et al. 2016), 373 K (Franck et al. 2006) and 420 K (O'Malley-James et al. 2013). Here we use 373 K as the upper temperature limit for habitability, which is higher than the moist greenhouse limit of 340 K, but not as optimistic as the 420 K limit, at which even extremophiles could not survive. For average surface temperatures higher than 373 K, liquid water may still be present at higher atmospheric pressures, and, in this case, hyperthermophilic extremophiles could survive. Alternatively, life could be confined to high-latitude or high-altitude cooler temperature refuges. However, we consider that these conditions are too harsh for life in general.

Models for the limits of the CHZ tend to focus on the atmosphere, considering its effects on the temperature and other conditions on the planet, but generally neglect the contribution of the geosphere for the habitability of a planet, which may constrain which atmosphere configurations are possible for a given planet.

For instance, a low mass terrestrial planet like Mars, when compared to a more massive terrestrial planet like Earth or a super-Earth, may not be able to sustain a high heat flux, geological activity and CO₂ flux from volcanoes for billions of years after its formation, making it hard to maintain an intense greenhouse atmosphere and its temperature above 0^°C even if it is inside the classical outer boundary of the CHZ. As another example, a very dry planet, but still with enough water to lubricate plate tectonics and maintain life, may be able to maintain part of its water and its temperature low for longer than an ‘aqua planet’ like Earth would if it is close to the classical inner edge of the CHZ (Abe et al. 2011). In a third example, consider that water can be lost via photolysis and hydrogen escape during the moist and runaway greenhouse, but the secular gaseous exchanges between surface and mantle can increase or decrease the surface reservoir volume of water and have an impact on the amount of water available for life, comparable to the loss due to photolysis (Bounama et al. 2001).

Considering that the CHZ gives a first approximation assessment of the habitability of a planet and corrections will be applied as more information is gathered about the planet and that it can be used as a first screening method to select the best candidates for future studies, its definition needs to be robust and simple. However, it seems that there are few models aiming at investigating the life span of the biosphere and the CHZ combining the co-evolution of the geosphere, atmosphere and biosphere in a proper Earth system treatment, and comparing its predictions in a more systematic way with the results of other studies. Specifically, the geosphere may play a greater role in planetary habitability than that considered in the literature until recently. In addition, as argued by Levenson (2011), we may not always need a multidimensional atmosphere model with detailed radiation codes to estimate a planet's surface temperature. Inclusion of the influence of the geosphere and the geophysical constraints on the definition and evolution of the CHZ may enrich discussion and narrow the uncertainties of the models. Therefore, in this work, we will pay special attention to geophysical processes as constraints for habitability.

In the current paper, we try to contribute to the understanding of the problem of the long-term existence of the biosphere and habitability of Earth by means of a model for the habitability of our planet or alike terrestrial planets considering the co-evolution of the geosphere, atmosphere and biosphere, taking into account the limiting factors such as temperature, CO₂ partial pressure lower limit for C3 and C4 plants and the amount of surface water. At the end, we compare our results with the literature.

Modelling the habitability

Our model combines the results of diverse models – for stellar, geophysical, atmospheric and biosphere evolution – from the literature. Namely, a convection model of the terrestrial mantle (McGovern and Schubert 1989; Franck and Bounama 1995; Franck 1998; Franck et al. 2002; Cowan and Abbot 2014), a quasi-grey atmosphere model (Lenton 2000; Chamberlain 1980), estimates of the variation of planetary albedo, effective solar flux and surface temperatures (Kasting et al. 1993; Williams and Kasting 1997), the loss of water via photolysis (Abe et al. 2011), a weathering model (Walker et al. 1981; Foley 2015) and a simple biosphere model (Volk 1987).

In this work, we try to assess the habitability of a planet in three ways. First, bysurface temperature boundaries. Surface water is present in the liquid state (provided an atmosphere) only for a narrow temperature range. In the same way, life as we know it, cannot endure too high temperatures. Second, by partial pressure of CO₂. Besides its role as a greenhouse gas, CO₂ is used in photosynthesis and may lead to the starvation of megaflora if its concentration in the atmosphere is very low. Third, by quantity of surface water. A minimum amount of water is necessary as it is essential for life, but too much water could hamper the planet habitability because there could be no dry land left for the occurrence of weathering, thus slowing down the carbon–silicate cycle.

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 1985):

(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 2007). 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 2007; Turcotte and Schubert 2014):

(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. (2007) and Arevalo Jr et al. (2009) 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. 2007; Arevalo Jr et al. 2009). See Table 1 for the initial concentrations and other parameters.

Table 1. Parameters for equation (2) using data from Arevalo Jr et al. (2009)

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 (1999a,b) and Korenaga (2006):

(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 (1999a): C ₁ = π^−0.5, C ₂ ≈ 13.5 (Korenaga 2006), and C ₃ ~ 2.5. This parameterization preserves the general shape of the formulation of McGovern and Schubert (1989) and Turcotte and Schubert (2014), 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 2006), 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 (1989) and Franck (1998):

(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. (2010) 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 1988; Ahrens 1989; Jambon and Zimmermann 1990; Smyth 1994; Dai and Karato 2009; Genda 2016; Kurokawa et al. 2018). Some parameterized models of mantle gas exchanges have used values around three modern oceans (McGovern and Schubert 1989; Franck and Bounama 1995, 1997; Franck 1998; Franck et al. 2000a, 2002). 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 (2014):

(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 (2014), 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 (2014) 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. (2011), 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 (1986), CH99a = Conrad and Hager (1999a), CH99b = Conrad and Hager (1999b), CA14 = Cowan and Abbot (2014), K06 = Korenaga (2006), K10 = Korenaga (2010), MS89 = McGovern and Schubert (1989), S80 = Schubert et al. (1980), S81 = Stacey (1981).

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. (2013), considering data of Baraffe et al. (1998), 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. (1993), with a larger increase in luminosity with time, and Gough (1981) and Franck et al. (1999), 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 (2000), Franck et al. (2002) and Levenson (2011), we consider a quasi-grey atmosphere (Chamberlain 1980; Chamberlain and Hunten 1990) 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 (2000), which used data of Kasting et al. (1993). However, we decided to make our own fit, because the original fits of Lenton (2000) 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 (1997), 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. (1993). The original fit in Williams and Kasting (1997) 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 (1997)'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. 2010). 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. 2011). The maximum diffusion-limited escape flux is given by Abe et al. (2011):

(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. 1993) 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. (1993), 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. 1981; Erkaev et al. 2007; Abe et al. 2011):

(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. 1981; Penz et al. 2008; Lammer et al. 2009; Kuramoto et al. 2013; Zahnle and Catling 2017).

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. 2005):

(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. (1981) 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. (2012) 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. 2010).

Following Walker et al. (1981), 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 1987; Colbourn et al. 2015). 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 (1983). According to the Gaia hypothesis of Lovelock and Margulis (1974), 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 1987) 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 1982; Caldeira and Kasting 1992; Lenton and von Bloh 2001) to 0.01 mbar in the case of C4 plants (Caldeira and Kasting 1992; Franck et al. 2000a,b, 2006). 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.

Results and discussion

The coupled equations of the model described above have been solved numerically iteratively until convergence below a preset value of ΔT _s between iterations. In our time variable, the present is set at t = 4.50 Gyr, which is less than the considered age of Earth of 4.54 Gyr (Dalrymple 2011) or the age of the Solar System, 4.57 Gyr (Amelin et al. 2002; Baker et al. 2005). This is because the model would not work well during the accretion period of Earth's formation and the Theia collision with proto-Earth and formation of the Moon, which probably happened during the first 50 to 100 million years of the Solar System formation (Halliday 2000). The solar luminosity function takes this into account through a short delay of 70 Myr: the solar luminosity at t = 0.0 refers to a star aged 70 Myr.

We estimate uncertainties by varying the present mean surface temperature and the pre-industrial CO₂ partial pressure in the range that both may have varied in the last 1–2 million years, which can be regarded as ‘now’ considering the time scales we are using. The upward variations are ~ 2^°C and ~0.020 mbar, and the downward variations are ~ 5^°C and ~0.110 mbar (Hansen et al. 2013; Foster et al. 2017).

First we will discuss the model results for the past of our planet in an attempt to test the model, comparing the results with the geochemical data and geophysical models. Then, we comment about results for today, t = 4.5 Gyr, and about the future of our planet and the fate of the biosphere.

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 1998; Sandu et al. 2011). 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 1986; Knauth et al. 2003; Robert and Chaussidon 2006) and geophysical (Franck et al. 2000a, 2006; Schwartzman 2002) 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. (1981) 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 (1989) 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. (1993) 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 2006; Hren et al. 2009; Blake et al. 2010; de Wit et al. 2016). Combining this with other possible phenomena that could alter Earth's temperature, for example, surface albedo (Rosing et al. 2010), cloud properties (Goldblatt and Zahnle 2010; Rondanelli and Lindzen 2012), total atmospheric pressure (Li et al. 2009) and, of course, other greenhouse gases that could be present at that time in different concentrations (Pavlov et al. 2000; Haqq-Misra et al. 2008; Wordsworth and Pierrehumbert 2013; Byrne and Goldblatt 2014) and with not only one mechanism, but a combination of several processes we could relieve this problem (Wolf and Toon 2014). 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. 2012), 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. 2008).

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 (1989). 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 2003).

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. (2011), B01 = Bounama et al. (2001), F98 = Franck (1998), F02 = Franck et al. (2002), FB97 = Franck and Bounama (1997), H09 = Hren et al. (2009), P12 = Pope et al. (2012).

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. 2012). 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. (2012) 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 1980; Schubert et al. 1980; Richter 1984; McGovern and Schubert 1989; Franck and Bounama 1995; Sandu et al. 2011), but consistent with the ~0.34 of the geochemical estimates (Arevalo Jr et al. 2009) and more recent parameterized models from the literature (Korenaga 2006; Lenardic et al. 2011; Korenaga 2012).

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. 2015), 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. 2013) or slower (Abe et al. 2011), 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 (1988) and Kasting et al. (1993).

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. 2013). 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 1992; McMahon et al. 2013). 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. (1993) and Kopparapu et al. (2013) 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 2002; Abe et al. 2011) predicts a much more remote collapse, at around 2.6 Gyr. Kopparapu et al. (2013) present an updated version of the model of Kasting et al. (1993), 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. (1993).

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 (2014), 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. (2011), 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.

Carbon dioxide concentration.

The fast and recent (in geological scale) anthropic increase in CO₂ atmospheric concentration will bring serious consequences for our species (and others), but the increase is expected to be a short-term phenomenon, and, in the long run, all the models predict a decrease in CO₂ atmospheric levels. In our model, the limit for C3 plants will be reached in 170( + 320, − 110) Myr, followed by the limit for C4 plants in 840( + 270, − 100) Myr. In the same way as for the temperature, these results should be viewed as conservative estimates, because the model provides only a first order approximation for the real levels. Adding to that, since CO₂ levels will drop slowly in the scale of hundreds of millions years, it is reasonable to expect that adaptations by plants and other photosynthetic organisms would allow them to survive some more time.

From Fig. 6, it is clear that the difference in estimates is in the range of hundreds of millions of years and our uncertainties are also of this order of magnitude. If the temperature evolution were close to the early rise in temperatures of Kopparapu et al. (2013) (Fig. 6(c)), the biosphere will have a very short life span of the order of one hundred million years and this should increase weathering and the CO₂ partial pressure fall, thus accelerating the CO₂ starvation of the biosphere.

In the same way as with the temperature evolution, it seems difficult to derive a common and convergent estimative for the CO₂ starvation curve. Nonetheless, the range of ~0.5 to 1.0 Gyr may be a conservative estimate for the C4 limit.

Magnetic field.

In our model we only considered the evolution of the atmosphere in terms of chemical composition (CO₂ and H₂O concentrations) and loss of water through hydrogen escape. Another important phenomenon that could influence over a planetary atmosphere is the erosion due to solar wind. Today, Earth's strong bipolar magnetic field ensures the maintenance of a dense atmosphere (Stadelmann et al. 2010). This was the case in the past, but may not be certain in the future.

The existence of a strong planetary dynamo-created magnetic field is dependent on the existence of a metal liquid core and correlated with the mass and rotational period of the planet (Zuluaga and Cuartas 2012; Zuluaga et al. 2013). Terrestrial planets more massive than Earth may have higher probabilities of maintaining a strong magnetic field and protect their atmospheres if their rotational periods are in the range of 1–3 days. With the same rotational period, terrestrial planets with Earth masses may easily develop a strong magnetic field, but its duration will limited to a couple of billion years. Only rapidly rotating terrestrial mass planets, with a rotational period around 1 day, have the best chances to develop long lived and intense magnetic fields (Zuluaga and Cuartas 2012).

The rotational period of Earth was smaller in the past, being ~20 h 0.90 Gyr ago and ~18 h 2.45 Gyr ago (Williams 2000). An extrapolation from this rate of increase in rotational period (~2 h per Gyr) into the future gives a period of ~28 h in 2.0 Gyr and one of ~35 h in 5.5 Gyr, when the Sun will exit the main sequence. These rotational periods are, at most, 50% longer than today and may still allow to maintain a strong magnetic field in the future. Other phenomena may influence and weaken the magnetic field in the future. For example, plate tectonics may be important in the maintenance of a strong planetary magnetic field (Nimmo 2002; Driscoll et al. 2013), but this phenomenon seems to be a second-order effect and is dependent on other variables, as the volume of surface water. Therefore, we may consider that the atmosphere stripping by solar wind in the future may be comparable to the present one, and so, negligible in a first approximation.

Future oceans.

We could not find any significant change in ocean mass or effective depth due to surface–mantle exchanges in the next ~1.5 Gyr, and so this factor will have little impact in the future habitability of our planet. This is in qualitative agreement with Bounama et al. (2001), but in discordance with their quantitative results and those of Franck et al. (2002), which predicted a decrease of up to a 27% in ocean mass in 1 Gyr. This small variation of water mass comes from the nearly balance of the in (5.8 × 10¹¹ kg yr⁻¹) and out (5.3 × 10¹¹ kg yr⁻¹) fluxes of water of the mantle, consistent with the ~ 10¹¹ kg yr⁻¹ estimates from the literature (Rea and Ruff 1996; Javoy 1998; Pineau et al. 1999; Bounama et al. 2001). Major losses of water would come only after that, when the higher temperatures of the future would displace water from the oceans into the atmosphere, effectively draining the oceans. This process would be accelerated by the hydrogen escape due to water photolysis and the reaction of oxygen with rocks. Lower water content on the surface would drive the recent in-and-out flux equilibrium out of balance and the mantle would start to lose water to the surface. The net flux direction of water would be: mantle→oceans→atmosphere→space. This process may be truncated when the surface reservoir become too low and the plates become less lubricated by surface water. Since we do not model water effects on plate velocity and plate friction, we cannot assess the magnitude of this effect during the period of loss of surface water.

Abe et al. (2011) have considered the possibility of a rapid loss of water and an exit from the trajectory into the runaway greenhouse so that the temperatures would not exceed 100^°C and the planet would remain habitable, even if dried out. A dry planet with a very low relative humidity could resist against higher solar fluxes and stay habitable for longer periods of time (Zsom et al. 2013). In our model, the loss of water is too slow to allow that scenario (Fig. 5), and the oceans would disappear in ~3.0 Gyr from now, close to the estimative of Abe et al. (2011). Increasing the efficiency factor ε for the heating of gas by the EUV from 0.2 to 1.0 accelerates the loss of water, but it still would take ~1.6 Gyr to deplete Earth's oceans, and temperatures could reach 420 K, with sterilizing effects on life. It seems unlikely that Earth will be able to lose water during the moist greenhouse stage fast enough to avoid the runaway greenhouse.

As a consequence, Earth will more probably follow the fate of Venus, a hot sterile world, and not that of Mars, a dry, but not so hot planet. Our world will still orbit the Sun for billions of years, barren of life, witnessing the collision of our Galaxy with the Large Magellanic Cloud in about 2.4 Gyr (Cautun et al. 2018) and the collision with the Andromeda galaxy in about 4.6 Gyr (van der Marel et al. 2019), before being engulfed by an enlarged red Sun in ~7.6 Gyr (Schröder and Smith 2008).

Conclusions

We have built a zero-dimensional model of the co-evolution of the geosphere, atmosphere and biosphere of our planet or alike terrestrial planets to estimate the life span of the biosphere considering as constraining factors: temperature, CO₂ partial pressure lower limits for C3 and C4 plants, and amount of surface water. Special attention was given to the contribution of the geosphere to planetary habitability and its role in constraining atmospheric configurations.

The end of the biosphere will happen long before the Sun goes red giant, with the biosphere facing increasingly more difficult conditions in the future until its collapse, a scenario consistent with that of Franck et al. (2006). The CO₂ partial pressure lower limit for C3 plants will be reached in 170 Myr, followed by the C4 plants limit in 840 Myr, limiting the carbon access for most of life as we know it today. The planet will reach 373 K in 1.63 Gyr, marking the thermal collapse of the biosphere. Some extremophiles may still live in refuges on higher latitudes or deep in the ground, but they will be nothing but ephemerous remnants of a previously exuberant biosphere. Water loss due to internal geophysical processes will not pose a problem for life, since there will be almost no variation in ocean mass or depth in the next one billion years or more. Loss of water due to hydrogen escape may be a greater threat, but a complete loss of the oceans will not occur before the sterilization of the surface due to high temperatures. The planet will go on, sterile, for billions of years, before being engulfed by our Sun.

We have found some convergence between the estimates of the planetary temperature increase trends of the models in the literature, indicating that ~1.6 Gyr is a reasonable conservative estimative for the life span of the biosphere and that the biosphere end will hardly happen earlier than ~1.5 Gyr. Nonetheless, it was not possible to derive a convergent estimate for the life span of the biosphere combining all the three criteria we intended to use. When comparing our predictions and those from the literature, the results vary within a range of the order of hundreds of millions years between different studies.

Regarding future work, one could include as input for the models more data on the past climate of Earth and the faint young sun problem. The integration of different phenomena may increase the reliability of the predictions of the evolution of the factors considered here in a more tightly-knit model. Three-dimensional climate models may be hard to work with, but their insights about the behaviour of the atmosphere regarding clouds and relative humidity may be integrated in zero-dimensional models providing more precise results for the co-evolution of the geosphere, atmosphere and biosphere of our planet, thus narrowing uncertainties in the predictions for the life span of the biosphere.

The role that the geosphere can have in planetary habitability deserves more attention since it can improve strategies for the selection of potentially habitable exoplanets. The parameterizations used in this work are intended to be extended from the Earth case, which is to be considered as a benchmark, to Earth-like exoplanets selected as targets for future surveys of biosphere signatures.