EDP Sciences
Free Access
Issue
A&A
Volume 605, September 2017
Article Number A71
Number of page(s) 21
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201730728
Published online 11 September 2017

© ESO, 2017

1. Introduction

One of the central issues regarding the potential habitability of extrasolar planets is the extent to which plate tectonics is required to maintain habitability (Southam et al. 2015). Plate tectonics is considered to be crucial for maintaining the activity of the carbon-silicate cycle over geological timescales. This helps to stabilize the climate and hence contributes to the habitability of Earth and possibly other planets (e.g. Walker et al. 1981; Kasting et al. 1993b). Indeed, the boundaries of the habitable zone (HZ), i.e. the region surrounding a star where liquid water can be stable on a planetary surface, are traditionally calculated under the assumption that plate tectonics operates effectively (Kasting et al. 1993b; Kopparapu et al. 2014).

On the one hand, even the Earth, where plate tectonics and a carbon-silicate cycle are active, could become uninhabitable by turning into a snowball if the degassing rate of CO2 from the interior becomes too low (Tajika 2007; Kadoya & Tajika 2014, 2015). On the other hand, the very way in which plate tectonics operates on Earth, when it started, and whether it is a stable or a transient feature in the tectonic history of our planet are all complex and still controversial matters (e.g. Tackley 2000; Bercovici 2003; Bercovici & Ricard 2003; van Hunen & van den Berg 2008; van Hunen & Moyen 2012; Gerya 2014; O’Neill et al. 2016). For the numerous super-Earths – large terrestrial extrasolar planets with masses between 1 and 10 M – which have been detected in the past few years (e.g. Batalha 2014), it has proven difficult to establish only on the basis of mass and radius whether plate tectonics is more or less likely to occur than on Earth. While some authors argue for a reduced tendency for plate tectonics to take place on these bodies (O’Neill & Lenardic 2007; Kite et al. 2009; Stamenković et al. 2012; Stein et al. 2013), others favour an increased tendency (Valencia et al. 2007; van Heck & Tackley 2011; O’Rourke & Korenaga 2012) or suggest that the tectonic behaviour of a rocky body can be strongly affected by the specific thermal conditions present after planetary formation and by the particular thermochemical history experienced by the interior (e.g. Noack & Breuer 2014; O’Neill et al. 2016).

Even if an Earth twin with the same mass and radius as our planet and at the same distance from a Sun-like star was detected, it would be very difficult to establish whether plate tectonics exists on such a body or not; finding rocky planets in the HZ is one of the major goals of exoplanetary missions such as the upcoming PLATO 2.0 (Rauer et al. 2014). In the absence of more detailed information and based on the limited evidence of the solar system, the stagnant lid may conservatively be considered as the most common tectonic mode under which terrestrial bodies operate. Because of the strong exponential dependence of mantle viscosity on temperature, the relatively cold upper layers of a rocky body naturally tend to be highly stiff and form a single, immobile plate, i.e. a stagnant lid. Contrary to tectonic plates, this lid does not participate in the dynamics of the mantle and does not allow surface materials to be directly recycled into the deep interior (e.g., Christensen 1984; Davaille & Jaupart 1993; Solomatov 1995). Mercury, Mars, and the Moon have been in a stagnant-lid mode for all or most of their history; Venus, also lacking evidence for plate tectonics at present, may have instead experienced one or more global resurfacing events in the past, possibly associated with episodes of subduction (e.g. Turcotte 1993).

Understanding to what extent a stagnant-lid or plate-tectonic planet may be habitable requires a joint effort involving the modelling of mantle melting and volcanic outgassing and the characteristics of the resulting atmosphere. Kite et al. (2009) investigated the evolution of melting and volcanism over the entire lifetime of the parent star for planets with Earth-like structures and compositions, masses between 0.25 and 25 M, in both a plate-tectonic and stagnant-lid mode of convection. These authors found that while plate tectonics generally allows for melt production and outgassing over the entire stellar evolution, for stagnant-lid bodies these processes tend to cease within a few billion years. Nevertheless, during this time, outgassing rates are predicted to be significantly higher for stagnant-lid than for plate-tectonic bodies and the planet’s mass plays a relatively minor role. O’Rourke & Korenaga (2012) focused specifically on modelling the interior evolution of stagnant-lid planets up to 10 M and argued that bodies more massive than the Earth might escape their stagnant-lid regime by producing a significant amount of negatively buoyant crust. Deep crustal layers of super-Earths, in fact, may lie well in the stability field of eclogite; because it is denser than the average peridotitic mantle, eclogite might promote some form of surface mobilization induced by crustal subduction. Noack et al. (2014) used simulations of finite-amplitude convection to investigate the evolution of CO2 outgassing on hypothetical stagnant-lid planets with Earth-like radius but different core sizes. These authors concluded that the steep slope of the melting temperature characterizing the mantle of planets with large cores (i.e. high average density) strongly prevents the generation of partial melt and, in turn, CO2 outgassing, thereby drastically reducing the potential habitability of this kind of terrestrial bodies.

However, on the basis of interior modelling only, without explicitly coupling the calculated outgassing rates of greenhouse volatiles with climate simulations, it is difficult to assess precisely the potential for habitability of a stagnant-lid (or plate-tectonic) body. Atmosphere modelling studies that calculate the boundaries of the HZ (such as Kasting et al. 1993b; Kopparapu et al. 2013; Wordsworth & Pierrehumbert 2013) usually need to assume a range of atmospheric pressures and concentrations that are plausible from solar system observations. These studies mainly evaluate the atmospheric processes that impact the planetary climate and may limit the habitability of the planet. To calculate the outer edge of the HZ, following Kasting et al. (1993b), it is generally assumed that the planet may provide the appropriate amount of CO2 needed to reach the maximum greenhouse limit via an Earth-like carbon-silicate cycle that stabilizes the planetary climate.

In this work we aim to address the questions of whether and how long a stagnant-lid planet could be habitable, i.e. host liquid water on its surface. To this end, we focus on modelling the coupled evolution of the interior and atmosphere of a hypothetical terrestrial planet with the same mass, radius, and composition as the Earth, but that is characterized by stagnant-lid tectonics throughout its history. We adopt a one-dimensional (1D) model of the thermal evolution of the crust, mantle, and core combined with detailed parameterizations of mantle melting and volatile extraction. The amount of outgassed volatiles (only H2O and CO2 are considered here) then provides the input for our 1D cloud-free radiative-convective atmosphere model that we use to calculate the evolution of the climate of the planet at 1 au from a Sun-like star along with the boundaries of the HZ. In Sect. 2 we present in detail our modelling framework, including the interior thermal evolution model (Sect. 2.1) with the parameterizations adopted to treat mantle melting (Sect. 2.2) and volatile outgassing (Sect. 2.3), the atmospheric model (Sect. 2.4), and the initial conditions and parameters (Sect. 2.5). Simulation results are presented in Sect. 3, where we discuss different aspects related to the evolution of the interior (Sect. 3.1), the outgassing of volatiles (Sect. 3.2), and the resulting atmosphere (Sect. 3.3). Discussion and conclusions follow in Sects. 4 and 5.

2. Theory and models

2.1. Thermal evolution of the interior

We employ a 1D model of parameterized stagnant-lid convection (e.g. Grott et al. 2011b; Morschhauser et al. 2011) to simulate the thermal evolution of the interior of an Earth-like planet over 4.5 Gyr starting from a post-accretion scenario when core formation and magma ocean solidification are completed. Although this parameterized approach cannot capture the complexity of the dynamics of the mantle, it compares well with 2D and 3D simulations of the evolution of stagnant-lid bodies such as Mars and Mercury, both in terms of thermal evolution (Plesa et al. 2015) and crust formation (Tosi et al. 2013).

Given an initial temperature profile for the entire planet (Fig. 1), we solve numerically the time-dependent energy balance equations for the core, mantle, and stagnant lid from which we obtain the evolution of the core-mantle boundary (CMB) temperature, mantle temperature beneath the lid, and thickness of the lid itself, respectively. Energy conservation for the core is given by (1)where t is the time, ρc, cc, and Vc are the density, heat capacity, and volume of the core, respectively, Tc and Ac are the temperature and surface area of the CMB, and qc is the heat flux out of the core into the mantle. For simplicity, we neglect in Eq. (1) the effects of core freezing, i.e. the accompanying release of latent heat of solidification and gravitational potential energy.

thumbnail Fig. 1

Diagram of the interior structure considered for the thermal evolution model and a schematic of the corresponding temperature profile. See text for the description of the various symbols.

Open with DEXTER

The mantle energy balance is governed by the following equation: (2)where ρm and cm are the density and heat capacity of the mantle; Vl and Al are the volume and outer area of the convecting part of the mantle, i.e. between the CMB and base of the stagnant lid (see Fig. 1); St is the Stefan number, which controls the release and consumption of latent heat upon mantle melting and solidification; Tm and Tl are the temperature of the upper mantle and base of the stagnant lid; ρcr, ccr, and Dcr, are the density, heat capacity, and thickness of the crust; is the latent heat of melting; ql and qc are the (parameterized) heat fluxes from the convecting mantle into the stagnant lid and from the convecting core into the mantle; and Qm is the mantle volumetric heating rate. On the right-hand side of Eq. (2), the term proportional to the crust growth rate dDcr/ dt accounts for the additional heat loss due to the transport of melt from the mantle source region to the surface, i.e. for the so-called heat-piping effect (Spohn 1990).

The evolution of the stagnant lid is obtained from the energy balance at its base (e.g. Spohn 1990), i.e. (3)where Dl is the thickness of the lid, km the thermal conductivity of the mantle, and ∂T/∂r the radial temperature gradient calculated at the base of the lid (i.e. at a radius r = Rl). The latter is obtained assuming steady-state heat conduction, i.e. (4)where r is the radial coordinate, kl the thermal conductivity, and Ql the heat production rate in the stagnant lid. Since this generally comprises both the crust and part of the mantle (Fig. 1), kl and Ql are replaced by the corresponding thermal conductivity and internal heating rate (kcr and Qcr for the crust or km and Qm for the mantle) as appropriate. Although neglecting the time dependence in Eq. (4) could affect the earliest transient phases of the evolution, this approximation is sufficiently accurate to capture the long-term thermochemical behaviour of the interior reliably, as demonstrated by comparisons of this approach with the outcomes of fully dynamic simulations (Tosi et al. 2013; Plesa et al. 2015).

The convective heat fluxes from the core into the mantle (qc) and mantle into the stagnant lid (ql) are obtained from boundary layer theory (e.g. Turcotte & Schubert 2002), which is used to determine the thickness of the two corresponding thermal boundary layers from scaling laws appropriate for stagnant-lid convection (Grasset & Parmentier 1998). In particular, the heat flux due to convection in the sublithospheric mantle is proportional to Ra1/3, where Ra is the thermal Rayleigh number defined as (5)with the coefficient of thermal expansion α, the gravitational acceleration g, and the mantle thermal diffusivity κm = km/ (ρmcm). The superadiabatic temperature difference ΔT that drives convection is given by the sum of the temperature drops across the upper and lower thermal boundary layers (6)where ΔTu = TmTl is the temperature drop across the top thermal boundary layer and ΔTd = TcTb that across the bottom boundary layer (Tb is the mantle temperature right above it). In the definition of Rayleigh number (Eq. (5)), the mantle viscosity η is calculated following Hirth & Kohlstedt (2003) assuming an Arrhenius law for wet diffusion creep as follows: (7)where A is a pre-exponential factor, R is the gas constant, is the water concentration in the mantle expressed in ppm (see Sect. 2.2), E and V are activation energy and activation volume for diffusion creep, and Pm is the pressure at the depth at which the upper mantle temperature Tm is calculated; see Table 1 for the values of these parameters and Fig. 2b for typical viscosity profiles calculated with Eq. (7).

