Carbon Cycling and Habitability of Earth-Sized Stagnant Lid Planets
Models of thermal evolution, crustal production, and CO2 cycling are used to constrain the prospects for habitability of rocky planets, with Earth-like size and composition, in the stagnant lid regime. Specifically, we determine the conditions under which such planets can maintain rates of CO2 degassing large enough to prevent global surface glaciation but small enough so as not to exceed the upper limit on weathering rates provided by the supply of fresh rock, a situation which would lead to runaway atmospheric CO2 accumulation and an inhospitably hot climate. The models show that stagnant lid planets with initial radiogenic heating rates of 100–250 TW, and with total CO2 budgets ranging from ∼10−2 to 1 times Earth's estimated CO2 budget, can maintain volcanic outgassing rates suitable for habitability for ≈1–5 Gyr; larger CO2 budgets result in uninhabitably hot climates, while smaller budgets result in global glaciation. High radiogenic heat production rates favor habitability by sustaining volcanism and CO2 outgassing longer. Thus, the results suggest that plate tectonics may not be required for establishing a long-term carbon cycle and maintaining a stable, habitable climate. The model is necessarily highly simplified, as the uncertainties with exoplanet thermal evolution and outgassing are large. Nevertheless, the results provide some first-order guidance for future exoplanet missions, by predicting the age at which habitability becomes unlikely for a stagnant lid planet as a function of initial radiogenic heat budget. This prediction is powerful because both planet heat budget and age can potentially be constrained from stellar observations. Key Words: Exoplanets—Habitability—Stagnant lid tectonics—Carbon cycle—Volcanism. Astrobiology 18, 873–896.
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 CO2 content increase, rates of chemical weathering increase in turn, thereby drawing CO2 out of the atmosphere and acting to cool the climate. Likewise, when temperature and atmospheric CO2 content decrease, weathering rates slow down and volcanic CO2 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 CO2 degassing by driving volcanism at mid-ocean ridges and island arcs, and by recycling surface CO2 into the mantle such that the mantle CO2 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 CO2 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 CO2 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 CO2 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, CO2 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 CO2 being released by magmatism and metamorphism, and thus allow a weathering feedback to operate; and (2) for how long, if at all, CO2 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 CO2 drawdown rate. In other words, the upper limit to CO2 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 CO2 drawdown rate only depends on the rate at which fresh rock is brought to the surface and no longer depends on temperature or atmospheric CO2 content. Thus weathering cannot adjust to balance CO2 outgassing, accumulation of atmospheric CO2 continues unabated, and an inhospitably hot climate develops, as long as volcanism lasts long enough to significantly increase atmospheric CO2 (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 CO2 outgassing rates can induce global glaciation (Kadoya and Tajika, 2014; Haqq-Misra et al., 2016). We assume that CO2 outgassing is a result of mantle volcanism, which releases mantle CO2, and metamorphic decarbonation of the crust, originally carbonated by weathering, as it is buried by subsequent lava flows. Thus the longevity and strength of CO2 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 CO2 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 CO2 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 CO2 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.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 CO2 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 Vman is the volume of convecting mantle beneath the stagnant lid, ρ is the bulk density of the mantle, cp is the heat capacity, Tp is the potential temperature of the upper mantle, t is time, Qman is the radiogenic heating rate in the mantle, Aman is the surface area of the top of the convecting mantle (base of stagnant lid), Fman is the heat flux from the convecting mantle into the base of the stagnant lid, fm is volumetric melt production, ρm is the density of mantle melt, Lm is the latent heat of fusion of the mantle, and ΔTm is the temperature difference between the melt erupted at the surface and the surface temperature. The volume of the convecting mantle is given by Vman = (4/3)π((Rp − δ)3 − Rc3), where Rp is the planetary radius, Rc is the core radius, and δ is the thickness of the stagnant lid. The surface area of the convecting mantle is then Aman = 4π(Rp − δ)2. In this study, we focus on Earth-sized planets and therefore use Rp equal to Earth's radius and Rc 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).
|ρ||Mantle density||4000 kg m−3||(1)|
|cp||Heat capacity||1250 J kg−1 K−1||(1)|
|ρm||Melt density||2800 kg m−3||(1)|
|Lm||Latent heat||600 × 103 J kg−1||(1)|
|Rp||Planetary radius||6378.1 km||below (1)|
|Rc||Core radius||3488.1 km||below (1)|
|k||Thermal conductivity||5 W m−1 K−1||(2)|
|c1||Scaling law constant for δ||0.5||(3)|
|d||Whole mantle thickness||2890 km||(3)|
|Ts||Surface temperature||273 K||below (3)|
|arh||Constant for rheological temperature scale||2.5||(4)|
|Ev||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)|
|κ||Thermal diffusivity||10−6 m2 s−1||below (4)|
|α||Thermal expansivity||3 × 10−5 K−1||below (4)|
|μn||Viscosity pre-exponential factor||4 × 1010 Pa s||below (4)|
|μr||Reference viscosity||2 × 1020 Pa s||below (4)|
|γmantle||Mantle adiabatic temperature gradient||10−8 K Pa−1||(9)|
|ρl||Lithosphere density||3300 kg m−3||(10)|
|(dφ/dP)S||Pressure derivative of melt fraction||1.5 × 10−10 Pa−1||(11)|
|γmelt||Adiabatic temperature gradient of melt||3.1 × 10−8 K Pa−1||below (12)|
|c2||Scaling law constant for v||0.05||(13)|
|c3||Scaling law constant for v||0.4||(14)|
|τrad||Radioactive decay constant||2.94 Gyr||(18)|
|A||Decarbonation temperature constant||3.125 × 10−3 K m−1||(20)|
|B||Decarbonation temperature constant||835.5 K||(20)|
|Distribution coefficient for CO2||10−4||(23)|
|χ||Fraction of reactable elements in crust||0.3||(29)|
|Molar mass of reactable elements in crust||0.055 kg mol−1||(29)|
|Tp0||Initial mantle potential temperature||2000 K||—|
|Tp||Mantle potential temperature (K)||(1)|
|Vman||Volume of convecting mantle (m3)||(1)|
|Aman||Surface area of top of convecting mantle (m2)||(1)|
|Qman||Radiogenic heat production in the mantle (W)||(1) and (19)|
|Fman||Mantle convective heat flux (W m−2)||(1)|
|fm||Volumetric melt production rate (m3 s−1)||(1)|
|ΔTm||Temperature difference between erupted melt and Ts (K)||(1) and below (12)|
|δ||Thickness of stagnant lid (m)||(2)|
|Tl||Temperature at base of stagnant lid (K)||(2) and (4)|
|θ||Frank-Kamenetskii parameter (unitless)||(3)|
|Rai||Internal Rayleigh number (unitless)||(3)|
|μi||Interior mantle viscosity (Pa s)||below (3)|
|x||Concentration of heat-producing elements (W m−3)||(5)|
|Vcrust||Volume of crust (m3)||below (5) and (15)|
|δc||Thickness of crust (m)||below (5) and (16)|
|Qcrust||Radiogenic heat production in the crust (W)||below (5) and (18)|
|Vlid||Volume of subcrustal stagnant lid (m3)||below (5)|
|xm||Concentration of heat-producing elements in the mantle (W m−3)||(6)|
|xc||Concentration of heat-producing elements in the crust (W m−3)||(7)|
|Tc||Temperature at the base of the crust (K)||(8)|
|Pi||Pressure where melting begins (Pa)||(9)|
|Pf||Pressure where melting ceases (Pa)||(10)|
|δ||Thickness of the lithosphere (m)||(10)|
|φ||Melt fraction (unitless)||(11)|
|dm||Depth where melting begins (m)||(12)|
|v||Convective velocity (m s−1)||(13) and (14)|
|Q0||Initial radiogenic heat production (W)||below (19)|
|Tdecarb||Temperature of metamorphic decarbonation (K)||(20)|
|δcarb||Decarbonation depth (m)||(22)|
|Rman||Mantle carbon reservoir (mol)||(23)|
|Fd||Mantle degassing flux (mol s−1)||below (23)|
|Rcrust||Crustal carbon reservoir (mol)||(25)|
|Fmeta||Metamorphic degassing flux (mol s−1)||(26)|
|Ctot||Total carbon budget (mol)||below (28)|
The thickness of the stagnant lid, δ, evolves with time based on an energy balance model (e.g., Schubert et al., 1979; Spohn, 1991),
where Tl 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, Fman, and lid base temperature, Tl, are calculated from scaling laws for mantle convection (Reese et al., 1998, 1999; Solomatov and Moresi, 2000; Korenaga, 2009):
Here, c1 is a constant, Ts is the surface temperature, d = Rp − Rc is the whole mantle thickness, θ is the Frank-Kamenetskii parameter, defined as θ = Ev(Tp − Ts)/(RTp2), which describes the strength of the temperature dependence of mantle viscosity, Ev is the activation energy for mantle viscosity, R is the gas constant, arh is a constant, and Rai is the internal Rayleigh number. We use c1 = 0.5 and arh = 2.5. Internal Rayleigh number is defined as Rai = ρgα(Tp − Ts)d3/(κμi), where g is gravitational acceleration, α is thermal expansivity, Ts is surface temperature, κ is thermal diffusivity, and μi is the interior mantle viscosity, defined at temperature Tp. Mantle viscosity is temperature-dependent following typical flow laws for diffusion creep in olivine as μi = μn exp(Ev/(RTp)), where μn is a constant. We use Ev ≈ 300 kJ mol−1 as the baseline value in this study (Karato and Wu, 1993); test cases using Ev = 200 kJ−1 mol and Ev = 400 kJ−1 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(Ev/(RTr)), where Tr = 1623 K is Earth's present-day mantle temperature. Using our baseline value of μn =4 × 1010 Pa s, μr ≈ 2 × 1020 Pa s.
Note that in (3) the temperature difference between the mantle interior and surface, Tp − Ts, and whole mantle thickness, d = Rp − Rc, are used instead of the temperature difference between the mantle interior and base of the lid, Tp − Tl, 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 Rp − δc ≤ z ≤ Rp, x = Qcrust/Vcrust ≡ xc, where δc is the thickness of the crust, Qcrust is the radiogenic heating rate within the crust, Vcrust is the volume of the crust, and xc is the crustal heat production rate per unit volume. Within the mantle lithosphere, that is, when Rp − δ ≤ z ≤ Rp − δc, x = Qman/(Vman + Vlid) ≡ xm, where xm is the mantle heat production rate per unit volume and Vlid = (4/3)π((Rp − δc)3 − (Rp − δ)3) 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−3 (see Section 4.3), which results in a downward velocity of crust of ≈6 × 10−12 m s−1 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−6 m2 s−1, 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 = Ts at z = Rp, T = Tl at z = Rp − δ, and by enforcing continuity of temperature and heat flux across the crust-mantle interface, that is, at z = Rp − δ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, Tc, 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)−1 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 Tc is within one degree kelvin of Tl, which is equivalent to switching when the difference between δc and δ is negligibly small.
2.2. Melting and crustal evolution
The melt production rate, fm, is calculated in the same manner as in Fraeman and Korenaga (2010). The pressure where melting begins, Pi, 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−8 K Pa−1). 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−3). The melt fraction, φ, is
where (dφ/dP)S ≈ 1.5 × 10−10 Pa−1.
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(dm − δ)πL, where v is the convective velocity, dm = Pi/(ρ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.45Rp 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 CO2 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 ΔTm = Tp − Piγ − Ts. Melt created where upwelling mantle crosses the solidus has a temperature Tp + Piγmantle. This melt cools by adiabatic decompression as it rises to the surface, by an amount Piγmelt, where γmelt is larger than γmantle. Subtracting the adiabatic cooling from the initial temperature gives the erupted melt temperature as Tp − Piγ, where γ = γmelt − γmantle = 6.766 × 10−4 K km−1 using an average value of γmelt ≈ 1 K km−1 (Driscoll and Bercovici, 2014). The erupted melt temperature subtracted by the surface temperature, Ts, then gives ΔTm.
where c2 and c3 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, Vcrust, 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 − δ)2min(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 xmelt is the concentration of heat-producing elements in the melt, x0 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, Qcrust, 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, fm. Associating the concentration of heat-producing elements in the mantle at a given time, xm, with x0 from (17), the rate of increase in Qcrust due to melting is given by (xm fm /φ)[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, 238U, 235U, 232Th, and 40K (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, Q0, that is entirely within the mantle. Q0 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 CO2 degassing rate decays with time, and whether CO2 degassing overwhelms CO2 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 CO2 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 CO2 drawdown via weathering is always equal to the rate of CO2 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 CO2 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 CO2, 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 % CO2) occurs across a narrow (≈100–200 K) temperature interval, chiefly by the breakdown of dolomite (CaMg(CO3)2). Magnesite is the principal host for CO2 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 % CO2 contour with a simple linear fit.
The temperature at which decarbonation occurs is then
where A = 3.125 × 10−3 K m−1, B = 835.5 K, and Tdecarb is temperature in kelvins. Pressure is converted into depth by assuming a lithospheric density of ρl = 3300 kg m−3. 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 = Rp − δcarb and is therefore given by
When δcarb < δ, crust decarbonates before being recycled into the mantle, and the evolution of the mantle carbon reservoir, Rman, is given by
where Rman 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 CO2 is outgassed from mantle to atmosphere, where weathering will deposit it in the crust. CO2 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 CO2 of = 10−4 (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 CO2 to the atmosphere. Assuming that only a fraction of the melt degasses its CO2 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, Rcrust, 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 Vcarb is the volume of carbonated crust, Vcarb = (4/3)π(Rp3 − (Rp − δcarb)3), and, as in (15), 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 goes to Rcrustfm/(Vcarb) when δc > δcarb and to 0 when δc < δcarb. Metamorphic CO2 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 Rcrust and Rman 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, Ctot = Rman + Rcrust, 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−1), ε 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−1. 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 CO2-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 Fd + Fmeta ≥ . Thus the threshold where supply-limited weathering begins is given by Fd + Fmeta = , and the critical mantle and crustal CO2 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 fmRman/(φ[Vman + Vlid]); 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 Rman >> Rcrust, (30) reduces to an expression for the critical mantle CO2 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 Vman + Vlid = (4/3)π(Rp3 − Rc3), ≈ 8 × 1022 to 8 × 1023 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 Vman + Vlid 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 × 1021 to 2 × 1022 mol of CO2 in Earth's mantle, so a planet would need to have a mantle CO2 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 CO2 budgets would lead to supply-limited weathering, as CO2 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 CO2 buildup. Thus the threshold for supply-limited weathering due to mantle volcanism in practice is the range of ≈8 × 1022 to 8 × 1023 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 Rcrust >> Rman): . With δcarb = 10–100 km, ≈ 1022 to 1023 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 CO2 budget when the CO2 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, Ctot, that resides in the crust. Thus Rcrust = cCtot, Rman = (1 − c)Ctot, and
where is the critical total planetary CO2 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 Rcrust, δ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 CO2 in the melt. In such situations, though, the total outgassing rate is low, on the order of 1010 to 1011 mol yr−1 or less in our models, and it would take billions of years for a thick CO2 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 Fd + Fmeta ≥ + Foffset. We set Foffset = 1014 mol yr−1, as a net CO2 input to the atmosphere of 1014 mol yr−1 will cause 100 bar of CO2 (≈ 1022 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 Foffset over the range Foffset = 1013 to 1015 mol yr−1. 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 CO2 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 CO2 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 CO2 degassing rates can be maintained at both 10% and 100% of Earth's present-day total degassing flux of ≈6 × 1012 mol yr−1 (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 CO2, 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 Tp0 = 2000 K, initial radiogenic heating rate of Q0 = 100 TW, and total CO2 content of Ctot = 1022 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; Qmantle 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 CO2 degassing flux is high throughout the lifetime of volcanism due primarily to metamorphic degassing of the carbonated crust. Metamorphic degassing dominates because CO2 is more concentrated in the crust than in the mantle and because the mantle CO2 reservoir is quickly depleted by volcanism and outgassing. CO2 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 CO2 degassing rates.
The planetary CO2 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 CO2 content allows high degassing rates to be maintained, primarily by metamorphic CO2 degassing, even as volcanism wanes (i.e., with a larger CO2 budget, more CO2 is deposited in the crust, such that degassing rates remain high even at low rates of melt production and crustal burial). The influence of Ctot and Q0 on the longevity of CO2 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 CO2 budget lead to longer-lived CO2 degassing. The lifetime of volcanism is not a function of planetary CO2 budget, and thus contours for the time when volcanism ends would form horizontal lines in Fig. 5A and 5B. Moving to large planetary CO2 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 CO2, becomes the dominant factor determining the longevity of adequate degassing. Moving toward lower CO2 budgets, the contours curve toward vertical lines, as the planetary CO2 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 CO2 inventory will quickly become globally frozen, even though volcanism can continue for billions of years.
Although large CO2 budgets are beneficial for promoting warm climates, there is an upper limit on the CO2 inventory for habitability provided by the transition to supply-limited weathering discussed in Section 3. The critical planetary CO2 inventory where weathering becomes supply-limited was calculated from the thermal evolution models as explained in Section 3; this is found to occur when Ctot ≈ 1023 mol at Q0 = 5 TW and Ctot ≈ 2 × 1022 mol at Q = 250 TW. The relationship between and Q0 is approximately linear. Larger heat production rates shift the boundary where weathering becomes supply-limited to lower values of Ctot (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 CO2 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 Q0 > 100 TW and Ctot ∼ 1020 to 1022 mol (corresponding to ≈7 × 10−5 to 7 × 10−3 mass % CO2), where habitable climates can potentially be maintained for up to 2–4 Gyr. Earth's total CO2 budget is not well constrained, but estimates range from ∼1022 to 1023 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 CO2 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 CO2 budgets similar to Earth's; larger CO2 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 Tp0, 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 Q0 = 250 TW degassing wanes in 3–4 Gyr. Moreover, a lower μr also shifts the onset of supply-limited weathering to lower Ctot, 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 Ctot–Q0 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 Q0 > 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 Ctot, 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 Ctot 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 Ctot. 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 Ctot, 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 Ctot and Q0 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 Rai1/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, fm, 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, Ev. 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 CO2 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−1 above and below our baseline value does not change the region of Ctot–Q0 space where potentially habitable climates are possible but does change the lifetime of sufficient CO2 degassing; with Ev = 200 kJ mol−1, degassing lasts longer than in the baseline model, while with Ev = 400 kJ−1 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 Ev is changed. For Ev = 200 kJ mol−1, μn = 7.3 × 1013 Pa s, and for Ev = 400 kJ mol−1, μn = 2.7 × 107 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 Ctot 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 Ctot, 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 CO2 supply largely remains in the mantle, and the large volume of the mantle means this carbon is diluted. Thus, the same Ctot 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 εfm with ε = 0.1 (Fig. 8). Eruption rates are typically < 100 km3 yr−1 for most of the period while volcanism is active, which results in ≈0.2 mm yr−1 of resurfacing when erupted lava is assumed to be distributed globally. For comparison, resurfacing rates at Kilauea Volcano in Hawaii are ≈5 mm yr−1 (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 104 km3 yr−1), and potentially uninhabitable, in our models, but such high volcanic rates are short lived (lasting on the order of 100 Myr).
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 CO2 budgets of ∼1020 to 1022 mol (≈7 × 10−5 to 7 × 10−3 mass % CO2), are potentially habitable for timescales of over 2 Gyr. The time when habitability becomes unlikely due to limited CO2 outgassing scales with initial radiogenic heat budget, with higher budgets leading to longer-lived outgassing. Thus planet heat budget, CO2 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 CO2 inventories. Accreted CO2 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). CO2 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 CO2 inventories. Ultimately, it is difficult to predict or constrain from stellar observations which exoplanet systems are likely to have planets with bulk CO2 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 CO2, 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 CO2 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, CO2 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 CO2 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, CO2 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 CO2 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 CO2 inventories up to ∼1024 to 1025 mol, much larger than the ∼1023 mol upper limit for stagnant lid planets when crustal decarbonation is active. With little or no land, the CO2 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 Ctot–Q0 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 CO2 locked in the crust, which would enhance CO2 outgassing during mantle volcanism. Melt liberation of crustal carbon is a potentially significant CO2 source on Earth (Lee et al., 2013), and including it would allow habitable conditions at lower total CO2 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 km3 yr−1 of melt flux (Crisp, 1984), which is similar to the melt fluxes seen in our models when volcanism is close to ceasing and CO2 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 CO2 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 Ctot–Q0 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 CO2 outgassing lasting longer than ∼105 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.
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 CO2 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 CO2 budgets of ∼1020 to 1022 mol (ranging from 10−2 to 1 times typical estimates for Earth's CO2 budget), can potentially maintain temperate climates for 1–5 Gyr. After this time, volcanism and CO2 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 CO2 budgets lower than ∼1020 mol, a planet's climate is estimated to be in a snowball state for its entire history, while above ∼1022 mol weathering would become supply-limited. With supply-limited weathering, CO2 outgassing overwhelms CO2 drawdown, such that an inhospitably hot, CO2-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 CO2 degassing, and that young planets are more likely to experience sufficiently high rates of CO2 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)
Clinopyroxene: Cpx(HP)—Holland and Powell (1996)
Epidote: Ep(HP)—Holland and Powell (1998)
CO2-H2O 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 CO2 in the melt is assumed to degas to the atmosphere. Moreover, all CO2 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 CO2 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 CO2 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 CO2 in this melt degasses; the rest is trapped within the crust and becomes part of the crustal CO2 reservoir. As a result, the volcanic degassing flux, Fd, is reduced by a factor ζ1:
We also allow for a fraction, ζ2, of the CO2 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 CO2 from the mantle to the crust is the same as in the main text, as all CO2 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 ζ1 = 0.1 and ζ2 = 0 and with ζ1 = 0.01 and ζ2 = 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 ζ2 > 0 is used; specifically as ζ2 approaches 1, model results converge to the case where crustal decarbonation is ignored (i.e.,Fig. 7B in the main text), as ζ2 = 1 implies no metamorphic CO2 outgassing and hence all the crustal carbon being recycled into the mantle by foundering. As a result, models with ζ1 = 0.1 and ζ2 = 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 CO2 budgets where long-lived degassing occurs shift to higher values. Thus, the range of Ctot 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 CO2 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 CO2 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 Ctot and Q0 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 ζ2 = 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 CO2 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 CO2 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 CO2 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.
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.