Carbon Cycling and Habitability of Earth-Sized Stagnant Lid Planets

1. Introduction

Thousands of exoplanets have been discovered to date, and many more are certain to follow (e.g., Batalha, 2014), raising the question of which, if any, of these planets might be suitable hosts for life. As a result, understanding the factors that allow a planet to be habitable, in particular for surface life that can be detected by remote observations, has become an increasingly important goal. While the question of what requirements a planet must meet in order to potentially harbor life has long been a topic of research, the recent boom in exoplanet discoveries brings more urgency to this question. Theoretical ideas about planetary habitability are now potentially testable with exoplanet observations, and theoretical considerations can be used to guide astronomers searching for life to those planets that are most likely to be habitable.

In order to host surface life, a planet must receive neither too much nor too little solar radiation for liquid water to be stable on the planet's surface. The focus for this study, and for most exoplanet astrobiology, is on surface life, because surface life can produce detectable signatures in a planet's atmosphere or on its surface. Life is certainly possible in oceans that lie beneath surface ice layers, such as those on Europa or Enceladus. However, such life would be difficult to detect by remote observations. The constraint that solar flux not preclude the existence of liquid surface water defines the “habitable zone,” the range of orbital distances around a star where liquid water may be stable (e.g., Kasting et al., 1993). However, lying within the habitable zone does not guarantee that a planet's surface conditions will be conducive to life. Extreme variations in greenhouse gas concentrations can cause temperatures far in excess of what Earth-like life can tolerate (the upper limit for known life is ≈400 K; Takai et al., 2008) and temperatures low enough to cause global glaciation (e.g., Walker et al., 1981; Sleep and Zahnle, 2001; Abbot et al., 2012; Kadoya and Tajika, 2014; Menou, 2015; Foley and Driscoll, 2016). On Earth, such variations are held in check by the long-term carbon cycle and the stabilizing effect it has on climate (e.g., Walker et al., 1981; Berner and Caldeira, 1997). As temperature and atmospheric CO₂ content increase, rates of chemical weathering increase in turn, thereby drawing CO₂ out of the atmosphere and acting to cool the climate. Likewise, when temperature and atmospheric CO₂ content decrease, weathering rates slow down and volcanic CO₂ accumulates in the atmosphere, acting to warm the climate (e.g., Kump et al., 2000; Berner, 2004).

The long-term carbon cycle on Earth is facilitated in many ways by plate tectonics (e.g., Kasting and Catling, 2003). Plate tectonics drives uplift and orogeny, which in turn spurs erosion and enhances weathering (e.g., West, 2012). Plate tectonics also leads to the creation of continents, which increase the area of exposed land and also enhance weathering (e.g., Foley, 2015). Finally, plate tectonics also promotes long-lived CO₂ degassing by driving volcanism at mid-ocean ridges and island arcs, and by recycling surface CO₂ into the mantle such that the mantle CO₂ reservoir is continuously resupplied. However, whether plate tectonics is required for the long-term carbon cycle to operate, and therefore stabilize a planet's climate, is not known (Foley and Driscoll, 2016; Lenardic et al., 2016).

A particularly interesting question is whether any form of climate regulation caused by chemical weathering can be maintained on a planet in the stagnant lid regime. In the stagnant lid regime, the lithosphere is rigid and nearly immobile, and no subduction occurs; convection only takes place beneath the rigid lid, in the form of driplike downwellings from the base of the lid (e.g., Ogawa et al., 1991; Davaille and Jaupart, 1993; Solomatov, 1995). As a result, recycling of surface material back into the mantle, in particular volatiles like carbon and water, is seemingly limited. The rigid, immobile surface also limits tectonic uplift, as the typical processes driving orogeny, that is, continental collision and subduction, are no longer active. However, stagnant lid planets still experience volcanism, which can release mantle CO₂ to the atmosphere, create fresh rock and topography at the surface, and potentially even allow for some surface recycling through burial by lava flows (e.g., Pollack et al., 1987; Hauck and Phillips, 2002; Hirschmann and Withers, 2008; Gillmann and Tackley, 2014; Lenardic et al., 2016).

This paper sets out to test whether the long-term carbon cycle can regulate atmospheric CO₂ levels, and thus maintain a stable climate over geological timescales, on a planet in the stagnant lid regime. Previous studies have considered the ability of stagnant lid planets to maintain temperate climates, primarily by focusing on CO₂ degassing rates on 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 exoplanets (Noack et al., 2017). However, a systematic study that considers both degassing and weathering limits to climate stability, for a range of mantle volatile and heat budgets, is lacking. We restrict our study to planets with an Earth-like size and composition. Little is known about mantle rheology, volcanism, and outgassing on planets with different bulk compositions, and hence different mantle mineralogies, than Earth, so it would be premature to model such planets. Moreover, the rheology of super-Earth mantles, even for an Earth-like composition, is also highly uncertain, as super-Earth mantles reach pressures much higher than experimental studies can currently attain (e.g., Karato, 2011). We thus focus on Earth-sized planets to avoid such uncertainties.