Table 1

Main parameters of the interior model.

In order to calculate the temperature difference in Eq. (6), the temperatures at the base of the lid (Tl) and at the base of the mantle (Tb) need to be determined. The latter is readily found by assuming that the mantle is vigorously convecting so that its radial thermal profile is adiabatic and using boundary layer theory to compute the thickness of the bottom thermal boundary layer to the top of which the adiabatic profile extends (see Fig. 1 and Eq. (12)). The lid temperature is obtained instead from scaling laws derived from numerical convection models with strongly temperature-dependent viscosity (e.g. Grasset & Parmentier 1998; Choblet & Sotin 2000). According to these models, Tl can be identified with the temperature at which the viscosity has grown by about one order of magnitude with respect to the viscosity of the convecting mantle. The lid temperature can then be expressed in terms of the mantle temperature and activation energy as (Grasset & Parmentier 1998)(8)with the factor Θ set to 2.9 to account for the effects of spherical geometry (Reese et al. 2005).

thumbnail Fig. 2

Panel a: Solidus temperatures of peridotite for different water contents from dry (grey) to water saturated (dark blue), liquidus temperature (red), and a typical initial temperature profile with Tm = 1700 K (black). Panel b: Viscosity profiles for different water concentrations calculated with Eq. (7) along with the temperature profile shown in panel a.

Open with DEXTER

The convective heat fluxes out of the mantle (ql) and core (qc) are (9)and (10)where dm and db are the thicknesses of the upper and lower thermal boundary layers, respectively. According to boundary layer theory, the first is given by (e.g. Turcotte & Schubert 2002) (11)where Racr is the critical Rayleigh number for the mantle. The second is given by (12)where fc is a factor accounting for the pressure dependence of the viscosity, ηc = η(Tb + Tc)/2 is the viscosity calculated at the average temperature attained in the lower thermal boundary layer, and Rai,cr is the local critical Rayleigh number. Following Deschamps & Sotin (2001), this is given by (13)where Rai is the thermal Rayleigh number for the entire mantle, i.e. (14)with ΔTi = (TmTs) + (TcTb).

With the above equations, the thermal evolution of the interior is obtained by advancing in time a radial temperature profile for the entire planet assuming that the temperature increases conductively in the stagnant lid (see Eq. (4)), linearly in the boundary layers, and adiabatically in the mantle and core.

It is important to note that the surface temperature (Ts) is held constant at 293 K throughout the evolution. Although this may seem inconsistent given that we use an atmospheric model to compute this quantity in response to the time-dependent outgassing of H2O and CO2 (Sect. 2.4), the effects of taking into account the evolution of Ts are negligible for the interior. Even at the highest surface temperatures obtained from the atmospheric model (Ts ≃ 430 K), the temperature-dependence of the viscosity (Eq. (7)) is sufficiently strong to guarantee that the planet never escapes from a stagnant-lid mode. Indeed, simulations conducted by keeping the surface temperature fixed at 450 K throughout the evolution differ by only 3–5% in the main output quantities (e.g. average temperature, crustal thickness, and pressure of outgassed volatiles).

2.2. Mantle melting, crust production, and element partitioning

The generation of partial melts and the accompanying production of crust are calculated by comparing the mantle temperature profile T(r) with the solidus temperature Tsol(r), which defines the temperature above which solid rocks begin to melt. At each radius r where the mantle temperature exceeds the solidus, we determine the local melt fraction φ(r) assuming that it increases linearly between the solidus and liquidus Tliq(r)(15)and compute the volume-averaged melt fraction as (16)where Vφ is the volume of the region where partial melting occurs. Since basaltic melts are expected to become denser than the mantle residue at around 8 GPa (e.g. Agee 2008), we neglect the extraction of partial melts produced at pressures greater than this value.

We consider a peridotitic solidus that we calculate following Katz et al. (2003) according to the assumed initial water concentration of the mantle (see Sect. 2.5). The presence of water depresses the solidus. Figure 2a shows solidus profiles for different water concentrations in the mantle from dry (grey line) to water saturated (dark blue line). For a given value of water concentration, the actual solidus takes the water-saturated solidus as its lower limit.

Upon melting and subsequent melt extraction, incompatible elements, such as H2O and radiogenic elements, are enriched in the crust and depleted in the mantle. In response to the extraction of water with the melt, the solidus tends to increase. We assume Tsol to increase linearly with depletion with a maximum change ΔTsol of 150 K corresponding to the solidus difference between peridotite and harzburgite (Maaløe 2004). The timescale for crustal growth depends on the rate at which undepleted mantle material can be supplied to the partial melt zone, which is a process that takes place according to the mantle convective flow velocity u. The crustal growth rate can be thus calculated as (17)where fp is a constant that describes the fraction of the surface covered by hot plumes in which partial melting takes place (Grott & Breuer 2010). The convective velocity is given by (18)where u0 is the characteristic mantle velocity scale (Spohn 1990). The factor fp in Eq. (17) accounts for the fact that partial melting can either take place in a global sublithospheric channel, where the mantle temperature exceeds the solidus everywhere; this situation corresponds to fp = 1 or in the head of isolated mantle plumes covering a limited portion of the surface, in which case fp< 1 (Grott et al. 2011b). Fully dynamic simulations of large stagnant-lid bodies such as Venus indicate that the planform of mantle convection is likely characterized by a certain number of hot upwellings (e.g. Li et al. 2007; Armann & Tackley 2012; Smrekar & Sotin 2012) with partial melting concentrated in plume heads. In the following, we thus consider a plume model with fp = 0.01 by adding a plume excess temperature ΔTd, corresponding to the temperature drop across the bottom thermal boundary layer, to the mantle temperature used to calculate melt fractions in Eq. (15).

As we show in Sect. 3, for most of the evolution our models are characterized by crust forming at a rate that is faster than the rate at which the stagnant lid thickens as a result of mantle cooling. In this case, we impose that the crust cannot grow thicker than the lid by setting Dcr = Dl, and that, in turn, it is recycled into the mantle by sublithospheric convection, a process that would also be facilitated by the basalt-eclogite transition (O’Rourke & Korenaga 2012).

As mentioned above, we take into account the partitioning of incompatible elements between crust and mantle caused by partial melting. In particular, we consider a model of accumulated fractional melting for the extraction of heat-producing elements (HPEs) and water from the mantle and their enrichment in the crust. The concentration Xliq of a given trace element in the liquid phase can be obtained from its bulk mantle concentration Xm as follows: (19)where δ is an appropriate partition coefficient.

As shown in Fig. 3, for partition coefficients smaller than one, partial melts and, in turn, the crust are strongly enriched in incompatible elements, and even further enriched at small melt fractions. For the long-lived HPEs uranium (U), thorium (Th), and potassium (K), we assumed δ = 0.001 (e.g. Blundy & Wood 2003), while for water, we used δ = 0.01 (e.g. Aubaud et al. 2004).

thumbnail Fig. 3

Ratio of the concentration of incompatible elements enriched in the liquid phase (Xliq) to the corresponding concentration in the solid mantle (Xm) as a function of melt fraction φ assuming fractional melting (Eq. (19)) and different partition coefficients δ.

Open with DEXTER

Knowing the depth-dependent melt fraction from Eq. (15), the average concentration of the various incompatible elements in the melt can be readily calculated as (20)The total mass of the extracted element (Mcr) that is enriched in the crust or, in the case of water, partly enriched in the crust and partly outgassed into the atmosphere (Sect. 2.3) is proportional to the crust production rate, (21)Finally, the concentration of the various incompatible elements in the residual mantle is reduced according to the extracted mass, (22)where M0 and Mm denote the initial and current mass of the silicate mantle, respectively. Also, Eq. (22) can be used to calculated the time-dependent volumetric mantle heating rate that is needed in Eq. (2) as follows: (23)where the index i refers to the four long-lived radioactive isotopes 235U, 238U, 232Th, and 40K, and λi and Hi are their respective decay constants and specific heat production rates that are chosen according to McDonough & Sun (1995). A similar expression holds for the heating rate of the crust, (24)where Xcr,i = Mcr,i/Mcr, i.e. the crustal concentration of the various elements is given by the extracted mass divided by the total mass of the crust.

2.3. Volatile outgassing

2.3.1. Outgassing of H2O

As mentioned in the previous section, the extraction of H2O from the interior is determined self-consistently according to a fractional melting model from which, using Eq. (20), we can calculate the average water concentration in the melt. The buoyant melt percolates from the source region through the lithosphere and crust via porous flow or forming dykes and sills, and eventually part of this melt is extruded at the surface. To calculate the actual amount of water that reaches the surface and can be potentially outgassed into the atmosphere, we thus need to assume a certain value for the ratio of intrusive-to-extrusive volcanism (rie). While this parameter affects the thermal history of the interior only marginally, it can have an important impact on the outgassing evolution as it affects the volume of melt available at the surface in a linear way. However, rie is difficult to constrain as it could vary with time according to the thickness of the lithosphere below which partial melt is generated and would be influenced by the porosity of the upper crust. For simplicity, we assume here rie = 2.5 (Grott et al. 2011b); this value is intermediate between values as low as 1, typical of some basaltic shields, and 5 or even more, characteristic of mid-ocean ridges and other volcanic complexes (White et al. 2006).

Whether the water contained in the extruded melts is outgassed into the atmosphere or retained in the solidifying melt depends on its solubility in surface lavas at the evolving pressure and temperature conditions of the atmosphere (Gaillard & Scaillet 2014), whereby the effect of pressure dominates. Figure 4 shows the solubility of H2O Fig. 4a and CO2 Fig. 4b in basaltic melts as a function of pressure from Newman & Lowenstern (2002). At each timestep of our simulations, we thus check whether the concentration of H2O and CO2 in surface lavas is large enough for these to be supersaturated in the two gases. In this case, the excess concentration is released into the atmosphere, which is then progressively built up, with the partial pressure of water calculated as (25)where is the mass of supersaturated water that can be effectively outgassed.

thumbnail Fig. 4

Saturation concentrations of H2O and CO2 in basaltic melts as a function of pressure. When surface melts are super-saturated, the excess pressure of the corresponding volatile is released into the atmosphere.

Open with DEXTER

The solubility of H2O is much larger than that of CO2 (by more than two orders of magnitude at atmosphere-relevant pressures below 100 bar). As we show in Sect. 3, it turns out that the outgassing of water can be significantly limited by its high solubility in the melts, while all extracted CO2 is easily released into the atmosphere.

2.3.2. Outgassing of CO2

Modelling of CO2 extraction and outgassing is complicated by the fact that carbon is not directly soluble in silicate minerals, but occurs in separate phases depending on pressure, temperature, and oxygen fugacity (e.g. Dasgupta & Hirschmann 2010). Under oxidizing conditions, carbon can be present as solid or molten carbonate, while under reducing conditions it occurs in one of its elemental high-pressure forms, i.e. graphite or diamond. Carbonate sediments form as a consequence of the interaction between atmospheric CO2 dissolved in rainwater and silicate rocks. Therefore, in contrast to the Earth where plate tectonics has been operating for billions of years, in stagnant-lid bodies where surface materials cannot be recycled into the interior via subduction, the mantle is likely to be characterized by relatively reducing conditions throughout its evolution. Indeed, there is evidence that the mantle of the Earth was more reduced for much of the Archean than it is now (e.g. Aulbach & Stagno 2016) and that it has become more and more oxidized in response to hydrogen escape (Catling et al. 2001) and to the subduction of oxidizing agents such as ferric iron, water, and carbonates (Kasting et al. 1993a; Lécuyer & Ricard 1999). In contrast, on Mars and the Moon, two bodies that have been likely characterized by a stagnant lid for most or all of their evolution, basalts are generally reduced with oxygen fugacities (fO2) ranging from about one log10-unit below the iron-wüstite buffer (IW) to one unit above it (e.g. Herd et al. 2002; Wadhwa 2008). Nevertheless, in the case of Mars, the composition of basaltic SNC meteorites, from which such reducing conditions are inferred, differs from that of old surface basalts found in the Gusev crater that require a more oxidized source instead (McSween et al. 2009). An explanation for this discrepancy is that meteorites and surface rocks formed from melting and crystallization of the same source but under different fO2 conditions (Tuff et al. 2013). The composition of the Gusev crater rocks could be then explained if oxidized surface material was recycled into the upper mantle early in the history of Mars, possibly through an episode of subduction (McSween et al. 2003). Our stagnant-lid models assume there is no effective mechanism for surface recycling, thus the oxidation state of the mantle should not change significantly over the evolution.

