© The Authors 2023

1 Introduction

The water inventory of a rocky planet originates from the time of formation, with water-bearing materials being delivered during accretion onto the protoplanet (O'Brien et al. 2018; Walsh et al. 2011). Water is brought to the surface and enters the atmosphere during the early magma ocean phase (Elkins-Tanton 2008; Hamano et al. 2013; Lebrun et al. 2013; Nikolaou et al. 2019) and later via volcanism throughout the lifetime of the planet (Tosi et al. 2017; Godolt et al. 2019). Water plays an important role in both interior and atmospheric processes. Its presence lowers the melting temperature of rocks (Katz et al. 2003) and their viscosity (Hirth & Kohlstedt 1996, 2004), with major implications for global-scale mantle dynamics (Nakagawa et al. 2015), planetary evolution (Tosi et al. 2017; Morschhauser et al. 2011), as well as for the emergence of plate tectonics (Peslier et al. 2017). In the atmosphere, water is involved in a number of feedback processes controlling the climate, with water vapor strongly contributing to greenhouse heating (Kasting 1988; Catling & Kasting 2017). Liquid water is also an essential component in carbonate-silicate weathering processes, which help stabilize the climate over geological timescales (Walker et al. 1981). Furthermore, liquid water is a crucial prerequisite for life (Westall & Brack 2018; Cockell et al. 2016).

Volcanic outgassing links the mantle and atmosphere and establishes feedback loops as water and other volatiles are removed from the mantle and brought into the atmosphere. The composition of volcanic gases, in turn, is in large parts determined by the composition and pressure of the atmosphere (Gaillard & Scaillet 2014). The planetary mass and interior structure play an additional important role in shaping the rate of volcanism and interior dynamics (Noack et al. 2014, 2017; Stamenković et al. 2012). Many previous interior-atmosphere studies of the habitability of rocky exoplanets over geological timescales have either considered only selected feedback processes, or investigated only planets with an Earth-like mass and structure (Noack et al. 2014, 2017; Tosi et al. 2017; Foley & Smye 2018; Höning et al. 2019; Godolt et al. 2019; Bower et al. 2019; Kite & Barnett 2020; Spaargaren et al. 2020; Liggins et al. 2022).

Most parameters driving the interior evolution of rocky planets, especially mantle parameters such as the redox state and initial water content, are difficult to constrain and often inaccessible to observation. The large size of the parameter space and the numerous feedback processes between interior and atmosphere require a statistical approach in order to obtain a thorough overview over which planets are the most likely to exhibit oceans. In this work, we present the results of more than 280 000 coupled interior-atmosphere evolutions of rocky planets with different initial conditions and interior structures (Sect. 2.1), where we investigate the emergence of surface water based on planet mass, age, mantle volatile content and redox state, and orbital distance to the host star. Our model combines a 1D parameterized convection model (Sect. 2.2) to simulate the mantle and core evolution of rocky planets up to 3 M_⊕, as well melting and volcanic outgassing, with a gray atmosphere model tracking the evolution of atmospheric composition, pressure, and temperature. We focus on stagnant-lid planets in order to establish a baseline of worlds that could sustain liquid water on their surface. These are less prone to sustain habitable conditions than plate-tectonics planets due to a strongly reduced recycling of volatiles into the mantle, which would otherwise favor long-term, temperate climates (Walker et al. 1981; Kasting et al. 1993). We include a comprehensive array of feedback processes:

We employ a chemical speciation model to self-consistently treat the outgassing of volatile C-O-H species from surface melts (Ortenzi et al. 2020; Sect. 2.3). The final composition of outgassed volatiles depends primarily on the redox state of the melt and the current composition and pressure of the atmosphere (Sect. 2.4). Under oxidizing conditions, the predominant outgassed species are H₂O and CO₂, while reducing conditions favor the outgassing of oxygen-poorer species, mainly H₂ and CO.

We include a water cycle through a simple scheme of surface water accumulation and evaporation (Sect. 2.5). We assume the atmosphere to be fully saturated in water. Any excess water condenses to form a surface ocean or surface ice.

To account for the removal of CO₂ from the atmosphere, we incorporate a stagnant-lid CO₂ weathering cycle (Höning et al. 2019; Sect. 2.6). On Earth, the long-term carbonate-silicate cycle is important for stabilizing the climate over geological timescales. This cycle is primarily driven by plate tectonics, which continuously recycles CO₂ into the mantle via subduction. On stagnant-lid planets, the cycle relies on the production of fresh crust through volcanism. In the presence of liquid water, CO₂ can be weathered and buried under subsequent volcanic eruptions (Foley & Smye 2018; Höning et al. 2019).

Finally, we track the evolution of H₂ in the atmosphere and the various processes it contributes to. According to the evolution of the stellar XUV flux (Sect. 2.7), H₂ can be lost due to atmospheric escape (Sect. 2.8) and replenished by volcanic outgassing, particularly under reducing conditions. Since the surface pressure and atmospheric composition play an important role in the outgassing of H₂O, a thick (potentially primordial) atmosphere can limit the outgassing of water, especially in the early, active phase of a planet's evolution. In addition, H₂ can act as a greenhouse gas via collision-induced absorption (Pierrehumbert & Gaidos 2011), which we also take into account (Sect. 2.4).

2 Methods

2.1 Planet interior structures

The diversity in densities of low-mass exoplanets hints at a high degree of variation in interior structures (Jontof-Hutter 2019). We modeled rocky planets between 0.5 and 3 Earth masses. While many potentially rocky exoplanets with larger masses have been observed, the large pressure gradient in planets more massive than 5–6 M_⊕ can prevent melt from reaching the surface and thus impede the outgassing of volatiles (Noack et al. 2017). Additionally, the mantle rheology and interior dynamics of massive super-Earths are poorly understood (Stamenković et al. 2011; Karato 2011). We modeled the interior structure of each planet with our interior structure code TATOOINE (Baumeister et al. 2020). Each planet consists of an iron-rich core and a silicate mantle of Earth-like composition. We considered planets with core radius fractions ranging from 0.3 (a "Moon-like" planet) up to 0.7 (a "Mercury-like" planet). From these modeled interior structures, we calculated the average density of the core (ρ_c) and mantle (ρ_m) for use in the parameterized convection model: $$\begin{array}{cc}{\rho }_{\text{c}}=\frac{M_{\text{c}}}{4/3\pi R_{\text{c}}^3},& {\rho }_{\text{m}}=\frac{M_{\text{p}}-M_{\text{c}}}{4/3\pi \left(R_{\text{p}}^3-R_{\text{c}}^3\right)}\end{array},$$(1)

where M_p and R_p are the mass and radius of the planet, and M_c and R_c are the core mass and radius, respectively. We assumed a constant gravitational acceleration g throughout the planetary mantle, with $$g=\frac{G{M}_{\text{p}}}{R_{\text{p}}^2},$$(2)

where G is the gravitational constant. Figure G.1 shows the range of investigated planets in the form of a mass-radius diagram.

2.2 1D parameterized convection model

We employed a one-dimensional (1D) parameterized convection model to simulate the thermal evolution of the mantle and core of stagnant-lid rocky planets, as well as the melting of mantle rocks and volcanic outgassing of volatiles (Stamenković et al. 2012; Grott et al. 2011; Tosi et al. 2017).

We focused on stagnant-lid planets, since the emergence of and transition into plate tectonics is still poorly understood, and (especially for exoplanets) the question of which planets are the most likely to have plate tectonics has proven controversial, with many studies giving contradicting results (O'Neill et al. 2016; Noack & Breuer 2014; Stein et al. 2013; Van Heck & Tackley 2011; Korenaga 2010; Valencia et al. 2007; O'Neill & Lenardic 2007). Furthermore, plate tectonics favors establishing stable, temperate climates due to the efficient recycling of volatiles into the mantle (Walker et al. 1981). On stagnant-lid planets, this recycling is strongly reduced. In this sense, our study provides a conservative baseline to assess whether a planet can sustain liquid water on its surface. In addition, except for Earth, the terrestrial planets of the Solar System are in a stagnant-lid regime at present day, and likely have been for a majority of their evolution (O'Neill et al. 2007; Tosi & Padovan 2021) (a direct comparison to present-day Venus can be found in Appendix C).

Partial melting occurs everywhere where the mantle temperature profile exceeds the solidus. The model accounts for the partitioning of incompatible trace elements (water, CO₂, and radiogenic elements) between mantle, crust, and, in the case of volatiles, the atmosphere via volcanic outgassing. In addition, the presence of water depresses the solidus temperature (Katz et al. 2003). We focused here on fully extrusive volcanism, where all the melt produced reaches the surface and is subject to outgassing. In Sect. 3.5.2, we discuss the impact of intrusive volcanism on our results. Once melt reaches the surface, dissolved volatiles can be outgassed into the atmosphere. This process is subject to a number of limiting factors, such as the solubility of volatiles in the melt and the evolving atmosphere composition and pressure. We modeled the composition of outgassed species with a speciation model based on the chemical equilibrium of volatiles between melt and atmosphere (Sect. 2.3).

We neglected an early magma ocean phase and focused here on the outgassing of water only via volcanism. This provides a conservative estimate of the total amount of water that can be expected to be outgassed. While magma ocean solidification can result in the formation of a thick steam atmosphere, on an Earth-like planet this would rapidly collapse to form an early ocean (Elkins-Tanton 2011; Lebrun et al. 2013). Additionally, young planets likely experience significant water loss shortly after formation due to high stellar XUV activity (Tian et al. 2018). In Appendix E, we discuss the influence of early steam and CO₂ atmospheres that may follow magma ocean solidification (Lebrun et al. 2013; Nikolaou et al. 2019). A detailed description of the convection model as well as the melting and outgassing scheme can be found in Appendices A and B.

2.3 Outgassing speciation model

We adapted the model by Ortenzi et al. (2020) to calculate the chemical speciation of volatiles within the C-O-H system during outgassing from surface melts, based on the amount of dissolved H₂O and CO₂. We followed the approach by Holloway (1998) and Grott et al. (2011) to calculate the concentration of CO₂ in melts. We assumed a sufficiently reduced mantle to allow for carbon to occur as graphite, with an oxygen fugacity $f_{{\text{O}}_2}$ ranging from −3 to +3 in log₁₀ units above and below the IW buffer. The (depth-dependent) concentration of CO₂ in the melt ($X_{\text{liq}}^{{\text{CO}}_2}$) is then given by the concentration of carbonate ($X_{\text{liq}}^{{\text{CO}}_3{}^{2-}}$) $$X_{\text{liq}}^{{\text{CO}}_2}\left(r\right)=\frac{b{X}_{\text{liq}}^{{\text{CO}}_3}{}^{{}^{2-}}\left(r\right)}{1+\left(b-1\right)X_{\text{liq}}^{{\text{CO}}_3}{}^{{}^{2-}}\left(r\right)},$$(3a) $$X_{\text{liq}}^{{\text{CO}}_3{}^{2-}}\left(r\right)=\frac{K_{\text{II}}{K}_{\text{I}}{f}_{{\text{O}}_2}}{1+K_{\text{II}}{K}_{\text{I}}{f}_{{\text{O}}_2}},$$(3b)

where K_II and K_I are equilibrium constants governing the reaction of forming carbonate and graphite from CO₂, respectively, and b is a constant. K_II, K_I, and b were all determined appropriately for Hawaiian basalts (Holloway 1998). We calculated the melt concentration of H₂O ($X_{\text{liq}}^{{\text{H}}_2\text{O}}$) based on a model of fractional melting as described in Appendix B, assuming a partition coefficient ${\delta }_{{\text{H}}_2\text{O}}=0.01$. The solubility of H₂O and CO₂ is governed by melt-gas equilibrium reactions according to Iacono-Marziano et al. (2012): $${\text{H}}_2{\text{O}}^{\left[\text{fluid}\right]}+{\text{O}}^{2-\left[\text{melt}\right]}\rightleftharpoons 2{\text{OH}}^{-\left[\text{melt}\right]},$$(4a) $${\text{CO}}_2{}^{\left[\text{fluid}\right]}+{\text{O}}^{2-\left[\text{melt}\right]}\rightleftharpoons {\text{CO}}_3{}^{2-\left[\text{melt}\right]}.$$(4b)

We assumed that all of the generated CO and H₂ are outgassed due to their low solubility in silicate melts. The final molar composition of outgassed species is then governed by the following gas-gas equilibria: $${\text{H}}_2{}^{\left[\text{fluid}\right]}+\frac{1}{2}{\text{O}}_2\rightleftharpoons {\text{H}}_2{\text{O}}^{\left[\text{fluid}\right]},$$(5a) $${\text{CO}}^{\left[\text{fluid}\right]}+\frac{1}{2}{\text{O}}_2\rightleftharpoons {\text{CO}}_2{}^{\left[\text{fluid}\right]},$$(5b)

which can be converted to weight fractions $X_{\text{outg}}^i$ based on the molar mass of the magma (see Table F.1), and ultimately into a mass rate (see Eq. (B.12)). We did not include CH₄ in the speciation model, which starts to become relevant only at litho-spheric pressures and colder temperatures that are not reached here. In none of the models investigated did the outgassed CH₄ concentration exceed 10⁻¹¹ ppm.

The rates of the gas-gas equilibrium reactions depend mainly on the melt temperature and oxygen fugacity (see e.g., Ortenzi et al. 2020; Gaillard & Scaillet 2014), with the oxygen fugacity being dependent on the degassing pressure and temperature as well. We assumed that all melt reaches the surface, where it is subject to the current atmospheric surface pressure. To determine the melt temperature, we calculated the volume-averaged temperature and pressure of the melt region and obtained the surface melt temperature by moving the melt adiabatically to the surface.

2.4 Atmosphere model

H₂O and CO₂ are potent greenhouse gases and can strongly modify the surface temperature of a planet. H₂, while not being a strongly absorbing molecule on its own, can also act as a greenhouse gas through collision-induced absorption (Pierrehumbert & Gaidos 2011), which can be especially relevant for planets with hydrogen-dominated atmospheres. To estimate the surface temperature for a wide range of atmospheric compositions, we adopted a simple two-stream radiative gray atmosphere model to calculate greenhouse heating at the surface (Catling & Kasting 2017). The surface temperature can be expressed in terms of the optical depth τ of the atmosphere and the equilibrium temperature T_eq of the planet: $$\begin{array}{cc}{T}_{\text{s}}=T_{\text{eq}}{\left(1+\frac{\tau }{2}\right)}^{1/4}& \text{with\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{T}_{\text{eq}}={\left(\frac{\left(1-A\right)S_{\odot }}{4\sigma }\right)}^{1/4}\end{array},$$(6)

where S_⊙ is the solar insolation at the top of the atmosphere, A is the bond albedo, and σ is the Stefan-Boltzmann constant. We assumed an Earth-like albedo of 0.3 for all planets considered in this study. Following Abe & Matsui (1985) and Pujol & North (2003), the optical depth of the atmosphere is given by $$\tau ={\displaystyle \sum _i{\tau }_i={\displaystyle \sum _i\frac{3k_i^{\prime }{P}_i}{2g}},}$$(7)

where g is the planet gravity, P_i is the partial pressure of a given atmosphere species i, and $k_i^{\prime }$ is the extinction coefficient relative to this pressure. $k_i^{\prime }$ can be expressed using the extinction coefficient k_0,i at standard atmospheric pressure P₀ $$k_i^{\prime }={\left(\frac{k_0{}_{,i}g}{3P_0}\right)}^{1/2}.$$(8)

In order to approximate the collision-induced absorption of H₂, we chose a value of $k_{0,{\text{H}}_2}=2\times {10}^{-2}{\text{m}}^2{\text{kg}}^{-1}$, which fits well to the results of Pierrehumbert & Gaidos (2011; see also Fig. G.2). A validation of our atmosphere model against a 3D climate model can be found in Höning et al. (2021).

2.5 Water condensation and evaporation

We assumed the atmosphere to be fully saturated in H₂O, with any excess outgassed water condensing into a surface ocean or forming ice. This provides an upper limit for the mass of water in the atmosphere, and consequently also for the contribution of water to greenhouse heating. We calculated the saturated partial pressure (in Pa) of H₂O from the saturation vapor pressure curve by Alduchov & Eskridge (1996): $$P_{\text{vapor}}=610.94\exp \left(\frac{17.625\left(T-273.15\right)}{T-30.1}\right)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}}T\ge 273.15K$$(9)

where T is the temperature in K. If the surface temperature dropped below the freezing point of water, we assumed that the surface of an existing ocean would freeze. In this case, we used a vapor pressure curve from Alduchov & Eskridge (1996) defined over a plane of ice: $$P_{\text{vapor}}^{\text{ice}}=611.21\exp \left(\frac{22.587\left(T-273.15\right)}{T+0.71}\right)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}}T<273.15\text{\hspace{0.17em}K.}$$(10)

