Habitability of Earth-like stagnant lid planets: Climate evolution and recovery from snowball states

TL;DR

This study models climate evolution on Earth-like stagnant lid planets, showing their potential for long-term habitability and recovery from snowball states depending on CO2 budgets and atmospheric exchange, with implications for exoplanet observations.

Contribution

It introduces a coupled model of mantle and climate evolution for stagnant lid planets, analyzing habitability and snowball recovery mechanisms, which is novel in planetary climate studies.

Findings

Habitability lasts 1-5 Gyrs across a range of CO2 budgets.

Recovery from snowball states depends on atmosphere-ocean exchange.

A minimum CO2 budget is required for snowball recovery.

Abstract

Coupled models of mantle thermal evolution, volcanism, outgassing, weathering, and climate evolution for Earth-like (in terms of size and composition) stagnant lid planets are used to assess their prospects for habitability. The results indicate that planetary CO$_2$ budgets ranging from $\approx 3$ orders of magnitude lower than Earth's to $\approx 1$ order of magnitude larger, and radiogenic heating budgets as large or larger than Earth's, allow for habitable climates lasting 1-5 Gyrs. The ability of stagnant lid planets to recover from potential snowball states is also explored; recovery is found to depend on whether atmosphere-ocean chemical exchange is possible. For a "hard" snowball with no exchange, recovery is unlikely, as most CO$_2$ outgassing takes place via metamorphic decarbonation of the crust, which occurs below the ice layer. However, for a "soft" snowball where there is exchange between atmosphere and ocean, planets can readily recover. For both hard and soft snowball states, there is a minimum CO$_2$ budget needed for recovery; below this limit any snowball state would be permanent. Thus there is the possibility for hysteresis in stagnant lid planet climate evolution, where planets with low CO$_2$ budgets that start off in a snowball climate will be permanently stuck in this state, while otherwise identical planets that start with a temperate climate will be capable of maintaining this climate for 1 Gyrs or more. Finally, the model results have important implications for future exoplanet missions, as they can guide observations to planets most likely to possess habitable climates.

Tables

Table 1. Table 1: Model Parameters and Assumed Values

Parameter	Meaning	Assumed baseline value	Equation
$\rho$	Mantle density	4000 kg m-3	(1)
$c_p$	Heat capacity	1250 J kg-1 K-1	(1)
${\rho }_m$	Melt density	2800 kg m-3	(1)
$L_m$	Latent heat	600 $\times {10}^3$ J kg-1	(1)
$R_p$	Planet radius	6378.1 km	below (1)
$R_c$	Core radius	3488.1 km	below (1)
$k$	Thermal conductivity	5 W m-1 K-1	(2)
$c_1$	Scaling law constant for $\delta$	0.5	(3)
$d$	Whole mantle thickness	2890 km	(3)
$T_s$	Surface temperature	273 K	below (3)
$a_{rh}$	Constant for rheological temperature scale	2.5	(4)
$E_v$	Activation energy for viscosity	300 kJ mol-1	(4)
$R$	Universal gas constant	8.314 J K-1 mol-1	(4)
$g$	Gravity	9.8 m s-2	below (4)
$\kappa$	Thermal diffusivity	${10}^{-6}$ m2 s-1	below (4)
$\alpha$	Thermal expansivity	$3\times {10}^{-5}$ K-1	below (4)
${\mu }_n$	Viscosity pre-exponential factor	$4\times {10}^{10}$ Pa s	below (4)
${\mu }_r$	Reference viscosity	$2\times {10}^{20}$ Pa s	below (4)
${\gamma }_{mantle}$	Mantle adiabatic temperature gradient	${10}^{-8}$ K Pa-1	(8)
${\rho }_l$	Lithosphere density	3300 kg m-3	(9)
${\left(d\varphi /dP\right)}_S$	Pressure derivative of melt fraction	$1.5\times {10}^{-10}$ Pa-1	(10)
$\gamma$	Effective adiabatic temperature gradient of melt	$6.766\times {10}^{-4}$ K km-1	below (11)
$c_2$	Scaling law constant for $v$	0.05	(12)
$D$	Distribution coefficient	0.002	(15)
${\tau }_{rad}$	Radioactive decay constant	2.94 Gyrs	(15)
$A$	Decarbonation temperature constant	$3.125\times {10}^{-3}$ K m-1	(17)
$B$	Decarbonation temperature constant	$835.5$ K	(17)
$D_{C{O}_2}$	Distribution coefficient for CO2	${10}^{-4}$	(19)
$E_{fsw}$	Seafloor weathering activation energy	$75000$ J mol-1	(21)
$\beta$	Sensitivity of seafloor weathering to eruption rate	$1/2$	(21)
$F_{fsw}^{\ast }$	Modern day Earth seafloor weathering rate	$0.5$ Tmol yr-1	(21)
$T_{bot}^{\ast }$	Modern day Earth ocean bottom temperature	$274$ K	(21)
${\eta }^{\ast }$	Modern day Earth eruption rate	$23\times {10}^9$ km3 yr-1	(21)
$ϵ$	Eruption efficiency	$0.1$	(21)
$\chi$	Weathering demand of basalt	5.8 mol kg-1	(22)
$T_s^{\ast }$	Present day Earth surface temperature	285 K	(24)
$p\mathrm{CO}_2^{}{}_{}{}^{\ast }$	Partial pressure of atmospheric CO2	30 Pa	(24)
${\overline{m}}_{C{O}_2}$	Molar mass of CO2	0.044 kg mol-1	below (24)
$A_{surf}$	Surface area of the planet	$5.11\times {10}^{14}$ m2	below (24)
$h_{\varphi }$	Atmosphere-ocean CO2 partitioning parameter	${10}^{20}$ mol	(25)
$T_i$	Initial mantle potential temperature	2000 K	§3.1

Table 2. Table 2: Model Variables

Variable	Meaning (units)	Equation
$T_p$	Mantle potential temperature (K)	(1)
$t$	Time (s)	(1)
$V_{man}$	Volume of convecting mantle (m3)	(1)
$A_{man}$	Surface are of top of convecting mantle (m2)	(1)
$Q_{man}$	Radiogenic heat production in the mantle (W)	(1) and (16)
$F_{man}$	Mantle convective heat flux (W m-2)	(1)
$f_m$	Volumetric melt production rate (m3 s-1)	(1)
$\Delta T_m$	Temperature difference between erupted melt and $T_s$ (K)	(1) and below (11)
$\delta$	Thickness of stagnant lid (m)	(2) and (9)
$z$	Height above planet’s center (km)	(2)
$T_l$	Temperature at base of stagnant lid (K)	(2) and (4)
$T_{bot}$	Temperature at the bottom of the ocean (K)	below (3) and (23)
$\theta$	Frank-Kamenetskii parameter (unit-less)	(3)
$R{a}_i$	Internal Rayleigh number (unit-less)	(3)
${\mu }_i$	Interior mantle viscosity (Pa s)	below (4)
$x_c$	Concentration of heat producing elements in the crust (W m-3)	below (4) and (6)
$Q_{crust}$	Radiogenic heat production in the crust (W)	below (4) and (15)
$V_{crust}$	Volume of crust (m3)	below (4) and (13)
$x_m$	Concentration of heat producing elements in the crust (W m-3)	below (4) and (5)
$V_{lid}$	Volume of subcrustal stagnant lid (m3)	below (4)
${\delta }_c$	Thickness of crust (m)	(7) and (14)
$T_c$	Temperature at the base of the crust (K)	(7)
$P_i$	Pressure where melting begins (Pa)	(8)
$P_f$	Pressure where melting ceases (Pa)	(9)
$\varphi$	Melt fraction (unit-less)	(10)
$d_m$	Depth where melting begins (m)	(11)
$v$	Convective velocity (m s-1)	(12)
$T_{decarb}$	Temperature of metamorphic decarbonation (K)	(17)
${\delta }_{carb}$	Decarbonation depth (m)	(18)
$R_{man}$	Mantle carbon reservoir (mol)	(19) and (27)
$F_d$	Mantle degassing flux (mol s-1)	(19)
$R_{crust}$	Crustal carbon reservoir (mol)	(20) and (28)
$F_{meta}$	Metamorphic degassing flux (mol s-1)	(20)
$V_{carb}$	Volume of carbonated crust (m3)	(20)
$F_{weather}$	Seafloor weathering flux (mol s-1)	(21)
$F_{sl}$	Supply limited weathering flux (mol s-1)	(22)
$T_s$	Surface temperature (K)	(24)
$p{\mathrm{CO}}_2$	Partial pressure of atmospheric CO2 (Pa)	(24)
${\varphi }_{aoc}$	Fraction of atmosphere-ocean CO2 in atmosphere (unit-less)	(25)
$R_{aoc}$	Atmosphere-ocean carbon reservoir (mol)	(26)
$C_{tot}$	Total carbon budget (mol)	§3.1
$Q_0$	Initial radiogenic heat production (W)	§3.1

Equations

$$V_{man}\rho c_p\frac{d{T}_p}{dt}=Q_{man}-A_{man}{F}_{man}-f_m{\rho }_m\left(c_p\Delta T_m+L_m\right),$$

$$\rho c_{p}(T_{p}-T_{l})\frac{d\delta}{dt}=-F_{man}-k\frac{\partial T}{\partial z}\Bigr{|}_{z=R_{p}-\delta},$$

$$F_{man}=\frac{c_{1}k(T_p-T_{bot})}{d}{\theta }^{-4/3}R{a}_i^{1/3}$$

$$T_l=T_p-\frac{a_{rh}R{T}_p^2}{E_v}.$$

$${-k\frac{\partial T}{\partial z}\Bigr{|}_{z=R_{p}-\delta}=\frac{k(T_{l}-T_{c})}{\delta-\delta_{c}}-\frac{x_{m}(\delta-\delta_{c})}{2}}$$

$$-k\frac{\partial T}{\partial z}\Bigr{|}_{z=R_{p}-\delta}=-\frac{x_{c}\delta_{c}}{2}+\frac{k(T_{l}-T_{bot})}{\delta_{c}}$$