Consistent with this picture, we make the assumption here that the mantle oxygen fugacity is sufficiently low for carbon to be available in its reduced form, which, at the pressure and temperature conditions where partial melting occurs in our models, is likely graphite, or possibly diamond (e.g. Stagno et al. 2013). Upon partial melting, some graphite dissolves into the melt in the form of carbonate ions and upon migration and extraction is subsequently outgassed as CO2 at the surface. We follow the approach applied by Grott et al. (2011b) to model the outgassing of CO2 on Mars, which in turn draws from the thermodynamic framework of redox melting introduced by Hirschmann & Withers (2008) to assess the solubility of CO2 in graphite-saturated magmas. The abundance of CO2 in the melt () depends on the concentration of carbonate () and on the assumed oxygen fugacity (fO2) (Grott et al. 2011b) as follows: (26)where b is a constant appropriate for Hawaiian basalts, and (27)where KII and KI, also appropriate for Hawaiian basalts, are equilibrium constants governing the reactions forming CO2 from graphite and oxygen and carbonate ions from CO2 (Holloway 1998).

As for H2O and heat producing elements, the average concentration of the CO2 extracted from the mantle () is obtained from Eq. (20) and its mass () by solving Eq. (21). The mass of CO2 that is actually outgassed () is then calculated by comparing with the saturation concentration of Fig. 4b and accounting for the extrusive-to-intrusive ratio of volcanism. Finally, the partial pressure of CO2 delivered to the atmosphere is given by (28)

2.4. Surface and atmospheric temperature

Taking into account the outgassed greenhouse gases H2O and CO2 from the interior, which build up a secondary atmosphere (assuming the planet lost its primordial atmosphere), we employ a 1D radiative-convective, cloud-free, stationary atmospheric model to calculate the resulting atmospheric temperature, pressure, and water content.

This radiative-convective atmospheric model is based on the atmospheric model of Kasting et al. (1984a) and Kasting et al. (1984b). This original model was further improved by Kasting (1988), Kasting (1991), Kasting et al. (1993b), Pavlov et al. (2000), and Segura et al. (2003). It has been adapted to account for H2O- and CO2-dominated planetary atmospheres over wide temperature and pressure ranges (see von Paris et al. 2008, 2010).

The model uses a variable, logarithmic-equidistant pressure grid based on the hydrostatic equilibrium, taking evaporation and condensation of water at the surface into account. The total atmospheric pressure P is given by P = Pback + PH2O, where PH2O is the partial pressure of H2O and Pback is the sum of the other partial pressures of molecules present in the atmosphere.

The energy transport that determines the temperature profile in the model is calculated via convection in the troposphere and radiative transfer in the stratosphere using a time-stepping approach. The convection is described by an adiabatic lapse rate based on Ingersoll (1969) accounting for evaporation and condensation of H2O and CO2 (see e.g. Kasting et al. 1993b; von Paris et al. 2010).

The radiative transfer is separated into two distinct frequency regimes: one for the shortwave radiation incident from the star (0.2376–4.545 μm subdivided into 38 spectral intervals) and one for the longwave radiation emitted by the planet (1–500 μm subdivided into 25 spectral intervals).

Only H2O and CO2 are considered as absorbing species in the thermal infrared regime. The spectral fluxes of the thermal radiation are obtained from the spectral intensity via a diffusivity approximation and assuming an isotropic upwelling intensity. The frequency dependence of the radiative transfer equation in this wavelength region is treated with the correlated-k method (e.g. Mlawer et al. 1997; Goody & Yung 1989). The temperature- and pressure-dependent absorption cross sections used to obtain the k-distributions are calculated with the line-by-line radiative transfer code SQuIRRL (Schreier & Schimpf 2001; Schreier & Böttger 2003) using the Hitemp 1995 database (Rothman et al. 1995). The k-distributions are calculated for temperatures from 100 to 700 K and for pressures from 10-5 to 103bar in nine equidistant logarithmic steps.

In addition to line absorption, collision-induced absorption (CIA) processes in the thermal radiative transfer scheme are also considered. Both self and foreign continua are taken into account. The description of the self and foreign continua of H2O and the foreign continuum of CO2 is based on the approach of Clough et al. (1989; CKD continuum). The implementation of the CKD continuum formulation is adopted from the line-by-line model SQuIRRL. The CO2 self-continuum formulation relies on the approach of Kasting et al. (1984b), which results from weak, pressure-induced transitions of the CO2 molecule near 7 μm and beyond 20 μm. Furthermore, we include the N2-N2 CIA as described in von Paris et al. (2013).

The radiative processes considered in the shortwave region include molecular absorption by H2O and CO2 as well as Rayleigh scattering by H2O, CO2, and N2. The numerical scheme is based on Kasting et al. (1984b) and Kasting (1988) and was improved and described in von Paris (2010) and Kitzmann et al. (2015). The equation of radiative transfer is solved by a δ-Eddington two-stream approximation (Toon et al. 1989). The frequency dependence of the radiative transfer equation in each spectral band is parameterized by a four-term correlated-k exponential sum in each interval (e.g. Wiscombe & Evans 1977). The absorption cross sections for the gaseous absorption of CO2 in the shortwave radiation are taken from Pavlov et al. (2000) using the HITRAN 1992 database (Rothman et al. 1992). Absorption coefficients for H2O in the shortwave part of the radiative transfer are updated based on the HITRAN 2008 database (Rothman et al. 2009).

The cross section σray,i(λ) for the Rayleigh scattering of a specific species i (N2, CO2, and H2O) for the complete spectral range are described by the approach of Vardavas & Carver (1984), which is also used by von Paris et al. (2010) and Kopparapu et al. (2013) as follows: (29)where λ is the wavelength in μm, the conversion factor 4.577 × 10-21 is taken from Allen (1973), Di is the depolarization factor, and r(λ) = [10-5 × Ai(1 + 10-3 × Bi/λ2)] 2 for CO2 and N2, where Ai and Bi are material parameters for the specific molecule i. The values for Di, Ai, and Bi for N2 and CO2 are taken from Vardavas & Carver (1984) and Allen (1973).

For H2O in Eq. (29), the depolarization factor DH2O is 0.17 (Marshall & Smith 1990). The refractivity r(λ) of water is determined by r(λ) = 0.85rdry(λ) (Edlén 1966) with the refractivity of dry air rdry approximated by Bucholtz (1995), (30)

2.5. Initial conditions and model parameters

We ran a series of simulations of the evolution of the interior by varying three main parameters that exert a first-order influence on the outgassing history of the planet. In particular, we varied the initial mantle temperature (Tm,0) between 1600 and 1800 K, the initial water concentration of the mantle () between 0 (corresponding to a dry mantle) and 2000 ppm, and the mantle oxygen fugacity (fO2) between one log10-unit below the IW buffer (IW-1) and two log10-units above it (IW+2).

The range of temperatures is chosen to limit the initial melt fractions to values below ~40%, above which the mantle would no longer deform via viscous creep, but rather exhibit a fluid-like behaviour (e.g. Costa et al. 2009) that the scaling laws we employ to model heat transfer via solid-state convection could not capture. The black line in Fig. 2a shows the initial temperature profile of the uppermost part of the mantle down to 8 GPa for Tm,0 = 1700 K. Figure 2b instead shows three initial viscosity profiles calculated with Eq. (7) along the temperature distribution of Fig. 2a for a dry mantle and for initial water concentrations of 500 and 1000 ppm.

The present-day upper mantle of the Earth is largely depleted in incompatible elements and, based on studies of mid-ocean ridge basalts, also relatively dry with a water concentration of 50–200 ppm (e.g. Saal et al. 2002). Estimates of the bulk Earth water abundance indicate a value within a relatively broad range between 550 and 1900 ppm (Jambon & Zimmermann 1990). Furthermore, planetary formation models predict that a very large amount of water (up to several Earth oceans) can be stored during the accretion of terrestrial planets that form near 1 au (Raymond et al. 2004). Indeed, the H2O storage capacity of nominally anhydrous minerals (olivine and pyroxene) can be as large as few thousand ppm at upper mantle pressures (Hirschmann et al. 2005). By choosing a range between 0 and 2000 ppm, we thus cover a broad parameter space of plausible bulk water contents.

The oxygen fugacity of meteoritic materials is generally low. For ordinary chondrites, for example, it lies about three log10-units below the IW buffer (IW-3). Shergottites, petrologically the most primitive of the Martian meteorites, instead have a slightly higher oxygen fugacity between IW and IW+1 (Righter et al. 2006). By assuming fO2 to vary between IW-1 and IW+2, we thus cover a broad range of redox conditions, from highly to moderately reducing, which leads to a similarly broad range of CO2 outgassing rates (see Sect. 3.2).

In all simulations we assumed that the core is initially super-heated with respect to the mantle by 200 K. The impact of this choice however is not very significant since for internally heated stagnant-lid bodies initial differences in the temperature drop across the bottom thermal boundary are rapidly eliminated by efficient extraction of heat from the core (e.g. Plesa et al. 2015).

The atmospheric model is used to calculate snapshots taking into account the outgassing of the greenhouse gases H2O and CO2 from the interior model and additionally the evolution of the luminosity of the Sun (see Gough 1981). Atmospheric snapshots are calculated in steps of 0.1 Gyr from 0.1 to 0.5 Gyr and in steps of 0.5 Gyr from 0.5 to 4.5 Gyr. The model considers molecular nitrogen (N2), CO2 and H2O as atmospheric gases. These are key component gases in the atmospheres of the terrestrial planets in our solar system. The concentration profile for N2 is an isoprofile of 1 bar. The concentration of CO2, which is also handled as an isoprofile in the model, results from the partial pressure of CO2 (PCO2) from the interior model divided by the background pressure Pback. The atmospheric profile of H2O is determined by assuming a completely saturated atmosphere (i.e. a relative humidity of 100%), such that the partial pressure of H2O (PH2O) is determined by the (temperature-dependent) saturation vapor pressure. Therefore, this approach is limited to temperatures below the critical point of water (647 K). The water reservoir is limited by PH2O outgassed from the interior. We assume an amount of N2 of 1 bar that is similar to the reservoirs found on Earth and Venus and, in addition, facilitates the comparison with other habitability studies, for example, by Kasting et al. (1993b).

We account for the solar evolution by increasing the luminosity of the Sun with time L(t) following Gough (1981), (31)when tt0. The value L0 is the present solar luminosity and t0 is the main-sequence lifetime of the Sun, which is 4.7 Ga. Additionally we vary the solar flux to determine the HZ boundaries. The solar input spectrum is based on a high-resolution spectrum of the Sun by Gueymard (2004). The mean solar zenith angle of 60° is used in the calculations.