The combination of water evaporation and greenhouse heating from water acts as a positive feedback process which can give rise to a runaway process commonly described by the runaway greenhouse scenario (Komabayasi 1967; Ingersoll 1969; Kasting 1988; Nakajima et al. 1992). While with this model we did observe these rapid transitions from temperate conditions to hot, dense atmospheres, caused by the runaway evaporation of water, we did not perform a full, rigorous description of the runaway greenhouse. To avoid ambiguities and to make clear that our process is not a complete representation of the runaway greenhouse scenario, we are calling these “hothouse” states here.

2.6 Carbon weathering cycle

On Earth, the long-term carbon-silicate cycle is an important process to stabilize the climate over geological time-scales (Walker et al. 1981). This cycle is primarily driven by plate tectonics, where CO₂ can be continuously recycled into the mantle with subducting plates, and fresh crust is constantly produced at mid-ocean ridges. Stagnant-lid planets, on the other hand, do not have subducting plates. A carbon cycle on these planets relies on the continuous production of new crust through hot-spot volcanism, which can be weathered in the presence of liquid water. The carbonated crust may then be buried by subsequent volcanic eruptions, sinking downward in the mantle. The rate of weathering is therefore closely coupled to the crust production rate. We followed the model by Höning et al. (2019), assuming that the rate of CO₂ weathering depends on the partial pressure of CO₂ in the atmosphere. We additionally assumed that all newly formed crust is subject to weathering. The weathering rate Φ_w can then be expressed as a function of the crustal growth rate $\frac{\text{d}{M}_{\text{cr}}}{\text{d}t}$ and the partial CO₂ pressure in the atmosphere $P_{{\text{CO}}_2}$, and scaled to the seafloor weathering rate on Earth (for more details, see Höning et al. 2019): $${\text{Φ}}_{\text{w}}=\frac{X_{\text{E}}{\xi }_{\text{E}}}{f_{\text{E}}{\varphi }_E}\left(\frac{\text{d}{M}_{\text{cr}}}{\text{d}t}\right){\left(\frac{P_{{\text{CO}}_2}}{P_{{\text{CO}}_2,\text{E}}}\right)}^{\alpha },$$(11)

where α = 0.23 is a scaling exponent, and $P_{{\text{CO}}_2,\text{E}}=4\times {10}^{-4}$ bar is the present-day partial pressure of CO₂ in Earth's atmosphere. The other parameters are factors scaling the weathering rate to the observed present-day seafloor weathering rate on Earth: X_e is the present-day Earth mid-ocean ridge CO₂ concentration in the melt, ξ_E is the proportion of seafloor weathering to the total weathering rate on Earth, f_E is the present-day fraction of carbonates that are recycled back into Earth's mantle, and ϕ_E is the fraction of carbonates that remain stable during subduction.

Carbonates are stable only up to a certain pressure-dependent temperature. Once the sinking carbonated crust reaches this temperature, it undergoes decarbonation and releases its CO₂, which will rise through cracks in the crust and eventually return to the atmosphere (Foley & Smye 2018; Höning et al. 2019). This means that, in contrast to plate tectonics, CO₂ is generally not recycled back into the mantle in the stagnant lid regime. We calculated the depth z_decarb at which decarbonation occurs following Höning et al. (2019): $$z_{\text{decarb}}=\frac{T_{\text{s}}-B_{\text{decarb}}}{A_{\text{decarb}}-\frac{T_{\text{m}}-T_{\text{s}}}{D_1+d_{\text{u}}}},$$(12)

where A_decarb = 3.125 × 10⁻³ K m⁻¹ and B_decarb = 835.5 K are constants related to the decarbonation temperature (Foley & Smye 2018), T_m is the mantle temperature, D₁ is the thickness of the stagnant lid, and d_u the thickness of the upper thermal boundary layer of the convecting mantle (see Sect. 2.2 for a detailed description of the convection model).

During the evolution of the planet, our model continuously tracks the depth of the previously weathered crustal layers. Once the decarbonation depth is reached, the carbon content of the layer is released as CO₂. To avoid numerical instabilities, the released CO₂ is first stored in a temporary volatile buffer, which releases 10% of its CO₂ content into the atmosphere at every time step (see also Höning et al. 2019).

Carbon weathering requires the presence of liquid water. Therefore, we set the weathering rate to zero if either no surface water is present, or if the surface temperature lies below the freezing point of water. For the latter, this means that any existing ocean in our model freezes over so that little exchange of CO₂ with the ocean is possible. Nevertheless, in both cases, the burying of previously formed carbonates and decarbonation continues as long as there is active volcanism.

2.7 Stellar evolution