$$T_c=\frac{T_{bot}(\delta -{\delta }_c)+T_l{\delta }_c}{\delta }+\frac{x_c{\delta }_c^2(\delta -{\delta }_c)+x_m{\delta }_c(\delta -{\delta }_c{)}^2}{2k\delta }.$$

$$P_i=\frac{T_p-1423}{120\times 10^{-9}-{\gamma }_{mantle}},$$

$$P_f={\rho }_{l}g\delta ,$$

$$\varphi =\frac{P_i-P_f}{2}{\left(\frac{d\varphi }{dP}\right)}_S,$$

$$f_m=17.8\pi R_{p}v(d_m-\delta )\varphi ,$$

$$v=c_2\frac{\kappa }{d}{\left(\frac{R{a}_i}{\theta }\right)}^{2/3},$$

$$\frac{d{V}_{crust}}{dt}=f_m-\left(f_m-4\pi (R_p-\delta {)}^2\text{min}\left(0,\frac{d\delta }{dt}\right)\right)(\tanh (({\delta }_c-\delta )20)+1).$$

$${\delta }_c=R_p-{\left(R_p^3-\frac{3V_{crust}}{4\pi }\right)}^{1/3},$$

$$\begin{array}{rl}\frac{d{Q}_{crust}}{dt}& =\frac{x_m{f}_m}{\varphi }[1-(1-\varphi {)}^{1/D}]-\\ & x_c\left(f_m-4\pi (R_p-\delta {)}^2\text{min}\left(0,\frac{d\delta }{dt}\right)\right)(\tanh (({\delta }_c-\delta )20)+1)-\frac{Q_{crust}}{{\tau }_{rad}}\end{array}$$

$$\begin{array}{rl}\frac{d{Q}_{man}}{dt}& =x_c\left(f_m-4\pi (R_p-\delta {)}^2\text{min}\left(0,\frac{d\delta }{dt}\right)\right)(\tanh (({\delta }_c-\delta )20)+1)-\\ & \frac{x_m{f}_m}{\varphi }[1-(1-\varphi {)}^{1/D}]-\frac{Q_{man}}{{\tau }_{rad}}.\end{array}$$

$$T_{decarb}=A(R_p-z)+B$$

$$\begin{array}{rl}{\delta }_{carb}=& \frac{{\delta }_c}{2}+\frac{k(T_c-T_{bot})}{{\delta }_c{x}_c}-\frac{Ak}{x_c}-\\ & \frac{k}{x_c}\sqrt{{\left(\frac{x_c{\delta }_c}{2k}+\frac{T_c-T_{bot}}{{\delta }_c}-A\right)}^2+\frac{2x_c}{k}(T_{bot}-B)}.\end{array}$$

$$F_d=\frac{f_m{R}_{man}[1-(1-\varphi {)}^{1/D_{C{O}_2}}]}{\varphi (V_{man}+V_{lid})},$$

$$F_{meta}=\frac{R_{crust}{f}_m}{2V_{carb}}(\tanh (({\delta }_c-{\delta }_{carb})20)+1)$$

$$F_{weather}=F_{sl}\left\{1-\exp \left[-\frac{F_{sfw}^{\ast }}{F_{sl}}\exp \left(-\left(\frac{E_{sfw}}{R}\right)\left(\frac{1}{T_{bot}}-\frac{1}{T_{bot}^{\ast }}\right)\right){\left(\frac{ϵf_m}{{\eta }^{\ast }}\right)}^{\beta }\right]\right\}$$

$$F_{sl}=ϵf_m\chi {\rho }_l,$$

$$T_{bot}=\max (270,1.02T_s-16.7);$$

$$T_s=T_s^{\ast }+5.6\left(\frac{\ln (p{\mathrm{CO}}_2/p{\mathrm{CO}}_2^{\ast })}{\ln (2)}\right),$$

$${\varphi }_{aoc}=\frac{0.78R_{aoc}}{R_{aoc}+h_{\varphi }}$$

$$\frac{d{R}_{aoc}}{dt}=F_d+F_{meta}-F_{weather},$$

$$\frac{d{R}_{man}}{dt}=-F_d,$$

$$\frac{d{R}_{crust}}{dt}=F_{weather}-F_{meta}.$$

$$\frac{d{R}_{aoc}}{dt}=F_d-F_{weather}$$

$$\begin{array}{rl}\frac{d{R}_{man}}{dt}& =\frac{R_{crust}}{V_{crust}}\left(f_m-4\pi (R_p-\delta {)}^2\text{min}\left(0,\frac{d\delta }{dt}\right)\right)(\tanh (({\delta }_c-\delta )20)+1)-F_d.\end{array}$$

Full text

astrobiology — planets and satellites: physical evolution — planets and satellites: terrestrial planets

1 Introduction

Determining the factors that allow long-lived, habitable climates to develop on rocky planets is a critical goal for astrobiology, especially in light of the large number of exoplanets that have been discovered in the last 20 years (e.g. Wright et al., 2011; Batalha, 2014; Burke et al., 2014). In particular, such knowledge can be used to guide future observing resources to those planets most likely to harbor life. In order to host surface life which can be remotely observed, it is thought that a planet must lie within the habitable zone, such that liquid water can exist on the planet’s surface (e.g. Kasting et al., 1993). However, even within the habitable zone, stabilization of climate via the carbonate-silicate cycle is necessary to prevent changes in solar luminosity, or variations in CO2 outgassing rate, from driving climate towards uninhabitably hot or cold states (e.g. Walker et al., 1981; Berner et al., 1983; Berner & Caldeira, 1997; Sleep & Zahnle, 2001; Abbot et al., 2012; Foley & Driscoll, 2016). Plate tectonics plays a critical role in facilitating the carbonate-silicate cycle on Earth, and has thus often been posited to be a necessary component for sustaining a long-lived habitable climate on any rocky planet (e.g. Kasting & Catling, 2003).

However, whether plate tectonics is actually required for the carbonate-silicate cycle to operate and stabilize a planet’s climate is unclear (e.g. Lenardic et al., 2016). Stagnant lid planets still experience volcanism, which can release a significant amount of CO2 to the atmosphere, as estimated in previous studies for Mars (e.g. Pollack et al., 1987; O’Neill et al., 2007; Grott et al., 2011), a hypothetical stagnant lid Earth (Tosi et al., 2017), and super-Earths (Noack et al., 2017; Dorn et al., 2018). Moreover, burial of carbonated crust under lava flows could lead to metamorphic decarbonation of this crust, providing an additional CO2 source to the atmosphere, or recycling of surface CO2 back into the mantle (e.g. Pollack et al., 1987; Elkins-Tanton et al., 2007). Foley & Smye (2018) modeled the evolution of CO2 outgassing rates for Earth-sized stagnant lid planets, considering both volcanic and metamorphic outgassing, and found that habitable climates could last for 2-5 Gyrs, depending on the amount of radiogenic elements the planet acquires during its formation, and the amount of CO2 in the mantle and surface reservoirs. Specifically there is an optimal range of CO2 budgets that allow for long-lived habitable climates. Too low a CO2 budget and even with significant volcanism, outgassing rates are too low to maintain a warm climate. In this case a snowball state would prevail, possibly with limit cycles where long-lived snowball states are punctuated by short-lived hothouse climates (Kadoya & Tajika, 2014; Menou, 2015; Haqq-Misra et al., 2016; Batalha et al., 2016). Too high a CO2 budget and outgassing rates overwhelm weathering, causing an uninhabitably hot climate to form. The higher the abundance of heat producing elements the planet acquires, the longer volcanism lasts, and, when CO2 budgets are within the optimum range for habitability, the longer habitable climates last as well. Valencia et al. (2018) obtained similar results considering planets that are tidally heated, in which case volcanism and outgassing could last as long as the planet’s orbit leads to significant tides.

While an important step in assessing the habitability of stagnant lid planets, Foley & Smye (2018) only tracks CO2 outgassing rates, and assumes that outgassing rates must be at least 10-100 % of Earth’s present day degassing rate in order to prevent a cold, snowball climate, and must be below the global weathering supply limit in order to prevent an extreme hothouse. Actual climate evolution is not tracked, and doing so might change the estimated range of CO2 budgets or radiogenic heating budgets that lead to long-lived habitability. Another issue not addressed in Foley & Smye (2018) is whether a stagnant lid planet could ever recover from a snowball climate, should one develop due to some perturbation in outgassing rate or weathering rate. The thermal evolution models used in this study, and in Foley & Smye (2018), track the evolution of global average outgassing rates over billion year timescales, but can not capture short timescale fluctuations in outgassing (on the order of 10s of millions of years) that could result from the time dependent nature of mantle convection. As such, perturbations that could trigger a snowball climate, even when long term average outgassing rates should be high enough to keep the climate warm, are possible. In this study a model for weathering and climate evolution is included with the thermal evolution and CO2 outgassing model of Foley & Smye (2018), to determine how climate actually evolves on stagnant lid planets, re-assess the conditions that lead to long-lived habitable climates, and assess whether snowball climates can be recovered from.

The models presented here make some important assumptions. First, only planets with a similar size and composition to Earth are considered, as the key material parameters governing mantle convection are highly uncertain for large super-Earth planets, where pressures in the mantle are far higher than experiments can reach (e.g. Karato, 2011; Stamenkovic et al., 2011), and for planets with significantly different compositions than Earth. I also assume all melt produced is buoyant and will rise to the surface, as the melt generation depth is shallower than typical estimates for the silicate melt density crossover, as shown in Foley & Smye (2018). Stagnant lid super-Earth planets, larger than 3 or 4 Earth masses, may actually experience limited outgassing, due to silicate melts being generated at such high pressures that the melts are denser than the surrounding mantle (Noack et al., 2017; Dorn et al., 2018). Another important assumption is that a constant solar luminosity, equal to the solar constant for Earth, is assumed. This is primarily done for simplicity, as the focus of this study is on how the thermal evolution of stagnant lid planets’ mantles controls their climate evolution and climate stability. Adding stellar evolution would include an additional independent variable in the model, overcomplicating the results and their interpretation.

2 Theory

2.1 Thermal evolution model