The planetary gravity acceleration is assumed to be the same as for Earth. Furthermore, a surface albedo of 0.22 is assumed, which is larger than the mean observed value of present Earth (0.13). Following the approach of Kasting (1988) and Segura et al. (2003), this high value is used to mimic the reflectivity of clouds in the planetary atmosphere. When assuming a relative humidity profile comparable to that of the Earth (Manabe & Wetherald 1967), this albedo leads to the global mean temperature of the Earth (288 K). Assumptions about the relative humidity and surface albedo can have a large impact on the planetary climate as shown, for example, by Godolt et al. (2016). We chose a relative humidity similar to those in Kasting et al. (1993b) and Kopparapu et al. (2013). This allows for a better comparison to these previous studies.

The orbital distance d of the planet to the star is calculated as (32)where S is the solar flux used in the model and S0 the solar constant (i.e. the solar flux at Earth’s orbit).

The important parameters used in the following computations are summarized in Table 2.

Table 2

Parameters of the atmospheric model.

3. Results

3.1. Thermochemical evolution of the interior

The thermochemical evolution of the mantle is summarized in Fig. 5, where we show the evolution of the mantle temperature (Fig. 5a), corresponding viscosity (Fig. 5b), crustal thickness, depth extent of the melt zone and melt fraction (Fig. 5c), and concentration of water in the mantle (Fig. 5d) for a series of simulations with different initial water concentrations and temperatures. The various colours refer to initial H2O concentrations from 250 ppm (grey) to 1000 ppm (dark blue), while no colour distinction is used to indicate the various initial temperatures. The thick black line in the four panels describes the evolution of a reference model characterized by intermediate values of the water content ( ppm) and initial mantle temperature (Tm,0 = 1700 K). In all models, the mantle oxygen fugacity is set to the IW buffer. The latter, however, albeit fundamental for the outgassing of both CO2 and H2O, has only a secondary effect on the evolution of the interior (see Sect. 3.2).

As shown in Fig. 5a, the thermal history is generally characterized by an initial heating phase during which convective cooling is not efficient enough to remove the internal heat generated by the decay of radioactive elements, a behaviour that is characteristic of the early evolution of the interior of stagnant-lid bodies (e.g. Morschhauser et al. 2011; Tosi et al. 2013). After this phase, which lasts between 500 and 1500 Myr depending on the model parameters, the mantle and the core (not shown) cool at a roughly constant rate of ~40 K/Gyr. The thermal history is largely controlled by the choice of the initial water concentration of the mantle; the higher is the latter, the lower the mantle viscosity (see Eq. (7) and Fig. 5b). A low viscosity causes convection, and hence heat loss, to be more efficient (Eq. (14)) with the consequence that models that have the highest initial water content – and hence the lowest reference viscosity – are also characterized by the shortest heating phase and, in turn, by the lowest temperature at the end of the evolution (see curves with ppm in Fig. 5a). The effect of the initial mantle temperature on the overall thermal evolution of the interior is minor compared to that of water. Because of the strong exponential dependence of the viscosity on temperature, in fact, an increase in the mantle temperature is accompanied by a viscosity reduction (Eq. (7)) that promotes convection and leads to a more rapid heat loss. On the contrary, upon cooling, the viscosity increases, slows convection down, and renders heat transfer less efficient. Relatively small changes in the initial temperature produce thus large variations in the heat flux. As a consequence, for a given choice of , the temperature is buffered at a nearly constant value, as expected according to the so-called Tozer effect (Tozer 1967). Models with different initial temperatures tend thus to converge rapidly and evolve in a similar fashion. Furthermore, the higher is, the earlier differences in the initial temperature tend to be removed.

thumbnail Fig. 5

Evolution of the interior for different initial mantle temperatures between 1600 and 1800 K, initial water concentrations () between 250 and 1000 ppm, and an oxygen fugacity corresponding to the IW buffer. Panel a: Mantle temperature is shown; b) corresponding mantle viscosity is shown; c) thickness of the crust (solid lines), the stagnant lid (dashed lines), and distribution of the melt zone and melt fraction (coloured area) are shown; and d) water concentration in the mantle is shown. The different colours, from grey to dark blue, indicate increasing initial water concentrations; the thick black lines in the four panels and the coloured melt zone in panel c indicate the reference model (ref.) with Tm,0 = 1700 K, ppm, and fO2 at the IW buffer. The different initial mantle temperatures, which can be identified in panel a, are not indicated with distinct colours.

Open with DEXTER

Since the solidus temperature strongly depends on the hydration state of the mantle (Fig. 2a), the initial H2O concentration also has a fundamental influence on the production of partial melt and on the formation of crust (and, in turn, on the outgassing history as we discuss in Sect. 3.2). As shown in Fig. 5c, the initial phase of mantle heating leads to the production of a large volume of partial melt, which causes the crust to grow the more rapidly the higher is. With the exception of some models with low Tm,0 and , in all cases the crust (solid lines in Fig. 5c) stops growing when it becomes as thick as the stagnant lid; there is a sharp bend in the crust curves when they cross the stagnant-lid curves indicated with dashed lines. When this happens, the erosion of the bottom part of the crust starts (Morschhauser et al. 2011) and continues until the rate at which the stagnant lid thickens because of mantle cooling overcomes the rate at which crust is produced. In our reference model, this phase lasts between ~800 and 3300 Myr (thick black lines in Fig. 5c). Nevertheless, melt and crust production continue over the entire evolution of the mantle, albeit at a decreasing rate. Neglecting the effects of the possible nucleation of an inner core could affect, at least in principle, the evolution of melt generation. The release of latent heat and of heat associated with the change in gravitational potential energy upon core freezing could slow down the cooling of the core with respect to the cooling of the mantle. The accompanying increase in the temperature drop across the bottom thermal boundary layer could thus favour the formation of plumes and the production of partial melt, which could potentially lengthen the degassing lifetime of the planet. Although the difference in the cooling rates of mantle and core caused by inner-core freezing is likely to be small (Grott et al. 2011a), this effect should be carefully quantified in the future.

Upon melting, incompatible components are removed from the mantle and extracted into the crust or outgassed at the surface. Figure 5d shows the evolution of the mantle water concentration. As expected for a mantle undergoing partial melting over its entire evolution, the concentration of water decreases continuously over time, reaching about 80% of its initial value after 4.5 Gyr. Similar to the crustal growth rate, the rate of water extraction diminishes at the time when subcrustal erosion starts as a consequence of the recycling of wet crust into the mantle; see e.g. the change in slope of the curve corresponding to the reference model at ~800 Myr in Fig. 5d.

3.2. Outgassing evolution

For a subset of the models discussed in the previous section (only those with Tm,0 = 1700 K), Fig. 6 shows the main features of the extraction and outgassing evolution of H2O and CO2. The possibility that the two volatiles are released into the atmosphere depends on whether their concentration in the melts that are extruded at the surface is higher than their saturation concentration at the evolving pressure conditions of the atmosphere (Sects. 2.3.1 and 2.3.2). Figure 6a shows the evolution of the concentration of H2O in surface melts (solid lines) and the evolution of the saturation concentration (black dashed line for the reference model only). The latter is greater than zero already at the beginning of the evolution because of the background pressure of 1 bar N2 used in the atmospheric model (see Sect. 2.4).

thumbnail Fig. 6

Extraction and outgassing evolution of H2O and CO2 for models with initial water concentrations () between 250 and 1000 ppm, an initial mantle temperature of 1700 K, and an oxygen fugacity at the IW buffer. Panel a: Concentration of H2O in the extracted melt (solid lines) and, for the reference model only, concentration of water at saturation (dashed line) are shown; b) average melt fraction is shown; c) partial pressure of outgassed H2O (also expressed as equivalent ocean thickness) is shown; d) partial pressure of outgassed CO2 is shown.

Open with DEXTER

Since the extraction of water is calculated according to a fractional melting model, the evolution of its surface concentration can be readily explained in terms of the evolution of the average melt fraction shown in Fig. 6b. In general, as the degree of partial melting increases because of mantle heating, the concentration of water that partitions into the melt decreases, and vice versa (Eq. (19) and Fig. 3). After a short initial transient phase during which the average melt fraction decreases, it rises over a time interval whose duration is controlled by the mantle viscosity, or, indirectly, by the assumed initial water concentration: the higher (i.e. the lower the mantle viscosity) is, the shorter is this interval (Fig. 6b). The water concentration in the surface melts evolves accordingly; after the initial transient phase, this water concentration first decreases and then rises as the mantle cools and melt fractions decline. Outgassing of H2O can only take place when the water concentration in the surface melts is larger than the saturation concentration, which increases monotonically with time according to the amount of volatiles that are progressively released into the atmosphere. In our reference model, water can be outgassed during the first billion years of evolution and later than about 2.5 Gyr (compare solid and dashed black lines in Fig. 6a). As a consequence, the partial pressure of H2O rises to ~6 bar during the first outgassing phase; this pressure remains constant until about 2.5 Gyr while surface melts are undersaturated in water and then rises again to reach ~9 bar at the end of the evolution, which would correspond to a global water layer of ~90 m (black line in Fig. 6c). Despite the high values of water concentration in the surface melt that are achieved during the second outgassing phase, the increase in the partial pressure of atmospheric water is relatively small compared to the first phase because the overall volume of partial melt produced declines significantly with time (Fig. 5c) and because the continuous increase of the total atmospheric pressure makes further degassing, of H2O in particular, increasingly difficult.

The outgassing evolution of CO2 is somewhat simpler than that of H2O for two reasons. On the one hand, CO2 outgassing is not calculated on the basis of fractional melting but depends on the assumed oxygen fugacity of the mantle (Sect. 2.3.2). On the other hand, since the saturation concentration of CO2 in surface melts is much lower than that of H2O (see Fig. 4), CO2 tends to be outgassed much more easily. Indeed, its surface concentration (which is not shown in Fig. 6) remains above the saturation level of the melt throughout the evolution. As a consequence, the partial pressure of outgassed CO2 rises monotonically as long as partial melt is produced, which, in these models, occurs over the entire evolution. In our reference model, in which we assumed an oxygen fugacity at the IW buffer, it reaches ~1.8 bar after 4.5 Gyr (Fig. 6d), although the increase in pressure after ~3 Gyr is minor owing to the small amount of partial melt produced during the last part of the evolution.

The effect of different initial water concentrations is significant for the outgassing history of both H2O and CO2. As expected, the higher is, the higher are the final partial pressures of the two gases. In the case of CO2, this is because a high water concentration in the mantle causes a strong decrease of the solidus temperature, which facilitates the production of large volumes of partial melt. In the case of H2O, in addition to the above reason, a high initial concentration clearly makes the amount of water available for outgassing high as well with the consequence that the final partial pressure of outgassed H2O ranges from 2.5 bar for ppm to 27 bar for ppm.

thumbnail Fig. 7

Evolution of the partial pressure of outgassed H2O for a model with Tm,0 = 1700 K, fO2 at the IW buffer, and two initial water concentrations ( and 1000 ppm) calculated by neglecting the evolving saturation concentration (dashed lines) and taking such a concentration into account (solid lines).

Open with DEXTER

As already discussed above and shown in Fig. 6c, the outgassing history of H2O is strongly influenced by the evolution of the concentration of water in surface melts. In general, we observed an initial phase during which these are largely supersaturated and H2O is efficiently outgassed. The duration of this phase is inversely proportional to the initial water concentration and varies between 800 and 1200 Myr for ranging from 1000 to 250 ppm. Afterwards, because of the increase of the saturation concentration from the increasing atmospheric pressure caused by the accumulation of both H2O and CO2, H2O outgassing slows down (for and 1000 ppm) or stops entirely (for and 500 ppm) for a period whose duration is also inversely proportional to . In Fig. 7 we show how significant the effect of evolving saturation conditions at the surface can be for the outgassing history of water. For two initial water concentrations of 500 and 1000 ppm, the figure illustrates the evolution of the partial pressure of H2O outgassed into the atmosphere obtained when the effect of the evolving saturation concentration is neglected (dashed lines) and taken into account (solid lines). In the latter case, the amount of outgassed water is significantly smaller throughout the evolution. As a result, after 4.5 Gyr, PH2O is only 43% (for ppm) and 32% (for ppm) of the partial pressure of H2O that would be achieved if all water extracted at the surface were outgassed into the atmosphere.