We focused on planets around G-type stars with one solar mass. We accounted for an increasing stellar insolation S_⊙ over the lifetime of the host star by using the parameterization by Gough (1981), $$S_{\odot }\left(t\right)=S_{\odot ,0}{\left(1+\frac{2}{5}\left(1-\frac{t}{t_0}\right)\right)}^{-1},$$(13)

where S_⊙,0 is the insolation at the planet at present day t₀ = 4.5 Gyr.

In order to model atmospheric escape processes, we followed Owen & Wu (2017) for a parametrization of the stellar XUV flux evolution: $$F_{\text{XUV}}\left(t\right)=\{\begin{array}{ll}\begin{array}{ll}{F}_{\text{sat}}\hfill & \text{for\hspace{0.17em}}t<t_{\text{sat}}\hfill \\ F_{\text{sat}}{\left(\frac{t}{t_{\text{sat}}}\right)}^{-1.5}\hfill & \text{for\hspace{0.17em}}t\ge t_{\text{sat}}\hfill \end{array}\hfill & \text{with\hspace{0.17em}}{F}_{\text{sat}}={10}^{-3.5}{S}_{\odot ,0},\hfill \end{array}$$(14)

with a saturation timescale of t_sat = 100 Myr.

2.8 Atmospheric escape

To model the transition from a primary H₂ to a secondary outgassed atmosphere and to treat the loss of later outgassed H₂, we included hydrodynamic escape of H₂. For hydrogen-dominated atmospheres, the maximum rate at which hydrogen can escape is limited by the amount of energy from XUV radiation that the atmosphere can absorb. The energy-limited mass-loss rate is given by $$\frac{\text{d}{M}_{\text{el}}}{\text{d}t}=\frac{\epsilon \pi R_{\text{p}}{R}_{\text{atm}}^2{F}_{\text{XUV}}}{G{M}_{\text{p}}},$$(15)

where ε is an efficiency factor we here take to be 0.15 following Kite & Barnett (2020), and R_atm is the planet radius at the top of the atmosphere, which we defined at 20 mbar.

In the case of hydrogen existing as a minor atmospheric component within a background of heavier species, the loss of hydrogen is limited by the rate at which it can be supplied from the lower parts of the atmosphere. This diffusion-limited escape provides an upper limit to hydrodynamic escape. The diffusion-limited mass loss rate can be expressed as $$\frac{\text{d}{M}_{\text{dl}}}{\text{d}t}=4\pi R_{\text{atm}}^2\frac{m_{{\text{H}}_2}}{N_A}{b}_{aj}{\chi }_{{\text{H}}_2}\left(\frac{1}{H_a}-\frac{1}{H_{{\text{H}}_2}}\right),$$(16)

where $m_{{\text{H}}_2}$ is the molar mass of molecular hydrogen, N_A is Avogadro's number, and ${{\rm X}}_{{\text{H}}_2}$ is the volume mixing ratio (VMR) of hydrogen in the atmosphere. $H_{{\text{H}}_2}$ and H_a are the unperturbed scale heights of H₂ and the background gas, respectively. b_aj is the binary diffusion coefficient between the escaping H₂ and the heavier background gas. In our case, the background gas consists of varying amounts of CO₂, CO, and H₂O. We calculated b_aj as the sum of the respective binary diffusion coefficients $b_{{\text{CO}}_2}=3\times {10}^{21}{\text{m}}^{-1}{\text{s}}^{-1}$, b_CO = 3 × 10²¹ m⁻¹ s⁻¹, and $b_{{\text{H}}_2\text{O}}=4.3\times {10}^{21}{\text{m}}^{-1}{\text{s}}^{-1}$ between H₂ and the background gases CO₂, CO, and H₂O, weighted by the respective relative volume mixing ratios ${\chi }_{{\text{CO}}_2}$, χ_CO, and ${\chi }_{{\text{H}}_2\text{O}}$ (see also Table F.2 for references): $$b_{aj}=b_{{\text{CO}}_2}{\chi }_{{\text{CO}}_2}+b_{\text{CO}}{\chi }_{\text{CO}}+b_{{\text{H}}_2\text{O}}{\chi }_{{\text{H}}_2\text{O}}.$$(17)

The transition from energy-limited to diffusion-limited escape, and thus from a hydrogen-dominated to a secondary atmosphere, is currently not well understood and requires the use of detailed hydrodynamical models (Owen 2019; Zahnle et al. 2019) which are out of the scope of this work, especially since volcanic outgassing continuously changes the atmospheric composition, which can make the hydrodynamical treatment challenging. Here, we opted for smoothly interpolating between energy-limited and diffusion-limited mass loss rates for intermediate H₂ fractions, with the interpolated mass loss rate given by $$\frac{\text{d}{M}_{\text{loss}}}{\text{d}t}=f_{\text{el}}\frac{\text{d}{M}_{\text{el}}}{\text{d}t}+\left(1-f_{\text{el}}\right)\frac{\text{d}{M}_{\text{dl}}}{\text{d}t},$$(18)

where the contribution of energy-limited escape f_el is given by a logistic function (see also e.g., Kite & Barnett 2020): $$f_{\text{el}}\left({\chi }_{{\text{H}}_2}\right)={\left(1+\exp \left(-\frac{{\chi }_{{\text{H}}_2}-{\chi }_0}{w}\right)\right)}^{-1},$$(19)

centered at a hydrogen VMR of χ₀ = 0.15 and a horizontal scaling w = 0.01. These parameters were chosen to account for a transition from purely energy-limited escape starting at a hydrogen VMR of 20% to a purely diffusion-limited escape at a mixing ratio of 10%.

We did not consider the photodissociation of water in the upper atmosphere and the subsequent loss of hydrogen. Significant water loss will occur once large amounts of water reach the stratosphere, which mainly occurs in planets undergoing a runaway greenhouse regime. On habitable planets, which are the main focus of this study, the tropopause "cold trap" prevents significant amounts of water vapor from reaching the stratosphere (Catling & Kasting 2017).

2.9 Investigated parameters and initial conditions

We adopted a Monte-Carlo sampling approach to model an entire population of planets, where the initial conditions for each planet are set to random values within given ranges (i.e., with a uniform distribution).

We computed the thermal evolution for a set of ≈280 000 initial conditions. Each evolution was run up to 8 Gyr to cover a wide range of potentially observable planets. We stopped the evolution earlier if the surface temperature exceeded 1500 K, at which point the surface rocks would be close to melting. In order to simulate the observation of planets with different ages, we selected snapshots of the evolution at up to five randomly chosen times after 100 Myr. For models which finished earlier (because the surface temperature had risen too high), we selected fewer snapshots accordingly to ensure a balanced sample of planet ages. This resulted in a final data set size of ≈ 1 000 000 planets.

We varied the initial water concentration in the mantle $X_{\text{m,0}}^{{\text{H}}_2\text{O}}$ between 100 and 1000 ppm, corresponding to relatively dry and wet conditions, respectively. We allowed the mantle oxygen fugacity to vary between three log₁₀ units below and above the iron-wüstite buffer (IW). Each planet in the parameter study was placed at a fixed distance to its (Sun-like) host star, ranging from a Venus-like orbit (0.723 au) to a Mars-like orbit (1.524 au). In addition, we set the mass of each planet between 0.5 and 3 M_⊕ and allowed for varied interior structures, where the radius of the core varied between 30 and 70% of the planet radius (Sect. 2.1). We fixed the initial mantle temperature T_m,0 at 1700 K and prescribed an initial temperature jump of 200 K at the CMB. The main model parameters used in our study are given in Table F.2.

We add a note on data presented in those figures below, where we show the surface state of water as a function of initial conditions (oxygen fugacity and initial mantle water content) for points from the data set described above. However, as multiple data points with different ages plot onto the same set of initial conditions, there is a risk to be biased toward later ages. To avoid this, for each figure, we shuffled the data set before plotting the points. We took advantage of the multiple data points per initial condition to map out regions of long-term behavior: We used a k-nearest neighbor algorithm on every point in the parameter space and colored the background according to the prevailing water phase at the surface that is most common among neighboring points.

3 Results

3.1 Characteristic planet evolutions

Figure 1 demonstrates the characteristic atmospheric evolution for two planets with one Earth mass and either an Earth-like interior structure (Figs. 1a, b, and c) or a large, Mercury-like core which makes up 70% of the interior (Figs. 1d, e, and f). Both planets are located at 1 au and the initial parameters are the same for both ($\log f_{{\text{O}}_2}=-0.05\text{\hspace{0.17em}IW}$, $X_{\text{m,0}}^{{\text{H}}_2\text{O}}=250\text{\hspace{0.17em}ppm}$). These parameters are in the range which enables prolonged habitable conditions for the Mercury-like planet, but puts the atmosphere of the Earth-like planet into a hothouse regime.

In both cases, an outgassed atmosphere of 0.1–1 bar is quickly built up within the first hundred million years. Due to the relatively reducing conditions of the mantle, this initial atmosphere consists mainly of CO, H₂, and CO₂ (Figs. 1b and 1e). The surface temperature is initially too low to allow for liquid water (Figs. 1a and d), but rises quickly due to the accumulation of greenhouse gases. Once the surface temperature exceeds the freezing point of water, the carbonate-silicate weathering cycle becomes active and is sufficiently strong to counteract the rate of CO₂ outgassing. This keeps the surface temperature just above the freezing point of water. The level of atmospheric H₂ remains relatively stable due to continuous outgassing supply and simultaneous loss from the top of the atmosphere.

Stagnant-lid planets lack an efficient long-term CO₂ sink. While the weathering cycle can temporarily remove CO₂ from the atmosphere, the slow sinking of carbonated crust causes decarbonation within the lid (Höning et al. 2019) and thus prevents recycling of carbonates into the mantle. As such, the combined mass of CO₂ in the crust and atmosphere reservoirs continuously grows through the influx of CO₂ from the mantle via volcanic outgassing. Besides the surface temperature, the mantle temperature and lid thickness determine the decarbonation depth (Eq. (12)). This gives rise to two scenarios for the evolution of the atmospheric CO₂ reservoir:

If the mantle is sufficiently hot, the crust continuously grows, burying carbonated layers. At some point, the bottom of the carbonated crust reaches the decarbonation depth and previously buried CO₂ degasses back into the atmosphere (see Sect. 2.6 and Höning et al. 2019). For the planet with Earth-like interior, this happens at around 0.8 Gyr (Fig. 1c). Decarbonation degassing quickly becomes the dominant source of atmospheric CO₂. A steady-state is reached between CO₂ ingassing caused by weathering of fresh rocks and CO₂ degassing from decar-bonation of old carbonated layers at the bottom of the crust, as both are tightly coupled to the crustal growth rate and thus to the rate of melt production. As the planet cools, volcanic activity slowly subsides. This reduces the weathering rate, which is tightly coupled to the crust production rate. At around 1.5 Gyr, crust production alone cannot sustain the weathering rates necessary to remove all outgassed CO₂. As the weathering rate is also proportional to the partial pressure of CO₂ (Eq. (11)), the amount of CO₂ in the atmosphere starts to grow until the weathering rate can counteract the increased CO₂ flux into the atmosphere (see also Fig. G.3). As long as decarbonation takes place, this renders CO₂ weathering largely ineffective as a climate control mechanism. Although decarbonation is closely linked to the crust production rate, the decarbonation flux remains relatively constant, even though the rates of volcanism and volcanic out-gassing decrease. This is the result of a hysteresis mechanism in the carbon cycle: Carbonated crustal layers reach the decarbonation depth around 200–300 Myr after being weathered at the surface (see Fig. G.4). The carbonated crust therefore retains a "memory" of past CO₂ weathering rates, when the rate of crust production was stronger. Upon decarbonation, relatively high amounts of CO₂ are returned to the atmosphere. The decarbonation flux is thus buffered against a decrease in crust production. In addition, the degassed CO₂ can be weathered again. This establishes a feedback loop, which stabilizes the decarbonation flux over geological timescales.

CO₂ starts accumulating in the atmosphere, driving the surface temperature upward over geological timescales (in addition to the increasing luminosity of the star), which in turn allows more water to enter the atmosphere. After a 2 Gyr-long habitable phase, the atmosphere transitions into a hothouse state, and the entire water reservoir evaporates to form a thick steam atmosphere with a surface pressure of around 45 bar. Once the surface heats up, decarbonation occurs in shallower and shallower layers. This positive feedback leads to a runaway decarbonation, where the entire crustal CO₂ reservoir is rapidly degassed back into the atmosphere (Fig. 1c). At this point, the evolution of this planet does not exhibit any qualitative change, since we do not consider water loss from steam atmospheres. The pressure is too high for much additional CO₂ to be outgassed, so little atmospheric evolution is possible at this point. Some part of the steam atmosphere would be lost via photodissociation, which would allow for more CO₂ outgassing. Such a planet would likely end up with a thick, Venus-like CO₂ atmosphere.

The evolution of the atmosphere is different in the case of the Mercury-like planet. The steep pressure gradient due to the higher gravity (compared to the Earth-like case) allows less melt to rise to the surface, resulting in less overall volcanic activity (Noack et al. 2017). Therefore, the crust grows more slowly. In addition, the thinner mantle cools faster. As a result, the growth of the decarbonation depth outpaces the bottom of the carbonated crust for most of the evolution of the planet, with very little decarbonation taking place throughout (Fig. 1f). Here, volcanic outgassing is the dominant source of CO₂. In contrast to the Earth-like planet, the crust acts as a long-term carbon reservoir, able to store any excess outgassed CO₂. Volcanism stops completely at around 4.5 Gyr as the core and mantle have cooled to temperatures that no longer make melting possible. Without volcanism, no volatiles are outgassed into the atmosphere, and the remaining H₂ is quickly lost due to atmospheric escape. Likewise, no CO₂ is removed by weathering since this depends on fresh basaltic rock delivered by volcanism, and no decarbonation occurs as the mantle is too cool and, due to the lack of crustal growth, carbonated layers are not transported toward the mantle. Over time, more water vapor enters into the atmosphere as a result of the rising surface temperature due to the increasing luminosity of the star. However, even at 8 Gyr, the planet remains habitable.

Fig. 1 Characteristic evolutions of the atmospheric temperature and pressure and CO2 flux in and out of the atmosphere for an Earth-mass planet at 1 au with an Earth-like (a–c) or a Mercury-like (d–f) interior structure. The colored background marks the state of water at the planet's surface. Both planets start with the same mantle water content of 250 ppm and an oxygen fugacity of 0.05 log10 units above the IW buffer. Panels (a–c) and (d–f) correspond respectively to the points 1 and 2 marked in Figs. 3 and 4.

3.2 Role of redox state on habitability

In this section, we focus on Earth-mass planets with Earthlike interior structures, with the effects of interior structure and planet mass being detailed in Sect. 3.3. We find that the mantle redox state and initial water content are the two main factors limiting the emergence of habitable conditions. Only a narrow range of these two parameters yields long-term stable habitable conditions (Fig. 2a). We can identify well-defined populations of planets with characteristic evolutions of surface habitability. At oxidizing conditions above the iron-wüstite (IW) buffer, the majority of planets are in a Venus-like hothouse regime with surface temperatures exceeding 400 K (Fig. 2b). At reducing conditions, one to two log10 units below the IW buffer, most planets with dry mantles have surface temperatures below the freezing point of water, whereas planets with wet mantles (initial mantle water content exceeding ~500 ppm) are in a hothouse state related to strong H₂ outgassing. Only over a narrow range of mildly reduced mantles are conditions suitable for liquid water to be stable on the surface over long time spans.

The evolution of the surface temperature and the partial pressure of water in the atmosphere are the main factors determining whether water can be outgassed and remain liquid on a planet's surface. The partial pressure of water in turn strongly depends on temperature. At cooler temperatures, most of the outgassed water is in the form of oceans or ice, with relatively small amounts in the atmosphere. At higher temperatures, more and more water vapor can be maintained in the atmosphere. Since water is a strong greenhouse gas, this acts as a positive feedback, where rising surface temperatures cause more water to evaporate (Catling & Kasting 2017). A rise in surface temperature can be caused by an increase in CO₂ or H₂ and, on longer timescales, an increase in solar luminosity. The outgassing of CO₂ and H₂ therefore predominantly shape the evolution of the surface temperature. Outgassing rates depend on the amount of melting in the interior. Water in the mantle decreases melting temperatures and viscosity, thus wetter mantles experience more melting and cause more volcanic activity.

The oxygen fugacity shapes the composition of outgassed species. At oxygen fugacities above the IW buffer, the main outgassed species are CO₂ and H₂O. At reducing mantle conditions, on the other hand, outgassing is dominated by the oxygen-poor species H₂ and CO. Increasing amounts of either CO₂ or H₂ in the atmosphere eventually push the planet into a hothouse regime, where any existing surface water quickly evaporates to form a H₂O-rich atmosphere, preventing further water outgassing. Strong CO₂ outgassing dominates on planets with an oxygen fugacity more than one log₁₀ unit above the IW buffer, which enter a hothouse state within the first billion year of their evolution (Fig. 3). Planets closer to the IW buffer may exhibit a short habitable phase. H₂ becomes the dominant outgassed greenhouse gas at an oxygen fugacity around one log10 unit below the IW buffer (see e.g., Fig. G.5 and Ortenzi et al. 2020). H₂ in the atmosphere is steadily lost due to atmospheric escape. Its presence in the atmosphere is therefore maintained only by a continuous replenishment from volcanic outgassing (Fig. 1b), and H₂ only builds up in the atmosphere if the rate of outgassing outweighs the rate of escape. We find that planets with significant H₂ outgassing may enter a hothouse state as well, specifically those with water-rich mantles for which outgassing rates of H₂ are high (Figs. 2a, b, Fig. 3). H₂ can build up faster on planets at high orbital distances, where atmospheric loss is less severe. These planets are more likely to end up in a H₂-induced hothouse as a result (Fig. 3).

Habitable conditions occur in the transition regime between CO₂- and H₂-dominated outgassing. Here, the combination of CO₂ and H₂ keeps the surface temperature above the freezing point of water, but a hothouse is prevented by weathering of excess CO₂ (Fig. 2d) and continuous loss of H₂ from the top of the atmosphere. These mechanisms favor a stabilization of the climate. Significant amounts of water can be outgassed under these conditions (Fig. 2c). Planets with very reduced, dry mantles can sustain habitable surface conditions for several billion years even at the orbit of Venus, albeit with very thin atmospheres.

3.3 Effect of interior structure and planetary mass

Planetary mass and the size of the iron core are additional important factors influencing the amount and lifetime of volcanism. Planets with large cores tend to undergo less partial melting due to larger hydrostatic pressure gradients in the mantle, which prevent melt from reaching the surface (Noack et al. 2017). Furthermore, the mantles of planets with large cores cool efficiently (Noack et al. 2017), which reduces the time a planet is volcanically active (as already shown in Figs. 1c and d). This prevents a hothouse state in many cases and enables the existence of habitable planets at a wider range of orbital distances, mantle water contents and oxygen fugacities (Fig. 4). Due to the lower outgassing rates, however, in many cases the surface remains frozen. By contrast, planets with small cores show strong, long-lasting volcanic activity, which limits the potential to develop habitable conditions (Fig. 4). The same applies to planets with higher mass, where the higher mantle volume supports long-lived volcanism. For very thick mantles however, the viscosity of the deepest mantle can become so large that a nonconvective, stagnant region is formed (Stamenković et al. 2012), shrinking the active, convective part of the mantle (see e.g., the case of 3 M_⊕, 30% core shown in Fig. G.6). With respect to outgassing, this resembles the behavior of a smaller planet with a larger core.

Fig. 2 States of a stagnant-lid Earth as a function of the mantle oxygen fugacity. Each point represents a snapshot of the planetary evolution at a randomly selected planet age between 100 Myr and 8 Gyr. The colored points in panel a show the prevailing water phase at the surface as a function of the oxygen fugacity ($f_{{\text{O}}_2}$) and initial water content of the mantle ($X_{{\text{H}}_2\text{O}}^0$). The color background shows the area where the majority of neighboring data points share the same surface conditions. In panels b-d, points are colored according to the greenhouse gas that contributes most to surface heating. Panel b shows the surface temperature, panel c shows the total mass of outgassed water, and panel d shows the total mass of CO2 removed from the atmosphere via weathering.

3.4 Planetary evolution pathways

We find that distinct pathways exist for the evolution of a rocky, stagnant-lid planet (Fig. 5), depending on the make-up of its mantle. If the planet's surface temperature remains low (e.g., for wide planetary orbits), outgassed water can start condensing and accumulating on the surface as ice or in liquid form. If CO₂ outgassing is strong early on (e.g., for planets with oxidized mantles), this will quickly result in atmospheric pressures unsuitable for water outgassing. The planet develops a thick, hot, CO₂-rich atmosphere. By contrast, with weak CO₂ outgassing, the planet can accumulate large amounts of water vapor, which will rapidly condense. Further CO₂ outgassing can eventually push these planets into a hothouse regime, where the entire ocean evaporates to form a hot steam atmosphere. In planets with wet mantles and reducing conditions, outgassing of H₂ can achieve the same effect. As discussed above, the weathering mechanisms on stagnant-lid planets do not permit a long-term removal of CO₂, preventing a return to habitable conditions even if all water vapor in the atmosphere were lost (a mechanism that we do not model here). Planets can stay in the habitable regime if the outgassing rates of greenhouse gases remain low over the volcanic lifetime of the planet, but high enough to keep the surface from freezing. This is the case for planets with oxygen fugacities around the IW buffer. Planets with dry, reduced mantles may never advance out of a frozen state due to limited outgassing of H₂ and CO₂.

The presence of primordial atmospheres can reduce the range of habitable conditions even further. Substantial H₂ atmospheres may not be lost quickly enough through atmospheric escape to allow the emergence of habitable surface conditions (Appendix D). While pure steam atmospheres collapse quickly into an ocean, the presence of enough CO₂ in a primordial atmosphere can strongly limit the range of interior conditions which yield habitable planets (Appendix E).

Fig. 3 Evolution of habitable conditions of a stagnant-lid Earth at different orbital distances. Each row depicts the evolution of surface conditions of a planet with Earth mass and core size at the orbit of Venus (0.723 au), Earth (1 au), and Mars (1.524 au), respectively. Each plot shows the prevailing surface conditions for water as a function of the oxygen fugacity and initial water content of the mantle, with the color background showing the area where the majority of neighboring data points share the same surface conditions. Similar to Fig. 2a, each point represents a snapshot of the evolution at a randomly selected planet age within the given age range. The marker in the middle row marks the evolution shown in detail in Figs. 1a, b, c.

3.5 Sensitivity to model parameters

We tested the model sensitivity with respect to changing a few key parameters to confirm the robustness of the results. Figure 6 summarizes this sensitivity analysis for planets with one Earth-mass.

3.5.1 Influence of CO₂ absorption coefficient

Although our adopted value of the CO₂ absorption coefficient $k_{0,{\text{CO}}_2}$ (0.05 m² kg⁻¹) was chosen to reproduce the present-day Earth climate sensitivity (Pujol & North 2003), which describes the response of Earth's climate to a doubling in CO₂, other values have been used in the literature (e.g., Elkins-Tanton 2008). Since the amount of greenhouse heating of CO₂ plays an important role for habitability, we tested the sensitivity of the model to a reduction of $k_{0,{\text{CO}}_2}$ to 0.001 m² kg⁻¹, which is generally used in magma ocean atmosphere modeling (Nikolaou et al. 2019; Elkins-Tanton 2008) and thus provides us with a lower bound on the greenhouse heating from CO₂. As shown in Fig. 3, this slightly extends the range of habitable planets toward more oxidizing conditions, as larger amounts of CO₂ are needed for a planet to be in a hothouse regime due to the lower efficiency of heating. Overall, the influence of the absorption coefficient is fairly minor, which can largely be attributed to the self-regulating nature of the feedback mechanisms.

3.5.2 Influence of ratio of intrusive-to-extrusive volcanism

In the models presented so far, we assumed that all melt produced in the mantle reaches the surface of the planet, where the supersaturated volatile species are outgassed into the atmosphere. However, in general a large part of the magma produced at depth is expected to be intrusive (White et al. 2006), where melt crystallizes within or at the base of existing crust, thus reducing the amount of volatiles that reach the surface. The degree of extrusive volcanism (f_extr) is difficult to constrain and can vary based on location, crust porosity and lithospheric thickness. To model the impact of reduced extrusive volcanism, we ran a model study for an intrusive-to-extrusive ratio of 2.5 (Tosi et al. 2017), corresponding to f_extr ≈ 0.286. The results are shown in Fig. 6c.

Similar to Sect. 3.5.1, the presence of intrusive volcanism extends the range of habitable planets to more oxidizing conditions to a small degree. With intrusive volcanism, the rate of outgassing is reduced. This primarily reduces the rate at which greenhouse gases, specifically CO₂, build up. Therefore, more CO₂ can be outgassed until the planet's atmosphere is too hot to allow for liquid water, which allows the planets to be habitable at more oxidizing conditions.

Fig. 4 As Fig. 3, but showing the prevailing surface conditions of water for different core sizes at different orbital distances. The two markers depict the initial conditions for the detailed evolutions in Figs. 1a–c (1), and d–f (2).

Fig. 5 Evolutionary tracks of water outgassing on stagnant-lid planets. The arrows illustrate potential pathways on which a planet may evolve over time. Only planets with a nonzero amount of water outgassing are shown. The points represent Earth-mass planets with cores sizes ranging from 30% to 70% of the planet's radius at randomly selected times in their evolution, with the color showing the redox state of the mantle.

Fig. 6 Influence of changing the CO2 absorption coefficient, introducing intrusive volcanism, and using temperature-dependent carbon weathering. The reference model shown here is the same as in Fig. 2.

3.5.3 Sensitivity to carbon weathering efficiency

The climate feedback from carbon weathering in the models presented here is fairly weak due to the sole dependence of the weathering model on the CO₂ partial pressure, which we chose to provide a conservative estimate of the emergence of habitable conditions. More efficient feedback mechanisms with a strong dependence on temperature have been suggested (Krissansen-Totton & Catling 2017), which would be able to remove much more CO₂ from the atmosphere and thus extend the range of habitable conditions. To gauge the effect of temperature-dependent weathering, we replaced the ${\left(\frac{P_{{\text{CO}}_2}}{P_{{\text{CO}}_2,\text{E}}}\right)}^{\alpha }$ term in Eq. (11) with a temperature-dependent term exp $\left[\frac{E_{\text{bas}}}{R}\left(\frac{1}{R{T}_{\text{p,E}}}-\frac{1}{R{T}_{\text{p}}}\right)\right]$ following Krissansen-Totton & Catling (2017) and Höning et al. (2019), where E_bas is the activation energy for basalt weathering, and T_p is the pore-space temperature which scales linearly with the surface temperature T_s as follows (Krissansen-Totton & Catling 2017): $$T_{\text{p}}=a_{\text{grad}}{T}_{\text{s}}+b_{\mathrm{int}}+9\text{\hspace{0.17em}K,}$$(20)

with constants a_grad = 1.02 and b_int = −16.7 K. T_p,E is the pore-space temperature based on the current mean surface temperature T_s,E = 285 K of Earth.

We find that the stronger feedback introduced by the temperature-dependence significantly increases the range of habitable planets toward more oxidizing mantle conditions (Fig. 6d), in particular for mantles with low water contents. Compared to the model with CO₂-dependent weathering, the mantle water content plays a larger role. This is mainly due to the lower rate of volcanism on planets with drier mantles. The temperaturedependent weathering feedback works to establish a CO₂ pressure in the atmosphere where the surface temperature is high enough so that the ingassed CO₂ flux from weathering is in steady-state with outgassing fluxes from volcanism and decarbonation (see also Höning et al. 2019). If the rate of volcanism is large (e.g., for high mantle water contents), the rates of outgassed and decarbonized CO₂ are also large (since the latter is directly tied to weathering rates). In this case, a steady-state is not possible, as the required steady-state temperature would be too high for liquid water. Weathering therefore cannot prevent a hothouse atmosphere. On the other hand, if the rate of volcanism is small, weathering can remove any excess CO₂ at temperature ranges which permit liquid water on the surface. A comparison of the surface temperature and pressure evolution for both weathering models is shown in Fig. G.7.

Overall, we find that the results and evolutionary paths shown above still hold well even with a stronger weathering feedback, although the exact range of $f_{{\text{O}}_2}$ yielding habitable planets in Fig. 6 is different. As a particular side effect of the temperature-dependence of weathering, the CO₂ content in the atmosphere is no longer driven solely by the rate of volcanism, but also by the increase in solar luminosity and other factors that increase the surface temperature of the planet, such as the water content of the atmosphere. We refer to Höning et al. (2019) for a more detailed discussion of the two weathering models for stagnant-lid planets.

4 Discussion and conclusions

The structure and composition of the interior are fundamental factors to determine whether a planet can be habitable or not. The mantle redox state in particular strongly constrains the space of potentially habitable planets. In general, it is difficult for planets with stagnant lids to remain habitable because of the limited number of pathways available to permanently remove outgassed CO₂ from the atmosphere. In fact, CO₂ tends to accumulate and heat the planet up until the atmosphere is too hot to allow liquid water. We modeled the (temporary) removal of CO₂ via silicate weathering, which is active only in the presence of liquid water, thus placing further limits on the removal of CO₂. Additional CO₂ could be lost through drag from escaping hydrogen (Tian 2015; Hunten et al. 1987) or through photodissociation in the upper atmosphere, processes we did not model here. Stagnantlid planets are common in the Solar System. Earth alone is in a plate tectonic regime. If this is also the case for rocky exoplanets, we would expect many of those to more closely resemble Venus than Earth, with hot, dense atmospheres even if they reside in the habitable zone of their host star.

Our simple atmosphere model cannot capture the full complexity of a planetary atmosphere. With more sophisticated atmospheric models (Scheucher et al. 2020; Wunderlich et al. 2020; Kaspi & Showman 2015; Schreier et al. 2014), the surface temperature would likely differ to some extent from the one calculated here, thus affecting both the habitability of planets and the critical amount of CO₂ or H₂ at which the planet transitions into a hot, nonhabitable regime. However, as discussed above, the main driver for the accumulation of surface water is the outgassing rate of CO₂ and H₂, which is mainly a function of the planetary interior. We note here that the outgassing rates are subject to the composition of the atmosphere, which may be different from the outgassed species due to atmospheric chemistry, which we do not take into account here. However, the solubilities of both CO₂ and H₂ in the melt are very low, and therefore these two species are least affected by partial pressures during outgassing. As such, while a difference in surface temperature could change the exact values of oxygen fugacities and water concentration in the mantle which yield habitable conditions, the general relations described above would still hold. As seen in Sect. 3.5.1, even a drastically reduced IR absorption coefficient of CO₂ has only a small effect upon the proportions of planets with surface water.

Here we considered the internal structure of a planet to be independent of its redox state. Yet, the latter and the size of the metallic core may well evolve jointly, based on the local oxidation level of the protoplanetary disk during formation, on the conditions of metal-silicate differentiation, and on the subsequent evolution of the magma ocean, complex processes whose mutual relations are still to be fully unraveled (Wade & Wood 2005; Frost & McCammon 2008; Zhang et al. 2017; Armstrong et al. 2019). Planets with deep magma oceans may develop rather oxidizing mantles (Deng et al. 2020), which would be less favored to develop long-term habitable conditions based on our results. In fact, our results indicate that a stagnant-lid Earth or Venus, having more oxidizing conditions, will always enter a hothouse state. In contrast, low-mass planets with large iron cores would likely have more reducing conditions due to more shallow magma oceans, which in turn would strongly favor long-term habitable conditions. Liggins et al. (2022) show that the mantle redox state imposes characteristic atmospheric compositions. The atmospheric composition of such planets is potentially detectable in exoplanets in the near-future via spectroscopic observations with instruments such as JWST (Greene et al. 2016). Our results suggest that planets which are habitable should generally have reduced atmospheric compositions. The unambiguous direct detection of an exoplanetary ocean is challenging, relying on direct imaging (Lustig-Yaeger et al. 2018), and is therefore likely to be limited to only a few planets in the foreseeable future. Promising first targets would therefore be small, dense exoplanets, which may offer the best chances in the pursuit to find planets capable of hosting liquid water on their surface.

Acknowledgements

We thank Brad Foley for his thorough and constructive review that helped to significantly improve an earlier version of this manuscript. PB and NT acknowledge support of the German Science Foundation (DFG) through the priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (TO 704/3-1) and the research unit FOR 2440 "Matter under planetary interior conditions" (PA 3689/1-1).

Appendix A 1D parameterized convection model

Fig. A.1 Schematic of the interior structure used in the thermal evolution model, alongside a diagram of the temperature profile (see the text for the explanation of the various symbols).

We modeled the thermal evolution of a planet's mantle by considering the energy balance between heat lost through the planetary surface and heat entering the mantle from the iron core, as well as the decay of radiogenic elements inside the mantle. Assuming the core to be fully liquid and convecting, its energy balance is given by $${\rho }_{\text{c}}{c}_{\text{c}}{V}_{\text{c}}\frac{\text{d}{\overline{T}}_{\text{c}}}{\text{d}t}=-q_{\text{c}}{A}_{\text{c}},$$(A.1)

where ${\overline{T}}_{\text{c}}$ is the average temperature in the core; ρ_c, c_c, and V_c are the density, specific heat capacity, and volume of the core; A_c is the area of the core-mantle boundary (CMB); and q_c is the heat flux out of the core.

The energy balance of the mantle is given by $$\begin{array}{l}{\rho }_{\text{m}}{c}_{\text{m}}{V}_{\text{m}}\left(1+St\right)\frac{\text{d}{\overline{T}}_{\text{m}}}{\text{d}t}=\hfill \\ -\left(q_1+\left[{\rho }_{\text{cr}}ℒ+{\rho }_{\text{cr}}{c}_{\text{cr}}\left(T_{\text{m}}-T_1\right)\right]\frac{\text{d}{D}_{\text{cr}}}{\text{d}t}\right)A_{\text{m}}+q_{\text{b}}{A}_{\text{b}}+Q_{\text{m}}{V}_{\text{m}},\hfill \end{array}$$(A.2)

where ${\overline{T}}_{\text{m}}$ is the average temperature of the convecting mantle; ρ_m and c_m are the density and specific heat capacity of the mantle; V_m and A_m are the volume and area of the convecting part of the mantle; St is the Stefan number, which describes the energy consumed and released upon mantle melting and solidification; q₁ and q_b are the heat fluxes out of and into the convecting part of the mantle respectively; T_m and T_l are the temperatures of the upper mantle and the bottom of the stagnant lid; ρ_cr, c_cr, and D_cr are the density, specific heat capacity, and thickness of the crust; ℒ is the latent heat of melting; and Q_m is the volumetric heating rate in the mantle from heat-producing elements. The main parameters used in this model are summarized in Table F.2.

The temperatures at the top of the core and mantle T_c and T_m are related to the volume-averaged temperatures ${\overline{T}}_{\text{c}}$ and ${\overline{T}}_{\text{m}}$ through scaling factors ɛ_c and ɛ_m $$\begin{array}{cc}{\overline{T}}_{\text{c}}={\epsilon }_{\text{c}}{T}_{\text{c}},& {\overline{T}}_{\text{m}}={\epsilon }_{\text{m}}{T}_{\text{m}}=\frac{1}{V_{\text{m}}}{\displaystyle \underset{V_{\text{m}}}{\int }{T}_{\text{ad}}}\left(r\right)\end{array}\text{d}V,$$(A.3)

where T_ad (r) is the adiabatic temperature profile in the mantle. For the core, we set ɛ_c = 1.2 following Stamenković et al. (2012), who found only a small dependence on planetary mass. For the mantle, ɛ_m is updated continuously during the mantle evolution based on the mantle temperature profile and on the varying thickness of the stagnant lid.

The evolution of the thickness of the stagnant lid (D₁) follows from the energy balance at the base of the lid (at radius R₁) $$\begin{array}{l}{\rho }_{\text{m}}{c}_{\text{m}}{V}_{\text{m}}\left(T_{\text{m}}-T_1\right)\frac{\text{d}{D}_1}{\text{d}t}=\hfill \\ -q_1+\left[{\rho }_{\text{cr}}ℒ+{\rho }_{\text{cr}}{c}_{\text{cr}}\left(T_{\text{m}}-T_1\right)\right]{\frac{\text{d}{D}_{\text{cr}}}{\text{d}t}-k_{\text{m}}\frac{\partial T}{\partial r}|}_{r=R_1}.\hfill \end{array}$$(A.4)

${\partial T/\partial r|}_{r=R_1}$ is the temperature gradient at the base of the lid, which we calculate assuming steady-state heat conduction: $$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{k}_1\frac{\partial T}{\partial r}\right)=-Q_1,$$(A.5)