The thermal evolution model follows after Foley & Smye (2018), which gives a thorough explanation of the model. As a result, only a brief overview of the model equations will be provided here. Pure internal heating is assumed, as this greatly simplifies the modeling and provides the most conservative estimate for the lifetime of outgassing and habitable climates. The thermal evolution models also include heat loss due to melting; all melt generated is assumed to contribute to mantle cooling and the growth of the crust, while only a fraction of the melt is assumed to erupt at the surface. These assumptions result in the following equation for the evolution of mantle temperature (e.g. Stevenson et al., 1983; Hauck & Phillips, 2002; Reese et al., 2007; Fraeman & Korenaga, 2010; Morschhauser et al., 2011; Driscoll & Bercovici, 2014; Foley & Smye, 2018):

[TABLE]

where $V_{man}$ is the volume of convecting mantle beneath the stagnant lid, $\rho$ is the bulk density of the mantle, $c_{p}$ is the heat capacity, $T_{p}$ is the potential temperature of the upper mantle, $t$ is time, $Q_{man}$ is the radiogenic heating rate in the mantle, $A_{man}$ is the surface area of the top of the convecting mantle, $F_{man}$ is the heat flux from the convecting mantle into the base of the stagnant lid, $f_{m}$ is the volumetric melt production rate, $\rho_{m}$ is the density of mantle melt, $L_{m}$ is the latent heat of fusion of the mantle, and $\Delta T_{m}$ is the temperature difference between the melt erupted at the surface and the surface temperature. The volume of the convecting mantle is $V_{man}=(4/3)\pi((R_{p}-\delta)^{3}-R_{c}^{3})$, where $R_{p}$ is the planetary radius, $R_{c}$ is the core radius, and $\delta$ is the thickness of the stagnant lid. The surface area of the convecting mantle is then $A_{man}=4\pi(R_{p}-\delta)^{2}$. As Earth-sized planets are considered in this study, $R_{p}$ is set equal to Earth’s radius and $R_{c}$ equal to Earth’s core radius.

The thickness of the stagnant lid, $\delta$, evolves with time as (e.g. Schubert et al., 1979; Spohn, 1991),

[TABLE]

where $T_{l}$ is the temperature at the base of the lid, $k$ is the thermal conductivity, and $z$ is height above the planet’s center. As volcanism is assumed to cause heat loss from the mantle directly to the surface, melt heat loss does not appear in (2). The mantle heat flux, $F_{man}$, and lid base temperature, $T_{l}$, are calculated from stagnant lid convection scaling laws as (Reese et al., 1998, 1999; Solomatov & Moresi, 2000; Korenaga, 2009):

[TABLE]

and

[TABLE]

In (3) & (4), $c_{1}$ and $a_{rh}$ are constants, with assumed values $c_{1}=0.5$ and $a_{rh}=2.5$. $T_{bot}$ is the temperature at the bottom of the oceans, which are assumed here to cover most of the planet’s surface and thus set the top boundary condition for the convecting mantle (see §2.3), $d=R_{p}-R_{c}$ is the whole mantle thickness, $\theta$ is the Frank-Kamenetskii parameter, defined as $\theta=E_{v}(T_{p}-T_{bot})/(RT_{p}^{2})$, where $E_{v}$ is the activation energy for mantle viscosity and $R$ is the gas constant, and $Ra_{i}$ is the internal Rayleigh number. The internal Rayleigh number is defined as $Ra_{i}=\rho g\alpha(T_{p}-T_{bot})d^{3}/(\kappa\mu_{i})$, where $g$ is gravitational acceleration, $\alpha$ is thermal expansivity, $\kappa$ is thermal diffusivity, and $\mu_{i}$ is the interior mantle viscosity, defined at temperature $T_{p}$. Mantle viscosity follows the temperature-dependent viscosity law $\mu_{i}=\mu_{n}\exp{(E_{v}/(RT_{p}))}$, where $\mu_{n}$ is a constant; $E_{v}=300$ kJ mol*-1* is the baseline activation energy used in this study (Karato & Wu, 1993), and test cases using $E_{v}=200$ kJ mol*-1* and $E_{v}=400$ kJ mol*-1* are shown in §4.4. A baseline reference viscosity of $\mu_{r}=2\times 10^{20}$ Pa s is assumed. The reference viscosity is defined as $\mu_{r}=\mu_{n}\exp{(E_{v}/(RT_{r}))}$, where $T_{r}=1623$ K is Earth’s present day mantle temperature; this then defines the constant $\mu_{n}$ as $\mu_{n}=\mu_{r}\exp{(-E_{v}/(RT_{r}))}$. With $E_{v}=300$ kJ mol*-1* $\mu_{n}=4\times 10^{10}$ Pa s (see Tables 1 & 2 for a complete list of parameters and variables in the model).

The temperature gradient at the base of the stagnant lid, as needed in (2), is calculated assuming one-dimensional heat conduction in steady-state, with constant radiogenic heating rates in the crust and mantle. The crustal radiogenic heating rate is defined as $x_{c}=Q_{crust}/V_{crust}$, where $x_{c}$ is the crustal heat production rate per unit volume, $Q_{crust}$ is the total radiogenic heating rate in the crust, and $V_{crust}$ is the volume of the crust. The mantle radiogenic heating rate $x_{m}=Q_{man}/(V_{man}+V_{lid})$, where $x_{m}$ is the mantle heat production rate per unit volume, $V_{lid}=(4/3)\pi((R_{p}-\delta_{c})^{3}-(R_{p}-\delta)^{3})$ is the volume of the sub-crustal stagnant lid, and $\delta_{c}$ is the thickness of the crust. The temperature profile, and hence lid base temperature gradient, can be solved for analytically in this case, giving

[TABLE]

when $\delta>\delta_{c}$, and

[TABLE]

when $\delta=\delta_{c}$. The temperature at the base of the crust, $T_{c}$, is given by

[TABLE]

The thermal conductivity, $k$, is assumed constant; the difference between crustal and mantle thermal conductivity is ignored for simplicity. As explained in Foley & Smye (2018), (6) is used when $T_{c}$ is within one degree Kelvin of $T_{l}$, which is equivalent to switching when the difference between $\delta_{c}$ and $\delta$ is negligibly small. The crustal temperature profile is also calculated assuming that advection is negligible. Foley & Smye (2018) showed that advection is small compared to diffusion, though not necessarily negligible. As a result of ignoring advection, the lid thickness may be slightly underestimated in our models, and hence conductive cooling overestimated. Advection may also suppress metamorphic decarbonation of the crust, so a case where decarbonation is assumed to not occur is presented in §4.4.

2.2 Melting and crustal evolution

The melt production rate, $f_{m}$, is calculated as in Foley & Smye (2018), which is based on Fraeman & Korenaga (2010). The pressure where melting begins, $P_{i}$, is

[TABLE]

where $\gamma_{mantle}$ is the average adiabatic temperature gradient in the upper mantle ($\approx 10^{-8}$ K Pa*-1*). Assuming melting stops at the base of the stagnant lid, the pressure where melting stops is

[TABLE]

where $\rho_{l}$ is the average density of the crust and lithosphere (assumed to be $\rho_{l}=3300$ kg m*-3*). The melt fraction, $\phi$, is

[TABLE]

where $(d\phi/dP)_{S}\approx 1.5\times 10^{-10}$ Pa*-1*.

The full melt production rate is then given by assuming a cylindrical region of upwelling mantle enters the melting zone beneath the stagnant lid (as in, e.g. Reese et al., 1998). A derivation is given in Foley & Smye (2018), which results in

[TABLE]

where $v$ is the convective velocity, $d_{m}$ is the depth where melting begins ($d_{m}=P_{i}/(\rho_{l}g$)), and assuming that $d\approx 0.45R_{p}$ as on Earth. The exact form of the melt production rate was found to not significantly impact the results in Foley & Smye (2018). The temperature difference between melt erupted at the surface and the surface temperature is $\Delta T_{m}=T_{p}-P_{i}\gamma-T_{s}$ where $T_{s}$ is the (atmospheric) surface temperature of the planet and $\gamma$ was found to be $\gamma=6.766\times 10^{-4}$ K km*-1* based on values in Driscoll & Bercovici (2014). The velocity, $v$, is calculated from scaling laws for internally heated stagnant lid convection (Reese et al., 1998, 1999; Solomatov & Moresi, 2000; Korenaga, 2009) as

[TABLE]

where $c_{2}$ is a constant.

All melt generated is assumed to contribute to crust growth, and any crust buried to depths below $\delta$ is assumed to founder into the mantle, such that the lithospheric thickness always follows (2). With the above assumptions, the evolution of crustal volume, $V_{crust}$, is

[TABLE]

The hyperbolic tangent function allows the crustal loss rate to go to zero when $\delta_{c}<\delta$, and equal the crustal growth rate plus the lid thinning rate when $\delta_{c}=\delta$. The term $4\pi(R_{p}-\delta)^{2}\textrm{min}(0,d\delta/dt)$ gives the rate at which the volume of stagnant lid is lost when $\delta$ is shrinking, and is 0 when $\delta$ is growing. The crustal thickness, $\delta_{c}$, is calculated as

[TABLE]

assuming that crustal thickness is constant spatially across the planet.

Heat producing elements are incompatible and partition preferentially into the melt. Accumulated fractional melting is assumed in this study (e.g. Fraeman & Korenaga, 2010; Morschhauser et al., 2011), resulting in the following equations for the evolution of the crustal heat production rate, $Q_{crust}$,

[TABLE]

and mantle heat production rate, $Q_{man}$,

[TABLE]

The distribution coefficient, $D$, is assumed to have a value $D=0.002$ (Hart & Dunn, 1993; Hauri et al., 1994). A constant decay constant for both the crust and mantle, $\tau_{rad}\approx 2.94$ Gyrs, is used based on an average of the four major radiogenic isotopes powering Earth’s mantle, 238U, 235U, 232Th, and 40K (Driscoll & Bercovici, 2014).

2.3 Carbon cycle and climate