Water outgassing is suppressed even more strongly as more and more oxidizing conditions are assumed for the mantle. A higher oxygen fugacity leads in fact to a higher amount of outgassed CO2 whose increasing pressure tends to prevent the concentration of H2O in surface melts from exceeding its saturation level. Figure 8a shows this effect for a model with Tm,0 = 1700 K, ppm, and fO2 ranging from IW-1 to IW+1. The corresponding evolutions of the outgassed CO2 are shown in Fig. 8b. For example, for fO2 at IW+1 (brown lines in Fig. 8), about 3 bar CO2 are outgassed within the first 750 Myr. This pressure is sufficient to completely stop the release of water whose partial pressure reaches ~4.7 bar at this time. Melts again become supersaturated (see also Fig. 6a) and additional 0.4 bar H2O can be outgassed only after ~4 Gyr.

thumbnail Fig. 8

Evolution of the partial pressure of outgassed H2O a) and CO2b) for a model with Tm,0 = 1700 K, ppm and oxygen fugacities ranging from one log10-unit below the IW buffer (IW-1) to one log10-unit above it (IW+1). The black line indicates the reference model.

Open with DEXTER

In order to better appreciate the effects of the model parameters on the outgassing of the two volatiles, we plotted in Fig. 9 the partial pressure of H2O (Fig. 9a) and CO2 (Fig. 9b) reached after 4.5 Gyr as a function of the initial water concentration between 100 and 2000 ppm and the oxygen fugacity between one log10-unit below and two log10-units above the IW buffer. As long as the mantle oxygen fugacity is relatively low, the amount of outgassed CO2 is very limited; for an oxygen fugacity not higher than the IW buffer, the maximum partial pressure of CO2 at the end of the evolution does not exceed 1 bar (Fig. 9b). Under these conditions, the maximum partial pressure of water increases linearly according to the initial water concentration in the mantle with ~1 bar outgassed for ppm and 65 bar for an extreme initial concentration of 2000 ppm (Fig. 9a). The partial pressure of outgassed CO2, however, increases roughly linearly with the mantle oxygen fugacity, with tens to hundreds of bars of CO2 at IW+1 and IW+2, respectively. Already for fO2 between IW and IW+1, corresponding to the presence of a few bars of CO2 in the atmosphere, PH2O increases with more slowly because of the increased solubility of water in surface lavas. For fO2 values higher than one log-unit above the IW buffer, the maximum H2O pressure is limited to 10–15 bar at most, even for the highest values of . On the other hand, as shown in Fig. 9b, the final partial pressure of CO2, which is much less soluble than water in basalts (Sect. 2.3.2), is largely determined by the choice of the oxygen fugacity. The initial concentration of water only affects significantly the amount of outgassed CO2 for values below ~400 ppm. In contrast to the enrichment of heat sources and water in partial melts, which is calculated in dependence of a partition coefficient (Eq. (19)), the concentration of CO2 is directly proportional to the melt fraction (Eq. (20)) and depends on the assumed oxygen fugacity (Eqs. (26) and (27)). At low water concentrations the solidus is relatively high and melt fractions are small (Fig. 2a). As a consequence, the effect of the latter on the extraction of CO2 is more significant than that of fO2. For ppm, instead, melt fractions become less relevant and CO2 outgassing is completely controlled by the oxygen fugacity.

thumbnail Fig. 9

Maximum partial pressures of outgassed H2O a) and CO2b) after 4.5 Gyr of evolution as a function of the initial water concentration in the mantle () and oxygen fugacity (fO2) for an initial mantle temperature of 1700 K.

Open with DEXTER

3.3. Atmospheric evolution

The interior evolution and resulting outgassing of CO2 and H2O into the planetary atmosphere lead to an evolution of the planetary climate. Since CO2 and H2O are important greenhouse gases, their abundance has a strong impact on the surface temperatures. Figure 10 shows the evolution of the surface temperatures of Earth-like stagnant-lid planets for different initial mantle concentrations of H2O (Fig. 10a) and oxygen fugacities (Fig. 10b). For these scenarios, the planet is located at an orbital distance of 1 au around the evolving Sun.

thumbnail Fig. 10

Evolution of the surface temperature for a: an oxygen fugacity at the IW buffer and different initial mantle concentrations of water (in different shades of blue), b: an initial mantle concentration of water of 500 ppm and different oxygen fugacities (in different shades of yellow and red) of an Earth-sized stagnant-lid planet located at 1 au around the evolving Sun. The reference scenario (Tm,0 = 1700 K, , and fO2 at the IW buffer) is indicated in black. The freezing point of water (273 K) is indicated by the dashed black line in both panels. The temperatures in panel a correspond to the outgassing evolutions of H2O and CO2 plotted in Figs. 6c and d, while the temperatures of panel b correspond to the outgassing evolutions of Figs. 8a and b.

Open with DEXTER

The stagnant-lid planets show habitable surface conditions over most of their history. Only during the early evolution up to 500 Myr, some scenarios – with relatively low initial mantle concentrations of water or low oxygen fugacities – show temperatures below 273 K. We use the term habitable for conditions in which the global mean surface temperature is above the freezing point of water (273 K) but still low enough to allow for liquid water on the surface of the planet.

The surface temperature evolution is controlled by the outgassing of CO2 and H2O from the interior and the increase in solar luminosity. The H2O outgassed from the interior serves as a water reservoir. The amount of water vapour in the atmosphere then depends on the surface temperature of the planet. As shown in Fig. 10a for an oxygen fugacity at the IW buffer, the initial mantle concentration of water has only a small effect on the atmospheric evolution since, at these temperatures, the water vapor concentrations in the atmosphere are lower than the outgassed water reservoir. For an initial mantle water concentration of 250 ppm, however, the outgassing of CO2 is much slower (see Fig. 6), which leads to a slower increase in surface temperatures. While at the lower temperatures of the atmosphere, which can be found during the early evolution, the CO2 abundance is important, it shows a negligible effect at later stages; the greenhouse effect of water vapor is so strong that the difference in CO2 of about 0.7 bar arising from the use of different initial values of the mantle water concentration (see Fig. 6d) does not exert a significant effect.

For different oxygen fugacities, the amount of CO2 outgassed into the atmosphere varies much more significantly among the various scenarios (Fig. 8). As shown in Fig. 10b, this leads to a larger impact upon the evolution of the surface temperatures with temperature differences up to 75 K after 4500 Myr. We performed atmospheric calculations for scenarios with oxygen fugacities up to IW+1, which leads to atmospheric CO2 partial pressures of about 20 bar. For higher partial pressures, the model boundary condition of a zero downwelling infrared radiation at the top-of-the-atmosphere becomes invalid. For these higher partial pressures we experienced difficulties in reaching a converged solution in radiative equilibrium.

Figure 11 shows the evolution of the partitioning of the water outgassed from the interior into the atmosphere and the ocean for the reference scenario with an initial mantle concentration of water of 500 ppm, an initial mantle temperature of 1700 K, and an oxygen fugacity at the IW buffer. Most of the water is stored in the ocean and the atmospheric abundance defined via the saturation vapor pressure is comparatively low at the surface temperatures obtained for an orbital distance of 1 au. For most of the evolution, the majority of the oceanic water can be liquid, except for the first 100 Myr, in which surface temperatures below 273 K are found and thus at least part of the oceanic water can be expected to be solid. After 4500 Myr, a water amount of about 9 bar, which corresponds to about 3% of an Earth ocean, would form an ocean on the stagnant-lid planet for this scenario, which corresponds to an equivalent depth of about 85 m.

thumbnail Fig. 11

Evolution of the partition of the outgassed H2O (black) between atmosphere (red) and ocean (blue) for the reference scenario with an initial mantle concentration of water of 500 ppm, an initial mantle temperature of 1700 K, and an oxygen fugacity corresponding to the IW buffer.

Open with DEXTER

thumbnail Fig. 12

Determination of the inner a) and outer b) HZ boundaries in solar irradiation S via evaporation of the entire water reservoir (PH2O,out), and the criteria of a minimum temperature of 273.15 K (dashed black line in panel b for an oxygen fugacity at the IW buffer and initial mantle water concentrations between 250 and 1000 ppm (indicated in different shades of blue). Owing to similar behaviours of the atmospheres, the model results are overlaid.

Open with DEXTER

Figure 12 shows how we determined the boundaries of the HZ for our stagnant-lid planets. For the inner boundary, the insolation has been varied to find the amount of irradiation at which the surface temperature is so high that the complete water reservoir outgassed from the interior is in its vapor phase according to the saturation vapor pressure of water (Fig. 12a). It can be inferred that although the amount of water to be evaporated increases when the initial mantle concentration of water is increased, the difference in insolation to evaporate the water reservoir completely is small. At these temperatures, the water vapor feedback is very efficient, leading to a strong increase in surface temperature and atmospheric water content for a small increase in solar irradiation. For the determination of the inner boundary of the HZ, we use the insolation where at least 90% of the outgassed amount of H2O is evaporated into the atmosphere. For the determination of the outer boundary of the HZ, the insolation was varied to find the irradiation that results in a global mean surface temperature of 273 K for the given CO2 amount outgassed from the interior. Figure 12b shows this procedure for different initial mantle water concentrations and an oxygen fugacity at the IW buffer. H2O in the atmosphere is determined via the surface temperature, which is set to 273 K, and therefore does not vary.

thumbnail Fig. 13

Evolution of the inner (in red) and outer (in blue) boundary of the HZ in a) unit of solar constants (SSun) and b) HZ boundary evolution in orbital distance (au) for the reference scenario. The impact of the increasing solar luminosity can be inferred from comparing the HZ distance evolution with (dashed) and without (solid) solar evolution.

Open with DEXTER

Figure 13 shows the evolution of the HZ boundaries for the reference scenario with an initial mantle water concentration of 500 ppm and mantle temperature of 1700 K. Figure 13a shows the evolution of the inner (red) and outer (blue) boundary with time in units of the solar constant. This figure shows that the outer edge of the HZ strongly varies in time in response to the outgassing from the interior, mainly of CO2, whose abundance increases during the evolution. The amount of water vapor in the atmosphere is constant since it is defined via the saturation vapor pressure at 273 K. The inner edge of the HZ, expressed in total solar irradiation, shows a small increase with time during the early evolution, which corresponds to an increase in water outgassed from the interior (see Fig. 7) since a larger water reservoir needs to be evaporated. At later evolutionary stages the inner boundary is found at a nearly constant insolation with time since the amount of water outgassed from the interior is approximately constant.

Figure 13b shows the orbital distances at which a planet would receive the insolation needed to reach the boundaries of the HZ. The solid lines show the boundaries when neglecting the solar evolution with time, while the dashed lines show the orbital distance including the fainter Sun during the early evolution (following Gough 1981). When neglecting the solar evolution, the orbital distances of the HZ boundaries show a similar evolution as the HZ boundaries in terms of solar irradiation (Fig. 13). Including the solar luminosity evolution leads to an increase in orbital distance of the inner HZ with time as the Sun brightens. At the outer edge of the HZ, the increase in orbital distance with time due to the increase in CO2 in the atmosphere is enhanced by the brightening Sun. However, the distances at the beginning of the evolution are shifted towards the Sun. The width of the HZ increases with time owing to the increase in atmospheric CO2 and is larger when accounting for the solar evolution because of the r-2 dependence of the stellar flux at the position of the planet.