The potential for habitability on stagnant lid planets is important to assess because they are common in the Solar System; only Earth has a mantle that convects in a plate-tectonic regime (e.g., Breuer and Moore, 2007). Moreover, predicting whether rocky exoplanets will have plate tectonics or stagnant lids is difficult due to our incomplete understanding of lithospheric rheology, exoplanet heat budgets, and material properties at extremely high pressure (O'Neill and Lenardic, 2007; Valencia et al., 2007; Kite et al., 2009; Korenaga, 2010b; Stamenkovic et al., 2011; van Heck and Tackley, 2011; Foley et al., 2012; Lenardic and Crowley, 2012; Noack and Breuer, 2014). Thus, if plate tectonics is not required for harboring life, then the search for habitable exoplanets is greatly simplified.

1.1. Modeling habitability of a stagnant lid planet

In order to test whether temperate climates can be sustained on Earth-like stagnant lid planets, we couple simple models for mantle thermal evolution, CO₂ outgassing, weathering, and crustal growth. The modeling focuses on two main issues: (1) if the supply of weatherable rock is sufficient for chemical weathering to balance the CO₂ being released by magmatism and metamorphism, and thus allow a weathering feedback to operate; and (2) for how long, if at all, CO₂ outgassing rates high enough to prevent a planet's entire surface from becoming frozen over can be maintained.

The first criterion is important to habitability, because the rate at which fresh, unweathered rock is brought to the near-surface environment, replenishing rock that has been altered via weathering, sets the ultimate upper limit on the CO₂ drawdown rate. In other words, the upper limit to CO₂ drawdown occurs when weathering immediately and completely alters all available fresh rock, as soon as it is supplied to the weathering zone (e.g., West et al., 2005; Mills et al., 2011; Foley, 2015), a situation known as “supply-limited” weathering (e.g., Stallard and Edmond, 1983; Kump et al., 2000; Riebe et al., 2004; West et al., 2005). On a planet where weathering becomes globally supply-limited, the CO₂ drawdown rate only depends on the rate at which fresh rock is brought to the surface and no longer depends on temperature or atmospheric CO₂ content. Thus weathering cannot adjust to balance CO₂ outgassing, accumulation of atmospheric CO₂ continues unabated, and an inhospitably hot climate develops, as long as volcanism lasts long enough to significantly increase atmospheric CO₂ (Foley, 2015).

In contrast, when weathering is not supply-limited, a feedback that stabilizes temperate surface temperatures is possible. The exact mechanism behind this feedback is unclear. In particular, the extent to which it depends on temperature as a result of weathering being limited by mineral dissolution reaction kinetics (e.g., Kump et al., 2000; Berner, 2004), precipitation rates as a result of limitation by pore water flow rates (Maher and Chamberlain, 2014), or some combination of the two, is still a topic of research. Given these large uncertainties involved in how silicate weathering operates, especially when extrapolated to exoplanets, we chose to simply determine when weathering is supply-limited and when it is not in our modeling. Calculating actual climate states is beyond this study's scope.

On Earth, physical erosion, which is greatly aided by large topographic gradients provided by mountain building and orogeny, is the primary mechanism for replenishing mineral supply for weathering on land (e.g., West et al., 2005; West, 2012), while mid-ocean ridge volcanism supplies fresh rock on the seafloor (e.g., Alt and Teagle, 1999). However, on a stagnant lid planet, mountain building is muted by the lack of plate collisions. We thus assume here that volcanism is the sole source of fresh, weatherable rock on a stagnant lid planet, and therefore associate volcanic eruption rates with the rate of mineral supply for weathering (as explained in more detail in Section 3).

The second criterion is important for habitability because even with an active weathering-climate feedback, low CO₂ outgassing rates can induce global glaciation (Kadoya and Tajika, 2014; Haqq-Misra et al., 2016). We assume that CO₂ outgassing is a result of mantle volcanism, which releases mantle CO₂, and metamorphic decarbonation of the crust, originally carbonated by weathering, as it is buried by subsequent lava flows. Thus the longevity and strength of CO₂ degassing is largely controlled by the longevity and rate of volcanic activity. The formulation for degassing, as well as for stagnant lid thermal evolution, is given in the next section. Note that our model results map out conditions where habitable climates are possible on stagnant lid planets, but cannot guarantee that such climates will exist in reality. Even when weathering is not supply-limited, and balance between weathering and degassing can be established, inhospitably hot climates are still possible due to very high degassing rates or weakening of the feedback between weathering and climate by factors such as a lack of topography (e.g., Maher and Chamberlain, 2014). Moreover, complexities not included in our modeling, such as additional greenhouse gases or atmospheric loss to space, could also cause uninhabitable climates to form even with sufficient CO₂ outgassing and fresh mineral supply rates.

Our model makes some important further assumptions. We assume model planets lie within the habitable zone and possess a large supply of surface water, such that weathering can occur (e.g., Berner, 2004; Maher and Chamberlain, 2014); that is, we assume that the availability of surface water and precipitation does not limit weathering. We assume that weathering can take place both on exposed, or subaerial, land and in crust on the seafloor (i.e., seafloor weathering), as seafloor weathering has been shown recently to be an important aspect of the carbon cycle on Earth (Gillis and Coogan, 2011; Coogan and Gillis, 2013; Coogan and Dosso, 2015; Krissansen-Totton and Catling, 2017). As a result, we do not make any explicit distinction between weathering on land or on the seafloor; all CO₂ released to the atmosphere is assumed to be weathered out and stored in the crust, so long as the upper bound on the weathering rate provided by the supply of weatherable rock is not exceeded. We also assume that the mantle has an oxidation state similar to Earth's mantle, such that CO₂ is a primary outgassing product of mantle volcanism. Such an assumption is reasonable for planets Earth-sized and larger, as mantle oxidation state is thought to largely be set by chemical reactions that occur in the lower mantle early in a planet's history, and are expected to occur on any planet large enough to have a mantle dominated by bridgmanite (Wade and Wood, 2005). Photolysis of water during a magma ocean stage can also contribute to mantle oxidation (Hamano et al., 2013).

2. Theory

2.1. Thermal evolution model

We model the thermal evolution of stagnant lid planets by assuming pure internal heating and that heat loss results from heat conducted across the lithosphere as controlled by mantle convection and from mantle melting and magmatism. Magmatism contributes to cooling the mantle because cooling of the magma from the temperature it is erupted or emplaced at to the surface temperature, and the release of latent heat upon solidification, both represent a net heat loss from the mantle (e.g., Richter, 1985; Bickle, 1986; Hauck and Phillips, 2002; Moore, 2003; Nakagawa and Tackley, 2012; Driscoll and Bercovici, 2014). This simple model captures the first-order controls on the thermal evolution of a rocky planet and greatly simplifies the model by allowing us to focus solely on mantle evolution. Including core heat flow would require a model for the thermal evolution of the core, adding more poorly constrained parameters to the problem. Finally, our assumptions result in a model that is conservative, in that it will give the fastest reasonable cooling history and hence earliest shutoff to large-scale volcanism and CO₂ degassing. Including heat flow into the mantle due to core cooling would only act to prolong volcanism (see Section 5 for further discussion). As a result, the model will produce minimum estimates for the lifetime of potentially habitable climates on stagnant lid planets.

For magmatic heat loss, we assume that all the melt generated rises to the surface and contributes to cooling the mantle and the growth of the crust. However, only a fraction of the melt is assumed to actually erupt on the surface (an important factor for weathering, as discussed in Section 3), with the rest solidifying at shallow depth; this assumption is analogous to mid-ocean ridge volcanism where only the layer of pillow lavas is erupted, with sheeted dikes and gabbros solidifying in the subsurface. As the whole melt layer is emplaced in the near-surface environment, we assume all of the melt contributes to mantle heat loss. This assumption is also conservative, as it will give the maximum mantle cooling rate and thus a minimum estimate of degassing lifetime. Melt migration through a thick stagnant lid is not well understood, however, and melt could stall during its ascent through the lid, or even at the lid base. Melt stalling would slow mantle cooling and prolong volcanism, but also reduce outgassing. Test cases where some melt is assumed to stall at the base of the stagnant lid are shown in Appendix B and give results consistent with those in the main text, where all melt rises to the surface.

Mantle energy balance under the preceding assumptions results in the following equation for the evolution of mantle temperature (e.g., Stevenson et al., 1983; Hauck and Phillips, 2002; Reese et al., 2007; Fraeman and Korenaga, 2010; Morschhauser et al., 2011; Driscoll and Bercovici, 2014):

where V_man is the volume of convecting mantle beneath the stagnant lid, ρ 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 (base of stagnant lid), F_man is the heat flux from the convecting mantle into the base of the stagnant lid, f_m is volumetric melt production, ρ_m is the density of mantle melt, L_m is the latent heat of fusion of the mantle, and ΔT_m is the temperature difference between the melt erupted at the surface and the surface temperature. The volume of the convecting mantle is given by V_man = (4/3)π((R_p − δ)³ − R_c³), where R_p is the planetary radius, R_c is the core radius, and δ is the thickness of the stagnant lid. The surface area of the convecting mantle is then A_man = 4π(R_p − δ)². In this study, we focus on Earth-sized planets and therefore use R_p equal to Earth's radius and R_c equal to Earth's core radius (see Tables 1 and 2 for all model parameters and variables, their assumed values, and units, and Fig. 1 for a sketch of the model setup).

The thickness of the stagnant lid, δ, evolves with time based on an energy balance model (e.g., Schubert et al., 1979; Spohn, 1991),

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. The melt heat loss terms that appear in (1) do not appear here because we assume that heat extracted from the mantle due to volcanism is entirely lost to the surface, rather than to the stagnant lid. The mantle heat flux, F_man, and lid base temperature, T_l, are calculated from scaling laws for mantle convection (Reese et al., 1998, 1999; Solomatov and Moresi, 2000; Korenaga, 2009):

and

Here, c₁ is a constant, T_s is the surface temperature, d = R_p − R_c is the whole mantle thickness, θ is the Frank-Kamenetskii parameter, defined as θ = E_v(T_p − T_s)/(RT_p²), which describes the strength of the temperature dependence of mantle viscosity, E_v is the activation energy for mantle viscosity, R is the gas constant, a_rh is a constant, and Ra_i is the internal Rayleigh number. We use c₁ = 0.5 and a_rh = 2.5. Internal Rayleigh number is defined as Ra_i = ρgα(T_p − T_s)d³/(κμ_i), where g is gravitational acceleration, α is thermal expansivity, T_s is surface temperature, κ is thermal diffusivity, and μ_i is the interior mantle viscosity, defined at temperature T_p. Mantle viscosity is temperature-dependent following typical flow laws for diffusion creep in olivine as μ_i = μ_n exp(E_v/(RT_p)), where μ_n is a constant. We use E_v ≈ 300 kJ mol⁻¹ as the baseline value in this study (Karato and Wu, 1993); test cases using E_v = 200 kJ⁻¹ mol and E_v = 400 kJ⁻¹ mol are shown in Section 4.2 and are consistent with our baseline model results. We also define the reference viscosity as μ_r = μ_n exp(E_v/(RT_r)), where T_r = 1623 K is Earth's present-day mantle temperature. Using our baseline value of μ_n =4 × 10¹⁰ Pa s, μ_r ≈ 2 × 10²⁰ Pa s.

Note that in (3) the temperature difference between the mantle interior and surface, T_p − T_s, and whole mantle thickness, d = R_p − R_c, are used instead of the temperature difference between the mantle interior and base of the lid, T_p − T_l, and thickness of the actively convecting mantle, d − δ. However, in the heat flux scaling law, both the mantle thickness and temperature difference cancel out, so the equation for the heat flux to the base of the lid is independent of the definition of mantle thickness or temperature difference.

We calculate the temperature profile in the lid using a one-dimensional diffusion equation in steady-state,

where x is the radiogenic heat production rate per unit volume. Within the crust, that is, when R_p − δ_c ≤ z ≤ R_p, x = Q_crust/V_crust ≡ x_c, where δ_c is the thickness of the crust, Q_crust is the radiogenic heating rate within the crust, V_crust is the volume of the crust, and x_c is the crustal heat production rate per unit volume. Within the mantle lithosphere, that is, when R_p − δ ≤ z ≤ R_p − δ_c, x = Q_man/(V_man + V_lid) ≡ x_m, where x_m is the mantle heat production rate per unit volume and V_lid = (4/3)π((R_p − δ_c)³ − (R_p − δ)³) is the volume of the subcrustal stagnant lid. Both the heat production, x, and thermal conductivity, k, are assumed constant in z; the difference between crustal and mantle thermal conductivity is ignored for simplicity.

Our temperature equation also ignores advective heat transport, which occurs when crust newly generated by volcanism pushes the preexisting crust down into the mantle. The Peclet number, Pe = Λu/κ, where Λ is a length scale, u is velocity, and κ is thermal diffusivity, describes the relative importance of heat diffusion to advective transport, such that at small Peclet numbers advection is negligible. Typical eruption rates in our models are on the order of ∼100 km yr⁻³ (see Section 4.3), which results in a downward velocity of crust of ≈6 × 10⁻¹² m s⁻¹ when crust generation is assumed to be uniform across the planet's surface. With a typical lid thickness of ∼100 km (see Section 4.1) and thermal diffusivity κ = 10⁻⁶ m² s⁻¹, Pe  ≈ 0.6. Thus diffusion is more important than advection, but advection is not entirely negligible. Advection of cold material from the surface downward through the lid will increase the magnitude of the lid-base temperature gradient above what (5) would give. Thus, from (2), the lid thickness will be underestimated by a small degree in our model as advection is neglected. As a result, conductive cooling of the planet is overestimated by our models, which is consistent with our other model simplifications that act to enhance mantle cooling. Very rapid downward advection of crust could also prevent metamorphic decarbonation (discussed further in Section 2.3), as it would suppress the geotherm at depth. In Section 4.2 we explore how a lack of metamorphic decarbonation influences our results.

Equation (5) is solved by using the boundary conditions T = T_s at z = R_p, T = T_l at z = R_p − δ, and by enforcing continuity of temperature and heat flux across the crust-mantle interface, that is, at z = R_p − δ_c. The lid temperature profile can be calculated analytically, meaning the temperature gradient at the base of the lid also has an analytical expression, and we can write

when δ > δ_c, and

when δ = δ_c. The temperature at the base of the crust, T_c, is given by

As described in more detail in Section 2.2, we assume all crust below the lid founders into the mantle, such that δ_c at all times is less than or equal to δ. As the crust grows due to volcanism and δ_c approaches δ, the (δ − δ_c)⁻¹ term in (6) forces the numerical solution to use very small time steps, greatly reducing computational efficiency. To prevent this issue from severely hampering code performance, we switch to (7) when T_c is within one degree kelvin of T_l, which is equivalent to switching when the difference between δ_c and δ is negligibly small.

2.2. Melting and crustal evolution

The melt production rate, f_m, is calculated in the same manner as in Fraeman and Korenaga (2010). The pressure where melting begins, P_i, is

using a parameterization based on the solidus for dry peridotite from Korenaga et al. (2002). Here, γ_mantle is the average adiabatic temperature gradient in the upper mantle (≈ 10⁻⁸ K Pa⁻¹). Melting is assumed to stop at the base of the rigid lid, giving the final pressure of melting as

where ρ_l is the average density of the crust and lithosphere (assumed to be ρ_l = 3300 kg m⁻³). The melt fraction, φ, is

where (dφ/dP)_S ≈ 1.5 × 10⁻¹⁰ Pa⁻¹.

The melt production rate is given by the product of the melt fraction and the flux of mantle material into the melting zone. Assuming a cylindrical region of upwelling mantle (as in, e.g., Reese et al., 1998), surrounded by downwelling, the volumetric flow rate of mantle into the melting region, which is equal to the flow rate of mantle out the side of the cylinder, is given by v(d_m − δ)πL, where v is the convective velocity, d_m = P_i/(ρ_lg) is the depth where melting begins, and L is the diameter of the cylindrical upwelling region. Dividing by the area of the circle with a radius equal to the horizontal length of convection cells, here assumed to be d, and multiplying by the surface area of the planet, the total melt production rate is

where we have assumed that d ≈ 0.45R_p as on Earth. The formulation for melt production rate assumes that downwellings are negligibly small, such that passive upwelling makes up approximately the entire area of a convection cell, as downwellings in internally heated stagnant lid convection are narrow (e.g., Solomatov and Moresi, 2000). One could also assume that the radius of the cylinder of upwelling mantle is only L/2 ≈ d/2, or half the radius of convection cells. Doing so would introduce a factor of 1/2 in (12), which has a negligible influence on our results (see Section 4.2). Note also that melt production, and hence crustal production and CO₂ outgassing, is continuous and global as long as volcanism is active, as melt is constantly generated by upwelling mantle into a global melt zone. From our thermal evolution models, the depth of melting is always < 200 km, which is above the depth range of ≈250–500 km where melt can become negatively buoyant (Reese et al., 2007). We therefore assume all melt is positively buoyant in our model.

The temperature difference between melt erupted at the surface and the surface temperature is ΔT_m = T_p − P_iγ − T_s. Melt created where upwelling mantle crosses the solidus has a temperature T_p + P_iγ_mantle. This melt cools by adiabatic decompression as it rises to the surface, by an amount P_iγ_melt, where γ_melt is larger than γ_mantle. Subtracting the adiabatic cooling from the initial temperature gives the erupted melt temperature as T_p − P_iγ, where γ = γ_melt − γ_mantle = 6.766 × 10⁻⁴ K km⁻¹ using an average value of γ_melt ≈ 1 K km⁻¹ (Driscoll and Bercovici, 2014). The erupted melt temperature subtracted by the surface temperature, T_s, then gives ΔT_m.

The velocity, v, is calculated from scaling laws for internally heated stagnant lid convection (Reese et al., 1998, 1999; Solomatov and Moresi, 2000; Korenaga, 2009), which give either

or

where c₂ and c₃ are empirical constants determined from mantle convection simulations. We use (13) as the scaling law for velocity in this study; test cases using (14) show no significant differences in terms of mantle degassing lifetime and the conditions for supply-limited weathering (see Section 4.2). Note that mantle viscosity is assumed to be Newtonian in this study, resulting in the above scaling law exponents.

Mantle melting and volcanism produce a crust at the surface that grows over time. If the crustal thickness, δ_c, approaches that of δ as calculated from (2), the crust can begin to influence convective heat flow and velocity in a way that is not captured by the scaling laws. However, crust also undergoes phase transitions as it is buried that influence its stability. In particular, basalt transforms into eclogite, which is dense relative to the underlying mantle, at a depth of ≈50 km on Earth (e.g., Hacker, 1996). The lower crust may therefore founder and sink into the mantle, preventing the crust from becoming so thick that it influences the thickness of the stagnant lid, and thus heat flow and velocity. The depth at which crustal foundering would occur is unknown. The basalt-to-eclogite transition will cause crust to become compositionally unstable at ≈50 km depth; however, with a thick stagnant lid the viscosity at this depth will be large, potentially preventing foundering. Although some studies have looked at stagnant lid convection with a basaltic crust and the eclogite transition (e.g., Johnson et al., 2014), heat flow scaling laws taking crustal foundering into account have not yet been developed. Lacking these scaling laws at the present time, we simply assume that crust delaminates when it reaches the base of the thermal lithosphere as dictated by convective instability, as in this region the viscosity contrast between lithosphere and mantle is low. As a result, we assume that lithospheric thickness follows (2), and thus all crust buried to depths below δ is recycled into the mantle. This assumption has been used in other thermal evolution studies (e.g., Morschhauser et al., 2011).

With the assumption that all the melt produced contributes to crustal growth, the evolution of crustal volume, V_crust, is

The hyperbolic tangent function is a mathematically convenient way to formulate a crustal loss rate that goes to zero when δ_c < δ and equals the crustal growth rate plus the lid thinning rate when δ_c = δ. The term 4π(R − δ)²min(0, dδ/dt) gives the rate at which the volume of stagnant lid is lost when δ is shrinking and is 0 when δ is growing. The crustal thickness, δ_c, is then calculated from the crustal volume as

assuming that crustal thickness is constant spatially (i.e., that melt eruption and crustal production are evenly distributed across the planet).

Incompatible elements, including heat-producing elements, are preferentially incorporated into the melt phase during mantle melting. As a result, the crust becomes enriched in heat-producing elements, and the mantle becomes depleted. We track the evolution of heat-producing elements in the mantle and crust by assuming accumulated fractional melting (e.g., Fraeman and Korenaga, 2010; Morschhauser et al., 2011),

where x_melt is the concentration of heat-producing elements in the melt, x₀ is the initial concentration in the solid, and D is the distribution coefficient (we use D = 0.002 [Hart and Dunn, 1993; Hauri et al., 1994]). The evolution of the crustal heat production rate, Q_crust, is therefore

and mantle heat production evolves as

The rate at which heat production is added to the crust due to melting is given by the product of the concentration of heat-producing elements in the liquid and the volumetric melt production rate, f_m. Associating the concentration of heat-producing elements in the mantle at a given time, x_m, with x₀ from (17), the rate of increase in Q_crust due to melting is given by (x_m f_m /φ)[1 − (1 − φ)^1/D], which is the first term on the right-hand side of (18). The second term on the right-hand side of (18) represents recycling of crustal heat-producing elements into the mantle due to crustal foundering, and the last term is the radiogenic decay of heat-producing elements. A constant decay constant, τ_rad ≈ 2.94 Gyr, is used based on an average of the four major radiogenic isotopes powering Earth's mantle, ²³⁸U, ²³⁵U, ²³²Th, and ⁴⁰K (Driscoll and Bercovici, 2014). The evolution of mantle heat production follows from (19); heat-producing elements are lost due to melting, gained from crustal foundering, and decay with time. All models start with a chosen initial heat production rate, Q₀, that is entirely within the mantle. Q₀ is varied over a range of 5–250 TW in our models. For reference, typical estimates for the Earth span ≈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).

2.3. Carbon cycle

In order to track how the CO₂ degassing rate decays with time, and whether CO₂ degassing overwhelms CO₂ drawdown via weathering, a simple model of carbon cycling between the surface and mantle is employed. Carbon is lost from the mantle during volcanism and deposited in the crust via weathering. As carbonated crust is buried by subsequent lava flows, temperature and pressure conditions are typically high enough for crustal decarbonation to occur (e.g., Kerrick and Connolly, 2001b). This metamorphic devolatilization results in an additional flux of CO₂ to the atmosphere, which also weathers out and is redeposited in the newly forming crust at the surface.

Formulating the carbon cycle model as described in the previous paragraph assumes that the rate of atmospheric CO₂ drawdown via weathering is always equal to the rate of CO₂ input via volcanic and metamorphic degassing, as long as the supply limit to weathering is not exceeded. This assumption is justifiable because the dependence of the weathering rate on surface temperature and atmospheric CO₂ content allows a balance between weathering and degassing to be quickly established (Sundquist, 1991; Berner and Caldeira, 1997). As the weathering rate increases with increasing temperature and partial pressure of atmospheric CO₂, a situation where the degassing flux is larger than the weathering flux will cause the weathering flux to increase until balance is re-established. An analogous process occurs for the situation when the weathering flux is greater than the degassing flux; cooling of the climate will cause the weathering rate to decrease until degassing and weathering come back into balance. Assuming that degassing and weathering are always in balance allows us to use a simple model that only tracks carbon in the mantle and the crust; no model for climate and atmospheric carbon evolution is required. By neglecting the atmosphere and ocean carbon reservoirs, the model used in this study will slightly overestimate the mantle and crustal carbon reservoirs. However, the atmosphere and ocean carbon reservoirs are small, so this overestimation is not significant (e.g., Sleep and Zahnle, 2001).

The depth at which crustal decarbonation occurs determines whether carbonated crust is recycled into the mantle by convective foundering or whether the crust degasses its carbon back to the atmosphere before it founders. The crustal decarbonation depth, δ_carb, is determined based on the topology of a P-T phase diagram section calculated for a carbonate-bearing oceanic metabasalt (e.g., Staudigel et al., 1989; Kerrick and Connolly, 2001b) (see Fig. 2). These calculations were performed using Perple_X 6.6.6 (Connolly, 2005), in combination with the Holland and Powell DS622 thermodynamic database (Holland and Powell, 2011). Details of the mineral activity-composition models used are provided in Appendix A. Decarbonation reactions have positive Clapeyron slopes across the depth interval relevant to our model; slopes become steeper at pressures above ≈30 kbar. The bulk of metamorphic decarbonation (>2.8 to <0.25 wt % CO₂) occurs across a narrow (≈100–200 K) temperature interval, chiefly by the breakdown of dolomite (CaMg(CO₃)₂). Magnesite is the principal host for CO₂ at temperatures less than 1000 K; calcite and aragonite have small stability fields. The effect of these metamorphic phase relations is incorporated into the model by approximating the pressure and temperature dependence of the 0.25 wt % CO₂ contour with a simple linear fit.

The temperature at which decarbonation occurs is then

where A = 3.125 × 10⁻³ K m⁻¹, B = 835.5 K, and T_decarb is temperature in kelvins. Pressure is converted into depth by assuming a lithospheric density of ρ_l = 3300 kg m⁻³. From (5), the temperature profile in the crust is

The height where decarbonation occurs is found by equating (20) and (21). The decarbonation depth, δ_carb, occurs at z = R_p − δ_carb and is therefore given by

When δ_carb < δ, crust decarbonates before being recycled into the mantle, and the evolution of the mantle carbon reservoir, R_man, is given by

where R_man has units of moles (e.g., Tajika and Matsui, 1992; Franck et al., 1999; Sleep and Zahnle, 2001; Driscoll and Bercovici, 2013; Foley, 2015). As crust decarbonates before being recycled into the mantle, the mantle only loses carbon over time when δ_carb < δ. The right-hand side of (23) is defined as the mantle degassing flux,

and gives the rate at which CO₂ is outgassed from mantle to atmosphere, where weathering will deposit it in the crust. CO₂ is incompatible and is therefore preferentially lost from the solid to the melt, analogous to the heat-producing elements. We use a distribution coefficient for CO₂ of  = 10⁻⁴ (e.g., Hauri et al., 2006). As we assume that all the melt produced travels to the surface or near surface, we also assume that all the melt produced contributes to degassing of CO₂ to the atmosphere. Assuming that only a fraction of the melt degasses its CO₂ to the atmosphere, with the rest being trapped in the crust, does not significantly change our overall results (see Appendix B).

The evolution of the crustal carbon reservoir, R_crust, then follows as

Metamorphic decarbonation results in no net change in the crustal carbon reservoir because carbon degassed from the base of the crust is immediately (on geological timescales) redeposited at the top of the crust via weathering. Volcanic degassing of the mantle therefore causes a net transfer of carbon from the mantle to the crust over the lifetime of a stagnant lid planet. Although not important for the evolution of the carbon crustal reservoir, the metamorphic degassing flux is important for climate and is given as

where V_carb is the volume of carbonated crust, V_carb = (4/3)π(R_p³ − (R_p − δ_carb)³), and, as in (15), a hyperbolic tangent function is used as a mathematically convenient way to parameterize a step-function-like change in crustal CO₂ around the decarbonation depth. Here a factor of 1/2 is included, so that the metamorphic degassing flux goes to R_crustf_m/(V_carb) when δ_c > δ_carb and to 0 when δ_c < δ_carb. Metamorphic CO₂ outgassing is assumed to be complete and to occur uniformly across the planet's surface, as we assume that crustal production is spatially uniform. The effect of incomplete metamorphic outgassing is explored in Appendix B, and a case where metamorphic decarbonation is ignored entirely is shown in Section 4.2.

When δ_c < δ_carb, the following equations for R_crust and R_man are used:

In this case, crust buried beneath the stagnant lid base will founder into the mantle, returning surface carbon to the mantle. In practice, δ_carb is usually found to be less than δ_c while volcanism is active, so models typically follow Equations 23–25 rather than Equations 27–28. The sum of the mantle and crustal carbon reservoirs is the total carbon budget, C_tot = R_man + R_crust, which is conserved. All models start with the entire carbon budget residing in the mantle and degassing to the surface over time.

3. The “Hot” Limit: Supply-Limited Weathering

As discussed in Section 1.1, the feedback between silicate weathering and climate that promotes habitable surface conditions only operates if the supply of fresh rock does not limit the weathering rate (e.g., West, 2012; Foley and Driscoll, 2016). To formulate the global supply limit to weathering on a stagnant lid planet, we assume that volcanism is the dominant source of weatherable rock at the surface and therefore write

where is the supply-limited weathering flux (in units of mol yr⁻¹), ε is the eruption efficiency, the fraction of melt produced that erupts at the surface, χ is the fraction of reactable elements in the crust, and is the average molar mass of the crust. Assuming a basaltic crust, χ = 0.3 and  = 0.055 kg mol⁻¹. The above formulation for the supply-limited weathering flux assumes complete carbonation of erupted basalt, which is the ultimate upper bound on silicate weathering. In practice, carbonation can limit the ability of pore water to reach fresh basalt and continue CO₂-rock reaction; pore space can be clogged by formation of carbonated minerals, and reaction products can coat fresh basalt surfaces (e.g., Kump et al., 2000; Kelemen et al., 2011; Beinlich et al., 2012). As a result, the true supply limit to weathering may be less than our estimate based on complete carbonation. In Section 4.2, we discuss the impact a lower supply limit to weathering has on our results.

Supply-limited weathering prevails when F_d + F_meta ≥ . Thus the threshold where supply-limited weathering begins is given by F_d + F_meta = , and the critical mantle and crustal CO₂ contents where weathering is supply limited are given by

While volcanism is active, melt fractions are typically on the order of 0.1. For  ≈ 1, and the mantle degassing flux reduces to f_mR_man/(φ[V_man + V_lid]); this simplified expression for the degassing flux is used in (30). There are two end-member scenarios: one where the planet's carbon is predominantly in the mantle and one where it is in the crust. If R_man >> R_crust, (30) reduces to an expression for the critical mantle CO₂ content for supply-limited weathering: . Assuming an eruption efficiency of 0.1 (Crisp, 1984; White et al., 2006) and that the volume of the crust is negligible such that V_man + V_lid = (4/3)π(R_p³ − R_c³),  ≈ 8 × 10²² to 8 × 10²³ mol for melt fractions ranging from φ  = 0.05 to 0.5. As crustal volume is negligible compared to mantle volume for a crust on the order of 100 km thick, the typical thickness produced by our thermal evolution models, we ignore crustal volume in calculating V_man + V_lid for the remainder of this section (crustal volume is not neglected in the thermal evolution models presented in Section 4).

Sleep and Zahnle (2001) estimate ≈7 × 10²¹ to 2 × 10²² mol of CO₂ in Earth's mantle, so a planet would need to have a mantle CO₂ budget ≈4–110 times larger than Earth's for supply-limited weathering to prevail due to mantle degassing. At melt fractions of φ < 0.05, even smaller mantle CO₂ budgets would lead to supply-limited weathering, as CO₂ becomes increasingly concentrated in the melt due to its incompatible nature. However, melt fractions are typically above 0.05 in our models while volcanism is active; lower melt fractions are only seen just before large-scale volcanism ceases. At this time, outgassing rates are too low and short-lived to lead to any significant atmospheric CO₂ buildup. Thus the threshold for supply-limited weathering due to mantle volcanism in practice is the range of ≈8 × 10²² to 8 × 10²³ mol given above. Metamorphic outgassing, on the other hand, can result in a lower threshold for supply-limited weathering. If carbon is concentrated in the crust (i.e., if R_crust >> R_man): . With δ_carb = 10–100 km,  ≈ 10²² to 10²³ mol. The same amount of carbon results in a higher degassing rate when it is concentrated in the crust versus in the mantle, for melt fractions of order 0.1, due to the much larger volume of the mantle. Thus supply-limited weathering will occur at a lower total planetary CO₂ budget when the CO₂ predominately resides in the crust.

The trade-offs between how carbon is distributed between crust and mantle, the depth of decarbonation, and the melt fraction, φ, can be illustrated by introducing c, the fraction of a planet's total carbon budget, C_tot, that resides in the crust. Thus R_crust = cC_tot, R_man = (1 − c)C_tot, and

where is the critical total planetary CO₂ budget above which weathering will be supply limited. Figure 3 shows that either increasing the fraction of carbon in the crust or decreasing the decarbonation depth lowers . Meanwhile, decreasing φ mutes the influence of on the fraction of carbon in the crust, because at lower melt fractions the threshold for supply-limited weathering due to volcanism is smaller and closer to that from metamorphic outgassing. During planetary evolution, c, δ_carb, and φ all change over time; c will increase as carbon is degassed from the mantle and deposited in the crust, δ_carb is inversely proportional to mantle temperature and will typically grow with time, and φ will decrease with time. A full thermal evolution and carbon cycle model is therefore needed to follow how R_crust, δ_carb, and φ evolve, and determine whether or not supply-limited weathering prevails; this is done in Section 4.

However, simply checking if the total degassing flux (i.e., the sum of the metamorphic and volcanic degassing fluxes) ever exceeds the supply-limited weathering flux can give misleading results about whether inhospitably hot climates will develop in practice. In particular, at low melt fractions mantle degassing can exceed the weathering supply limit as a result of the high concentration of CO₂ in the melt. In such situations, though, the total outgassing rate is low, on the order of 10¹⁰ to 10¹¹ mol yr⁻¹ or less in our models, and it would take billions of years for a thick CO₂ atmosphere to form. Mantle volcanism with a low melt fraction does not persist for such long timescales, because it only occurs when volcanism is close to ceasing entirely. We thus define the threshold for supply-limited weathering to form an uninhabitable hothouse climate as F_d + F_meta ≥  + F_offset. We set F_offset = 10¹⁴ mol yr⁻¹, as a net CO₂ input to the atmosphere of 10¹⁴ mol yr⁻¹ will cause 100 bar of CO₂ (≈ 10²² mol) to accumulate in 100 Myr. The boundary for the onset of supply-limited weathering shown in our model results below (Section 4) was found to be insensitive to F_offset over the range F_offset = 10¹³ to 10¹⁵ mol yr⁻¹. When supply-limited weathering prevails, either as a result of mantle or metamorphic outgassing, it typically persists for the lifetime of volcanism. In principle as the mantle cools and δ_carb grows, weathering could transition back out of the supply-limited regime; however, this is not seen in our models, as δ_carb typically stays low (order of 10 km) during the lifetime of volcanism.

4. The “Cold” Limit: Longevity of CO₂ Degassing

As discussed in Section 1.1, the other requirement for habitability we consider is whether degassing rates are high enough to prevent global surface glaciation; this threshold is calculated using climate models and by assuming that atmospheric CO₂ is set by a balance between degassing and weathering (e.g., Kadoya and Tajika, 2014; Menou, 2015; Batalha et al., 2016). For an Earth-like planet, Kadoya and Tajika (2014) estimate that a degassing flux approximately equal to Earth's present-day degassing flux is needed to keep the surface temperature above freezing for solar luminosities of ≈0.7–0.8S* (where S* is the present-day solar luminosity), while a degassing flux of ∼10% the present day is sufficient for solar luminosities of ≈S*. Haqq-Misra et al. (2016) find that a degassing flux of ∼10% present-day Earth's allows for warm climates even for Archean solar luminosities. To be conservative, we take the estimates from Kadoya and Tajika (2014) as reasonable limits on the degassing flux necessary for preventing a snowball climate, and calculate how long CO₂ degassing rates can be maintained at both 10% and 100% of Earth's present-day total degassing flux of ≈6 × 10¹² mol yr⁻¹ (Marty and Tolstikhin, 1998).

Snowball climate states triggered by low degassing rates may not be permanent. It is possible that the climate will oscillate between globally frozen and ice-free states, as weathering shuts down during a snowball phase. With weathering shut down, volcanism replenishes atmospheric CO₂, triggering melting and a short-lived warm climate phase. During the ice-free warm climate phase, precipitation and weathering ensue, plunging the climate back into a snowball state (e.g., Menou, 2015). However, the period of time spent under glaciated conditions is much longer than the time spent at warm conditions during these oscillations (Kadoya and Tajika, 2014; Menou, 2015; Batalha et al., 2016). For simplicity, we will consider that snowball climates prevail below the degassing thresholds set in the previous paragraph, even though warm climates may exist for brief durations.

4.1. Thermal evolution model results

A typical model evolution is shown in Fig. 4. The model uses baseline parameter values (see Table 1), an initial mantle temperature of T_p0 = 2000 K, initial radiogenic heating rate of Q₀ = 100 TW, and total CO₂ content of C_tot = 10²² mol. The initial conditions result in a low convective heat loss, less than the heat production rate, and high magmatic heat loss, far in excess of internal heat production. As a result, the mantle rapidly cools until magmatic and convective heat loss approximately balance mantle heat production. Note that without magmatic heat loss, the mantle would have initially warmed until convective heat loss became large enough to balance mantle heat production on its own. Including magmatic heat loss therefore results in lower mantle temperatures because it increases total heat loss. Also during this early phase of magmatic-dominated cooling, a thick crust forms and becomes enriched in heat-producing elements; Q_mantle falls from 100 TW to ≈50 TW within ≈200 Myr, primarily due to partitioning of heat-producing elements into the crust. Eventually the radiogenic heating rate decays to the point where convective heat loss on its own can balance heat production. As a result, the mantle cools more rapidly, and volcanism quickly ceases. Convective heat loss then tracks mantle radiogenic heating for the remainder of the evolution, and the more rapid mantle cooling rate is sustained. The CO₂ degassing flux is high throughout the lifetime of volcanism due primarily to metamorphic degassing of the carbonated crust. Metamorphic degassing dominates because CO₂ is more concentrated in the crust than in the mantle and because the mantle CO₂ reservoir is quickly depleted by volcanism and outgassing. CO₂ degassing is larger than Earth's present-day degassing rate for ≈2 Gyr, nearly the entire lifetime of volcanism. Transferring carbon from the mantle to the crust via volcanism, and continually burying this crust through the metamorphic decarbonation depth, is thus an effective means of sustaining high CO₂ degassing rates.

The planetary CO₂ budget and radiogenic heating rate are the two most important factors governing how long substantial degassing rates can be maintained on a stagnant lid planet. High radiogenic heating budgets promote long-lived volcanism and crustal burial through the decarbonation depth, while a large initial CO₂ content allows high degassing rates to be maintained, primarily by metamorphic CO₂ degassing, even as volcanism wanes (i.e., with a larger CO₂ budget, more CO₂ is deposited in the crust, such that degassing rates remain high even at low rates of melt production and crustal burial). The influence of C_tot and Q₀ on the longevity of CO₂ degassing is summarized in Fig. 5A, which shows the time when the total degassing flux drops below Earth's present-day degassing flux, and Fig. 5B, which shows when the total degassing flux drops below 10% of Earth's present-day degassing flux. As expected, both a larger radiogenic heating rate and a larger planetary CO₂ budget lead to longer-lived CO₂ degassing. The lifetime of volcanism is not a function of planetary CO₂ budget, and thus contours for the time when volcanism ends would form horizontal lines in Fig. 5A and 5B. Moving to large planetary CO₂ budgets (i.e., toward the right edge of Fig. 5A and 5B), the contours asymptotically approach horizontal lines because the end of volcanism, rather than the supply of CO₂, becomes the dominant factor determining the longevity of adequate degassing. Moving toward lower CO₂ budgets, the contours curve toward vertical lines, as the planetary CO₂ content becomes the primary factor controlling when degassing rates above our designated thresholds end. Thus, even with a high radiogenic heating budget, a planet with a low CO₂ inventory will quickly become globally frozen, even though volcanism can continue for billions of years.

Although large CO₂ budgets are beneficial for promoting warm climates, there is an upper limit on the CO₂ inventory for habitability provided by the transition to supply-limited weathering discussed in Section 3. The critical planetary CO₂ inventory where weathering becomes supply-limited was calculated from the thermal evolution models as explained in Section 3; this is found to occur when C_tot ≈ 10²³ mol at Q₀ = 5 TW and C_tot ≈ 2 × 10²² mol at Q = 250 TW. The relationship between and Q₀ is approximately linear. Larger heat production rates shift the boundary where weathering becomes supply-limited to lower values of C_tot (i.e., making it easier for a planet to experience supply-limited weathering) because higher heating rates cause more volcanism, and thus a higher concentration of planetary carbon in the crust, and higher temperatures, which lower the decarbonation depth (see Section 3).

Thus a large swath of the most favorable region of parameter space for long-lived CO₂ degassing would lead to inhospitably hot climates due to supply-limited weathering and is therefore uninhabitable (Fig. 5A and 5B). Nevertheless, there is still a sizeable region of parameter space, at Q₀ > 100 TW and C_tot ∼ 10²⁰ to 10²² mol (corresponding to ≈7 × 10⁻⁵ to 7 × 10⁻³ mass % CO₂), where habitable climates can potentially be maintained for up to 2–4 Gyr. Earth's total CO₂ budget is not well constrained, but estimates range from ∼10²² to 10²³ mol (Marty and Tolstikhin, 1998; Sleep and Zahnle, 2001; Dasgupta, 2013), while the radiogenic heating budget is estimated at ≈70–100 TW of Hadean heat production (e.g., Korenaga, 2006). A planet with an Earth-like CO₂ and radiogenic heating budget would therefore be able to sustain a temperate climate for up to 1–2 Gyr. For longer maintenance of habitable conditions, planets with Earth's size and composition would need larger radiogenic heating budgets than Earth and total CO₂ budgets similar to Earth's; larger CO₂ budgets would cause supply-limited weathering.

4.2. Sensitivity of results to parameters and model formulation

Two key parameters that are not well constrained are μ_r, the reference viscosity, and T_p0, the initial mantle potential temperature. Varying these parameters over a reasonable uncertainty range leads to some noticeable differences but does not significantly change the overall results (Fig. 6). A lower μ_r (μ_r is changed by changing the pre-exponential constant, μ_n) shortens the lifetime of degassing, because the mantle convects more vigorously and cools more rapidly (Fig. 6A); planets with initial radiogenic heating budgets less than 100 TW see significant degassing and volcanism end in less than 1 Gyr, while even at Q₀ = 250 TW degassing wanes in 3–4 Gyr. Moreover, a lower μ_r also shifts the onset of supply-limited weathering to lower C_tot, because more vigorous convection leads to a larger concentration of carbon in the crust, and a thinner lithosphere, and hence shallower decarbonation depth. Thus planets with low-viscosity mantles are less likely to sustain habitable climates. A larger μ_r slows mantle cooling and allows degassing to last longer (Fig. 6B) but does not significantly expand the range of C_tot–Q₀ parameter space where long-lived degassing is possible as compared to the baseline model (Fig. 5). Lowering the initial mantle temperature suppresses volcanism at low radiogenic heating rates below ≈60 TW and therefore shrinks the region where habitable conditions are possible (Fig. 6C). However, for Q₀  > 100 TW there is no significant change to the degassing lifetime. Finally, a larger initial mantle temperature does not change the calculated degassing longevities but does cause the onset of supply-limited weathering to occur at lower C_tot, as a hotter mantle leads to more outgassing and therefore concentration of carbon in the crust, and a lower decarbonation depth (Fig. 6D).

Another key parameter is the eruption efficiency, ε, or the fraction of melt produced that erupts at the surface. The eruption efficiency dictates the rate at which weatherable rock is supplied to the surface and is therefore crucial for determining whether weathering becomes supply-limited or not. The boundary for the onset of supply-limited weathering scales linearly with ε. Thus setting ε = 1 shifts the transition to supply-limited weathering to values of C_tot one order of magnitude higher than those shown in Figs. 5–7, where ε  = 0.1, and gives the most optimistic possible estimate for when supply-limited weathering would prevail on a stagnant lid planet. Choosing lower ε would shift the weathering regime boundary to lower C_tot. Moreover, as discussed in Section 3, the supply limit to weathering could be lower than we estimate if erupted basalt can not be completely carbonated. Incomplete carbonation would also shift the supply-limited weathering regime to lower values of C_tot, in the same manner as lowering ε. If only 1/100 or less of erupted basalt can be carbonated, then supply-limited weathering will prevail over nearly the entire range of C_tot and Q₀ where sufficient degassing rates take place; very inefficient crustal carbonation could thus prevent habitable climates on stagnant lid planets entirely.

Figure 7 shows the results of changing additional key aspects of the model. Changing the convective velocity scaling law from Equation 13 to Equation 14 (i.e., from a scaling exponent of 2/3 to 1/2) causes degassing lifetimes to be ≈500 Myr longer overall than our baseline model results where (13) is used (Fig. 7A). The reason degassing and volcanism lasts longer with a Ra_i^1/2 scaling law is that convective velocities are lower, resulting in lower melt production rates and thus less depletion of heat-producing elements in the mantle. However, the difference between the two velocity scaling laws on degassing lifetime and the transition to supply-limited weathering is minor compared to other sources of uncertainty in the modeling. We also test the effect of reducing the melt production rate, f_m, by a factor of 2 as compared to our standard formulation given in (12), as explained in Section 2.2. The results are nearly identical to our baseline model and actually show longer-lived degassing as a result of greater retention of heat-producing elements in the mantle (Fig. 7C).

We also test the effect of the distribution coefficient, D, which controls the degree to which heat-producing elements are preferentially partitioned into the crust during melting, and the activation energy for viscosity, E_v. The upper limit for the distribution coefficient is D ≈ 0.01 for a mantle dominated by pyroxene (Hart and Dunn, 1993; Hauri et al., 1994). Using D = 0.01 produces no discernible differences in the lifetime of CO₂ degassing (Fig. 7D). Distribution coefficients of D  >∼ 0.1–1 are needed to significantly impede enrichment of heat-producing elements in the crust, and produce longer degassing lifetimes than calculated in our baseline models. Varying the activation energy for viscosity by 100 kJ mol⁻¹ above and below our baseline value does not change the region of C_tot–Q₀ space where potentially habitable climates are possible but does change the lifetime of sufficient CO₂ degassing; with E_v = 200 kJ mol⁻¹, degassing lasts longer than in the baseline model, while with E_v = 400 kJ⁻¹ mol it ends sooner. The stronger the dependence of viscosity on temperature, the more lid thickness grows with decreasing mantle temperature. Thus a higher activation energy causes the lid to grow more rapidly as the mantle cools and shut off volcanism, and outgassing, sooner. Note that when we change activation energy, the reference viscosity, μ_r, is held constant; this is accomplished by changing μ_n when E_v is changed. For E_v = 200 kJ mol⁻¹, μ_n = 7.3 × 10¹³ Pa s, and for E_v = 400 kJ mol⁻¹, μ_n = 2.7 × 10⁷ Pa s.

Finally, we test how our results are affected if metamorphic decarbonation does not occur. Considering a case where crust does not decarbonate as it is buried is useful for two main reasons: (1) Exoplanets can have a wide range of chemical compositions, and crustal decarbonation may not occur at the temperature and pressure conditions prevalent in the crust on all planets. We note, however, that our model assumes an olivine-dominated mantle and basaltic crust and is thus not directly applicable to planets with different mineralogies. (2) Rapid burial of the crust could suppress temperatures in the lid by advecting cold rocks from the surface downward faster than thermal diffusion can heat them to the steady-state geotherm. For melt production rates typically seen in our models, advection would only have a minor influence on the temperature profile in the crust as Pe < 1 (Section 2.1). However, if volcanism is confined to only a handful of places on the surface, local crustal burial rates would be higher, and crust could potentially founder into the mantle before it heats enough to undergo decarbonation. Figure 7B shows that the range of C_tot that can support long-lived degassing shifts to higher values when crustal decarbonation is no longer active. However, the supply limit to weathering also shifts to larger C_tot, meaning that the size of the region of parameter space where long-lived habitable climates are possible is approximately the same, regardless of whether metamorphic decarbonation occurs or not. Without metamorphic decarbonation, the onset of supply-limited weathering shifts to higher values because the planet's CO₂ supply largely remains in the mantle, and the large volume of the mantle means this carbon is diluted. Thus, the same C_tot results in an overall lower degassing rate when only mantle volcanism is active than when both mantle volcanism and crustal decarbonation can take place.

4.3. Erupted melt fluxes and habitability

Although active volcanism is beneficial for habitability to the extent that it promotes a temperate climate, extreme volcanic resurfacing rates could prevent life from developing by rapidly sterilizing the surface. The resurfacing rate at which life would be detrimentally impacted by lava flows is not well known. Despite this uncertainty, we can at least compare our eruption rates to those seen on Earth. Eruption rates are calculated from our thermal evolution models as εf_m with ε = 0.1 (Fig. 8). Eruption rates are typically < 100 km³ yr⁻¹ for most of the period while volcanism is active, which results in ≈0.2 mm yr⁻¹ of resurfacing when erupted lava is assumed to be distributed globally. For comparison, resurfacing rates at Kilauea Volcano in Hawaii are ≈5 mm yr⁻¹ (e.g., Quane et al., 2000), and life is present across the volcano. Thus eruptions are unlikely to be a major impediment to life on stagnant lid planets for most of the period when volcanism is active. Initial volcanic rates are often very high (on the order of 10⁴ km³ yr⁻¹), and potentially uninhabitable, in our models, but such high volcanic rates are short lived (lasting on the order of 100 Myr).

5. Discussion

5.1. Implications for exoplanet observations and target selection

The model results indicate that Earth-sized stagnant lid planets with extrapolated initial radiogenic heating rates greater than Earth's, and with total CO₂ budgets of ∼10²⁰ to 10²² mol (≈7 × 10⁻⁵ to 7 × 10⁻³ mass % CO₂), are potentially habitable for timescales of over 2 Gyr. The time when habitability becomes unlikely due to limited CO₂ outgassing scales with initial radiogenic heat budget, with higher budgets leading to longer-lived outgassing. Thus planet heat budget, CO₂ budget, and age are critical factors for habitability of stagnant lid planets. Whether planets falling within these limits are common or not is not well constrained, in particular regarding CO₂ inventories. Accreted CO₂ inventories are likely strongly dependent on the size and composition of the proto-planetary disk and the ability of volatile-rich planetesimals to be scattered inward toward the terrestrial planet-formation zone, likely by interaction with giant planets like Jupiter (Raymond et al., 2004). When volatile-rich planetesimals are scattered inward, significant radial mixing occurs, with water and other volatiles delivered to terrestrial planets in a highly stochastic manner. In models of solar system formation, Raymond et al. (2004) found planets formed with water concentrations ranging from 2 orders of magnitude lower than Earth's to 2 orders of magnitude larger, with the majority of planets falling within a range of ∼1–100 ocean masses of water (i.e., ∼0.1–1 to ∼10–100 times Earth's total water abundance). CO₂ delivery is likely similar, as it would also be caused by inward scattering of volatile-rich planetesimals, meaning stochastic processes probably lead to a similar range of final planetary CO₂ inventories. Ultimately, it is difficult to predict or constrain from stellar observations which exoplanet systems are likely to have planets with bulk CO₂ budgets within the optimum range for habitability indicated by our models.

Radiogenic heating budget and age, on the other hand, can be constrained by stellar observations, and both stars with sufficiently high radionuclide abundances and those young enough to potentially host habitable stagnant lid planets are apparently common. Unterborn et al. (2015) showed that the thorium abundance of Sun-like stars varies from 59% to 251% the solar value, for stars with abundances of major rock-forming elements that would produce planet mineralogies similar to Earth. These thorium variations would lead to radiogenic heating budgets ranging from ≈0.5–2.5 times Earth's. Moreover, most stars measured in Unterborn et al. (2015) had higher thorium abundances than the Sun. Planet age can be inferred from star age, measurable by asteroseismology (e.g., Silva Aguirre et al., 2015; Burgasser and Mamajek, 2017), and results from Kepler stars show many planet-hosting stars with ages of 2–4 Gyr (Silva Aguirre et al., 2015). Thus, the likelihood of planets meeting the CO₂, radiogenic heating budget, and age constraints for maintaining habitable climates in the stagnant lid regime appears to be relatively high. As a result, plate tectonics may not be required for habitability, although as we discuss in Section 5.2, it still has some major benefits compared to stagnant lid convection. Furthermore, in order to avoid a fate similar to that of Venus, stagnant lid planets would still need to lie within the habitable zone, and in order to avoid becoming barren, frozen worlds like Mars, would need to be large enough for the planet's gravity to prevent rapid escape of atmospheric greenhouse gases.

Our results can also naturally be used to guide target selection for future missions. In general, the model results indicate that young planets and/or those that orbit stars with high U and Th abundances are good targets for searching for biosignatures, as these planets are potentially habitable even in a stagnant lid regime. The modeling presented here could even be used in future applications to directly constrain the prospects for habitability of individual planets as they are detected. Whether high heat production rates also enhance a planet's prospects for plate tectonics is not clear, however, as how mantle temperature and heat production rate influence a planet's propensity for plate tectonics is debated (O'Neill et al., 2007; Korenaga, 2010a; Foley et al., 2014; Foley and Driscoll, 2016). Another potential way to use the results for target selection is to estimate the amount of tidal heating that would be expected on exoplanets based on their orbits. Planets experiencing significant tidal heating will be able to maintain long-lived CO₂ degassing, even if the abundance of heat-producing elements in the mantle is low. Tidal heating is potentially important for rocky planets in the habitable zone of M dwarfs (Jackson et al., 2008), or planets or moons undergoing orbital interactions with other planets, similar to Io (Moore, 2003). In particular, stagnant lid planets orbiting M dwarfs at the outer edge of the habitable zone can see significant tidal heating last for ∼10 Gyr, before their orbits circularize (Driscoll and Barnes, 2015); circularization is much faster for planets near the inner edge of the habitable zone.

5.2. Comparison to plate-tectonic planets

On a plate-tectonic planet, CO₂ degassing is easier to maintain. Ridges bring hot mantle all the way to the surface, meaning that volcanism at divergent plate margins will continue until the mantle potential temperature drops below the zero pressure solidus of ≈1400 K (Herzberg et al., 2000). Plate tectonics also allows for continual metamorphic degassing through subduction. Slab carbon loss is most likely driven by dissolution in fluids released by metamorphic dehydration reactions (Ague and Nicolescu, 2014), as temperature and pressure conditions along modern slab top geotherms do not cause significant CO₂ release (Kerrick and Connolly, 2001a, 2001b). Thus, arc degassing will continue as long as significant slab dewatering and fluid-flux melting of the mantle wedge occur. Hydrous wedge melting takes place at low temperatures of ≈900–1100 K (Hacker et al., 2003). For mantles cooler than this, all slab carbon will be recycled into the mantle; but without slab or ridge degassing, CO₂ release to the atmosphere will be negligible, and the surface will likely be frozen. How long it takes the mantle of a plate-tectonic planet to cool below these temperature thresholds depends on the details of an individual planet's evolution. However, Kite et al. (2009) provide a reasonable estimate based on simple models that planetary mantles will stay warm enough for ridge degassing for over 10 Gyr. Slab degassing can be maintained at lower temperatures and would thus last even longer.

Plate tectonics also allows a stabilizing weathering feedback, and thus temperate climates, to be maintained at larger total CO₂ budgets than stagnant lid planets when the area of exposed land is >∼ 1% of the plate tectonic planet's surface (Foley, 2015). A large land area, combined with the high erosion rates typical of mountainous and tectonically active regions, produces a high supply limit to weathering, which requires a very large total degassing flux in order to be reached. Thus, a plate-tectonic planet can potentially maintain a temperate climate at CO₂ inventories up to ∼10²⁴ to 10²⁵ mol, much larger than the ∼10²³ mol upper limit for stagnant lid planets when crustal decarbonation is active. With little or no land, the CO₂ budget where weathering becomes supply-limited on a plate-tectonic planet is approximately equal to the estimates shown here for stagnant lid planets with crustal decarbonation.

5.3. Uncertainties in model formulation

Although some differences can be seen, the overall results, in particular the size of the region of C_tot–Q₀ parameter space where habitable climates are possible, are consistent over a wide range of parameters and to changes in the formulation for convective velocity, metamorphic decarbonation of the crust, melt transport, and volcanic outgassing. Nevertheless, there are effects not included in this study that could influence the results. The most obvious is size, as many rocky exoplanets more massive than Earth have been discovered (e.g., Batalha, 2014). Larger planets have a smaller surface-area-to-volume ratio and thus retain their heat for longer, all else being equal. Therefore, increasing planet size is expected to prolong mantle melting and volcanism. However, melt will be formed at higher pressures due to the larger pressure gradient on a massive planet and could end up being too dense to erupt. Noack et al. (2017) argue that planets larger than ≈2–3 Earth masses will experience no outgassing because dense melt will prevent eruption. Future studies of melting and melt density at super-Earth conditions are needed to confirm this result, but it is possible that there is a critical planet size above which a stagnant lid planet would not be habitable.

Considering only Earth-sized planets, there are additional important uncertainties to highlight. Our assumptions about melt heat loss maximize this effect and thus hasten the end of volcanism. A more realistic model, particularly regarding the degree to which melt solidifying at shallow depth contributes to net mantle cooling, may allow degassing to be sustained longer, and over a wider range of parameter space, than presented here. We also ignore melt liberation of CO₂ locked in the crust, which would enhance CO₂ outgassing during mantle volcanism. Melt liberation of crustal carbon is a potentially significant CO₂ source on Earth (Lee et al., 2013), and including it would allow habitable conditions at lower total CO₂ budgets. Thus, the true prospects for habitability on stagnant lid planets could be better than our results indicate. Core cooling would also lead to the formation of plumes, which could sustain volcanism even after pressure release melting due to passive mantle upwelling has ceased. However, plumes alone may not be a sufficient source of volcanism to sustain a temperate climate. Plumes on the modern Earth result in ∼1–3 km³ yr⁻¹ of melt flux (Crisp, 1984), which is similar to the melt fluxes seen in our models when volcanism is close to ceasing and CO₂ outgassing levels are too low to prevent glaciation. The effect of mantle volatile content on viscosity and the solidus is also not included. The viscosity of wet olivine is 1–2 orders of magnitude lower than dry olivine (Hirth and Kohlstedt, 2003), meaning that the mantle could become more viscous over time as water is outgassed to the surface. At the same time, the solidus will shift from lower to higher temperatures as mantle water content declines (e.g., Katz et al., 2003). As a result, the mantle cooling rate would decline as water is expelled from the mantle, acting to keep the mantle hotter, but the solidus temperature will also increase. New models would be needed to determine which of these competing effects of mantle hydration state on CO₂ outgassing are dominant.

There are also factors that could make habitability less likely than our results indicate. First, as discussed in Sections 3 and 4.2, we assume a supply limit to weathering based on complete carbonation of basalt. In reality, the amount of basalt that can be carbonated will be controlled by the ability of pore fluids to reach fresh mineral surfaces. The competition between reaction-driven cracking, acting to open up new pore space, versus permeability loss due to mineral precipitation, will determine the amount of basalt that can be carbonated (e.g., Rudge et al., 2010). As a result, only a fraction of surface basalt may actually be available for weathering. If only ∼1/100 of surface basalt can be carbonated, then nearly all the parameter space where long-lived degassing is possible will result in supply-limited weathering, and habitability of stagnant lid planets will be unlikely. Metamorphic degassing may also be less significant than we model, or even entirely absent, if eruption is concentrated in a few localized areas such that crustal burial at these sites is rapid. However, we show that even without metamorphic decarbonation there is still a sizable region of C_tot–Q₀ parameter space where habitable climates are possible (see Section 4.2).

We also assume complete outgassing of all melt produced to the atmosphere. However, a significant fraction of melt generated could stall in the crust or stagnant lid and not outgas to the atmosphere, potentially making it more difficult than our models predict to sustain degassing rates high enough to avoid a snowball climate. Simple test cases designed to capture the effect of melt stalling and incomplete degassing are presented in Appendix B and show no significant difference from our models where all the melt degasses. Finally, we assume that volcanism and outgassing are continuous in time, as melt is constantly being created by mantle upwelling into the melting region. Continuous volcanism is important because pauses in CO₂ outgassing lasting longer than ∼10⁵ yr could trigger snowball states (Tajika, 2008). While melt is continuously being generated in the mantle in our models, it is possible that melt transport through the lid could induce episodicity in eruptions that our simple thermal evolution model cannot capture. Ultimately, future work involving convection models and models of melt transport through a stagnant lid will be needed to test whether significant pauses in stagnant lid planet volcanism are likely, and the effect they would have on habitability.

6. Conclusions

Models of the thermal, magmatic, and degassing history of rocky planets with Earth-like size and composition demonstrate that a carbon cycle capable of regulating atmospheric CO₂ content, and stabilizing climate to temperate surface temperatures, can potentially operate on geological timescales on planets in the stagnant lid regime. Plate tectonics may not be required for habitability, at least in regard to sustaining a stable, temperate climate on a planet. Specifically, stagnant lid planets with radiogenic heating budgets similar to or larger than Earth's, and with total CO₂ budgets of ∼10²⁰ to 10²² mol (ranging from 10⁻² to 1 times typical estimates for Earth's CO₂ budget), can potentially maintain temperate climates for 1–5 Gyr. After this time, volcanism and CO₂ outgassing largely cease, and climate would likely cool below the water freezing point and induce global surface glaciation (though plumes, which are neglected in our modeling, could prolong volcanism and temperate climates). At CO₂ budgets lower than ∼10²⁰ mol, a planet's climate is estimated to be in a snowball state for its entire history, while above ∼10²² mol weathering would become supply-limited. With supply-limited weathering, CO₂ outgassing overwhelms CO₂ drawdown, such that an inhospitably hot, CO₂-rich atmosphere forms. Thus, the amount of carbon accreted to a planet during formation is critical for whether it can sustain habitable surface conditions in a stagnant lid regime. Our results also indicate that high internal heating rates favor long-term habitability by prolonging volcanism and CO₂ degassing, and that young planets are more likely to experience sufficiently high rates of CO₂ outgassing than older planets. Thus, planets orbiting high thorium or uranium abundance stars or young stars are good targets for searching for biosignatures, because they are more likely to be habitable, even if the mode of surface tectonics is stagnant lid rather than plate tectonics.

A. Appendix: Details of Modeling Phase Equilibria

The following solution models were used in the computation of phase relations (Fig. 2); all other phases considered were treated as ideal:

• Chlorite: Chl(HP)—Holland et al. (1998)

• Biotite: Bio(HP)—Powell and Holland (1999)

• Garnet: Gt(HP)—Holland and Powell (1998)

• Amphibole: GlTrTsMg, GlTrTsPg—Wei and Powell (2003), White et al. (2001, 2003)

• Clinopyroxene: Cpx(HP)—Holland and Powell (1996)

• Epidote: Ep(HP)—Holland and Powell (1998)

• CO₂-H₂O fluid: Holland and Powell (1991)

• Magnesite: M(HP)—Holland and Powell (1998)

• Orthopyroxene: Opx(HP)—Holland and Powell (1996)

• Dolomite: Do(HP)—Holland and Powell (1998)

• Phengite: Pheng(HP)—Holland and Powell (1998)

• Spinel: Sp(HP)—Holland and Powell (1998)

B. Appendix: Effect of Melt Stalling and Incomplete Degassing

In the thermal evolution and outgassing models presented in the main text, all the melt produced is assumed to percolate to the surface or near surface. As such, all the CO₂ in the melt is assumed to degas to the atmosphere. Moreover, all CO₂ liberated from the crust by metamorphic reactions is also assumed to degas to the atmosphere. In this appendix, we present test cases where degassing of CO₂ from the melt and from metamorphic decarbonation is incomplete and where melt is allowed to stall at the base of the stagnant lid and not contribute to surface volcanism and outgassing. These tests assess the sensitivity of our results to the assumption of complete melt percolation to the surface and complete outgassing of CO₂ from the melt and metamorphic reactions.

We first consider a case where all the melt is still assumed to rapidly rise to the surface, but only a fraction of the CO₂ in this melt degasses; the rest is trapped within the crust and becomes part of the crustal CO₂ reservoir. As a result, the volcanic degassing flux, F_d, is reduced by a factor ζ₁:

We also allow for a fraction, ζ₂, of the CO₂ released by metamorphic reactions to remain trapped in the crust and eventually recycle into the mantle. As a result, the metamorphic degassing flux is

and the evolution of the mantle and crustal carbon reservoirs follow

when δ_carb < δ_c, and (28) and (27) otherwise. Note that transfer of CO₂ from the mantle to the crust is the same as in the main text, as all CO₂ that is stripped from the mantle by melting is assumed to be deposited in the crust regardless of whether it outgases into the atmosphere first or is trapped in the crust during solidification. The rest of the thermal evolution and carbon cycle model is the same as presented in the main text.

The results of test cases with ζ₁ = 0.1 and ζ₂ = 0 and with ζ₁ = 0.01 and ζ₂ = 0 are identical to the results of the baseline model in the main text (Fig. B1A and B1B). As it is the metamorphic decarbonation that dominates outgassing for most of the lifetime of volcanism in our models, incomplete degassing of the melt does not significantly change the results. The results are changed when ζ₂ > 0 is used; specifically as ζ₂ approaches 1, model results converge to the case where crustal decarbonation is ignored (i.e.,Fig. 7B in the main text), as ζ₂ = 1 implies no metamorphic CO₂ outgassing and hence all the crustal carbon being recycled into the mantle by foundering. As a result, models with ζ₁ = 0.1 and ζ₂ = 0.9 (i.e., where only 10% of the carbon that is metamorphically lost from the crust degasses to the atmosphere) produce very similar degassing lifetimes to the no-crustal-decarbonation case in the main text (Fig. B1C). However, because volcanic outgassing is also assumed to be incomplete, overall degassing rates are lower, and both the transition to supply-limited weathering and the range of total CO₂ budgets where long-lived degassing occurs shift to higher values. Thus, the range of C_tot where habitable climates are possible is approximately the same size as in the models presented in the main text, but shifted to higher values as a result of incomplete outgassing of both melt and CO₂ released by metamorphism.

Another important scenario to test is one where some fraction of melt produced in the mantle cannot rise through the stagnant lid and reach the surface or near surface. We therefore perform models where only a fraction, ξ, of the melt produced reaches the surface and contributes to mantle cooling, crustal growth, and CO₂ outgassing. The rest is assumed to be trapped at the base of the lid and therefore remains in the mantle. The model is thus modified as follows. The mantle energy balance becomes

Heat-producing elements in the crust and mantle evolve as

and the mantle and crustal carbon reservoirs follow

when δ_carb < δ_c, and

when δ_carb > δ_c.

With melt stalling at the base of the lid, the model results are different than those presented in the main text in detail, but the range of C_tot and Q₀ that allows for potentially habitable climates is nearly identical (Figs. B1D and B1E). In fact, the main difference when melt is assumed to stall at the base of the lid is that degassing lifetimes are everywhere longer; the mantle cools more slowly with less efficient melt heat loss and fewer radionuclides partitioning into the crust. When metamorphic decarbonation is also assumed to be limited by setting ζ₂  = 0.9, we obtain results similar to the case in the main text where decarbonation is excluded (Fig. B1F). As explained above, assuming that the majority of metamorphically liberated CO₂ is unable to outgas is effectively the same as assuming the crust does not undergo decarbonation. Again, degassing lifetimes are found to be longer with melt stalling because the mantle loses less heat from volcanism and fewer heat-producing elements to the crust. Ultimately, our test cases demonstrate that the results in the main text are robust to the effects of incomplete outgassing of CO₂ and melt stalling in the lower lithosphere, at least as far as these affects can be incorporated into our simple thermal evolution models. Ultimately, more sophisticated modeling involving two-phase flow of melt and CO₂ vapor through the solid lithosphere is needed to more thoroughly assess how much melt produced in the mantle is actually able to degas at the surface, and the effect this has on long-term climate evolution.

Acknowledgments

We thank Ed Kite, Lena Noack, and two anonymous reviewers whose comments helped us significantly improve the manuscript.

Author Disclosure Statement

No competing financial interests exist.