Cycling of CO2 is tracked in a similar fashion to Foley & Smye (2018), though with some important differences: in this study CO2 in the atmosphere and oceans is explicitly tracked, surface temperature is calculated from a simple climate model, and weathering rates are calculated to link outgassing, atmospheric CO2, and climate. Carbon is assumed to outgas to the atmosphere-ocean system via both mantle volcanism and metamorphic decarbonation of buried crust. Once in the atmosphere-ocean system weathering removes CO2 to the crust, where subsequent volcanism buries the crust until it either experiences metamorphic decarbonation, or the carbonated crust founders back into the mantle. The atmosphere-ocean system is assumed to always be in equilibrium, and CO2 is partitioned between the atmosphere and ocean accordingly (see (25)).

The temperature-pressure conditions where carbonate-bearing metabasalt experiences decarbonation are determined from the P-T phase diagram shown in Foley & Smye (2018) (see also Staudigel et al., 1989; Kerrick & Connolly, 2001). With a simple linear fit to the 0.25 wt. % CO2 contour, the temperature at which decarbonation occurs is

[TABLE]

where $A=3.125\times 10^{-3}$ K m*-1*, $B=835.5$ K, $T_{decarb}$ is temperature in Kelvins, and pressure is converted into depth by assuming a lithospheric density of $\rho_{l}=3300$ kg m*-3*. The depth where decarbonation occurs, $\delta_{carb}$, is then given by (see Foley & Smye, 2018, for details)

[TABLE]

There are three carbon fluxes in the model: the mantle degassing flux due to volcanism, $F_{d}$, the metamorphic degassing flux due to decarbonation of the crust, $F_{meta}$, and the weathering flux, $F_{weather}$. The mantle degassing flux is given by the product of the concentration of CO2 in mantle melt and the mantle melt production rate,

[TABLE]

where $R_{man}$ is the mantle CO2 reservoir, in mols of CO2. The formulation for $F_{d}$ takes into account the fact that CO2 is an incompatible element, like heat producing elements, with a distribution coefficient of $D_{CO_{2}}=10^{-4}$ (e.g. Hauri et al., 2006). The metamorphic degassing flux is given by the product of the concentration of CO2 in the crust and the volumetric melt production rate, as

[TABLE]

where $R_{crust}$ is the crustal CO2 reservoir in mols, $V_{carb}$ is the volume of carbonated crust, $V_{carb}=(4/3)\pi(R_{p}^{3}-(R_{p}-\delta_{carb})^{3})$, and, as in (13), a hyperbolic tangent function is used as a mathematically convenient way to parameterize a step-function like change in crustal CO2 around the decarbonation depth. Here a factor of 1/2 is included, so that the metamorphic degassing flux is $R_{crust}f_{m}/(V_{carb})$ when $\delta_{c}>\delta_{carb}$ and zero when $\delta_{c}<\delta_{carb}$. Note that in this formulation metamorphic outgassing does not begin until the crust first grows deep enough to reach the decarbonation depth. Metamorphic CO2 outgassing is assumed to be complete and to occur uniformly across the planet’s surface, as crustal production is assumed to be spatially uniform. Incomplete metamorphic outgassing was found to not significantly impact the results in Foley & Smye (2018). The above formulation for the metamorphic outgassing flux also assumes that the crust acts as a single reservoir with instantaneous mixing, when in reality mixing between different crustal layers will not occur. However, assuming a mixed crustal reservoir does not significantly impact the results, as shown in Appendix A, and greatly simplifies the model.

The rate of CO2 drawdown from the atmosphere-ocean system is calculated by assuming that seafloor weathering is the dominant weathering process. Seafloor weathering is assumed because stagnant lid planets are likely to have crusts that are predominantly basaltic (Breuer & Moore, 2007), as there is no subduction to produce continental crust as on Earth. Furthermore, without continents, planets will not have well defined ocean basins, so oceans may be more wide spread across the planets’ surface. Of course no actual constraints on topography and the size of oceans for real stagnant lid exoplanets are available, so the assumption here of seafloor weathering being the primary weathering process may not be correct. However, it is a simple first order assumption to make based on our knowledge of stagnant lid versus plate tectonic planets. The weathering flux includes a supply limit, as in Foley (2015) & Foley & Smye (2018), and is temperature-dependent as in Krissansen-Totton & Catling (2017). The plate speed dependence of seafloor weathering from Krissansen-Totton & Catling (2017) is replaced with a dependence on the melt eruption flux, as it is ultimately the crustal production rate that matters for seafloor weathering. However the dependence of seafloor weathering on ocean pH is ignored, as it is significantly weaker than the temperature dependence.

The CO2 drawdown rate is thus

[TABLE]

where $F_{sl}$ is the supply limit to seafloor weathering, $F_{sfw}^{*}$ is the modern day Earth seafloor weathering rate and is a constant, $E_{sfw}$ is the activation energy for seafloor weathering, $T_{bot}^{*}$ is the modern day Earth ocean bottom temperature, $\epsilon$ is the eruption efficiency, or the fraction of melt produced that erupts at the surface (and thus $\epsilon f_{m}$ is the eruption rate), and $\eta^{*}$ is the eruption rate at modern Earth mid-ocean ridges. To prevent numerical convergence problems associated with very small atmosphere-ocean carbon reservoir sizes, $F_{weather}$ is set to 0 when surface temperatures fall below 260 K. In reality weathering could continue in a snowball state, in particular if exchange between atmosphere and sub-ice ocean occurs (e.g. Tajika, 2008), so climate could cool even beyond the 260 K cutoff imposed here. However, surface temperatures typically only fall to 260 K as volcanism and outgassing are waning, at which point habitable climates are assumed to end anyway in this study (see §4.1). The baseline values of $F_{sfw}^{*}$ and $E_{sfw}$ are taken from (Krissansen-Totton & Catling, 2017), the baseline value for $\epsilon$ is $\epsilon=0.1$ (Foley & Smye, 2018), and for $\eta^{*}$ is $\eta^{*}=23\times 10^{9}$ km3 yr*-1* (Crisp, 1984) (see Table 1).

The supply limit to weathering is (Foley, 2015; Foley & Smye, 2018)

[TABLE]

where $\chi$ is the weathering demand for complete carbonation of basalt, analogous to the formulation in Kump (2018) for continental crust. The weathering demand is calculated from the mass fraction of reactable oxides in basalt: CaO, MgO, FeO, K2O, and Na2O. Mass fractions from Gale et al. (2013) are used, which gives 5.8 moles of oxide per kg of basalt. Hence $\chi=5.8$ mol kg*-1* of net CO2 drawdown by completely carbonating basalt. The temperature at the base of the ocean is related to the surface (atmospheric) temperature of the planet as (Krissansen-Totton & Catling, 2017)

[TABLE]

the temperature is not allowed to fall below 270 K at the bottom of the ocean because even during a cold climate the bottom of the ocean is not expected to freeze. $T_{bot}^{*}$ is $T_{bot}^{*}=1.02T_{s}^{*}-16.7$, where $T_{s}^{*}=285$ K is Earth’s present day surface temperature. The surface temperature is calculated from the simple parameterization given in Krissansen-Totton & Catling (2017),

[TABLE]

where $p\mathrm{CO}_{2}$ is the partial pressure of atmospheric CO2, and $p\mathrm{CO}_{2}^{*}$ is the pre-industrial $p\mathrm{CO}_{2}$ for the Earth ($p\mathrm{CO}_{2}^{*}=30$ Pa). The partial pressure of atmospheric CO2 is defined as $p\mathrm{CO}_{2}=R_{atm}\bar{m}_{CO_{2}}g/A_{surf}$, where $R_{atm}$ is the size of the atmospheric CO2 reservoir in mol, $\bar{m}_{CO_{2}}$ is the molar mass of CO2, and $A_{surf}$ is the surface area of the planet. Finally, $R_{atm}$ is related to the ocean plus atmosphere CO2 reservoir, $R_{aoc}$, as $R_{atm}=\phi_{aoc}R_{aoc}$, where

[TABLE]

as in Mills et al. (2011); here $h_{\phi}=10^{20}$ mol. Note that the climate parameterization used in this study can not capture a snowball, as this would significantly change the parameterization for surface temperature as a function of atmospheric CO2. However, the ability of planets to recover from snowball states is assessed in this study based on outgassing rates (see §3.2).

When $\delta_{carb}<\delta$, crust decarbonates before being recycled into the mantle. The evolution of the combined atmosphere-ocean CO2 reservoir is then given by (e.g. Tajika & Matsui, 1992; Franck et al., 1999; Sleep & Zahnle, 2001; Driscoll & Bercovici, 2013; Foley, 2015)

[TABLE]

the evolution of the mantle carbon reservoir, $R_{man}$, by

[TABLE]

and the evolution of the crustal carbon reservoir, $R_{crust}$, by

[TABLE]

As CO2 is not returned to the mantle when crustal rocks decarbonate before reaching the base of the stagnant lid, volcanic degassing of the mantle causes a net transfer of carbon from the mantle to the crust over the lifetime of a stagnant lid planet in this case.

When $\delta_{c}<\delta_{carb}$, the following equations for $R_{aoc}$, $R_{man}$, and $R_{crust}$ are used:

[TABLE]

[TABLE]

[TABLE]

In this case, crust buried beneath the stagnant lid base will founder into the mantle, returning surface carbon to the mantle. However, in practice $\delta_{carb}$ is usually found to be less than $\delta_{c}$ while volcanism is active, so models typically experience crustal decarbonation and do not recycle CO2 from the crust back into the mantle.

3 Model setup

3.1 Initial conditions

The above equations for the evolution of mantle temperature, crust volume, CO2 cycling, and surface temperature are solved simultaneously as a set of coupled ordinary differential equations. In order to solve the system initial conditions and some additional constraints must be specified. First, two key input quantities that are varied in the models are the initial radiogenic heating rate, $Q_{0}$, and the total CO2 budget of the surface and mantle reservoirs, $C_{tot}$; as CO2 among these reservoirs is conserved, $C_{tot}=R_{crust}+R_{man}+R_{aoc}$. All models start with a chosen initial heat production rate, $Q_{0}$, that is entirely within the mantle. $Q_{0}$ is varied over a range of 5-250 TW, encompassing typical estimates for the Earth’s initial mantle radiogenic heating rate, which span $\approx 70-100$ TW when present day heating rates in the mantle and continental crust are combined and extrapolated back in time (e.g. Korenaga, 2006). Models also start with an initial crust volume of zero, that is no crust is present; as volcanism causes the crust to grow, heat producing elements are then partitioned into the crust. Models are also run for a wide range of $C_{tot}=3\times 10^{18}-3\times 10^{24}$ mol, bracketing the Earth’s estimated value of $\sim 10^{22}$ mol (Sleep & Zahnle, 2001).