Figure 14 shows the HZ boundaries after 4.5 Gyr of evolution for different interior scenarios in comparison to the HZ boundaries as determined by Kasting et al. (1993b) and Kopparapu et al. (2013). For low oxygen fugacities, the inner edge of the HZ is located at similar insolations, with a slight shift towards the Sun until fO2 approaches a value of IW+0.5. At IW+1, and for low initial water concentrations, the inner boundary is located at smaller insolations. While the amount of CO2 outgassed from the interior increases with the oxygen fugacity, the amount of water decreases because of its increased solubility in surface lavas. This behaviour is therefore the result of the competing influences of H2O and CO2 as radiative gases, which are present in different abundances.

thumbnail Fig. 14

Habitable zone boundaries for different oxygen fugacities (relative to the IW buffer) and initial water concentrations in the mantle (indicated in different shades of blue) in comparison to other HZ boundaries (Kasting et al. 1993b, in red; and Kopparapu et al. 2013, in green). Solid and dashed lines refer to inner and outer HZ boundaries, respectively.

Open with DEXTER

The increase in insolation needed to reach the inner boundary of the HZ when moving from low to medium oxygen fugacities (from IW-1 to about IW or IW+0.5 depending on the initial water concentration in the mantle) is caused by an increase in planetary albedo as the amount of CO2 increases, while the greenhouse effect in the atmosphere owing to CO2 and H2O is nearly constant (not shown). At an oxygen fugacity of IW+1, the planetary albedo still increases, but also the greenhouse effect of CO2 becomes more efficient, as CO2 rises from about 4.5 to 14.2 bar for an initial mantle concentration of water of 250 ppm upon increasing the oxygen fugacity from IW+0.5 to IW+1. This leads to a noticeable shift of the inner HZ away from the Sun because, already at less stellar irradiation, higher surface temperatures are reached as a result of the high amount of CO2 in the atmosphere. In addition, the amount of water released from the interior decreases from about 2.1 to 1.6 bar for the same scenarios (for an initial water mantle concentration of 250 ppm when moving form IW+0.5 to IW+1), hence less water is available that can be evaporated.

At the outer boundary of the HZ, the amount of solar irradiation needed to keep the global mean surface temperature at the freezing point of water decreases with increasing oxygen fugacity, hence CO2 outgassed from the interior, owing to the increasing greenhouse effect. At oxygen fugacities between IW+0.5 and IW+1, however, the decrease is much smaller, showing that about the same insolation is needed to obtain 273 K for atmospheres containing about 10 bar more CO2. This is because not only the greenhouse effect is increased for increasing amounts of CO2 but also scattering becomes more efficient.

4. Discussion

4.1. Stagnant-lid versus plate-tectonic planets: interior and outgassing evolution

There are a few studies on the long-term evolution of the habitability of the Earth and other planets in which calculations of CO2 outgassing based on actual models of the thermal and melting history of the mantle have been used in combination with atmospheric models of varying complexity to predict the resulting climate (e.g. Kadoya & Tajika 2014, 2015; Foley & Driscoll 2016). From the point of view of modelling of the interior, the most important difference between our approach and those above is that we assume our planet to operate in a stagnant-lid mode of convection rather than in the mobile-lid (or active-lid) mode realized through plate tectonics. On the one hand, in stagnant-lid bodies the presence of a thick, immobile lithosphere limits the production of partial melt, making this possible only at large depths where the solidus is relatively high. As the mantle cools and the lithosphere thickens, melting tends first to wane (Fig. 5c) and then to vanish completely within a few billion years (Kite et al. 2009). On the other hand, in a planet with an active lid as the Earth where thin tectonic plates form the lithosphere, melting can take place over a wide pressure range up to shallow depths of few kilometers beneath mid-ocean ridges with the consequence that outgassing can typically occur over the entire lifetime of the host star despite the cooling of the interior (Kite et al. 2009).

In contrast to a plate-tectonic body, for a stagnant-lid body, volatiles cannot be cycled as efficiently as they are cycled, for example, in the water and carbon cycles of the Earth, even though our model allows for the possibility of some reintroduction of crustal – and hence volatile-rich – material into the mantle and therefore is not strictly a one-way street, at least not for water. Still, as Fig. 5d shows, there is a secular loss of water from the deep interior. The absence of subduction in our models has also an indirect influence on the chemistry of the interior that may feed back into the degassing behaviour and, in turn, the atmosphere. The recycling of water, which in the Earth is considered an oxidizing agent reintroduced into the mantle chiefly at subduction zones (e.g. McCammon 2005) is comparatively small in our models. Therefore, the evolution from a reducing environment towards a more oxidizing environment, which has been postulated at least for the upper mantle of the Earth (Frost & McCammon 2008), would be strongly inhibited in a stagnant-lid planet. This justifies a posteriori our assumption that the oxygen fugacity can be held fixed at a rather low value. In the context of plate tectonics and subduction, this assumption would probably have to be dropped, although how and to which extent these processes would alter the oxygen fugacity of the mantle is poorly known. At any rate, an evolution towards more oxidizing conditions would imply that carbon would change from its reduced, bound form into the volatile form of CO2 that can be degassed, as long as any carbon is available. In the absence of significant sinks of carbon at the surface, an enormous CO2 pressure would build up in the atmosphere. This is also suggested by Fig. 9b, although that diagram was not calculated for such a scenario.

Our model of redox melting considerably simplifies the treatment of CO2 extraction and outgassing, which could become extremely complex in a plate-tectonic setting. In thermal evolution models based on plate tectonics, the amount of outgassed CO2 is typically treated independently of the mantle composition and simply related to the seafloor spreading rate and depth of melt generation (e.g. Kadoya & Tajika 2015; Foley & Driscoll 2016). Yet carbon in the upper mantle of the Earth is stored in a number of accessory phases rather than being bound to silicates as in the case of water (Dasgupta & Hirschmann 2010). Together with the fact that under oxidizing conditions, the presence of carbon can dramatically reduce the solidus of domains of carbonated mantle (Dasgupta & Hirschmann 2006), it is clear that modelling the melting and extraction of carbon under these conditions would be a challenging task, particularly when using 1D models that cannot account for lateral variations in the volatile distribution.

Thermal models based on plate-tectonic convection generally assume this mode to operate throughout the planetary evolution. Although the surface rock record of the Earth and numerical models suggest that some form of surface mobilization and recycling may have already been present during the Archean, subduction in the early Earth, if existent at all, was probably characterized by an episodic, intermittent behaviour (see van Hunen & Moyen 2012, and references therein). It is only during the past two to three billion years that modern-style subduction, consisting of a continuous, regular creation of new seafloor at mid-ocean ridges and of recycling of cold plates at subduction zones, is likely to have been active (e.g. Condie 2016). The possibility that plate tectonics may have not operated uniformly over the evolution of the Earth, in particular that it could have followed an initial stagnant-lid state (e.g. Korenaga 2013; Moore & Webb 2013) and could even be a transient rather than an end-member state (O’Neill et al. 2016) makes the modelling of long-term mantle melting and outgassing particularly difficult and subject to major uncertainties. Furthermore, the present understanding of the basic physical mechanisms that are ultimately responsible for the generation of plate tectonics is still incomplete (Bercovici 2003). Since the solar system provides us with examples of bodies that do not show any evidence for present or past plate tectonics, our simplifying assumption of considering a stagnant-lid planet appears thus capable of delivering robust results.

Owing to the similarity of Earth and Venus in terms of size and composition and to the fact that the latter does not show evidence of plate tectonics today, it is also tempting to compare our outgassing results with the composition of the atmosphere of Venus; since we have limited our atmosphere calculations to temperatures below the critical point of water, we cannot reproduce Venus-like climates. At present, the atmosphere of Venus contains 92 bar of CO2 and, with only 20 ppm of water vapour, is essentially dry; this is a figure that is actually not very surprising, at least as far as the amount of CO2 is concerned. Earth and Venus are thought to have similar compositions (e.g. Fegley 2005). But while on Venus CO2 resides mostly in the atmosphere, on Earth it is largely contained in carbonate rocks. Indeed, as much as the equivalent of 60 bar of CO2 is thought to be stored in carbonates on Earth (Kasting 1988), whose formation is strongly facilitated by the presence of water, which Venus lost during its early evolution (see e.g. Kasting et al. 1984a; Taylor & Grinspoon 2009). As shown in Fig. 9b, by choosing relatively oxidizing conditions for the mantle (fO2 between IW+1.5 and IW+2), it is not difficult to obtain a 90-bar present-day atmosphere, nearly independently of the initial water concentration of the mantle. The build-up of such a thick atmosphere would also prevent water from being released from surface lavas and would thus lead to a strongly water-depleted atmosphere. Our models generally predict CO2 to be continuously outgassed during the evolution, albeit at a decreasing rate (Fig. 8b). The abundance of CO2 in the atmosphere of the Earth is buffered over geological timescales by surface weathering and the carbon-silicate cycle. Whether a comparable process able to buffer the CO2 concentration over long timescales is also at play on Venus is unclear (Taylor & Grinspoon 2009). Yet understanding whether such a process actually operates would help to clarify whether the surface pressure of Venus was acquired early and remained nearly stable over the evolution (Gillmann & Tackley 2014) or evolved substantially over time, which is the scenario predicted by our models. A direct comparison with Venus is also not straightforward because we assume that the stagnant-lid regime persists throughout the evolution. Venus, in contrast, with its young surface, is thought to have experienced a large-scale resurfacing event about 1 Gyr ago (Romeo & Turcotte 2010). This event, which replenished the ancient surface with new volcanic material, may also have been accompanied by a significant injection of volatiles into the atmosphere, whose extent, however, is difficult to constrain.

The question as to whether the atmospheric pressure of Venus has a primordial origin or evolved over time is also related to a fundamental assumption of our models, namely that the atmospheres that we considered are generated solely by secondary outgassing of H2O and CO2 from the interior and that any primordial atmosphere, whether accreted from the nebula or degassed by a magma ocean has been lost, except for the presence of 1 bar of N2. The assumption of 1 bar N2 can be motivated by a comparison of Venus and Earth, both of which hold about the same amount of N2. Furthermore, this assumption allows for a better comparison with other habitability studies (Sect. 4.2) such as those performed by Kasting et al. (1993b) and Kopparapu et al. (2013). At any rate, the impact of this choice is not very significant. We tested the effect of a lower amount of N2 on the surface temperatures for our reference scenario ( ppm and fO2 at the IW buffer) after 4.5 Gyr and found that with a N2 pressure of 0.1 bar we obtain a surface temperature of 357.45 K, while with only 0.01 bar we obtain 356.40 K, i.e. less than 10 K lower than in the case with 1 bar N2.

The neglect of a thick primordial atmosphere such as one in excess of tens to hundreds of bars of H2O and CO2 that could be generated by catastrophic outgassing of a magma ocean (e.g. Elkins-Tanton 2008; Lebrun et al. 2013) has a potentially much more significant influence on the outcomes of our models. In contrast to previous studies, no matter whether based on plate-tectonic or stagnant-lid convection, our models account for the effect of growing atmospheric pressure on the outgassing of H2O and CO2 (see Figs. 7 and 8). Although with the kind of atmospheres that we have considered the outgassing of CO2 is hardly influenced by this effect because of its low solubility (Fig. 4b), the outgassing of H2O, which on the contrary is much more soluble in basaltic lavas (Fig. 4a), can be strongly suppressed, particularly when the mantle oxygen fugacity is large enough for significant quantities of CO2 to be released into the atmosphere. Indeed, as shown in Fig. 9a, depending on the actual amount of outgassed CO2, this poses strict limits on the amount of water that can be brought to the surface and atmosphere via volcanism. Therefore, depending on the pressure of a primordial atmosphere and on its resilience to escape, volcanic outgassing of H2O could be easily suppressed from the very beginning of the evolution, and subsequent outgassing of CO2 could also be hindered.