where k₁ is the thermal conductivity and Q₁ the heat production rate in the lid. Parts of the stagnant lid can be comprised of crustal material, so when solving Eq. (A.5), we set k₁ and Q₁ to crustal values (k_cr and Q_cr) from the surface to the base of the crust or mantle values (k_m and Q_m) from the base of the crust to the base of the stagnant lid as appropriate. The boundary conditions for Eq. A.5 are the given surface temperature T_s and the temperature at the stagnant lid base T₁.

Numerical convection models suggest that the viscosity contrast across the upper thermal boundary layer is typically about one order of magnitude (Grasset & Parmentier 1998). Based on this, the lid base temperature T₁ can be calculated from the mantle temperature and activation energy (Grasset & Parmentier 1998) $$T_1=T_{\text{m}}-2.9\frac{R{T}_{\text{m}}^2}{E^{*}},$$(A.6)

where the factor of 2.9 accounts for the effects of spherical geometry (Reese et al. 2005).

The convective heat flux out of the mantle q₁, assuming that the upper thermal boundary layer is small so that the radial temperature profile is close to linear, can then be expressed as $$q_1=k_{\text{m}}\frac{T_{\text{m}}-T_1}{d_{\text{u}}},$$(A.7)

where the thickness of the upper TBL d_u can be calculated from boundary layer theory $$d_{\text{u}}=D_{\text{m}}{\left(\frac{R{a}_{\text{crit}}}{Ra}\right)}^{1/3},$$(A.8)