In the case of CO2, initial conditions are important as the distribution of CO2 between the different reservoirs determines the planet’s climate. Two end member limits are considered here: an initially hot climate (or “hot start” scenario) and an initially cool climate (or “cold start” scenario). In the initially hot climate, 99.99 % of the CO2 in the planet’s surface and mantle reservoirs is placed in the atmosphere-ocean system, with the rest residing in the mantle. This initial condition is a possible result of magma ocean solidification, which is predicted to result in significant outgassing from the mantle to the surface, and hence an initially CO2 rich atmosphere (Zahnle et al., 2007; Elkins-Tanton, 2008). In the cold start scenario, the initial surface temperature is held fixed at 273 K, and CO2 partitioned between reservoirs accordingly. That is, the amount of CO2 required to produce $T_{s}=273$ K is placed in the atmosphere-ocean system, and the remainder is placed in the mantle. Colder initial surface temperatures are not considered, since the climate model in this study can not explicitly capture snowball climates. The initial temperature of the mantle, $T_{i}$, is typically set at $T_{i}=2000$ K, a reasonable value for a recently solidified silicate mantle (e.g. Abe, 1997; Solomatov, 2000; Foley et al., 2014). However, a range of initial mantle temperatures from $T_{i}=1500-2200$ K is tested (see §4.4). Finally, the initial stagnant lid thickness is calculated from (2) in steady-state at the specified initial mantle temperature, surface temperature, and heat production rate.

3.2 Calculating recovery from snowball states

A major goal of this paper is determining whether stagnant lid planets could recover from snowball episodes. To determine this, two snowball scenarios are considered: a “hard” snowball, where there is no exchange between atmosphere and ocean and no weathering during the snowball episode, and a “soft” snowball where atmosphere and ocean are assumed to maintain equilibrium. In the soft snowball, the partitioning fraction of CO2 between atmosphere and ocean, $\phi_{aoc}$, is modified by setting $h_{\phi}=5.3\times 10^{20}$ mol, as in Mills et al. (2011), which results in more CO2 being partitioned into the ocean, and less in the atmosphere, during the snowball phase than during a non-snowball climate. I assume that in order to recover from the snowball phase, the partial pressure of atmospheric CO2 must reach $2.5\times 10^{4}$ Pa of CO2, as in Mills et al. (2011) and consistent with estimates from pervious studies (e.g. Caldeira & Kasting, 1992; Hoffman & Schrag, 2002). For a soft snowball, this constraint on $p\mathrm{CO}_{2}$ means $\approx 1.6\times 10^{20}$ mol of CO2 are needed in the atmosphere-ocean reservoir to melt the ice layer, and for a hard snowball $\approx 3\times 10^{19}$ mol of CO2 are needed in the atmosphere.

For a hard snowball, only volcanism releasing CO2 above the ice layer can contribute to warming the climate. Metamorphic outgassing is likely confined to low lying regions that are beneath the ice layer, so I assume that all CO2 released by crustal decarbonation is released to the oceans, and therefore does not contribute to warming the climate during a hard snowball state. Mantle volcanism, however, is assumed to occur at volcanoes that reach above the ice layer, and thus release of mantle CO2 is assumed to contribute to warming the climate. In reality some of the CO2 released by mantle volcanism may also be confined to the sub-ice oceans, so the assumptions here present an optimistic scenario for recovering from a hard snowball state. In the soft snowball, both metamorphic and volcanic outgassing contribute to warming the climate, as equilibrium between atmosphere and ocean is assumed.

The thermal evolution models produce a time history of outgassing and volcanism; this history is then used to determine if a hypothetical snowball state could be recovered from. For a hard snowball, the outgassing rate due to mantle volcanism, $F_{d}$, is integrated over time as:

[TABLE]

where $t^{\prime}$ is the time when a hypothetical hard snowball state is assumed to occur, and $t_{end}$ is 10 Ga, the time when the thermal evolution models are stopped. Once the integrated value of CO2 released reaches the threshold needed to melt the ice layer, $p\mathrm{CO}_{2}=2.5\times 10^{4}$ Pa, then the snowball is considered to be recovered from. The integration is performed at each timestep of the model output, and the latest time during model evolution at which a hard snowball state can be recovered from is recorded. That is, for any time after this there is not enough CO2 that can be outgassed from the mantle, before volcanism dies out, to melt the ice layer. Thus, the age at which the planet can no longer recover from a hard snowball due to insufficient degassing, is determined.

As with the hard snowball case, in the soft snowball case I determine the planet age after which a soft snowball state can not be recovered from. For each timestep in the model output, corresponding to a time $t^{\prime}$, I calculate whether metamorphic and volcanic outgassing can bring the system out of a hypothetical soft snowball state, by increasing the atmosphere-ocean CO2 reservoir size to at least $1.6\times 10^{20}$ mol. In the case of a soft snowball, carbon cycling between the crust, mantle, and atmosphere needs to be tracked as weathering will still be active, drawing CO2 out of the combined ocean-atmosphere system, and the crustal CO2 reservoir will be changing over time. Thus at each timestep, $t^{\prime}$, equations (26) & (28) are solved forward in time to see if enough CO2 can accumulate in the atmosphere to melt the ice. The degassing flux, $F_{d}(t)$, and the volcanism rate, $f_{m}(t)$, are taken from the thermal evolution model outputs. The initial conditions for the atmosphere-ocean and crustal CO2 reservoirs for the hypothetical snowball state, here labeled $R_{aoc}^{*}(t^{\prime})$ and $R_{crust}^{*}(t^{\prime})$, are determined as follows. I assume that in order to trigger the snowball state, atmospheric CO2 levels must have fallen to 15 Pa (e.g. Donnadieu et al., 2004; Mills et al., 2011). The amount of CO2 in the ocean is then determined by (25) with $h_{\phi}=5.3\times 10^{20}$, hence also setting $R_{aoc}^{*}(t^{\prime})$. The crustal CO2 abundance at the start of the snowball state, $R_{crust}^{*}(t^{\prime})$, is given by $R_{crust}^{*}(t^{\prime})=R_{crust}(t^{\prime})+R_{aoc}(t^{\prime})-R_{aoc}^{*}(t^{\prime})$, where $R_{crust}(t^{\prime})$ and $R_{aoc}(t^{\prime})$ are the amounts of CO2 in the crust and atmosphere-ocean, respectively, at time $t^{\prime}$ of the thermal evolution model output (i.e. before the hypothetical soft snowball climate develops). In other words, I assume that enough CO2 has been removed from the atmosphere-ocean system in order to trigger a snowball state, and that this CO2 has been weathered out into the crust.

Some models are run where metamorphic decarbonation of the crust is assumed to not occur, and all crustal carbon recycles into the mantle. Calculating snowball recovery in this case is slightly different than in the cases where decarbonation is allowed to occur. For a soft snowball state, the same procedure as described in the preceding paragraph is followed, except equations (29)-(31) are integrated, with $f_{m}(t)$ as output from the thermal evolution model results used in the integration. Initial conditions for the crust, mantle, and atmosphere-ocean carbon reservoirs are the same as described above. For a hard snowball, (30) & (31) are integrated forward in time, with the seafloor weathering flux set to zero, so that depletion of the crustal and mantle carbon reservoirs are taken into account. That is, as there is no exchange of CO2 from the atmosphere to the ocean and hence back into the crust during a hard snowball, the crustal CO2 reservoir will be depleted by foundering into the mantle, and the mantle CO2 reservoir will be depleted by volcanic degassing.

4 Results

4.1 Baseline results and influence of initial conditions

Example thermal evolution model results are given in Figures 1 & 2. Specifically two models with identical model parameters, but different initial distributions of CO2 between mantle and atmosphere-ocean, are shown to compare how these different initial conditions influence the overall model evolution. The parameter values are as given in Table 1, which are considered the baseline model parameters for this study. The different initial climate states have no influence on the evolution of mantle temperature or mantle heat loss (Figure 1A & E, D & H), because the difference in surface temperature between the two cases does not significantly influence stagnant lid convection, and climate states for the two cases converge early in planetary evolution (Figure 2). However, the volcanic CO2 outgassing rate and evolution of the mantle CO2 reservoir are strongly affected by the initial condition (Figure 1B & F, C & G). With a hot start, the mantle is depleted in CO2, and CO2 is never recycled back into the mantle over time due to metamorphic decarbonation of the crust. As a result, $R_{man}$ and $F_{d}$ are very low throughout model evolution. With an initially cold climate, mantle outgassing is initially significant, but the mantle is quickly depleted of CO2 and metamorphic outgassing becomes the primary source of atmospheric CO2. The total outgassing rate and amount of CO2 in the atmosphere-ocean reservoirs follow nearly identical evolutionary paths for either initial condition. As a result, the climate evolution is not strongly affected by the initial climate state (Figure 2). For both a hot and cold start, atmospheric CO2 and surface temperature evolutionary paths are nearly identical, after a short initial adjustment phase. During this adjustment phase outgassing leads to a rapid warming for the initially cold climate case, and weathering leads to rapid cooling of the initially hot climate. The initial climate condition has no influence on when mantle volcanism ends, at which point outgassing ceases and a cold climate is expected to develop, as explained below. The insensitivity of climate evolution to the two end member initial conditions used in this study is also confirmed by the more extensive set of results presented below.