4.2. Habitability

The term “habitability” denotes in a broad sense conditions that are hospitable to life. Life as we know it has three basic requirements, namely (i) an energy source; (ii) a solvent; and (iii) the potential for molecular complexity (see Horneck 2006). The classical HZ (Huang 1960; Kasting et al. 1993b) refers to the annulus region around a star where a planet could support liquid water as a solvent on its surface. The wider, non-classical HZ however (see e.g. Lammer et al. 2009) extends further out to include, for example icy moons beyond the snowline, which are widely accepted to have subsurface oceans.

In exoplanetary science, habitability is commonly defined using the classical HZ concept although clearly the existence of liquid water is not the only criterion. Numerous other potentially important factors have been discussed in the literature. These include (i) the availability of nutrients, for example in the form of the elements CHNOPS (see e.g. Horneck et al. 2016); (ii) the formation and maintenance of an atmosphere that enables the presence of liquid water (Grenfell et al. 2010); (iii) the existence of a magnetic field that could protect the surface from high-energy particles (Lundin et al. 2007); (iv) climate stabilization, possibly via the presence of a large moon (Laskar et al. 1993); (v) environmental diversity and climate-stabilizing feedbacks, for example due to the presence of plate tectonics which favours the carbonate–silicate cycle (e.g. Korenaga 2011; Höning & Spohn 2016); (vi) planetary protection from impacts due to the presence of gas giants (Horner & Jones 2008); and (vii) the properties of the central star (Beech 2011). It should be noted that complex molecules such as proteins have a limited range of temperature over which they are stable. Most proteins can exist up to 50 °C, while some can exist up to 120 °C (see e.g. Lineweaver & Chopra 2012). A detailed review of the concept of habitability can be found in Cockell et al. (2016).

In this study we focused on the availability of liquid water as a criterion for habitability as used in the concept of the classical HZ. We used a mean surface temperature of 273 K as the lowest habitable temperature. However, recent studies have shown that liquid water on a planetary surface is also possible for global mean temperatures below 273 K with temperatures as low as ~235 K since, locally, the surface temperatures are above 273 K and may still allow for liquid water (e.g. Charnay et al. 2013; Wolf & Toon 2013; Kunze et al. 2014; Shields et al. 2014; Godolt et al. 2016).

At high temperatures, we have not defined a limit of habitability in terms of surface temperature, but rather in terms of water reservoir available from outgassing. We defined the upper limit for habitability at the stellar irradiation at which the entire water reservoir outgassed from the interior would be in its vapor phase according to the assumption of phase equilibrium at the surface, as previously done also by Pollack (1971). Therefore this limit occurs at different surface temperatures. This approach differs from most of the definitions proposed in the literature. A hard limit used by Kasting et al. (1993b) and Kopparapu et al. (2013) is the runaway greenhouse limit. This limit defines the stellar irradiation at which a planet with a water reservoir of one Earth ocean and an Earth-like atmosphere does not allow for climate states with liquid water on the planetary surface as the water feedback tends to increase the surface temperatures and pressures above the critical point of water. Kasting et al. (1993b) also suggested another habitability limit, the water loss limit. It is defined at a stellar irradiation at which a water reservoir of one Earth ocean (270 bar) may be lost to space within 4.5 Gyr. Water loss is thought to be efficient if stratospheric water vapor volume mixing ratios exceed a critical value of about 3 × 10-3 (Kasting et al. 1993b). This limit is estimated to occur at temperatures of about 340 to 350 K (e.g. Selsis et al. 2007; Wolf & Toon 2015). Some of our stagnant-lid scenarios show higher global mean surface temperatures at later stages of the evolution and may therefore lose their water reservoir, which is even smaller than one Earth ocean, via photodissociation of H2O molecules and subsequent loss of hydrogen to space.

The comparison of our model results with the HZ boundaries in Fig. 14 shows that the inner boundary of the HZ determined for the stagnant-lid scenarios lies almost for all cases in between the boundaries determined by Kasting et al. (1993b) and Kopparapu et al. (2013). In Fig. 14 we show the inner boundary of the HZ that these two studies determined as the runaway greenhouse limit for an Earth-like planet with a water reservoir of one Earth ocean (270 bar). The main difference between these two studies lies in the spectral databases and continua used to derive the radiative fluxes in the planetary atmosphere. It has been shown that especially the treatment of the radiative properties of water vapor can lead to different results at the inner edge of the HZ (Yang et al. 2016). Here we used yet another spectral database (HITEMP1995) and continuum assumptions different from those of Kasting et al. (1993b) and Kopparapu et al. (2013). Furthermore, we used a forward climate modelling approach, where we specify the stellar irradiation, the atmospheric pressure and composition (except for water vapor, which is calculated) and calculate the surface temperatures self-consistently, while Kasting et al. (1993b) and Kopparapu et al. (2013) apply inverse climate modelling. In this approach, the temperature profile is fixed by specifying a surface temperature, following a dry or moist adiabat until a defined stratospheric temperature of, for example 200 K is met. For this temperature, the stellar irradiation needed to balance the outgoing infrared radiation is then derived. The 1D climate model used in this study would result in an insolation of about 1.32 SSun for the inner boundary of the HZ, when using the inverse climate modelling approach, specifying the surface temperature at about 647 K (the critical point of water), a stratospheric temperature of 200 K, a water reservoir of one Earth ocean, and deriving the insolation needed via the assumption of radiative energy balance at top of the atmosphere. Hence, the classical inner HZ boundary determined with our 1D model would lie in between the results of Kasting et al. (1993b; at 1.41 SSun) and Kopparapu et al. (2013; at 1.06 SSun). This result arises because we use a different database and continua. In comparison to this boundary considering an Earth ocean as water reservoir, the inner boundary of the HZ of the stagnant-lid planets generally lies further away from the star the smaller the water reservoirs obtained via outgassing from the interior are.

The assumption of a saturated atmosphere as made in our study and in Kasting et al. (1993b) and Kopparapu et al. (2013) has been questioned by recent 3D modelling studies (see e.g. Leconte et al. 2013a). The impact of different relative humidities upon 1D climate modelling results and their comparison to 3D model results has been investigated by Godolt et al. (2016). The assumption of a fully saturated atmosphere overestimates the surface temperatures at higher temperatures. Because of this the inner boundaries determined with this assumption are pessimistic. For lower relative humidities, the inner boundary of the HZ would move towards the star. For even smaller water reservoirs than those found here, the HZ could move inwards even further, as shown by, for example Abe et al. (2011), Zsom et al. (2013), or Leconte et al. (2013b), since the greenhouse effect of water would be very low at the resulting low atmospheric volume mixing ratios.

We found that the inner edge of the HZ may indeed depend on the amount of CO2 in the atmosphere. Kasting et al. (1993b) found that the runaway greenhouse limit for a planet with a water reservoir of one Earth ocean does not depend on the amount of CO2 in the atmosphere, while the water loss limit occurs at smaller orbital distances if the atmosphere contains more CO2. A study by Popp et al. (2016) also found that the inner boundary of the HZ as determined by water loss can be strongly influenced by the CO2 content of the atmosphere. Since we obtain a much smaller water reservoir than one Earth ocean by outgassing from the interior of a stagnant-lid planet, the atmospheres at the inner edge of the HZ are not necessarily dominated by water vapor as is the case in the runaway greenhouse scenarios as modelled by Kasting et al. (1993b). Therefore, increasing the amount of CO2 to a few tens of bars completely changes the mean composition of the atmospheres and can therefore also have an impact on the planetary climate.

The outer boundaries of the HZ determined for the stagnant-lid planets lie partly within and partly outside the HZ as determined by Kasting et al. (1993b) and Kopparapu et al. (2013) depending mainly on the oxygen fugacity. Our scenarios, which lie within the HZ, have lower partial pressures of CO2 than those assumed by the other two studies. For higher oxygen fugacities, the outer boundary of the HZ lies outside the outer boundary determined by Kasting et al. (1993b) and Kopparapu et al. (2013), which is surprising. We cannot infer a maximum greenhouse effect from our model calculations, which range up to CO2 partial pressures of 22 bar. At high partial pressures of CO2, we found the outer edge of the HZ to be fairly constant. In contrast, when applying the inverse climate modelling approach used by Kasting et al. (1993b) and Kopparapu et al. (2013) and assuming a stratospheric temperature of 150 K, we obtain a maximum greenhouse effect at CO2 partial pressures of about 4 bar. With this approach the outer boundary of the HZ would be located at 0.33 SSun, which is at a slightly lower insolation as for Kasting et al. (1993b; at 0.36 SSun) and Kopparapu et al. (2013; at 0.35 SSun) because of the effect of different databases and continua. Exploring the effect of inverse versus forward climate modelling upon the boundaries of the HZ in detail is beyond the scope of this paper. At the outer edge of the HZ, the assumption of a saturated atmosphere does not show a large impact as the water concentrations are still very low at 273 K. It has been shown in other studies (e.g. Pierrehumbert & Gaidos 2011) that other greenhouse gases, such as molecular hydrogen, which may still be present from a primordial atmosphere, may expand the outer boundary of the HZ far beyond the boundary defined here.

In the present study we only considered a stagnant-lid planet around the Sun. The surface temperature and habitable zone evolution of Earth-like stagnant-lid planets around other types of central stars can be expected to be different for two reasons. Firstly, the luminosity evolution is different for other stellar types, and secondly, the spectral stellar flux distribution is known to have a large impact on the planetary climate due to the wavelength-dependent absorption and scattering properties of the atmospheric compounds.

4.3. Influence of surface processes

One of the limitations of the models that we presented is the lack of consideration of surface processes that may alter the volatile budget predicted on the base of interior outgassing.

Liquid water on the surface, whether in the form of a global or regional ocean, could provide a sink for the atmospheric CO2. Such a sink is temperature dependent with a higher amount of CO2 trapped in the ocean for colder temperatures (see e.g. Pierrehumbert 2010; Kitzmann et al. 2015). Since we assume that all CO2 outgassed from the interior resides within the atmosphere, at the inner boundary of the HZ we obtain the maximum heating, hence the smallest insolation, because at these temperatures some CO2 could still be bound in the ocean. At the outer HZ boundary, part of the CO2 would be dissolved in the ocean (more than at the inner edge due to the temperature dependence of the solubility), which would lead to lower CO2 concentration in the atmosphere and hence probably a lower greenhouse effect. However, the ocean is smaller than the Earth ocean, hence less CO2 would be dissolved.

Apart from the ocean, there are other sinks for carbon and water that would have to be considered in a more detailed model. Among the most complex ones is probably the sequestration of carbon dioxides by carbonate minerals such as calcite (CaCO3) (e.g. Walker et al. 1981) and of water in various hydrated silicates. While in the modern Earth the biogenic production of CaCO3 is most familiar and has been producing carbonate sediments for hundreds of millions of years, there are also abiogenic processes that may also operate in a planet that does not, or does not yet, host life. One such process of this type is the low-temperature (<60 °C) alteration of the uppermost part of the ocean floor in the modern Earth, in which carbonates form from basaltic rock, but its quantitative significance is controversial (e.g. Alt & Teagle 1999). It is nonetheless worthwhile considering because it does not require plate tectonics; a carbon dioxide source and the existence of some sort of hydrothermal activity in basalt, which is a widespread rock type on other terrestrial planets without plate tectonics, should be sufficient. For instance, van Berk et al. (2012) modelled Mg–Fe-rich carbonates in the Comanche outcrop in Gusev crater on Mars, which was analysed by the Mars Exploration Rover Spirit, and found that one possible scenario for their formation involves a water-rich environment at temperatures around 280 K and 0.5–2 bar PCO2, as may have prevailed during the Noachian. Alt & Teagle (1999) analysed borehole cores from various ocean drilling sites that probed young oceanic crust and observed a decrease of bound carbon with depth into the extrusive rock; most of the CO2 is stored in the upper 600 m. They found an average content of 0.214 wt.% CO2 in the bulk rock, while Staudigel (2014) arrived at a total carbonate uptake of 0.355 wt.% and an uptake of 0.45 wt.% of crystal-bound water. A synopsis of data from the literature suggests that although young volcanic crustal rock may not take up a major part of this CO2 and water, the oceanic crust remains an active sink as it ages and continues to react for some tens of millions of years (e.g. Staudigel 2014).