where D_m = R₁ – R_b is the thickness of the convecting part of the mantle, Ra_crit is the critical Rayleigh number of the mantle, and Ra is the Rayleigh number for the entire mantle $$Ra=\frac{{\alpha }_{\text{m}}\left(P_{\text{m}},T_{\text{m}}\right){\rho }_{\text{m}}g\text{Δ}T{D}_{\text{m}}^3}{{\kappa }_{\text{m}}{\eta }_{\text{m}}},$$(A.9)

with the mantle viscosity η_m, mantle thermal diffusivity K_m = k_m/(ρ_mc_m), the pressure- and temperature-dependent coefficient of thermal expansion α_m with the pressure at the top of the mantle P_m = ρ_mg(D₁ + d_u), and ΔT = (T_m – T₁) + (T_c − T_b), which is the sum of temperature differences across both boundary layers.

The coefficient of thermal expansion α, which has a strong influence on the heat transport in the convecting mantle, is often assumed to be constant, although it is known from experimental data that this parameter can vary considerably with both pressure and temperature (Fei & Ahrens 1995). This becomes especially important if one considers the modeling of super-Earths (Miyagoshi et al. 2018). We used here the temperature- and pressure-dependent parameterization of α from Tosi et al. (2013), $$\alpha \left(P,T\right)=\left(a_0+a_{1}T-a_2{T}^{-2}\right)\exp \left(-a_{3}P\right),$$(A.10)