The thermal evolution models are used to calculate how long temperate climates, with surface temperatures between 273 K and 400 K, last, for a range of CO2 budgets ($C_{tot}$) and initial radiogenic heating rates ($Q_{0}$). The longevity of habitable climates is only calculated for the time period when volcanism is active, as climate evolution without volcanism is uncertain. In the model it is assumed that weathering is entirely inactive when volcanism is inactive, as without volcanism there is no supply of fresh rock to the surface. However, even with no volcanism weathering can still take place, as long as there are still fresh, unaltered minerals in the weathering zone that groundwater can reach. Thus climate would likely cool over time without volcanism, until either the climate has cooled to a permanent snowball state or weathering ceases due to complete carbonation of the crust. Note that in Figure 2 climate is still temperate when volcanism ends; however in reality climate would likely cool further due to continued weathering and the absence of outgassing. In this study I assume that planets will always cool to a snowball state when volcanism ends, which is a conservative assumption; habitable climates may persist after volcanism has ended (or exist before it begins, in cases where the mantle is initially cool), so the estimates presented here give the minimum lifetimes for habitable climates. Habitable lifetime is thus calculated as the minimum of the time when volcanism ends and the time when temperatures fall outside of the 273-400 K range, minus the maximum of the time when volcanism begins and surface temperatures first enter the 273-400 K habitable range.

The lifetimes of habitable climates as a function of $Q_{0}$ and $C_{tot}$ are shown for both hot and cold start scenarios (Figure 3). At very low CO2 budgets of $\sim 10^{19}$ mol or less, habitable climates are not possible, as outgassing rates are not sufficient to prevent surface temperatures from immediately dropping below the freezing point of water. That is, even when most of the planet’s CO2 has been outgassed into the atmosphere-ocean reservoir, there is insufficient warming to keep the planet above freezing. With increasing $C_{tot}$, habitable climate lifetime increases, until CO2 budgets are so high that planets would enter the supply limited weathering regime. Here habitable climates are still possible for a narrow range of $C_{tot}$, where even with supply limited weathering the resulting surface temperature is still less than $400$ K. However, as $C_{tot}$ continues to increase, planets reach a point where climates are always too hot to support life as we know it. Thus there is an optimum range of $C_{tot}$ for habitability, from $~{}10^{19}-10^{23}$ mol.

Initial mantle heat production rate also plays a critical role in habitable climate lifetime, by dictating how long volcanism and outgassing can last. For the baseline model (Figures 3A & D), $Q_{0}\gtrapprox 50$ TW is needed for habitable climates to last at least 1 Ga. With increasing $Q_{0}$, the habitable lifetime increases in turn, as long as $C_{tot}$ is within the range where habitable climates are possible; this is because higher rates of radiogenic heat production prolong volcanism and outgassing. With a lower initial mantle temperature of $T_{i}=1500$ K (3B & E), the mantle is initially too cold for significant volcanism; thus radiogenic heating must warm the mantle up before volcanism and outgassing can begin. As a result, a higher $Q_{0}$ is needed before habitable climates lasting 1 Ga or more can be maintained, and habitable lifetimes are generally lower than in the baseline model, where the initial mantle temperature is larger. When the mantle reference viscosity is larger than the baseline value (3C & F), habitable lifetimes are longer because the mantle cools more slowly, and hence outgassing lasts longer.

A key feature of Figure 3 is that the habitable lifetimes for hot and cold start scenarios closely match. For the baseline model, there is almost no difference between models with either an initially hot or initially cold climate. With $T_{i}=1500$ K or $\mu_{r}=5\times 10^{21}$ Pa s there are some noticeable differences between the two initial conditions at low values of $Q_{0}$, but overall habitable lifetimes are still remarkably consistent. The differences are caused by the short lived nature of volcanism at low $Q_{0}$. For the hot climate initial condition, volcanism does not last long enough for metamorphic decarbonation to begin or for weathering and outgassing to come into balance; the band of $\sim$ 1 Ga habitable lifetimes at $Q_{0}\approx 25$ TW for $\mu_{r}=5\times 10^{21}$ Pa s and at $Q_{0}\approx 50$ TW for $T_{i}=1500$ K are caused by volcanism beginning, but not lasting long enough to draw down the initially CO2 rich atmosphere. A small window of temperate climates thus results. However in the cold start cases, higher rates of volcanism are needed to keep the climate warm, in light of weathering and the much lower initial atmospheric CO2 content. Thus higher heat production rates are required before even short lived habitable climates can develop with a cold start initial condition. Outside of these limited regions of parameter space, however, the initial condition does not have a significant influence on long term evolution of climate, outgassing, and mantle thermal evolution, as long as the initial climate state is not a snowball (see §4.2 & 4.3), or so hot that liquid water cannot exist and weathering cannot operate. A climate that starts off initially hot cools rapidly due to weathering, while one that starts cold warms due to outgassing, until both models forget their initial conditions and follow the same evolution.

4.2 Recovery from snowball states

However, different initial conditions do influence whether snowball states could be recovered from. As explained above in §3.2, for each model, I calculate the planet age after which both a hard snowball and soft snowball state can no longer be recovered from. As outgassing is required to recover from a snowball state, these ages are related to the lifetime of volcanism; without CO2 outgassing any snowball state would be permanent. The planet ages after which recovery is impossible for both hard and soft snowball states, with both hot and cold start initial conditions, are shown in Figure 4 for models with baseline parameter values. To recover from a hard snowball state, mantle volcanism must release at least $3\times 10^{19}$ mol of CO2 to the atmosphere. Thus, the mantle CO2 content is critical for whether a hard snowball state can be recovered from. With an initially hot climate, most of the planet’s CO2 resides in the atmosphere, before being weathered out into the crust, and then cycling through the crust and returning to the atmosphere via metamorphic decarbonation. The mantle starts with only a small fraction of the planet’s CO2, and does not gain CO2 from the surface reservoirs as decarbonation reactions do not allow for recycling of CO2 back into the mantle. As a result, planets with an initially hot climate can not recover from a hard snowball state at any point during their evolution, unless the CO2 budget is very large ($\gtrapprox 3\times 10^{23}$ mol). Such large values of $C_{tot}$ are expected to result in supply limited weathering and uninhabitably hot climates, so planets within the range of $C_{tot}$ that allow for habitable climates would not be able to recover from hard snowball states, if the climate is initially hot.

However, with an initially cold climate, more CO2 initially resides in the planet’s mantle, and recovery from a hard snowball state is possible. Even here, snowball recovery is only possible early in the planet’s evolution, as the mantle quickly becomes depleted in CO2 via outgassing. For example, with $C_{tot}=10^{21}$ mol, recovery from a hard snowball state is only possible for the first 500 million years of a planet’s history. Increasing $C_{tot}$ generally prolongs the time period over which recovery from a hard snowball is possible (Figure 4B). Within the range of $C_{tot}$ where habitable climates are possible, planets can recover from hard snowball states when they are younger than up to 1-1.5 Gyrs. The initial radiogenic heating rate is also important and has two competing effects; increasing $Q_{0}$ allows volcanism to last longer, and mantle volcanism is essential for recovery from a hard snowball state. However, increasing $Q_{0}$ also increases the rate of volcanism, and hence depletion of mantle CO2. Thus at low values of $Q_{0}$, increasing $Q_{0}$ allows snowball recovery to occur at later planet ages, because the end of volcanism is the limiting factor for whether a snowball state can be recovered from in this case. However, at even larger values of $Q_{0}$, increasing $Q_{0}$ lowers the age when snowball recovery becomes impossible, as the high rates of internal heating deplete the mantle of CO2 more rapidly.

An important end-member case is when a hard snowball state can never be recovered from, even if it occurs as a planet’s initial surface condition. This case is important because it marks the point where a planet’s initial climate state would dictate its later climate evolution; that is a planet that starts in a snowball would be stuck in this snowball for it’s entire history, while a planet that starts with a hot or temperate climate would be able to avoid a snowball and potentially develop a habitable climate for at least some portion of its evolution. As $2.5\times 10^{4}$ Pa of CO2 in the atmosphere, corresponding to $3\times 10^{19}$ mol, are required to recover from a snowball, a simple limit to when a hard snowball can never be recovered from is when $C_{tot}<3\times 10^{19}$ mol. This limit is plotted on Figure 4B as a dot-dashed white line, and corresponds well to the model results where initial snowball states can not be recovered from. Thus, planets with $C_{tot}<3\times 10^{19}$ mol would be permanently stuck in a snowball state if they are initially frozen over, unless the luminosity of the star they orbit increases to such an extent that it melts the ice layer.

Recovery from a soft snowball state is generally much easier than recovery from a hard snowball (Figures 4C & D). Soft snowball states can be recovered from for nearly the entire time period when volcanism is active. Thus, planets that manage to fall into a soft snowball state, due to some perturbation to their climate and carbon cycle, should be able to recover from this state throughout the time period when the models predict that habitable climates are possible. There is also no discernible difference between an initially hot climate state and initially cold climate state, in terms of their ability to recover from a soft snowball. As outgassing that leads to snowball recovery is dominated by metamorphic outgassing, the initial mantle CO2 content is no longer important. However, even in the soft snowball case there is still a limit where an initially frozen climate would remain frozen throughout the planet’s history, due to an inability to recover from the snowball climate. For a soft snowball, $1.6\times 10^{20}$ mol of CO2 must reside in the atmosphere-ocean system in order to melt the ice. Thus, planets with $C_{tot}<1.6\times 10^{20}$ mol will be unable to escape an initial snowball climate, and such a state would be permanent. The carbon budget below which snowball states are always permanent is actually higher for a soft snowball than hard snowball; this is because a significant portion of the outgassed CO2 will be sequestered in the ocean during a soft snowball, where it does not warm the climate, and because seafloor weathering will remain active, removing CO2 from the combined ocean-atmosphere system.

4.3 Hysteresis in stagnant lid planet habitability