An important difference between the modern Earth and the stagnant-lid Earth-like planet considered in this study is that in the latter, there is no mechanism that transports carbonated or hydrated rock back into the mantle at a steady rate, thus completing the carbon and water cycles, although periodic crustal delamination may still operate. Moreover, the mode of crust production is also different, as mid-ocean ridges would also be absent: the crust would be built by intrusions and, in a top-down manner rather than lateral addition, by extrusives. If the supply of fresh basalt ceases, for instance, because the crust becomes too thick to allow the ascent of lava to the surface or because the mantle becomes too cool, the crustal carbon and water sinks are not renewed and therefore are saturated after a certain period. In the atmospheric evolution, such a situation would result in a delay of the atmospheric accumulation of that volatile, as opposed to the suppressed or reduced accumulation expected in evolutions without a sink. In a planet with plate tectonics, the availability of carbon sinks would be even more important to prevent the build-up of a massive CO2 atmosphere because the increase in oxygen fugacity may reinforce CO2 degassing, as discussed above (see Sect. 4.1).

The weathering of impact ejecta is closely related to the sequestration of carbon and water in normal crust. Early in the history of the terrestrial planets, the much higher impact flux has produced abundant fresh rock surfaces of mostly mafic to ultramafic chemistry that were available for weathering (Koster Van Groos 1988; Zahnle & Sleep 2002). The latter authors have carried out model calculations for CO2 that suggest that ejecta weathering could bind enough of it to suppress the greenhouse effect entirely during the first few hundred millions of years. The possible effects of impacts go beyond the production of reactive surfaces, however. Depending on the sizes, velocities, and compositions of impactors, the impact processes themselves may result either in the erosion or the replenishment of the atmosphere of the target planet (Ahrens 1993; Pham et al. 2011; de Niem et al. 2012). Moreover, very large impacts also affect the deeper interior and may trigger volcanic activity, which in turn affects the atmosphere (Marchi et al. 2016). The effects especially of rare, large impacts depend strongly on the parameters of the specific event and cannot be predicted in a unique general manner.

5. Conclusions

The upcoming generation of exoplanetary missions holds promise for detecting small rocky planets orbiting Sun-like stars in the HZ (Rauer et al. 2014). The measurement of mass and radius alone does not well constrain whether these bodies have plate tectonics. We have thus studied the interior and climate evolution of an Earth-like planet operating in the stagnant-lid mode of convection, which presently characterizes all solid bodies of the solar system other than the Earth.

We analysed the evolution of the surface temperature and boundaries of the HZ based solely on volcanic degassing of H2O and CO2 from the interior. The partial pressures of H2O and CO2 that accumulate over time in the atmosphere are controlled by the initial water concentration and by the oxygen fugacity of the mantle, respectively. As the total atmospheric pressure grows, the release of H2O becomes more and more difficult because of its high solubility in surface lavas. After 4.5 Gyr, no more than a few tens of bars H2O can be outgassed, even when assuming initial mantle concentrations in excess of 1000 to 2000 ppm. As a consequence, an Earth-sized ocean cannot be built up by secondary outgassing from the interior. On the other hand, CO2, because it is much less soluble than water, can be outgassed throughout the evolution with partial pressures that depend on the assumed redox state of the mantle. While at reducing conditions (fO2 from IW-1 to IW), only up to 1 bar CO2 is outgassed, between 10 and 200 bar can be outgassed if the mantle is more oxidizing (fO2 from IW+1 to IW+2).

At 1 au, for all cases that we analysed (i.e. up to fO2 at IW+1), the amount of H2O and CO2 outgassed from the interior leads to surface temperatures that allow for liquid water on the surface over nearly the entire evolution. A stagnant-lid Earth could then be habitable over geological timescales even though the obtained surface temperatures (~350–420 K) may eventually lead to water loss to space and are hardly compatible with the limits for complex terrestrial life.

The width of the HZ after 4.5 Gyr is controlled by the amount of outgassed CO2, which in turn is determined by the mantle oxygen fugacity. For fO2 at IW + 1 and assuming a relatively high mantle water concentration (e.g. 1000 ppm), for example, the HZ requires an insolation between ~0.2 and 1.2 SSun. In contrast, for fO2 at IW-1, the HZ is considerably thinner and an insolation between 0.7 and 1.15 SSun is necessary if the same H2O concentration is assumed.

The outer edge of the HZ is mostly influenced by the outgassed CO2, and is the farther away from the Sun the higher fO2. The inner edge is characterized instead by a complex shape dependent on the amount of CO2 and on that of H2O. When the partial pressure of H2O becomes too low in response to a CO2 increase, in fact, the limited water reservoir can easily evaporate so that the inner edge of the HZ is significantly shifted away from the Sun.

In conclusion, the joint modelling of interior evolution, volcanic outgassing, and accompanying climate is crucial for a robust determination of the habitability conditions on rocky exoplanets.

Acknowledgments

We are grateful to Craig O’Neill for his constructive comments on an earlier version of this work. N. Tosi acknowledges support from the Helmholtz Association (project VH-NG-1017) and from the German Research Foundation (DFG) through the Priority Programme 1833 “Building a habitable Earth” (grant TO 704/2-1). A.-C. Plesa acknowledges support from the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office through the Planet TOPERS alliance and from the DFG (SFB-TRR 170). T. Ruedas was supported by the DFG (grant Ru 1839/1-1 and SFB-TRR 170). This is TRR 170 Publication No. 29.

References

All Tables

Table 1

Main parameters of the interior model.

Table 2

Parameters of the atmospheric model.

All Figures

thumbnail Fig. 1

Diagram of the interior structure considered for the thermal evolution model and a schematic of the corresponding temperature profile. See text for the description of the various symbols.

Open with DEXTER
In the text
thumbnail Fig. 2

Panel a: Solidus temperatures of peridotite for different water contents from dry (grey) to water saturated (dark blue), liquidus temperature (red), and a typical initial temperature profile with Tm = 1700 K (black). Panel b: Viscosity profiles for different water concentrations calculated with Eq. (7) along with the temperature profile shown in panel a.

Open with DEXTER
In the text
thumbnail Fig. 3

Ratio of the concentration of incompatible elements enriched in the liquid phase (Xliq) to the corresponding concentration in the solid mantle (Xm) as a function of melt fraction φ assuming fractional melting (Eq. (19)) and different partition coefficients δ.

Open with DEXTER
In the text
thumbnail Fig. 4

Saturation concentrations of H2O and CO2 in basaltic melts as a function of pressure. When surface melts are super-saturated, the excess pressure of the corresponding volatile is released into the atmosphere.

Open with DEXTER
In the text
thumbnail Fig. 5

Evolution of the interior for different initial mantle temperatures between 1600 and 1800 K, initial water concentrations () between 250 and 1000 ppm, and an oxygen fugacity corresponding to the IW buffer. Panel a: Mantle temperature is shown; b) corresponding mantle viscosity is shown; c) thickness of the crust (solid lines), the stagnant lid (dashed lines), and distribution of the melt zone and melt fraction (coloured area) are shown; and d) water concentration in the mantle is shown. The different colours, from grey to dark blue, indicate increasing initial water concentrations; the thick black lines in the four panels and the coloured melt zone in panel c indicate the reference model (ref.) with Tm,0 = 1700 K, ppm, and fO2 at the IW buffer. The different initial mantle temperatures, which can be identified in panel a, are not indicated with distinct colours.

Open with DEXTER
In the text
thumbnail Fig. 6

Extraction and outgassing evolution of H2O and CO2 for models with initial water concentrations () between 250 and 1000 ppm, an initial mantle temperature of 1700 K, and an oxygen fugacity at the IW buffer. Panel a: Concentration of H2O in the extracted melt (solid lines) and, for the reference model only, concentration of water at saturation (dashed line) are shown; b) average melt fraction is shown; c) partial pressure of outgassed H2O (also expressed as equivalent ocean thickness) is shown; d) partial pressure of outgassed CO2 is shown.

Open with DEXTER
In the text
thumbnail Fig. 7

Evolution of the partial pressure of outgassed H2O for a model with Tm,0 = 1700 K, fO2 at the IW buffer, and two initial water concentrations ( and 1000 ppm) calculated by neglecting the evolving saturation concentration (dashed lines) and taking such a concentration into account (solid lines).

Open with DEXTER
In the text
thumbnail Fig. 8

Evolution of the partial pressure of outgassed H2O a) and CO2b) for a model with Tm,0 = 1700 K, ppm and oxygen fugacities ranging from one log10-unit below the IW buffer (IW-1) to one log10-unit above it (IW+1). The black line indicates the reference model.

Open with DEXTER
In the text
thumbnail Fig. 9

Maximum partial pressures of outgassed H2O a) and CO2b) after 4.5 Gyr of evolution as a function of the initial water concentration in the mantle () and oxygen fugacity (fO2) for an initial mantle temperature of 1700 K.

Open with DEXTER
In the text
thumbnail Fig. 10

Evolution of the surface temperature for a: an oxygen fugacity at the IW buffer and different initial mantle concentrations of water (in different shades of blue), b: an initial mantle concentration of water of 500 ppm and different oxygen fugacities (in different shades of yellow and red) of an Earth-sized stagnant-lid planet located at 1 au around the evolving Sun. The reference scenario (Tm,0 = 1700 K, , and fO2 at the IW buffer) is indicated in black. The freezing point of water (273 K) is indicated by the dashed black line in both panels. The temperatures in panel a correspond to the outgassing evolutions of H2O and CO2 plotted in Figs. 6c and d, while the temperatures of panel b correspond to the outgassing evolutions of Figs. 8a and b.

Open with DEXTER
In the text
thumbnail Fig. 11

Evolution of the partition of the outgassed H2O (black) between atmosphere (red) and ocean (blue) for the reference scenario with an initial mantle concentration of water of 500 ppm, an initial mantle temperature of 1700 K, and an oxygen fugacity corresponding to the IW buffer.

Open with DEXTER
In the text
thumbnail Fig. 12

Determination of the inner a) and outer b) HZ boundaries in solar irradiation S via evaporation of the entire water reservoir (PH2O,out), and the criteria of a minimum temperature of 273.15 K (dashed black line in panel b for an oxygen fugacity at the IW buffer and initial mantle water concentrations between 250 and 1000 ppm (indicated in different shades of blue). Owing to similar behaviours of the atmospheres, the model results are overlaid.

Open with DEXTER
In the text
thumbnail Fig. 13

Evolution of the inner (in red) and outer (in blue) boundary of the HZ in a) unit of solar constants (SSun) and b) HZ boundary evolution in orbital distance (au) for the reference scenario. The impact of the increasing solar luminosity can be inferred from comparing the HZ distance evolution with (dashed) and without (solid) solar evolution.

Open with DEXTER
In the text
thumbnail Fig. 14

Habitable zone boundaries for different oxygen fugacities (relative to the IW buffer) and initial water concentrations in the mantle (indicated in different shades of blue) in comparison to other HZ boundaries (Kasting et al. 1993b, in red; and Kopparapu et al. 2013, in green). Solid and dashed lines refer to inner and outer HZ boundaries, respectively.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.