where the coefficients a₀ - a₃ are chosen assuming a lower mantle composition of 80% perovskite/20% periclase (see Table F.2 for parameter values). The temperature profile in the convecting part of the mantle is assumed to be adiabatic: $$\frac{\text{d}T}{\text{d}P}=\frac{\alpha \left(P,T\right)}{{\rho }_{\text{m}}{c}_{\text{m}}}T.$$(A.11)

To account for a potentially nonconvecting zone near the CMB due to the effect of high pressures on the mantle viscosity, we used the parameterization from Stamenković et al. (2012). Assuming that this conductive layer is close to convective stability, the thickness can be approximated from boundary layer theory using a critical Rayleigh number $R{a}_{\text{crit}}^{\text{CMB}}$ based on the viscosity contrast Δ_η across the layer: $$R{a}_{\text{crit}}^{\text{CMB}}\left(\text{Δ}\eta \right)=\max \left\{R{a}_{\text{crit}},11.74\cdot \ln {\left(\text{Δ}\eta \right)}^4\right\},$$(A.12)

with $$\text{Δ}\eta \text{\hspace{0.17em}=\hspace{0.17em}max\hspace{0.17em}}\left\{\frac{\eta \left(R_{\text{c}}\right)}{\eta \left(R_{\text{b}}\right)},\frac{\eta \left(R_{\text{b}}\right)}{\eta \left(R_{\text{c}}\right)}\right\}.$$(A.13)

Here, R_b is the radius at the top of the conductive layer, and R_c is the radius of the CMB.

The thickness of this layer is then given by: $$d_{\text{b}}={\left(R{a}_{\text{crit}}^{\text{CMB}}\left(\text{Δ}\eta \right)\frac{{\kappa }_{\text{m}}\min \left\{\eta \left(R_{\text{b}}\right),\eta \left(R_{\text{c}}\right)\right\}}{{\alpha }_{\text{m}}\left(P_{\text{b}},T_{\text{b}}\right){\rho }_{\text{m}}g\left|T_{\text{c}}-T_{\text{b}}\right|}\right)}^{1/3}.$$(A.14)

Especially for planets more massive than Earth, the conductive CMB layer can make up a significant part of the mantle. Therefore, we treated the heat fluxes from the CMB lid into the mantle and from the core into the CMB lid separately, and assume time-dependent thermal conduction across the layer. The heat fluxes are given by the temperature gradients at the top and bottom of the conductive CMB layer: $$q_{\text{b}}={-k_{\text{m}}\frac{\partial T}{\partial r}|}_{r=R_{\text{b}}},$$(A.15) $$q_{\text{c}}={-k_{\text{m}}\frac{\partial T}{\partial r}|}_{r=R_{\text{c}}},$$(A.16)

which we determined by solving the time-dependent heat equation across the CMB lid: $$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{k}_{\text{m}}\frac{\partial T}{\partial r}\right)=-Q_{\text{m}}+{\rho }_{\text{m}}{c}_{\text{m}}\frac{\partial T}{\partial t}.$$(A.17)