The inability of planets with low CO2 budgets to recover from an initial snowball state means hysteresis is possible. Planets with an initially frozen climate may be stuck in this state permanently, and hence uninhabitable for surface life, while planets that start with an initially temperate or hot climate will be able to sustain habitable conditions, for at least a limited amount of time. Figure 5 shows the baseline model predictions for the lifetime of habitable climates for planets with an initially hot climate, and denotes the regions of parameter space where initial snowball states can not be recovered from. Specifically, with $C_{tot}<3\times 10^{19}$ mol an initial hard snowball state could not be recovered from. However, at this $C_{tot}$ limit, most planets will not be capable of sustaining habitable climates regardless of the initial condition; when $C_{tot}<10^{19}$ mol habitable climates can not be maintained for 1 Gyrs, even with $Q_{0}=250$ TW. Thus hysteresis based on initial condition is not a significant effect for a hard initial snowball. However, for a soft snowball hysteresis is more important, as a higher CO2 budget of $1.6\times 10^{20}$ mol is required to recover from this state. The $C_{tot}<1.6\times 10^{20}$ mol limit encompasses a significant portion of the region of $Q_{0}-C_{tot}$ space where habitable climates can be maintained. Another important aspect of hysteresis in stagnant lid planet climate evolution is that planets that fall into a snowball state at some later point during their evolution could become stuck in these states, as shown in Figure 4. In particular planets that enter a hard snowball state are unlikely to recover, unless significant quantities of CO2 still reside in the mantle, meaning planets one might otherwise predict to be habitable could end up stuck in an unrecoverable snowball, should some climate or carbon cycle perturbation kick them into this state during their evolution. Soft snowball states can generally be recovered from while volcanism is active, as long as $C_{tot}>1.6\times 10^{20}$ mol, so hysteresis during later planetary evolution is less likely in this case.

4.4 Sensitivity of model results

The overall model results are consistent when varying the initial mantle temperature, reference viscosity, viscosity activation energy, and when crustal decarbonation is ignored. However, there are some key differences worth discussing. Habitable lifetimes for a hot start initial condition are shown in Figure 6; only models with an initially hot climate are shown, as those with an initially cold climate show very similar habitable lifetimes. Decreasing the initial mantle temperature lowers the habitable climate lifetimes for all models, because the lower initial mantle temperature decreases the longevity of volcanism. Furthermore, higher values of $Q_{0}$ are required for habitable climates to last at least 1 Gyr, because the initial mantle temperature is not high enough for volcanism, and thus the mantle must heat before volcanism can commence (Figure 6B). With a higher initial mantle temperature of $T_{i}=2200$ K, there is no discernible difference compared to the baseline model, as volcanism can begin immediately for all models with either $T_{i}=2000$ K, as in the baseline models, or with $T_{i}=2200$ K (Figure 6C). Changing the reference viscosity has an important influence on habitable climate lifetime; with a larger reference viscosity, habitable climates last longer due to a slower cooling rate of the mantle, while with a lower reference viscosity the opposite is true, and habitable climate lifetimes are everywhere shorter (Figure 6D & E). The activation energy for viscosity influences habitable lifetimes in a similar manner. Lower values of $E_{v}$ prolong temperate climates by causing the mantle to cool more slowly, while a higher $E_{v}$ causes the mantle to cool more rapidly, and hence shortens the lifetime of temperate climates (Figure 6F & G). Finally, when decarbonation of the crust is assumed to not occur, such that all crustal CO2 recycles back into the mantle, higher total CO2 budgets are needed for long-lived temperate climates (Figure 6H). Specifically, $C_{tot}>10^{20}$ mol is needed for a habitable climate to last at least 1 Gyr, compared to $C_{tot}>10^{19}$ mol when decarbonation of the crust is considered. However, temperate climates extend to larger $C_{tot}$ in the no decarbonation case, before outgassing overwhelms weathering and leads to uninhabitability hot climates. The reason for this behavior is the much larger volume of the mantle compared to the crust. The same total CO2 budget will lead to overall lower outgassing rates when outgassing is primarily due to mantle volcanism rather than crustal decarbonation, as the large volume of the mantle dilutes the CO2 concentration of rocks experiencing melting and outgassing.

Similar effects of the key model parameters, discussed above, are seen when looking at planets’ ability to recover from soft snowball states (Figure 7). With a lower initial mantle temperature, planets with initial radiogenic heat production rates lower than $\approx 75$ TW would never be able to recover from a soft snowball state, due to a lack of volcanism; at higher $Q_{0}$ the age where snowball recovery becomes impossible is similar to the baseline model (Figure 7B). With $T_{i}=2200$ K, results are indistinguishable from the baseline model (Figure 7C). Consistent with the results for habitable climate lifetime, a higher reference viscosity extends the age range where snowball states can be recovered from, while a lower reference viscosity shortens this range (Figure 7D & E), as a result of how mantle viscosity influences mantle cooling rate. A lower activation energy for viscosity increases the planet age where snowball recovery becomes impossible, while a higher activation energy has the opposite effect, again due to the way activation energy influences mantle cooling rate (Figure 7F & G). Finally, when crustal decarbonation is assumed to not occur, the planet age when recovery from a soft snowball state is no longer possible is similar to the baseline model, that includes decarbonation (Figure 7H).

The planet age when recovery from a hard snowball state is no longer possible is shown in Figure 8, in this case for an initially cold climate, as with an initially hot climate recovery from a hard snowball is nearly always impossible. The results are similar to those for soft snowball recovery, with the same effects of initial temperature and viscosity. One notable case is when decarbonation of the crust is assumed to not occur; here the age when hard snowball recovery becomes impossible shows a very similar pattern in $C_{tot}-Q_{0}$ space to the soft snowball recovery case, but with lower ages throughout the parameter space (Figure 8H). The reason for this is that without decarbonation of the crust, the mantle CO2 reservoir is constantly resupplied by foundering of carbonated crust, and thus not quickly depleted as in the other models.

Finally, the effect of uncertainties in the seafloor weathering function is explored (Figure 9). Only the influence of the seafloor weathering parameters on the habitable climate lifetime is shown, as the maximum age when planets can recover from either a soft or hard snowball climate state is almost entirely unaffected by varying these parameters. Changing the seafloor weathering activation energy, $E_{sfw}$, does not significantly change the longevity of habitable climates, though with a larger $E_{sfw}$, habitable climates can last longer at low $C_{tot}\sim 10^{19}$ mol, as a result of more effective climate buffering. The reference seafloor weathering rate, $F_{sfw}^{*}$ has a similar effect; when $F_{sfw}^{*}$ is lower habitable climates can last longer at low values of $C_{tot}$ then when $F_{sfw}^{*}$ is higher, as the lower overall weathering rate means more CO2 can build up in the atmosphere and keep the climate warm. The parameter with the biggest influence is the exponent describing the dependence of seafloor weathering on volcanic eruption rate, $\beta$. With $\beta=1$, i.e. a strong dependence of seafloor weathering rate on volcanic eruption rate, planets with lower initial radiogenic heating rates and lower total CO2 budgets ($C_{tot}<\sim 10^{20}$ mol) can maintain habitable climates for longer than when $\beta=1/2$ (the baseline value) or $\beta=0$. The reason for this is that a strong dependence of seafloor weathering on eruption rate causes the seafloor weathering rate to decline as volcanism wanes; as a result, lower outgassing rates still allow enough CO2 to remain in the atmosphere for a temperate climate to exist.

5 Discussion

The results are generally consistent with Foley & Smye (2018), which constrained the lifetime of potentially habitable climates on stagnant lid planets by calculating outgassing rates and estimating when supply limited weathering was expected to prevail. However, this study improves on the results of Foley & Smye (2018) by explicitly modeling atmospheric CO2 abundance and climate evolution, to directly constrain when habitable climates can exist. As a result habitable climates are found to exist over a wider range of $C_{tot}$ than in Foley & Smye (2018); $C_{tot}\sim 10^{19}-10^{23}$ mol, compared to $C_{tot}\sim 10^{20}-10^{22}$ mol. The reason for the difference is two fold. First, seafloor weathering results in generally lower CO2 drawdown rates than continental weathering; the present day Earth seafloor weathering rate, $F_{sfw}^{*}$, determines the overall scale of the weathering rate, and $F_{sfw}^{*}$ is approximately a factor of 10 lower than the modern day continental weathering rate (Krissansen-Totton & Catling, 2017). As a result, lower outgassing rates can still support a temperate climate, as the low CO2 drawdown rate allows more CO2 to remain in the atmosphere. A consequence of the slower rate of seafloor weathering is that this study finds outgassing rates need only be $\sim 10$ % Earth’s present day degassing rate to sustain temperate climates, consistent with the lower bound explored in Foley & Smye (2018) and with Haqq-Misra et al. (2016). Second, while supply limited weathering begins at $C_{tot}\sim 10^{22}$ mol in this study, consistent with Foley & Smye (2018), uninhabitabily hot climates do not develop until $C_{tot}\sim 10^{23}$ mol. That is, planets must be well within the supply limited weathering regime before sufficient CO2 builds up in the atmosphere to result in extremely hot climates. Thus, the prospects for habitability of stagnant lid planets are even better than previously estimated.

There are important caveats, however. As highlighted in §4.3, at low CO2 budgets, $C_{tot}<\sim 10^{20}$ mol, planets would be incapable of ever recovering from a snowball state, including just after their formation. Thus these planets would experience drastically different climate evolutions based on their initial conditions; if they start in a snowball state, this snowball will be permanent and the planet will never be habitable for surface life, while if they start with a temperate or hot climate, then a habitable climate can persist for $\approx 1$ Gyr or more, depending on the radiogenic heating rate. Planets with these low CO2 budgets also would be unable to recover from a snowball climate state if it developed later in the planet’s history as well. However, planets with larger CO2 budgets, $C_{tot}>\sim 10^{20}$, mol can recover from initial snowball states. Such planets are found to be capable of recovering from soft snowball states, should one develop, throughout nearly their entire lifetime spent with significant outgassing and a temperate climate. However, hard snowball states are much more difficult to recover from, and planets developing such a climate will often be stuck in this state for the remainder of their lifetime.

The initial distribution of carbon between mantle and atmosphere is thus important for planetary climate evolution. This distribution is highly uncertain, but most studies of magma ocean solidification indicate that this process depletes the mantle, and enriches the atmosphere, in CO2 (e.g. Zahnle et al., 2007; Elkins-Tanton, 2008; Hamano et al., 2013; Salvador et al., 2017). Thus planets that have initially hot climates are probably more common, as high energy collisions during planet formation are expected to lead to magma oceans (Elkins-Tanton, 2012). As a result, initial snowball states that become permanent are probably rare for the “typical” stagnant lid planet, though “typical” stagnant lid planets would be unable to recover from a hard snowball state, should such a state develop later during the planet’s evolution. An important implication of the hysteresis in climate evolution is that some planets that would be predicted to be habitable and possess temperate climates, based on their CO2 and heat budgets, may in fact lie in permanent snowball states.