Appendix B Melting, trace element partitioning, and volatile outgassing

We computed the distribution of partial melt in the mantle by comparing the local mantle temperature profile T(r) against the solidus T_sol(r) and liquidus T_liq(r) temperatures. We assumed the amount of partial melt to vary linearly between the solidus and the liquidus: $$\varphi \left(r\right)=\frac{T\left(r\right)-T_{\text{sol}}\left(r\right)}{T_{\text{liq}}\left(r\right)-T_{\text{sol}}\left(r\right)}.$$(B.1)

We did not consider melting above a pressure of 8 GPa. Under these conditions, melt becomes denser than the surrounding mantle rocks and cannot reach the surface (Agee 2008).

The presence of water in the mantle depresses the solidus and liquidus curves. We calculated wet solidus and liquidus curves following a parameterization by Katz et al. (2003). The volumeaveraged, extractable melt fraction $\overline{\varphi }$ in the mantle is then given by $$\overline{\varphi }=\frac{1}{V_{\varphi }}{\displaystyle \underset{V_{\varphi }}{\int }\varphi \left(r\right)\text{\hspace{0.17em}d}V},$$(B.2)

where V_φ is the total volume of the melt zone (i.e., where the temperature lies above the solidus).

Knowing the volume of melt produced, we can calculate the evolution of the crustal thickness D_cr. We adopted the plume model description by Grott et al. (2011). Partial melting in the mantle is generally restricted to localized upwelling plumes. We assumed a plume covering fraction of f_p = 0.01, and added the temperature difference across the bottom thermal boundary layer to the local temperature profile when evaluating the melt fraction in Eq. (B.1). In addition to the amount of available melt, the crustal growth rate depends on the rate at which fresh mantle material can be supplied to the partial melt zone, which is governed by the convective velocity u of the mantle. The crustal growth rate is given by $$\frac{\text{d}{D}_{\text{cr}}}{\text{d}t}=f_{\text{p}}u\overline{\varphi }\frac{V_{\varphi }}{4\pi R_{\text{p}}^3},$$(B.3)

where the convective velocity is $$u=u_0{\left(\frac{Ra}{R{a}_{\text{crit}}}\right)}^{2/3},$$(B.4)

where u₀ is a characteristic mantle convective velocity scale. We imposed the additional constraint that the crust cannot grow larger than the lid. Once the crust reaches the bottom of the lid, any excess crust is recycled back into the mantle, and the crustal growth rate is set to be equal to the lid growth rate. During crustal formation, we treated the release and consumption of latent heat during mantle melting and crystallization via the Stefan number, which is recalculated at every time step (See also Eq. A.2) $$St=\frac{ℒ}{c_{\text{m}}}\frac{V_{\varphi }}{V_{\text{m}}}\frac{\text{d}\overline{\varphi }}{\text{d}{T}_{\text{m}}}.$$(B.5)

We considered the partitioning of radiogenic elements and water between crust and mantle due to melt production and crust formation, and the subsequent enrichment of the crust in these elements We considered a model of fractional melting to calculate the partitioning of trace elements between melt and mantle rocks. The concentration in the melt $X_{\text{liq}}^i$ of a given trace element i at radius r is then given by $$X_{\text{liq}}^i\left(r\right)=\frac{X_{\text{m}}^i}{\varphi \left(r\right)}\left[1-{\left(1-\varphi \left(r\right)\right)}^{1/{\delta }_i}\right],$$(B.6)

where $X_{\text{m}}^i$ is the corresponding bulk concentration in the mantle and δ_i a trace-element-specific partition coefficient. We assumed δ_i = 0.001 for heat-producing elements (Blundy & Wood 2003), and δ_i = 0.01 for water (Aubaud 2004). The average concentration in the melt then follows as $${\overline{X}}_{\text{liq}}^i=\frac{1}{\overline{\varphi }{V}_{\varphi }}{\displaystyle \underset{V_{\varphi }}{\int }{X}_{\text{liq}}^i\left(r\right)\varphi \left(r\right)\text{d}V}.$$(B.7)

Enriched melt is transported to the surface and forms a crust. The total mass of an incompatible element $M_{\text{cr}}^i$ in the crust is given by the crust production rate and the average concentration in the melt $$\frac{\text{d}{M}_{\text{cr}}^i}{\text{d}t}=4\pi R_{\text{p}}^2{\rho }_{\text{cr}}{\overline{X}}_{\text{liq}}^i\frac{\text{d}{D}_{\text{cr}}}{\text{d}t}.$$(B.8)

At the same time, the mantle will be depleted in trace elements accordingly, with the concentration of the trace elements in the mantle and crust given by $$X_{\text{m}}^i=\frac{X_{\text{m,0}}^i{M}_{\text{m,0}}-M_{\text{cr}}^i}{M_{\text{m}}},X_{\text{cr}}^i=\frac{M_{\text{cr}}^i}{M_{\text{cr}}},$$(B.9)

where M_m,0 and M_m are the initial and current mass of the mantle, respectively, and $X_{\text{m,0}}^i$ is the initial mantle concentration of the respective trace element.

The enrichment of heat-producing elements in the crust and depletion in the mantle leads to different volumetric heating rates in crust and mantle, which can be calculated as follows $$Q_{\text{m}}\left(t\right)={\rho }_{\text{m}}{\displaystyle \sum _i{X}_{\text{m}}^i\left(t\right)H_i\exp \left(\frac{\ln 2\cdot \left(4.5\text{\hspace{0.17em}Gyr}-t\right)}{{\tau }_{1/2}^i}\right)},$$(B.10) $$Q_{\text{cr}}\left(t\right)={\rho }_{\text{cr}}{\displaystyle \sum _i{X}_{\text{cr}}^i\left(t\right)H_i\exp \left(\frac{\ln 2\cdot \left(4.5\text{\hspace{0.17em}Gyr}-t\right)}{{\tau }_{1/2}^i}\right)},$$(B.11)

where i specifies one of the four long-lived radioisotopes ⁴⁰K, ²³²Th, ²³⁵U and ²³⁸U, with corresponding specific heat production rates H_i and half-lives ${\tau }_{1/2}^i$ based on bulk-silicate-Earth abundances from McDonough & Sun (1995).

Not all melt produced in the mantle is able to reach the surface, but is instead intruded into the lid, solidifies there, and is therefore unavailable for outgassing. To model this, we need to assume a fraction of extrusive volcanism f_extr, which we set to 1 in this study (i.e., all melt reaches the surface). While this is not fully realistic for an Earth-like planet, it provides an upper limit for the outgassed species. In Sect. 3.5.2 in the main text, we also tested the model results with a more realistic value (Tosi et al. 2017) of f_extr = 0.286 and show that this ultimately serves to further increase the number of habitable planets.

Volatile species will be partially outgassed into the atmosphere once the melt reaches the surface. To outgas a volatile species, the melt needs to be supersaturated with respect to the atmosphere, and any excess concentration can be released into the atmosphere. The redox state of the mantle plays a large role regarding which volatile species are outgassed, with an oxidized mantle favoring H₂O and CO₂, while a reduced mantle is dominated by H₂ and CO outgassing (Ortenzi et al. 2020). The chemical outgassing model based on Ortenzi et al. (2020) (as described in Sect. 2.3) calculates the outgassed mass fraction $X_{\text{outg}}^i$ of volatile species based on the chemical equilibrium between melt and atmosphere, taking into account the solubility of each species. The outgassed mass $M_{\text{outg}}^i$ of volatile species is then given by $$\frac{\text{d}{M}_{\text{outg}}^i}{\text{d}t}=f_{\text{extr}}{X}_{\text{outg}}^i\frac{\text{d}{M}_{\text{cr}}^i}{\text{d}t}.$$(B.12)

We can then calculate the partial pressure of each species according to its atmospheric mass and the presence of other species (e.g., Bower et al. 2019). The actual mass enriched in the crustin Eq. (B.8) is then reduced by the outgassed mass. We also assumed that all surface volcanism takes place above any potential ocean surface, since the pressure at the bottom of the ocean would limit outgassing. Ocean coverage and depth are difficult to estimate as they are dependent on surface topography. Similar to the assumption of fully extrusive volcanism, this assumption provides an upper limit to volatile outgassing. Likewise, we did not take into account so-called water worlds, which are planets with several tens of kilometers of water oceans. Due to the high pressures at the ocean bottom, volatile outgassing would likely be severely limited (Noack et al. 2016; Krissansen-Totton et al. 2021).

Appendix C Venus-like planets

Venus is the quintessential hothouse planet. In order to test the ability of our model to reproduce a Venus-like scenario, we simulated the evolution of 5000 planets with Venus-like interior structures and orbital distance, while varying the mantle water content and oxygen fugacity as described in the methods section 2.9. We find that at present day (4.5 Gyr) all modeled planets are in an extreme greenhouse state. No habitable planets are present (Fig. C.1a), and surface pressures in planets with mantle oxygen fugacities above the IW buffer are comparable to those of present-day Venus (Fig. C.1b). Our model tends to overestimate the surface temperature compared to actual Venus. This stems from the fact that we do not consider the loss of water through photodissociation, which leaves water in the atmosphere as a potent greenhouse gas. Many of the water-rich atmospheres shown in Fig. C.1 would evolve into thick, dry, CO₂-dominated atmospheres. Furthermore, Venus' high albedo due to its global cloud cover is not represented in the model, which contributes to explaining the higher surface temperatures. A more in-depth discussion of the evolution of Venus using the outgassing model discussed here can be found in Höning et al. (2021).

Fig. C.1 Thermal (a) and pressure (b) state of Venus-like planets at present day, after 4.5 Gyr of evolution. The color map indicates the oxygen fugacity of the mantle. The red dashed lines mark the present-day surface temperature and atmospheric pressure of Venus.

Appendix D Influence of primordial H₂ atmospheres

So far, we have assumed that any primordial atmosphere is lost at the point when our evolution models start. This provides us with an upper limit to any outgassed secondary atmosphere. However, while the life-time of primordial atmospheres can be very short, on the order of tens of millions of years, especially for close-in, low-mass planets (Kite & Barnett 2020; Lammer et al. 2014), thick H₂ atmospheres may survive magma ocean solidification. We investigated the influence of a primordial H₂ atmosphere by running a number of evolution models of Earth-like planets with initial atmospheric pressures of up to 300 bar. This upper limit is motivated by the amount of hydrogen an Earth-like planet may accrete on formation (Lammer et al. 2014) and by the maximum amount it can lose over time so that most planets we consider here are still rocky (Howe et al. 2020). This results in a reasonable range of different planet evolutions, but different choices could change the proportions of planets with different atmospheres. Low-mass planets in particular may not be able to accrete a hydrogen atmosphere of that extent during formation. We find that the presence of primordial H₂ atmospheres with pressures above ≈ 50 – 100 bar can significantly reduce the amount of habitable planets (Fig. D.1). This is due to two main factors: First, sufficiently thick H₂ atmospheres may not be lost completely, and thus the planet never reaches surface temperatures that are low enough to allow for liquid water. Second, it can easily take several hundred million years for extensive primordial H₂ atmospheres to be completely lost. At this early stage, a planet is volcanically very active, but the outgassed CO₂ is not removed from the atmosphere because the existing H₂ atmosphere inhibits the carbon-silicate cycle that requires liquid water. In addition, additional outgassing of H₂ can further prolong the lifetime of an H₂ atmosphere. As a result, even though the H₂ atmosphere may be eventually lost, too much CO₂ has been outgassed during that time to allow habitable conditions. The planets which may avoid a hothouse in these cases are those with very low amounts of CO₂ outgassing, that is, planets with dry mantles and oxygen fugacities well below the IW buffer.