Another important caveat is that changes in solar luminosity are not considered in this study. By using a climate parameterization based on the present day Earth, a present day solar luminosity is implicitly assumed. In reality, stars are dimmer earlier in their history, and luminosity increases over time (e.g. Gough, 1981). This change in luminosity could impact both the lifetime of habitable climates predicted by the thermal evolution models, and the ability of stagnant lid planets to recover from snowball climates. In particular, planets could be colder than our models predict early in their history, and warmer later on. However, seafloor weathering does include a strong negative feedback that acts to stabilize climate in the face of varying solar luminosity (Brady & Gíslason, 1997; Coogan & Gillis, 2013; Coogan & Dosso, 2015; Krissansen-Totton et al., 2018). In fact, seafloor weathering can buffer climate against changes in luminosity more effectively than continental weathering, because seafloor weathering contains only a very weak direct dependence on atmospheric CO2, through changes in ocean pH, and is instead much more sensitive to temperature; this results in stronger buffering with solar luminosity (Krissansen-Totton et al., 2018). The main effect changes in solar luminosity would have on the habitable lifetimes calculated here is that at low luminosity, larger $p\mathrm{CO}_{2}$ is needed in order to keep the climate temperate. Thus when $C_{tot}$ is low enough that there is not enough carbon available to meet the required $p\mathrm{CO}_{2}$ for a temperate climate, habitability will not be possible until luminosity increases. This threshold value of $C_{tot}$ is expected to be higher with a lower luminosity than the $C_{tot}\sim 10^{19}$ threshold we find in the models presented here. Likewise, if luminosity is higher, then lower values of $C_{tot}$ are expected to result in uninhabitability hot climates, than what is presented in this study.

Solar luminosity also has an important influence on recovery from snowball climates. With a lower luminosity, more CO2 in the atmosphere is required to melt the ice layer, while with a higher luminosity less CO2 is required. Thus, recovery from both a hard and a soft snowball state will be harder when luminosity is lower, and easier when luminosity is higher. In particular, the value of $C_{tot}$ required to recover from either a hard or soft snowball climate will increase when luminosity is lower, and decrease when luminosity is higher. Fortunately, the effect of luminosity on snowball recovery can easily be captured with the models presented here, as it was found that the limit where a snowball state can never be recovered from is simply when the CO2 budget is lower than what is needed to melt an ice layer. Thus, one could determine different limits on $C_{tot}$ needed to recover from snowball states at different solar luminosities, and the models in this study would indicate that planets with CO2 budgets below this limit will never be able to recover from a snowball.

As seen in §4.4, the model results are largely consistent over a wide range of parameter values. Moreover, Foley & Smye (2018) performed even more extensive testing, in particular looking at the influence of incomplete degassing of both the crust and mantle, and found that these effects also do not significantly impact the results. I note that as the thermal evolution model and supply limited weathering formulation used in this study are similar to Foley & Smye (2018), the same limitations of this model discussed in Foley & Smye (2018) apply here. There are interesting implications of how some of the key parameters influence the prospects for habitability of stagnant lid planets, however. In particular, the higher the mantle reference viscosity, or lower the activation energy, the longer volcanism lasts and the longer habitable climates persist. Thus planets with mantle mineralogies that lead to larger overall viscosities, or lower activation energies, compared to Earth would be better homes for life. Ultimately the models of stagnant lid planet habitability presented here could be used to constrain the likelihood of finding temperate climates on exoplanets as they are discovered. In particular if abundances of radionuclides and major rock forming elements in the star can be measured, and used to constrain the composition of the planet, direct estimates of the likelihood of habitability can be made, providing a valuable tool in the search for extra-terrestrial life.

6 Conclusions

Models coupling mantle thermal evolution, volcanism, outgassing, CO2 cycling, and climate evolution for Earth-sized stagnant lid planets constrain the conditions that lead to long-lived temperate climates on such planets. Specifically, the lifetime of habitable climates is calculated as a function of the initial radiogenic heating budget of the mantle and total CO2 budget of the mantle and surface reservoirs. Planets with CO2 budgets of $C_{tot}\sim 10^{19}-10^{23}$ mol are found to sustain habitable climates for 1-5 Gyrs. When $C_{tot}$ is lower than this limit planets will remain frozen over for their entire lifetime, and above this limit planets will remain uninhabitably hot. The initial radiogenic heating rate is also important; heating rates of $\approx 50$ TW are needed for habitable climates to last for at least 1 Gyr, and increasing the heating rate increases the lifetime of volcanism, and hence temperate climates. Stagnant lid planets are assumed to be dominated by seafloor weathering in this study, and the strong climate buffering ability of seafloor weathering, along with its lower overall weathering rate as compared to continental weathering on Earth, allow a wide range of CO2 budgets to support long-lived habitable climates.

The ability of stagnant lid planets to recover from potential snowball episodes is considered. If a hard snowball, where there is no exchange between atmosphere and ocean, develops, recovery is very difficult as most CO2 outgassing occurs via metamorphic decarbonation of the crust, below the ice layer. With a soft snowball, where there is atmosphere-ocean exchange, recovery is possible under most conditions. However, both snowball scenarios place limits on the CO2 budget a planet must possess in order to ever recover from a snowball state. Planets with very small CO2 budgets, lower than $\sim 10^{19}-10^{20}$ mol, would be unable to recover from a snowball climate, even if this occurs at the start of the planet’s evolution. Thus such planets could exhibit hysteresis in their evolution, where planets that start off with a snowball climate, or develop one later, will be stuck in these states permanently, while planets that start with a warm or hot climate will be able to sustain habitable climates for 1 Gyr or more. Hysteresis means that planets one might otherwise expect to be habitable, based on their radiogenic heat and CO2 budgets, may instead be found in a permanent snowball state. However, it is not clear that planets are likely to become stuck in these snowball states, as magma ocean solidification tends to release mantle CO2 to the atmosphere, and produce an initially warm climate. Overall, the results of this study can be used to estimate the planetary factors that allow for long-lived habitability, aiding in target selection for future missions searching for biosignatures on exoplanets.

I thank an anonymous reviewer for comments that helped to significantly improve the manuscript.

Appendix A Crustal carbon reservoir treatment

The models presented in the main text treat the crustal carbon reservoir as a single reservoir with instantaneous mixing. However, in reality there is no mixing in the crust, as it is part of the rigid lid, and crust decarbonating had its CO2 concentration set as that parcel of crust formed at the surface. Such a simplification greatly helps with formulating the model, and here I show that treating the crustal carbon reservoir as a single, well mixed reservoir, does not significantly affect the results presented in the main text. The metamorphic outgassing flux calculated as in the main text is nearly identical, throughout most of a stagnant lid planet’s thermal evolution, to the metamorphic outgassing flux that results from assuming that the CO2 concentration of decarbonating crust is set when that parcel of crust forms.

The concentration of CO2 in crust formed at time $t$ is $C_{crust}=F_{weather}(t)/f_{m}(t)$, while the velocity of downgoing crust is $w_{crust}(t)=f_{m}(t)/A_{surf}$. The age of crust that is decarbonating at any time $t$ can be determined by finding the crust age, $t_{crust}$, that satisfies

[TABLE]

Thus the metamorphic outgassing rate where crustal CO2 concentration varies with depth, is given by

[TABLE]

The metamorphic outgassing rate as formulated in the main text, $F_{meta}$, and $F_{meta}^{\prime}$ are compared in Figure 10 for a few representative models, while the decarbonation depth and velocity of downgoing crust, $w_{crust}$, is shown in Figure 11 for the same models. The two metamorphic outgassing rates are nearly identical for most of the period over which volcanism is active. There are differences in outgassing rate very early in planet evolution, in particular for the cases where the mantle viscosity is high or initial mantle temperature is low; in these cases the mantle must heat up before volcanism can begin, and hence the start of CO2 outgassing is delayed. In either formulation, metamorphic outgassing begins at the same time, as even with the formulation in the main text where instantaneous mixing in the crustal CO2 reservoir is assumed, metamorphic outgassing can not take place until the crust first extends deep enough to reach the decarbonation depth. Also note that the crust is always at least somewhat carbonated, as all crust that forms experiences weathering as it is deposited at the surface.

Very early in model evolution there can be significant differences between the two outgassing formulations as the weathering rate is either initially rapid (for a hot start initial condition) or sluggish (for a cold start). In particular, with $T_{i}=1500$ K or $\mu_{r}=5\times 10^{21}$ Pa s, there is a rapid draw down of CO2 that causes a cold climate, which then warms as metamorphic outgassing ramps up. During this cold climate phase, the CO2 concentration in the crust is reduced, and hence $F_{meta}^{\prime}$ shows a short-lived dip. However, at later times for all models, the majority of the planet’s CO2 has been outgassed from the mantle and is circulating through the crust. Over long time scales, the weathering rate balances the outgassing rate, meaning that the concentration of CO2 in newly created crust is $F_{meta}^{\prime}/f_{m}$. With $F_{meta}(t)^{\prime}=C_{crust}(t_{crust})f_{m}(t)$, the concentration of CO2 in newly created crust is the same as the concentration of CO2 in the crust at the decarbonation depth. Thus once metamorphic outgassing becomes the dominant CO2 source to the atmosphere, the concentration of CO2 in the crust is constant, and the outgassing rate simply declines as the melt production rate declines. While the simple formulation of metamorphic outgassing used in the main text might miss some interesting changes in outgassing, and hence climate, early in a planet’s history, it is a good approximation for the metamorphic outgassing rate for most a planet’s history. The overall results of this study, which focus on the longevity of habitable climates on stagnant lid planets, are therefore unlikely to be significantly affected by the assumption of an instantaneously mixed crustal carbon reservoir.