Fig. D.1 Influence of a primordial H2 atmosphere on the emergence of habitable surface conditions. Each column corresponds to different initial pressures of H2, with each row marking planets at different orbital distances to their host star.

Appendix E Influence of primordial steam and CO₂ atmospheres

It is likely that a magma ocean would form a thick steam atmosphere, although these may collapse shortly after magma ocean solidification to form an early ocean (Elkins-Tanton 2011). To determine the influence of an early post magma ocean atmosphere, we investigated two end member compositions: Pure steam atmospheres up to 200 bar, and pure CO₂ atmosphere up to 5 bars. Pure steam atmospheres have little impact on the long-term evolution of the planets (Fig. E.1). Due to the positive climate feedback of water vapor, these atmospheres are not stable, and other sources of heating, such as from other greenhouse gases, are needed to sustain a steam atmosphere. Shortly after the start of the evolution, these atmospheres therefore collapse and rain out to form an ocean. They do not contribute further to the warming of the surface, but merely provide a reservoir of water. On the other hand, even small amounts of initial CO₂ atmospheres can severely limit the occurrence of habitable conditions (Fig. E.2). In our model, CO₂ is only removed via silicate weathering, which requires the presence of liquid water. If the initial amount of CO₂ in the atmosphere is too high to permit liquid water, there exists no pathway for a planet to lose CO₂ and become habitable. Therefore, the CO₂ pressures given in Fig. E.2 represent a conservative estimate on the amount of CO₂ which can still yield habitable conditions.

Fig. E.1 Influence of a primordial steam atmosphere on the emergence of habitable surface conditions. Each column of plots corresponds to different initial pressures of H2O, with each row marking planets at different orbital distances to their host star.

Fig. E.2 Influence of a primordial CO2 atmosphere on the emergence of habitable surface conditions. Each column of plots corresponds to different initial pressures of CO2, with each row marking planets at different orbital distances to their host star.

Appendix F Tables

Appendix G Additional figures

Fig. G.1 Mass-radius range of terrestrial planets considered in this study. The color map indicates the thickness of the iron core relative to the planet radius.

Fig. G.2 Amount of atmospheric H2 needed to keep the surface temperature at 280 K as a function of orbital distance (after Pierrehumbert & Gaidos 2011). Dashed lines show the results from Pierrehumbert & Gaidos (2011), solid lines show the results of our atmosphere model, assuming a hydrogen absorption coefficient of $k_{0,{\text{H}}_2}=2\times {10}^{-2}{\text{m}}^2{\text{kg}}^{-1}$.

Fig. G.3 Example of limited weathering due to limited crustal growth. The blue dashed line shows the partial pressure of CO2 which would be necessary for weathering to be in equilibrium with the influx of CO2 from outgassing and decarbonation, based on the current crustal growth rate. The orange line shows the actual CO2 pressure in the atmosphere. The blue shaded region marks where the system is in disequilibrium, where the influx of C02 is larger than can be removed by weathering. This example correspond to the Earth-like model in Fig. 1a–c.

Fig. G.4 Example of the "memory effect" present in the stagnant-lid carbon cycle. The blue line shows the age of the crustal layers upon reaching the decarbonation depth zdecarb, and therefore the time that has passed since the weathering rate was "recorded" into the crust. The shaded region marks where the crust is not yet old enough to reach the decarbonation depth. At 2.3 Gyr, the surface begins to heat up and thus younger and younger crust decarbonates as the decarbonation depth decreases. This example correspond to the Earth-like model shown in Fig. 1a–c.

Fig. G.5 Example of the composition of outgassed species as a function of oxygen fugacity. Assumed here is a 1 bar buffer atmosphere, and a volatile content in the melt of 1000 ppm H2O and CO2, respectively, at a melt temperature of 1500 K.

Fig. G.6 As Fig. 3, but showing the prevailing surface conditions for water for different planet masses and core sizes.

Fig. G.7 As in Fig. 1, but showing a comparison between a CO2-dependent (a, b, c) and a temperature-dependent weathering model (d, e, f). All other model parameters are the same as in Fig. 1a–c.

References

All Tables

All Figures

Fig. 1 Characteristic evolutions of the atmospheric temperature and pressure and CO2 flux in and out of the atmosphere for an Earth-mass planet at 1 au with an Earth-like (a–c) or a Mercury-like (d–f) interior structure. The colored background marks the state of water at the planet's surface. Both planets start with the same mantle water content of 250 ppm and an oxygen fugacity of 0.05 log10 units above the IW buffer. Panels (a–c) and (d–f) correspond respectively to the points 1 and 2 marked in Figs. 3 and 4.

In the text

Fig. 2 States of a stagnant-lid Earth as a function of the mantle oxygen fugacity. Each point represents a snapshot of the planetary evolution at a randomly selected planet age between 100 Myr and 8 Gyr. The colored points in panel a show the prevailing water phase at the surface as a function of the oxygen fugacity ($f_{{\text{O}}_2}$) and initial water content of the mantle ($X_{{\text{H}}_2\text{O}}^0$). The color background shows the area where the majority of neighboring data points share the same surface conditions. In panels b-d, points are colored according to the greenhouse gas that contributes most to surface heating. Panel b shows the surface temperature, panel c shows the total mass of outgassed water, and panel d shows the total mass of CO2 removed from the atmosphere via weathering.

In the text

Fig. 3 Evolution of habitable conditions of a stagnant-lid Earth at different orbital distances. Each row depicts the evolution of surface conditions of a planet with Earth mass and core size at the orbit of Venus (0.723 au), Earth (1 au), and Mars (1.524 au), respectively. Each plot shows the prevailing surface conditions for water as a function of the oxygen fugacity and initial water content of the mantle, with the color background showing the area where the majority of neighboring data points share the same surface conditions. Similar to Fig. 2a, each point represents a snapshot of the evolution at a randomly selected planet age within the given age range. The marker in the middle row marks the evolution shown in detail in Figs. 1a, b, c.

In the text

	Fig. 4 As Fig. 3, but showing the prevailing surface conditions of water for different core sizes at different orbital distances. The two markers depict the initial conditions for the detailed evolutions in Figs. 1a–c (1), and d–f (2).
In the text

Fig. 5 Evolutionary tracks of water outgassing on stagnant-lid planets. The arrows illustrate potential pathways on which a planet may evolve over time. Only planets with a nonzero amount of water outgassing are shown. The points represent Earth-mass planets with cores sizes ranging from 30% to 70% of the planet's radius at randomly selected times in their evolution, with the color showing the redox state of the mantle.

In the text

	Fig. 6 Influence of changing the CO2 absorption coefficient, introducing intrusive volcanism, and using temperature-dependent carbon weathering. The reference model shown here is the same as in Fig. 2.
In the text

	Fig. A.1 Schematic of the interior structure used in the thermal evolution model, alongside a diagram of the temperature profile (see the text for the explanation of the various symbols).
In the text

	Fig. C.1 Thermal (a) and pressure (b) state of Venus-like planets at present day, after 4.5 Gyr of evolution. The color map indicates the oxygen fugacity of the mantle. The red dashed lines mark the present-day surface temperature and atmospheric pressure of Venus.
In the text

	Fig. D.1 Influence of a primordial H2 atmosphere on the emergence of habitable surface conditions. Each column corresponds to different initial pressures of H2, with each row marking planets at different orbital distances to their host star.
In the text

	Fig. E.1 Influence of a primordial steam atmosphere on the emergence of habitable surface conditions. Each column of plots corresponds to different initial pressures of H2O, with each row marking planets at different orbital distances to their host star.
In the text

	Fig. E.2 Influence of a primordial CO2 atmosphere on the emergence of habitable surface conditions. Each column of plots corresponds to different initial pressures of CO2, with each row marking planets at different orbital distances to their host star.
In the text

	Fig. G.1 Mass-radius range of terrestrial planets considered in this study. The color map indicates the thickness of the iron core relative to the planet radius.
In the text

	Fig. G.2 Amount of atmospheric H2 needed to keep the surface temperature at 280 K as a function of orbital distance (after Pierrehumbert & Gaidos 2011). Dashed lines show the results from Pierrehumbert & Gaidos (2011), solid lines show the results of our atmosphere model, assuming a hydrogen absorption coefficient of $k_{0,{\text{H}}_2}=2\times {10}^{-2}{\text{m}}^2{\text{kg}}^{-1}$.
In the text

Fig. G.3 Example of limited weathering due to limited crustal growth. The blue dashed line shows the partial pressure of CO2 which would be necessary for weathering to be in equilibrium with the influx of CO2 from outgassing and decarbonation, based on the current crustal growth rate. The orange line shows the actual CO2 pressure in the atmosphere. The blue shaded region marks where the system is in disequilibrium, where the influx of C02 is larger than can be removed by weathering. This example correspond to the Earth-like model in Fig. 1a–c.

In the text

Fig. G.4 Example of the "memory effect" present in the stagnant-lid carbon cycle. The blue line shows the age of the crustal layers upon reaching the decarbonation depth zdecarb, and therefore the time that has passed since the weathering rate was "recorded" into the crust. The shaded region marks where the crust is not yet old enough to reach the decarbonation depth. At 2.3 Gyr, the surface begins to heat up and thus younger and younger crust decarbonates as the decarbonation depth decreases. This example correspond to the Earth-like model shown in Fig. 1a–c.

In the text

	Fig. G.5 Example of the composition of outgassed species as a function of oxygen fugacity. Assumed here is a 1 bar buffer atmosphere, and a volatile content in the melt of 1000 ppm H2O and CO2, respectively, at a melt temperature of 1500 K.
In the text

	Fig. G.6 As Fig. 3, but showing the prevailing surface conditions for water for different planet masses and core sizes.
In the text

	Fig. G.7 As in Fig. 1, but showing a comparison between a CO2-dependent (a, b, c) and a temperature-dependent weathering model (d, e, f). All other model parameters are the same as in Fig. 1a–c.
In the text