Mantle Degassing Lifetimes through Galactic Time and the Maximum Age Stagnant-lid Rocky Exoplanets Can Support Temperate Climates

Foley & Smye (2018) and Foley (2019) gave a thorough description of our model, and all of the key governing equations are listed in Appendix A.2. Here we will only highlight the differences between the model of Foley & Smye (2018) and the model presented in this paper. For an Earth-like CMF = 0.33, we use the scaling laws from Valencia et al. (2006, 2007) to determine the average mantle density, planet radius, mantle thickness, and surface gravity as a function of planet mass:

$$\begin{eqnarray}&&\rho ={\rho }_{\oplus }{\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{0.2},\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&{R}_{p}={R}_{\oplus }{\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{0.27},\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&d={d}_{\oplus }{\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{0.28},\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&g=\displaystyle \frac{{{GM}}_{p}}{{R}_{p}^{2}},\end{eqnarray} \tag{ A4 }$$

where M_p is the planet mass, and G is the gravitational constant. Reference Earth values are ρ_⊕ = 4450 kg m⁻³, R_⊕ = 6378 km, and d_⊕ = 2890 km.

To vary the CMF, we use the scaling laws from Noack et al. (2016; see also Foley et al. 2020), which assume that all iron resides in the core. The resulting equations for R_p , R_c , and ρ, as a function of CMF and planet mass, are

$$\begin{eqnarray}&&{R}_{p}=(7121-2021\times \mathrm{CMF}){\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{0.265}\mathrm{km},\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}\begin{array}{rcl}{R}_{c} & = & (19200\times \mathrm{CMF}-31760\times {\mathrm{CMF}}^{2}\\ & & +\,18100\times {\mathrm{CMF}}^{3}){\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{0.252}\mathrm{km},\end{array}\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&\rho =\displaystyle \frac{(1-\mathrm{CMF}){M}_{p}}{\tfrac{4}{3}\pi ({R}_{p}^{3}-{R}_{c}^{3})}.\end{eqnarray} \tag{ A7 }$$

Gravity is still calculated from Equation (A4).

We assume that all other material properties are independent of planet size and CMF. This is a simplification because thermal conductivity, expansivity, and viscosity are all functions of pressure, and the larger the planet or core, the higher the pressures reached in the mantle. However, robust parameterizations for how to incorporate these pressure effects are currently lacking. There is still significant uncertainty about how key material properties change at the extreme temperature and pressure conditions of super-Earth lower mantles, which are not yet experimentally accessible (e.g., Duffy et al. 2015). Moreover, how significant variation of key material properties with pressure, and hence depth in a super-Earth mantle, modifies the dynamics of the convecting mantle has not been extensively studied; scaling laws for convective heat flux and velocity that take these pressure effects into account have not yet been developed.

Lacking robust parameterizations for viscosity, thermal expansivity, and thermal conductivity pressure effects, we instead randomly vary the mantle reference viscosity, μ_ref, in our models over a range of 4 orders of magnitude. We focus on viscosity because, of the key mantle material properties, it shows the strongest dependence on pressure and mantle composition, as it can vary by orders of magnitude (e.g., Karato & Wu 1993; Hirth & Kohlstedt 2003). We use a standard Arrhenius temperature-dependent viscosity law in our models,

$$\begin{eqnarray}&&{\mu }_{i}={\mu }_{n}\exp \left(\displaystyle \frac{{E}_{v}}{{{RT}}_{p}}\right),\end{eqnarray} \tag{ A8 }$$

where μ_i is the mantle interior viscosity, T_p is the mantle potential temperature, E_v = 300 kJ mol⁻¹ (e.g., Karato & Wu 1993) is the activation energy, and R is the universal gas constant. The reference viscosity is defined at the reference temperature T_ref = 1623 K, or approximately Earth's present-day mantle potential temperature. The constant μ_n is then adjusted in each model run to match the chosen reference viscosity. Our results therefore explicitly demonstrate how the degassing lifetime depends on the mantle reference viscosity, which itself may vary with planet size or composition.

The equations for calculating the rate of metamorphic degassing due to crustal burial from Foley & Smye (2018) are reformulated in terms of pressure, rather than depth (see Appendix A.2.3). As such, they can then be applied to planets with variable sizes and CMFs and hence different surface gravities. Finally, we have improved the melting model to include depletion of the mantle and a subsequent increase in the melting temperature. We follow the method of Tosi et al. (2017) and assume that the solidus can increase by up to 150 K upon full depletion of the mantle, the difference in the solidi of harzburgite and peridotite. The degree of mantle depletion is calculated based on the volume of crust present at each time step (see Appendix A.2.1). The model also tracks each of the four major HPEs separately, rather than treating them together with an average decay constant as in Foley & Smye (2018). Here we are interested in observationally constrained variations in each of the four major HPEs, so we naturally must treat each HPE separately in the model (see Appendix A.2.1). The remainder of the model equations are general and can be applied to planets with different masses and CMFs.

Here we give a complete description of the coupled thermal evolution and volatile cycling model. Assuming pure internal heating, and that all melt produced contributes to the cooling of the mantle, the mantle thermal evolution is given by

$$\begin{eqnarray}&&{V}_{\mathrm{man}}\rho {c}_{p}\displaystyle \frac{{{dT}}_{p}}{{dt}}={Q}_{\mathrm{man}}-{A}_{\mathrm{man}}{F}_{\mathrm{man}}-{f}_{m}{\rho }_{m}\left({c}_{p}{\rm{\Delta }}{T}_{m}+{L}_{m}\right),\end{eqnarray} \tag{ A9 }$$

where V_man is the volume of the actively convecting mantle, ρ is the average density of the mantle, c_p is the heat capacity, T_P is the potential temperature of the mantle, t is time, Q_man is the total radiogenic heat production rate of the mantle, A_man is the surface area of the top of the convecting mantle (the base of the stagnant lid), F_man is the heat flux from the mantle, f_m is the volumetric melt production rate, ρ_m is the density of melt, ΔT_m is the temperature difference between the erupted melt and the surface temperature, and L_m is the latent heat of the mantle (e.g., Stevenson et al. 1983; Hauck & Phillips 2002; Reese et al. 2007; Fraeman & Korenaga 2010; Morschhauser et al. 2011; Driscoll & Bercovici 2014; Foley & Smye 2018; Foley 2019). The volume of the convecting mantle is ${V}_{\mathrm{man}}={(4/3)\pi (({R}_{p}-\delta )}^{3}-{R}_{c}^{3})$, where R_p is the planet radius, R_c is the core radius, and δ is the thickness of the stagnant lid. The surface area of the top of the convecting mantle is then ${A}_{\mathrm{man}}=4\pi {({R}_{p}-\delta )}^{2}$. Finally, the temperature difference between the erupted melt and the surface temperature is ΔT_m = T_p − P_i γ − T_s , where γ is the average adiabatic temperature gradient of the mantle melt, estimated as γ ≈ 2 × 10⁻⁸ K Pa⁻¹ in Foley & Smye (2018), and _Pi is the pressure where melting begins (derived below in Section A.2.1).

The thickness of the stagnant lid, δ, is then (e.g., Schubert et al. 1979; Spohn 1991)

$$\begin{eqnarray}&&\rho {c}_{p}({T}_{p}-{T}_{l})\displaystyle \frac{d\delta }{{dt}}=-{F}_{\mathrm{man}}-k{\left.\displaystyle \frac{\partial T}{\partial z}\right|}_{z={R}_{p}-\delta },\end{eqnarray} \tag{ A10 }$$

where T_l is the temperature at the base of the stagnant lid, k is the thermal conductivity (assumed to be the same throughout the crust and mantle for simplicity), and z is the height above the planet's center. The mantle heat flux, F_man, and lid base temperature, T_l , are calculated from the following scaling laws for stagnant-lid convection (Reese et al. 1998, 1999; Solomatov & Moresi 2000; Korenaga 2009):

$$\begin{eqnarray}&&{F}_{\mathrm{man}}=\displaystyle \frac{{c}_{1}k({T}_{p}-{T}_{s})}{d}{{\rm{\Theta }}}^{-4/3}{{Ra}}_{i}^{1/3}\end{eqnarray} \tag{ A11 }$$

and

$$\begin{eqnarray}&&{T}_{l}={T}_{p}-\displaystyle \frac{{a}_{{rh}}{{RT}}_{p}^{2}}{{E}_{v}},\end{eqnarray} \tag{ A12 }$$

where c₁ and a_rh are constants (assumed to be c₁ = 0.5 and a_rh = 2.5), and T_s is the surface temperature, here fixed to 273 K, as surface temperature fluctuations of order 100 K that could result from changes in atmospheric CO₂ do not significantly impact the evolution of the underlying mantle in these models (Foley 2019). The mantle thickness is d = R_p − R_c , and Θ is the Frank–Kamenetskii parameter, ${\rm{\Theta }}={E}_{v}({T}_{p}-{T}_{s})/({{RT}}_{p}^{2})$, where E_v is the activation energy for mantle viscosity, and R is the universal gas constant. The internal Rayleigh number, Ra_i , is defined as Ra_i = ρ g α(T_p − T_s )d³/(κ μ_i ), where g is gravity, α is the thermal expansivity of the mantle, κ is the thermal diffusivity of the mantle, and μ_i is the viscosity at the mantle potential temperature of T_p (see Equation (A8)).

A.2.1. Melting and Crustal Evolution

Large-scale mantle melting and subsequent volcanism take place when passively upwelling mantle is hot enough to cross the solidus beneath the lid. As in, e.g., Fraeman & Korenaga (2010) and Foley & Smye (2018), the pressure at which melting begins, P_i , is calculated from the intersection of the mantle adiabat, with adiabatic gradient γ_mantle, and the dry peridotite solidus from Takahashi & Kushiro (1983),

$$\begin{eqnarray}&&{P}_{i}=\displaystyle \frac{{T}_{p}-{T}_{{\rm{0}}}^{{\rm{sol}}}}{120\times {10}^{-9}-{\gamma }_{{\rm{m}}{\rm{a}}{\rm{n}}{\rm{t}}{\rm{l}}{\rm{e}}}},\end{eqnarray} \tag{ A13 }$$

where T^sol ₀ is the melting solidus temperature at zero pressure, and the adiabatic gradient in the mantle is γ_mantle ≈ 2 × 10⁻⁸ K Pa⁻¹. The solidus temperature at zero pressure depends on the degree of mantle depletion, which we estimate based on the volume of crust present, as explained below. We also do not allow the pressure where melting begins to exceed 10 GPa, the approximate pressure where silicate melts become denser than solids (e.g., Dorn et al. 2018). Melting stops at the base of the lid, which occurs at pressure

$$\begin{eqnarray}&&{P}_{f}={\rho }_{l}g\delta ,\end{eqnarray} \tag{ A14 }$$

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

$$\begin{eqnarray}&&\phi =\displaystyle \frac{{P}_{i}-{P}_{f}}{2}{\left(\displaystyle \frac{d\phi }{{dP}}\right)}_{S},\end{eqnarray} \tag{ A15 }$$

where (d ϕ/dP)_S ≈ 1.5 × 10⁻¹⁰ Pa⁻¹. The melt production rate, f_m , is calculated as (see Foley & Smye 2018 for a derivation)

$$\begin{eqnarray}&&{f}_{m}=17.8\pi {R}_{p}v({d}_{\mathrm{melt}}-\delta )\phi ,\end{eqnarray} \tag{ A16 }$$

where v is the characteristic convecting mantle velocity and d_melt = P_i /(ρ_l g) is the depth where melting begins. The convecting mantle velocity is (Reese et al. 1998, 1999; Solomatov & Moresi 2000; Korenaga 2009)

$$\begin{eqnarray}&&v={c}_{2}\displaystyle \frac{\kappa }{d}{\left(\displaystyle \frac{{{Ra}}_{i}}{{\rm{\Theta }}}\right)}^{2/3},\end{eqnarray} \tag{ A17 }$$

where c₂ is a constant.

Melting produces a crust whose thickness, δ_c , and volume, V_crust, evolve over time. To calculate the evolution of the V_crust, we assume that all melt produced contributes to the growth of the crust, and that all crust buried to depths below the lithospheric thickness, δ, founders into the mantle. The resulting equation is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dV}}_{\mathrm{crust}}}{{dt}} & = & {f}_{m}-\left({f}_{m}-4\pi {({R}_{p}-\delta )}^{2}{\rm{\min }}\left(0,\displaystyle \frac{d\delta }{{dt}}\right)\right)\\ & & \times \,(\tanh (({\delta }_{c}-\delta )20)+1).\end{array}\end{eqnarray} \tag{ A18 }$$

The second term on the right-hand side of Equation (A18) describes the rate of crust loss due to foundering of the crust; the hyperbolic tangent function formulation allows this crust loss rate to go to zero when δ_c < δ. The term $4\pi {({R}_{p}-\delta )}^{2}{\rm{\min }}(0,d\delta /{dt})$ captures the loss of crust when the lid thickness is decreasing and δ_c = δ and is zero otherwise (that is, if either the lid thickness is growing or the crust ends before the base of the stagnant lid). The crustal thickness is calculated from the volume of crust as

$$\begin{eqnarray}&&{\delta }_{c}={R}_{p}-{\left({R}_{p}^{3}-\displaystyle \frac{3{V}_{\mathrm{crust}}}{4\pi }\right)}^{1/3}.\end{eqnarray} \tag{ A19 }$$

To incorporate how depletion of the mantle influences the solidus, and thus later melt production, we increase T^sol ₀ linearly with crust thickness following (Tosi et al. 2017)

$$\begin{eqnarray}&&{T}_{{\rm{0}}}^{{\rm{sol}}}=1423+{\rm{\Delta }}{T}_{{\rm{s}}{\rm{o}}{\rm{l}}}\times {\rm{\max }}\left(1,\displaystyle \frac{{\delta }_{c}}{{\delta }_{ref}}\right).\end{eqnarray} \tag{ A20 }$$

Here 1423 K is the dry peridotite solidus temperature at zero pressure from Takahashi & Kushiro (1983), ΔT_sol = 150 K is the increase in the solidus upon full depletion (which is set here to the difference in the zero-pressure solidus temperatures for peridotite and harzburgite), and ${\delta }_{\mathrm{ref}}=0.2{V}_{\mathrm{man}}^{0}/{A}_{\mathrm{surf}}$ is the reference crust thickness produced upon full depletion of the mantle. Here ${V}_{\mathrm{man}}^{0}=(4/3)\pi ({R}_{p}^{3}-{R}_{c}^{3})$, and ${A}_{\mathrm{surf}}=4\pi {R}_{p}^{2}$ is the surface area of the planet. The solidus can thus increase by up to 150 K due to mantle depletion. We explore other compositional effects that affect the solidus (e.g., water) in the main text.

The HPEs are preferentially partitioned into the crust during mantle melting due to their incompatible nature. We track this partitioning for all four long-lived HPEs assuming accumulated fractional melting. The evolution of the crustal heat production rate for a given HPE is thus

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dQ}}_{c,i}}{{dt}} & = & \displaystyle \frac{{x}_{m,i}{f}_{m}}{\phi }[1-{(1-\phi )}^{1/{D}_{i}}]\\ & & -\,{x}_{c,i}\left({f}_{m}-4\pi {({R}_{p}-\delta )}^{2}{\rm{\min }}\left(0,\displaystyle \frac{d\delta }{{dt}}\right)\right)\\ & & \times \,(\tanh (({\delta }_{c}-\delta )20)+1)-\displaystyle \frac{{Q}_{c,i}}{{\tau }_{i}}.\end{array}\end{eqnarray} \tag{ A21 }$$

Here Q_c,i is the heat production rate in the crust resulting from one of the four HPEs tracked in the model. The total crustal heat production rate Q_crust = Q_c,U238 + Q_c,U235 + Q_c,Th + Q_c,K. Each HPE has a specified decay constant, τ_i ; distribution coefficient, D_i ; and crustal and mantle heat production rate per unit volume, x_c,i and x_m,i , respectively (e.g., τ_U238 is the decay constant for ²³⁸U, D_U238 is the distribution coefficient for ²³⁸U, and x_c,U238 is the heat production per unit volume in the crust due to ²³⁸U). The well-known half-lives taken from Turcotte & Schubert (2002) are used to calculate the decay constants, τ_i . The chosen distribution coefficients of D_U238 = D_U235 = 0.0012 and D_Th =0.0029 are from Beattie (1993), and D_K = 0.0011 is from Hart & Brooks (1974), assuming 60% olivine and 40% pyroxene in the mantle. The evolution of mantle heat production from a given radioactive isotope is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dQ}}_{m,i}}{{dt}}={x}_{c,i}\left({f}_{m}-4\pi {({R}_{p}-\delta )}^{2}{\rm{\min }}\left(0,\displaystyle \frac{d\delta }{{dt}}\right)\right)\\ \,\times \left(\tanh (({\delta }_{c}-\delta )20)+1\right)-\displaystyle \frac{{x}_{m,i}{f}_{m}}{\phi }\left[1-{\left(1-\phi \right)}^{1/{D}_{i}}\right]\\ \,-\displaystyle \frac{{Q}_{m,i}}{{\tau }_{i}}.\end{array}\end{eqnarray} \tag{ A22 }$$

As before, the total heat production rate of the mantle is Q_man = Q_m,U238 + Q_m,U235 + Q_m,Th + Q_m,K. The total heat production rates per unit volume in the crust and mantle are also sums of the heat production rates of the four HPEs: x_c = x_c,U238 + x_c,U235 + x_c,Th + x_c,K and x_m = x_m,U238 + x_m,U235 +x_m,Th + x_m,K.

A.2.2. Crustal Geotherm

Equation (A10) requires as an input the conductive heat flux at the base of the lid. The temperature profile for steady-state one-dimensional heat conduction with constant heat production rates in the crust and mantle, neglecting advection, is used to determine this heat flux (for details, see Foley & Smye 2018):

$$\begin{eqnarray}&&{\left.-k\displaystyle \frac{\partial T}{\partial z}\right|}_{z={R}_{p}-\delta }=\displaystyle \frac{k({T}_{l}-{T}_{c})}{\delta -{\delta }_{c}}-\displaystyle \frac{{x}_{m}(\delta -{\delta }_{c})}{2}\end{eqnarray} \tag{ A23 }$$

when δ > δ_c , and

$$\begin{eqnarray}&&{\left.-k\displaystyle \frac{\partial T}{\partial z}\right|}_{z={R}_{p}-\delta }=-\displaystyle \frac{{x}_{c}{\delta }_{c}}{2}+\displaystyle \frac{k({T}_{l}-{T}_{s})}{{\delta }_{c}}\end{eqnarray} \tag{ A24 }$$

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

$$\begin{eqnarray}&&{T}_{c}=\displaystyle \frac{{T}_{s}(\delta -{\delta }_{c})+{T}_{l}{\delta }_{c}}{\delta }+\displaystyle \frac{{x}_{c}{\delta }_{c}^{2}(\delta -{\delta }_{c})+{x}_{m}{\delta }_{c}{\left(\delta -{\delta }_{c}\right)}^{2}}{2k\delta }.\end{eqnarray} \tag{ A25 }$$

The mantle radiogenic heating rate per unit volume, x_m =Q_man/(V_man + V_lid), where ${V}_{\mathrm{lid}}={(4/3)\pi {(({R}_{p}-{\delta }_{c})}^{3}-({R}_{p}-\delta )}^{3})$, is the volume of the subcrustal stagnant lid. The crustal radiogenic heat rate per unit volume is x_c = Q_crust/V_crust, where Q_crust is the total radiogenic heating rate in the crust, and V_crust is the volume of the crust.

A.2.3. CO₂ Outgassing Rates

The CO₂ is outgassed to the atmosphere due to mantle melting, subsequent volcanism, and metamorphic breakdown of carbonated minerals as the crust is buried (Foley & Smye 2018; Foley 2019). Foley & Smye (2018) showed that the temperature as a function of depth where metamorphic decarbonation occurs can be approximated as a simple linear relationship,

$$\begin{eqnarray}&&{T}_{\mathrm{decarb}}=A{\rho }_{l}g({R}_{p}-z)+B,\end{eqnarray} \tag{ A26 }$$

where A = 9.66 × 10⁻⁸ K Pa⁻¹, B = 835.5 K, T_decarb is the temperature in Kelvin, and ρ_l = 3300 kg m⁻³ is the lithosphere density. The depth where decarbonation occurs, δ_carb, is (for details, see Foley & Smye 2018)

$$\begin{eqnarray}\begin{array}{rcl}{\delta }_{\mathrm{carb}} & = & \displaystyle \frac{{\delta }_{c}}{2}+\displaystyle \frac{k({T}_{c}-{T}_{s})}{{\delta }_{c}{x}_{c}}-\displaystyle \frac{A{\rho }_{l}{gk}}{{x}_{c}}\\ & & -\,\displaystyle \frac{k}{{x}_{c}}\sqrt{{\left(\displaystyle \frac{{x}_{c}{\delta }_{c}}{2k}+\displaystyle \frac{{T}_{c}-{T}_{s}}{{\delta }_{c}}-A{\rho }_{l}g\right)}^{2}+\displaystyle \frac{2{x}_{c}}{k}({T}_{s}-B)}.\end{array}\end{eqnarray} \tag{ A27 }$$

With the decarbonation depth determined by Equation (A27), the metamorphic outgassing flux is given by

$$\begin{eqnarray}&&{F}_{\mathrm{meta}}=\displaystyle \frac{{R}_{\mathrm{crust}}{f}_{m}}{2{V}_{\mathrm{carb}}}(\tanh (({\delta }_{c}-{\delta }_{\mathrm{carb}})20)+1),\end{eqnarray} \tag{ A28 }$$

where R_crust is the size of the crustal CO₂ reservoir (in moles), and the hyperbolic tangent function allows F_meta to go to zero when δ_c < δ_carb and R_crust f_m /(V_carb) when δ_c > δ_carb. The outgassing rate due to mantle melting is

$$\begin{eqnarray}&&{F}_{d}=\displaystyle \frac{{f}_{m}{R}_{\mathrm{man}}[1-{\left(1-\phi \right)}^{1/{D}_{{{\rm{CO}}}_{2}}}]}{\phi ({V}_{\mathrm{man}}+{V}_{\mathrm{lid}})},\end{eqnarray} \tag{ A29 }$$

where R_man is the size of the mantle CO₂ reservoir, and ${D}_{{{CO}}_{2}}={10}^{-4}$ is the distribution coefficient for CO₂ (Hauri et al. 2006). The evolution of the mantle and crustal reservoirs when crustal decarbonation is active (e.g., δ_carb < δ_c ) follows

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}_{\mathrm{man}}}{{dt}}=-{F}_{d},\end{eqnarray} \tag{ A30 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}_{\mathrm{crust}}}{{dt}}={F}_{d},\end{eqnarray} \tag{ A31 }$$

while when crustal decarbonation is inactive (δ_c < δ_carb),

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dR}}_{\mathrm{crust}}}{{dt}} & = & {F}_{d}-\displaystyle \frac{{R}_{\mathrm{crust}}}{{V}_{\mathrm{crust}}}\left({f}_{m}-4\pi {\left({R}_{p}-\delta \right)}^{2}{\rm{\min }}\left(0,\displaystyle \frac{d\delta }{{dt}}\right)\right)\\ & & \times \,(\tanh (({\delta }_{c}-\delta )20)+1),\end{array}\end{eqnarray} \tag{ A32 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dR}}_{\mathrm{man}}}{{dt}} & = & -{F}_{d}+\displaystyle \frac{{R}_{\mathrm{crust}}}{{V}_{\mathrm{crust}}}\left({f}_{m}-4\pi {\left({R}_{p}-\delta \right)}^{2}{\rm{\min }}\left(0,\displaystyle \frac{d\delta }{{dt}}\right)\right)\\ & & \times \,(\tanh (({\delta }_{c}-\delta )20)+1).\end{array}\end{eqnarray} \tag{ A33 }$$

The total CO₂ budget of the mantle and crust, C_tot, is conserved, such that C_tot = R_man + R_crust. As in Foley & Smye (2018), the planet starts with all of the CO₂ residing in the mantle, and CO₂ is outgassed over time to the surface. The entire allotment of each HPE also initially resides in the mantle, as before planetary evolution begins, we assume that there is no crust present (crust formation does not take place until mantle convection, and subsequent volcanism, begins). The abundance of each HPE is linked to the stellar observed abundances, as explained in Appendix A.3.

Finally, an important limit for the carbon cycle and habitability is the global weathering supply limit. This limit is the upper bound on the weathering rate and set by the rate at which CO₂ drawdown will occur if all available fresh surface rock is completely carbonated as soon as it is brought to the surface. On a stagnant-lid planet, this limit is assumed to be set by volcanism. Using the average composition of basalt on Earth, the total amount of CO₂ that can be drawn down by crustal carbonation is χ = 5.8 mol kg⁻¹ of basalt (Foley 2019). We assume the crust will have a basaltic composition, as it is a result of primary mantle melting of a peridotitic mantle. Exoplanets, however, could have different mantle bulk compositions that would lead to different crustal compositions upon melting. However, the weathering demand, χ, only increases by about a factor of 2 for the extreme ultramafic composition end-member, peridotite, and only decreases by about a factor of 2 if the crust is felsic, like Earth's continental crust. Average compositions from Taylor & McLennan (1985) for continental crust and Warren (2016) for peridotite were used in making these estimates. The total variation in weathering demand is thus approximately a factor of 4 from ultramafic to felsic end-members. This would not significantly change the planetary conditions (e.g., mantle CO₂ budget, size, internal heat production rate, etc.) that control when planets enter a supply-limited weathering regime and thus develop hothouse climates (Foley & Smye 2018; Foley 2019).

The weathering supply limit, F_sl, assuming all erupted basalt is available for weathering, is

$$\begin{eqnarray}&&{F}_{{\rm{sl}}}=\epsilon {f}_{m}\chi {\rho }_{l}\end{eqnarray} \tag{ A34 }$$

with units of mol yr⁻¹; is the fraction of mantle melt produced that erupts at the surface ( = 0.1 is assumed). In the models, supply-limited weathering is assumed to lead to an inhospitable hothouse climate when F_d + F_meta > F_sl + 10¹⁴ mol yr⁻¹ at any point during planetary evolution, as in Foley & Smye (2018). The total outgassing rate exceeding the weathering supply limit by 10¹⁴ mol yr⁻¹ means that hot climates, with T_s ≥ 400 K, would form in ∼100 Myr, well within the typical degassing lifetimes of modeled planets.

We define the initial total heat production rate of the mantle, Q₀, as a function of the specific power (P_X ) produced by each HPE (X), their concentration within the planet (C_X ), and the mass of the mantle (m_man). We calculate the initial HPE abundance as Q₀ = P_X × C_X × m_man. In order to quantify the range of HPE concentrations in rocky exoplanets (${C}_{X}^{\mathrm{planet}}$), we randomly sample within the observationally constrained stellar abundance distributions for each HPE (Figure 1 of main text) and apply the scaling relationship to convert from stellar abundance to mantle concentration as a function of planet formation time, t, after the birth of the Milky Way:

$$\begin{eqnarray}&&\displaystyle \frac{{C}_{X}^{\mathrm{planet}}(t)}{{C}_{X}^{\mathrm{Earth}-\mathrm{like}}(t)}=\left(\displaystyle \frac{{f}_{X}^{\mathrm{planet}}}{{f}_{X}^{\mathrm{Earth}}}\right)\times \displaystyle \frac{{C}_{X}^{\mathrm{star}}}{{C}_{X}^{\mathrm{Sun}}},\end{eqnarray} \tag{ A35 }$$

where ${C}_{X}^{\mathrm{star}}$ and ${C}_{X}^{\mathrm{Sun}}$ are the concentrations of the element in the randomly selected stellar abundance and the Sun, respectively, and f_X is a correction for volatility effects during planet formation. Here ${C}_{X}^{\mathrm{Earth}-\mathrm{like}}(t)$ represents the predicted initial concentration of an HPE if it was a "cosmochemically Earth-like" planet (Frank et al. 2014) forming at some t (Figure A1). We rescale the results of Frank et al. (2014) such that the Earth-like abundance of each HPE at the time the Earth formed (t = 8 Gyr) is that of Palme & O'Neill (2003). This lowers our choice of ${C}_{X}^{\mathrm{Earth}-\mathrm{like}}$ by ∼25% for each HPE compared to those presented in Frank et al. (2014), who adopted the initial Earth HPE abundances of Turcotte & Schubert (2002).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We define the concentration of an HPE in both stars and the Sun as its molar ratio with Mg (e.g., X/Mg). We normalize relative to Mg because it is more likely to remain in the mantle, as opposed to Si, which may partition into the core (Hirose et al. 2013). Equation (A35) then becomes

$$\begin{eqnarray}&&{C}_{X}^{{\rm{p}}{\rm{l}}{\rm{a}}{\rm{n}}{\rm{e}}{\rm{t}}}(t)=\left(\displaystyle \frac{{f}_{X}^{{\rm{p}}{\rm{l}}{\rm{a}}{\rm{n}}{\rm{e}}{\rm{t}}}}{{f}_{X}^{{\rm{E}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{h}}-{\rm{l}}{\rm{i}}{\rm{k}}{\rm{e}}}}\right)\ast \displaystyle \frac{{\left(X/{\rm{M}}{\rm{g}}\right)}_{{\rm{s}}{\rm{t}}{\rm{a}}{\rm{r}}}}{{\left(X/{\rm{M}}{\rm{g}}\right)}_{{\rm{S}}{\rm{u}}{\rm{n}}}}\ast {C}_{X}^{{\rm{E}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{h}}-{\rm{l}}{\rm{i}}{\rm{k}}{\rm{e}}}(t).\end{eqnarray} \tag{ A36 }$$

We define f_X as the fraction that an element is enriched or depleted relative to Mg during planet formation relative to the host star:

$$\begin{eqnarray}&&{f}_{X}=\displaystyle \frac{{\left(X/\mathrm{Mg}\right)}_{\mathrm{planet}}}{{\left(X/\mathrm{Mg}\right)}_{\mathrm{star}}}.\end{eqnarray} \tag{ A37 }$$

For those elements that fractionate relative to Mg during planet formation, f_X will be greater than 1 in the case of enrichment and less than 1 in the case of depletion. Of the HPEs, U and Th are both refractory and not expected to fractionate between the star and planet relative to Mg. That is to say, f_U and f_Th are ∼1 for both a rocky exoplanet and its host star, as well as the Earth and Sun. For the current bulk silicate Earth (mantle + crust), the present-day molar ratios of Th/Mg and U/Mg are 3.4 × 10⁻⁸ and 9 × 10⁻⁹, respectively (McDonough 2003). Comparatively, the Sun's current composition (Lodders et al. 2009), which is usually defined to be equal to the abundances in CI chondrites, has molar values of Th/Mg = 3.5 × 10⁻⁸ and U/Mg = 9.8 × 10⁻⁹. Among all chondrites, our best proxies for planet-forming materials, their whole-rock Th/Mg and U/Mg abundance ratios, are within <10% of the CI value, except for CV chondrites, which are 1.4 times the solar value (Wasson et al. 1988). The model of Desch et al. (2018), which computes how refractory elements redistribute themselves in protoplanetary disks, predicts small deviations (<20%) of the molar ratios of elements at least as refractory as Mg in planetary materials (50% condensation temperature ${T}_{c}^{\mathrm{Mg}}$ = 1336 K; Lodders 2003). The Th and U have 50% condensation temperatures of 1659 and 1610 K, respectively (Lodders 2003). Based on these models, along with chondrite abundances and the Earth's abundances, we expect that the Th/Mg and U/Mg ratios in a rocky exoplanet should match within tens of percent of the ratios in the star. We therefore set both the Earth-like and planet values of f_Th and f_U to 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast to U and Th, K is substantially depleted in the Earth relative to Sun and CI chondrites, with ${f}_{{\rm{K}}}^{\mathrm{Earth}}$ = 0.19 (McDonough 2003; Lodders et al. 2009). This is part of the well-known planetary volatility trend observed in the compositions of the Earth and other planets: elements less refractory than Mg (${T}_{c}\lt {T}_{c}^{\mathrm{Mg}}$) are depleted in the Earth, relative to Mg and CI chondrites, by amounts that increase with decreasing condensation temperature (McDonough 2003). The 50% condensation temperature of K is ${T}_{{\rm{c}}}^{{\rm{K}}}=1006$ K (Lodders 2003), and K is thus "moderately volatile" in comparison with Mg, Th, and U. The relative volatility of K is reflected in the range of K/Mg among chondrites, which is wider than the spread in Th/Mg or U/Mg. The solar value of K/Mg is closest to the CI value of 3.6 × 10⁻³; however, it can vary from as low as 1.3 × 10⁻³ in CV chondrites to as high as 4.7 × 10⁻³ in EH chondrites and 3.2 × 10⁻³ in EL chondrites, i.e., from 0.4 times CI in CV to 1.3 and 0.9 times CI in EH and EL chondrites. This demonstrates that fractionation of K occurs among planetary materials, although the Earth's depletion by a factor of 5 remains unexplained (see, however, Desch et al. 2020). Without knowledge of the mechanism that depletes K relative to Mg during planet formation, we cannot constrain ${f}_{{\rm{K}}}^{\mathrm{planet}}$. Changes in ${f}_{{\rm{K}}}^{\mathrm{planet}}$ relative to ${f}_{{\rm{K}}}^{\mathrm{Earth}}$ will directly change ${C}_{{\rm{K}}}^{\mathrm{planet}}$ by an equal amount. Initially, we set ${f}_{{\rm{K}}}^{\mathrm{planet}}={f}_{{\rm{K}}}^{\mathrm{Earth}}=0.19$ and discuss the consequences of variable ${f}_{{\rm{K}}}^{\mathrm{planet}}$ in the main text. The ⁴⁰K is not expected to fractionate relative to K in any planet formation scenario; therefore, it will be depleted by the same amount as bulk K between a star and planet.

Zoom In Zoom Out Reset image size

For our models, then, ${C}_{X}^{\mathrm{planet}}$ is simply a function of the ratio of (X/Mg) between the star and the Sun:

$$\begin{eqnarray}&&{C}_{X}^{{\rm{p}}{\rm{l}}{\rm{a}}{\rm{n}}{\rm{e}}{\rm{t}}}(t)=\displaystyle \frac{{\left(X/{\rm{Mg}}\right)}_{{\rm{s}}{\rm{t}}{\rm{a}}{\rm{r}}}}{{\left({\rm{X}}/{\rm{Mg}}\right)}_{{\rm{S}}{\rm{u}}{\rm{n}}}}\ast {{\rm{C}}}_{X}^{{\rm{E}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{h}}-{\rm{l}}{\rm{i}}{\rm{k}}{\rm{e}}}(t).\end{eqnarray} \tag{ A38 }$$

This method is similar to that used by Frank et al. (2014) in their determination of a cosmochemically Earth-like planet; however, our model is able to account for the system-to-system variation in HPE concentrations due to inefficient mixing within the Galaxy. For all HPE concentrations relative to the Sun, we adopt the solar composition model of Lodders et al. (2009).

Zoom In Zoom Out Reset image size

We compile measured stellar Th abundances from Unterborn et al. (2015) and Botelho et al. (2019) for a sample of 72 solar twins and analogs (Figure 1, left). Solar twins and analogs are stars of similar metallicity (that is, iron abundance), mass, and surface temperature to the Sun. Because of these similarities to the Sun, systematic uncertainties in abundance measurements due to the assumed stellar atmosphere model are minimized, particularly of trace elements like Th. The reported Th abundances do not include abundance information for Mg, instead providing only Si. To correct for this, we assume a constant Si/Mg molar ratio equal to that of the Sun (Si/Mg = 0.95; Lodders et al. 2009) for each of these stars. The distribution of stellar Si/Mg abundances is 0.9 ± 0.2 by mole (Hinkel & Unterborn 2018); thus, the uncertainty introduced in our conversion from Th/Si to Th/Mg does not drastically increase the range of planetary Th concentrations that we explore across our Monte Carlo thermal models. Assuming these abundances follow a lognormal distribution, we calculate an average current-day value (t = 12.5 Gyr after the birth of the Milky Way) of ${\left(\mathrm{Th}/\mathrm{Mg}\right)}_{\mathrm{star}}=1.21$ times our chosen solar value (Table A3) with the 95% confidence interval between 0.77 and 1.88 times solar (Figure 1, left).

Zoom In Zoom Out Reset image size

Unlike thorium, uranium has yet to be measured in young Sun-like stars. Additionally, the isotopic ratio of ²³⁵U/²³⁸U has not been measured in any system outside of the solar system. In the absence of direct observational constraints on U, we adopt Eu as its proxy. The ratios of r-process elements (i.e., Eu, Th, and U) in stars are remarkably well correlated for extremely old, metal-poor stars with r-process enhancements (Beers & Christlieb 2005; Roederer et al. 2009; Barbuy et al. 2011; Hansen et al. 2017, and references therein). The nucleosynthetic origins of third-peak r-process elements are observationally and theoretically correlated (Goriely & Arnould 2001; Frebel et al. 2007). For the Sun, log(U/Eu) ≈ −1 (by mole); that is, U is depleted by a factor of 10 relative to Eu. Given that ultra- and hyper-metal-poor stars in particular may reflect element production and enrichment from single or a few events, we would expect a similarly narrow variation in both U and Eu, with little change in their abundances relative to each other due to their coproduction. The concentration of europium, defined as Eu/Mg, is then a viable proxy for predicting the system-to-system variation in U abundances. Europium is also as refractory as U, with a 50% condensation temperature of 1356 K (Lodders 2003). Therefore, the concentration of Eu as determined from Eu/Mg is, like U, not expected to vary between host star and exoplanet by more than tens of percent, i.e., less than the total observed range. For comparison, the bulk silicate Earth has Eu/Mg by mole of 10.5 × 10⁻⁸ (McDonough 2003), while the Sun's value is nearly the same at 9.6 × 10⁻⁸ (Lodders et al. 2009). We therefore adopt Eu/Mg as a proxy for the distribution of U stellar abundances and set f_Eu to 1 for both the planet and Earth-like values. Both ²³⁵U and ²³⁸U isotope abundances are sampled independently from the Eu distribution.

Though still difficult to measure, Eu is observed in stars with reasonable frequency. We adopt a data set of 2040 FGK stars with measured Eu and Mg abundances in the Hypatia catalog (Figure 1, middle; Hinkel et al. 2014). We assume this data set to be an upper limit of the range of Eu abundances relative to the Sun, as these measurements are inherently less precise than abundance determinations from solar twins and analogs. Assuming that these abundances follow a lognormal distribution and that U/Eu is constant throughout the Galaxy (U/Eu = 1/10), we calculate an average current-day (Eu/Mg)_star =(U/Mg)_star = 0.93 times the solar ratio (Table A3), with 95% of of our sample falling between 0.45 and 1.92 times solar (Figure 1, middle).

Bulk K/Mg ratios show a larger range of variation than Th/Mg and Eu/Mg ratios (Figure 1, right). There are 179 FGK stars with both K and Mg reported in the Hypatia catalog (Hinkel et al. 2014). Assuming that these abundance ratios follow a lognormal distribution, we calculate an average current-day (K/Mg)_star of 1.13 times the solar ratio (Table A3), with 95% of all data falling between 0.35 and 3.63 times solar (Figure 1, right). Only the single isotope ⁴⁰K is radioactive, but no data exist for ⁴⁰K/K ratios outside of the solar system. For this model setup, a variation in a planet's ⁴⁰K/K will effectively act as an increase or decrease in bulk K, similar to our discussion of f_K above. We discuss the consequences of the variable volatility of K on our model results in the main text and Appendix B.

By taking our data from the Hypatia catalog, we are implicitly combining abundances from different sources with different measurement uncertainties. Because we adopt the median abundance value if multiple sources are available for the same star, we likely overestimate the range of any abundance ratio in Figure 1.

Our Monte Carlo models described below randomly sample within each distribution of Figure 1 independently. There is observational evidence from solar twins that stellar Th/Eu is roughly constant through Galactic time (Botelho et al. 2019), suggesting that these abundances are correlated. It is not observationally known, however, whether this extends to U/Th in metal-rich stars. Both U and Th are produced via the r-process, and thus their abundances in stars do correlate somewhat. However, ⁴⁰K is produced via the s-process, and its correlation with the r-process elements is not known. Bulk K is produced via explosive oxygen burning (Shimansky et al. 2003) and has been observationally shown to correlate with the α-elements (e.g., Mg; Zhang et al. 2006); it is therefore unlikely to correlate with Eu, Th, or U through Galactic time. The models of Frank et al. (2014) capture each of these behaviors; thus, our treatment of ${C}_{X}^{\mathrm{planet}}(t)$ somewhat captures the correlation between Eu, U, and Th (Figure A1). Our treatment of the system-to-system variations implied from Figure 1 causes our model to explore some areas of U and Th parameter space that are unlikely. Given that we find these elements to have a minor effect on the longevity of mantle degassing, however, these inclusions of Eu, Th, and U correlations are not likely to substantially change our determinations of Age${}_{\max }$.

For this work, we adopt a Monte Carlo method for determining the degassing lifetime of planets with a fixed size, CMF, and formation age (in terms of time after Galaxy formation). In each Monte Carlo suite, ∼10⁴–10⁵ models are run. In each run, C_X ^planet is determined using Equation (A38) by independently sampling from the lognormal distributions of each HPE as shown in Figure 1 and multiplying by the ${C}_{X}^{\mathrm{Earth}-\mathrm{like}}$ corresponding to its formation time after the birth of the Milky Way (t = 0; Figure A1). This method allows us to simulate both the GCE of the HPEs over Galactic history and the system-to-system variation of these HPEs due to their stochastic distribution due to inefficient mixing of these elements throughout the Galaxy.

Unlike the HPE abundances, which can be constrained by stellar observations, other factors that influence a planet's thermal evolution cannot be constrained observationally. We account for variation in two of the most important of these factors, the initial mantle temperature and mantle reference viscosity, in our models by randomly sampling from uniform distributions across reasonable uncertainty ranges. For the initial mantle temperature, we use a range of 1700–2000 K, and for the mantle reference viscosity, we use a range of 10¹⁹–10²³ Pa s. The uniform distribution for the reference viscosity is sampled in log space, or, in other words, we sample from a uniform distribution of ${\mathrm{log}}_{10}({\mu }_{\mathrm{ref}})=19-23$. The range of mantle viscosities considered covers 2 orders of magnitude above and below the typical estimates for the average viscosity of Earth's mantle based on postglacial rebound studies (e.g., Mitrovica & Forte 2004). The initial mantle temperature is not well constrained for Earth or any planet, but our assumed range covers the spread typically used in studies of solar system planets (e.g., Hauck & Phillips 2002; Fraeman & Korenaga 2010; Morschhauser et al. 2011; Breuer & Moore 2015; Foley et al. 2020). With the reference viscosity, initial mantle temperature, and radiogenic heat production rate set, the initial stagnant-lid thickness is calculated assuming that the conductive heat flux at the base of the lid matches the advective heat flux supplied by the convecting mantle. The last remaining parameter to set is the total carbon budget of the mantle and surface reservoirs, C_tot.

For our suites of Monte Carlo models, the concentration of CO₂ in the mantle at the model start time is set to ∼5 × 10⁻³ wt%, regardless of planet mass or other assumed properties. We examined the effects of higher or lower mantle C contents and found that the lifetime of mantle degassing was insensitive to the total mantle C concentration above ∼10⁻⁵ wt%, in agreement with Foley & Smye (2018). Below this threshold, the mantle is more likely to exhaust its entire C budget, causing degassing to end, even if there is a sufficient HPE budget to support volcanism for longer. Above ∼10⁻⁵ wt%, there is sufficient mantle CO₂ to support degassing until the planet's dwindling heat budget no longer supports mantle melting and volcanism. We also assume no CO₂ in the atmosphere at the model start time and, since the models start with no crust present, no CO₂ in the crust either. Whether a planet's C would initially lie entirely in the mantle (as we assume), the atmosphere, or a combination of the two is not known. However, Foley (2019) found that the initial distribution of C between surface and interior does not significantly affect subsequent climate evolution, at least when liquid water is present and silicate weathering can occur. For an Earth-sized planet, ≈5 × 10⁻³ wt% CO₂ scales to a total C budget of 5 × 10²¹ mol; this value is about a factor of 2 lower than the estimate for the Earth given by Sleep & Zahnle (2001).

Since models are presented in terms of a CO₂ concentration in wt%, larger planets will have larger total C budgets than smaller planets for the same CO₂ concentration. In calculating rates of CO₂ outgassing, our models assume that the concentration of CO₂ in the source region to mantle melts is set by the (time-evolving) bulk mantle CO₂ concentration. We use melt-solid partition coefficients to then calculate the CO₂ concentration in the resulting melt and hence the rate of CO₂ outgassing. In doing so, we implicitly assume that the mantle oxidation state of the planets we model is the same as the present-day Earth. A more reduced mantle would favor production of more reduced gases at the expense of CO₂ and lower the CO₂ outgassing rate (e.g., Tosi et al. 2017). More reduced conditions also reduce the solubility of C species in mantle melt (e.g., Grewal et al. 2020). A more reduced mantle would thus be equivalent to lower mantle CO₂ concentrations in our models.

All parameters varied in our model are shown in Table A1.

Table A2. Sample of Likely Solid Rocky Planets

	Mass	Density	Age	Teq	% Probability	% Probability
Planet	(M⊕)	(g cm−3)	(Gyr)	(K)	Age ≤ Age${}_{\max }^{\mathrm{Avg}}$	Age ≤ Age${}_{\max }^{\mathrm{UL}95 \% \mathrm{CI}}$
K2-36 b	3.9 ± 1.1	7.3 ± 2.1	1.4 ± 0.45	1224	100	100
Kepler-80 e	2.6 ± 0.75	6.5 ± 1.9	2 ± 1	629	70	99
Kepler-65 d	4.14 ± 0.8	6.5 ± 1.2	2.9 ± 0.7	1117	56	100
Kepler-138 c	1.97 ± 1.5	6.3 ± 4.9	4.68 ± 4.17	402	26	44
Kepler-105 c	4.6 ± 0.9	11.2 ± 2.2	3.17 ± 0.6	997	45	100
Kepler-345 c	2.2 ± 0.9	7.0 ± 2.9	2.75±1.7	575	42	83
Kepler-197 c	5.3 ± 3.1	15.6 ± 9.2	5.37 ± 3.1	930	2	63
TRAPPIST-1 b	1.374 ± 0.07	5.4 ± 0.3	7.6 ± 2.2	401	1	5
TRAPPIST-1 g	1.321 ± 0.038	5.0 ± 0.1	7.6 ± 2.2	199	1	5
TRAPPIST-1 c	1.308 ± 0.056	5.4 ± 0.1	7.6 ± 2.2	341	1	5
TRAPPIST-1 f	1.039 ± 0.031	5.0 ± 0.1	7.6 ± 2.2	219	0	4
Kepler-36 b	3.83 ± 0.1	6.3 ± 0.2	4.79 ± 0.65	978	0	85
HD 219134 b	4.74 ± 0.19	6.3 ± 0.03	11 ± 2.2	1015	0	1
HD 219134 c	4.36 ± 0.22	6.9 ± 0.4	11 ± 2.2	782	0	1
Kepler 93 b	4.54 ± 0.85	6.5 ± 1.2	6.6 ± 0.9	1037	0	17
Kepler-68 c	2.04 ± 1.75	14.4 ± 12.3	6.31 ± 0.82	941	0	6
HD 136352 b	4.62 ± 0.45	7.2 ± 0.7	8.2 ± 3.2	911	5	22

Download table as:  ASCIITypeset image

As the Galaxy aged, HPEs were produced via supernova processes, while those already existing decayed. This is all while the nonradioactive rocky planet–building elements were continuously being produced, thus diluting the concentration of HPEs as the Galaxy aged (Figure A1). This, in turn, caused the average Q₀, and thus degassing lifetimes, to also decrease as the Galaxy aged (Figure A1, middle). The exact makeup of a planet's HPE budget also affected the degassing lifetimes. Those planets older than ∼8 Gyr had initial heat budgets composed primarily of ²³⁵U. For planets younger than 8 Gyr, ⁴⁰K instead became the dominant component of its initial HPE budget. The shorter half-life of ²³⁵U (704 Myr) compared to ⁴⁰K (1.25 Gyr) causes older planets to release their internal heat more rapidly. As a result, the oldest, ²³⁵U-dominant planets degas for only ∼1 Gyr longer than the ⁴⁰K-dominant planets formed at the same time as the Earth, despite having nearly double the Q₀ and five times as much ²³⁵U.

These results suggest that as the Galaxy evolves, newer planets may form with a lower HPE budget than older ones (Figure A1). While we do not extrapolate our results into the future, this potential decrease in HPE abundance would cause degassing lifetimes for stagnant-lid exoplanets to also decrease. In order to counteract the decay and subsequent dilution of radionuclides as the Galaxy evolves, the likelihood of a rocky exoplanet supporting a temperate climate may be more reliant on the local emplacement of HPE-rich material from supernova sources to boost its radiogenic heat budget above the average from GCE. Exploring these effects is beyond the scope of this paper, but our results point to the need for more work exploring the spatial and temporal distribution of the long-lived HPEs due to GCE.

We first consider possible sources of error or bias in our estimated probability distributions for the mantle abundances of the four major HPEs in rocky exoplanets as derived from stellar chemistry observations. Of the four HPEs we considered, a rocky exoplanet's mantle concentration of ⁴⁰K is dependent on both the disk's bulk abundance of ⁴⁰K and depletion of K relative to the refractory rock-forming elements during planetary formation. Unlike Th and U, K is moderately volatile, meaning that there could be significant differences between the measured stellar abundances and the resulting planetary abundances. In our solar system, accretion of chondrites themselves cannot explain the depletion of K by the factor of 5 seen in Earth relative to the Sun. Moderately volatile elements within chondrules (igneous melt spherules), however, appear to be depleted with respect to the matrix grains by factors comparable to the requisite factor of 5 depletion of K relative to Mg between the Earth and Sun. Perhaps the preferential accretion of chondrules (rather than chondrites) by Earth may explain its depletion (Desch et al. 2020); however, no first-principle models exist that quantitatively predict the degree of depletion of the volatile elements in planetary materials (e.g., see review by Alexander 2005). Whether preferential melting, vaporization, or differentiation play a role in the degree of K volatilization during planet formation is an open but critical question in understanding whether a rocky exoplanet is likely to be temperate today. Until such time as a theory is developed to quantify the range of possible fractionation, all that can be said is that unlike the Th/Mg and U/ Mg ratios, the K/Mg ratio of a planet cannot be expected to match its host star, and a volatility/fractionation factor must be assumed for exoplanetary systems. If planets were composed of purely chondritic material, stellar abundances matching planetary abundances might be the case, but based on our solar system, it is more likely that the planet has a lower K/Mg than its host star. This degree of volatilization then remains unknown. Given its effect on the degassing lifetime of rocky exoplanets, work by both the meteoritics and planetary formation communities is critically needed to quantify this parameter in order to better constrain the evolution of rocky exoplanets.

Additionally, the ⁴⁰K/K ratio for other stars, and hence planetary systems, is unconstrained. In our previous models, we simply assumed the same depletion factor between the Earth and Sun for all planets and hence the same volatility and ⁴⁰K/K ratio as the Earth. However, system-to-system variation in these properties would lead to an even wider distribution of ⁴⁰K compared to those in Figure 1. The probability of a planet acquiring a larger ⁴⁰K concentration could therefore be higher than we previously estimated.

To assess how changes in volatility and ⁴⁰K/K affect Age${}_{\max }$ (Figure A5), we reran our models for determining the degassing lifetime as a function of planet formation time and mass (Figure 2) and arbitrarily changed the abundance of ⁴⁰K (f_x ^K; Equation (A37)) by factors of 0.5 and 2. Halving a planet's ⁴⁰K abundance decreases Age${}_{\max }^{\mathrm{Avg}}$ by ∼500–700 Myr and Age${}_{\max }^{\mathrm{UL}95 \% \mathrm{CI}}$ by ∼600–750 Myr, with this difference increasing as the mass increases (Figure B1). Doubling the K abundance yields increases in Age${}_{\max }^{\mathrm{Avg}}$ by ∼700–850 Myr and Age${}_{\max }^{\mathrm{UL}95 \% \mathrm{CI}}$ by 0.7–1 Gyr. While observations of a host star's K/Mg abundance can allow us to estimate the degassing lifetime of a rocky exoplanet, differences in the volatilization and/or the system's ⁴⁰K/K abundances compared to Earth will have a direct effect on the location of Age${}_{\max }$. Next, we will explore the degree of potential K volatility and the likelihood of systems having non-Earth ⁴⁰K/K ratios for comparison to these estimates of ⁴⁰K-dependent Age${}_{\max }$.

Zoom In Zoom Out Reset image size

An efficient method of ⁴⁰K enrichment might be supernova injection of material into the protoplanetary disk. Supernova remnants frequently exhibit large-scale bipolar or unipolar morphologies (e.g., Leonard et al. 2000; Ouellette et al. 2007; Larsson et al. 2013; Moranchel-Basurto et al. 2017; Reilly et al. 2017), as well as filaments and clumpy "bullets" of ejected material on smaller scales (e.g., Willingale et al. 2002; Fesen & Milisavljevic 2016; Zhou et al. 2016; García et al. 2017). These sorts of small-scale structures are seen in hydrodynamic simulations (e.g., Hammer et al. 2010; Ellinger et al. 2012; Wongwathanarat et al. 2015). As a result, star- and planet-forming structures near the supernova can be enriched by bullets of material that represent the abundances in a specific part of the explosion, rather than the average yields. This material can be injected as gas into a cloud core or as dust into a protoplanetary disk (Young et al. 2009; Ellinger et al. 2012; Pan et al. 2012). It is possible to estimate the maximum enhancement that such an injection can yield by comparing the amount of ⁴⁰K in an ejecta clump to the total inventory in the protoplanetary material.

Based on ALMA observations, the dust mass in M dwarf disks is of order 10⁻⁵ M_⊙, and the total disk mass for a dust-to-gas ratio of 1:100 is 10⁻³ M_⊙ (Ward-Duong et al. 2018). If a 10⁻³ M_⊙ disk initially contains a solar K concentration of 3.7 ppm (3.7 μg K g^–1 of disk) and the primordial solar ⁴⁰K/K (0.14%) of Lodders (2003), the masses of the bulk K and ⁴⁰K in the disk are 3.7 × 10⁻⁹ and 5.17 × 10⁻¹² M_⊙, respectively. For an individual 15 M_⊙ supernova progenitor (Figure B2; Ellinger et al. 2012), we calculate a peak mass fraction of bulk K and ⁴⁰K potentially injected into the disk of ∼10^−3.6 and ∼10^−4.7, respectively (Figure B3). This corresponds to ⁴⁰K/K ∼ 10^−1.1 (8%). The average mass of a bullet in our model is 10⁻⁵ M_⊙. For a 10⁻⁵ M_⊙ bullet mass, then, we estimate that ∼2.5 × 10⁻¹¹ M_⊙ of K and 2 × 10⁻¹² M_⊙ of ⁴⁰K are delivered to the disk via the injection of supernova material. Injecting a single 10⁻⁵ M_⊙ bullet into our disk raises the total mass of K in the disk by <1% and ⁴⁰K by ∼39%. This would raise the disk's primordial ⁴⁰K/K to ∼0.19%. The injection of this superenriched supernova-derived ⁴⁰K material will account for an increase in the protoplanetary disk's ⁴⁰K/K by only ∼36%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the disk mass increases above 10⁻³ M_⊙, the relative mass contribution of high-⁴⁰K/K material from these injections to the protoplanetary disk begins to diminish. For the minimum-mass solar nebula of 0.013 M_⊕ (Hayashi 1981; Kuchner 2004), the contribution from a single injection increases the disk's ⁴⁰K/K by <1%. Thus, while ⁴⁰K/K can increase somewhat due to the direct injection of supernova material, this mechanism is likely to only affect the radiogenic heat budgets of planets forming in low-mass, M dwarf disks, and only marginally, leading only to a <500 Myr increase in Age${}_{\max }$ (Figure B3). Larger bullet masses on the order of 10⁻⁴ M_⊙, however, are possible (Young et al. 2009; Ellinger et al. 2012), which would allow for a significant enrichment for higher-mass disks. The relative occurrence rates of these high-mass bullets are unconstrained and beyond the scope of this paper. Thus, we offer only a conservative estimate of ⁴⁰K enrichment via supernova injection. We argue then that supernova events are not likely to increase a planet's radiogenic heat budget and degassing lifetime enough to drastically change Age${}_{\max }$ for all but those planets orbiting low-mass stars. Instead, a system's bulk K abundance and the degree of K volatilization are the most likely first-order sources for increasing a rocky exoplanet's radiogenic heat budget and degassing lifetime. In this sense then, observations of the host star's K/Mg coupled with modeling of the range of volatility are the most likely avenues for quantifying whether a planet contains sufficient K to be actively degassing today.

We estimate how increasing the mantle water content would influence Age${}_{\mathrm{Max}}$, taking into account the competing effects of water content on mantle viscosity and solidus. Viscosity has been found to scale as ${(1/{\chi }_{m})}^{r}$, where χ_m is the mantle water concentration, expressed here as the water mass fraction, and r is an experimentally determined exponent, r = 0.7–1.2, according to the compilation in Hirth & Kohlstedt (2003). Drier mantles would thus have higher viscosities, up to a point (there is a threshold water content where the material is effectively dry, and further decreases in water content will not affect viscosity), and wetter mantles would have lower viscosities. At the same time, water content also affects the mantle solidus, with melting temperatures decreasing as mantle water content increases (e.g., Kushiro et al. 1968). These two effects counteract each other: a drier mantle has a higher viscosity, acting to extend the degassing lifetime and Age_max, but also a higher solidus temperature, which inhibits melting and acts to lower the degassing lifetime and Age_max.

Ultimately, to fully capture these competing effects, new models incorporating the combined effects of water on viscosity and solidus, and how mantle water content evolves over time due to in- and outgassing from the mantle (e.g., McGovern & Schubert 1989; Crowley et al. 2011; Spaargaren et al. 2020), would be needed. However, we can provide a first-order estimate using two additional model suites where the solidus temperature is first lowered by 100 K everywhere, then by 200 K, and the mantle reference viscosity is systematically varied as before. From the results of these model suites, we find that, on average, the degassing lifetime increases by ≈0.43 Gyr per 100 K decrease in solidus temperature for a 1 M_⊕ planet and ≈0.60 Gyr for 6 M_⊕ (Figure B4). We then use the diffusion creep viscosity flow laws from Hirth & Kohlstedt (2003) and the solidus parameterization from Katz et al. (2003) to estimate how Age${}_{\max }$ would change with increasing water content. We start from the water content where the wet viscosity flow law is equal to the dry viscosity flow law and increase the water content from this point up to ≈0.3 wt% (Figure B5). For comparison, the Earth's mantle water content is estimated to be 10⁻² wt% (Ohtani 2019). We do this for a range of water viscosity dependency exponents, r.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For both 1 and 6 M_⊕ planets, increasing the water content first decreases Age${}_{\mathrm{Max}}^{\mathrm{Avg}}$ in comparison to a dry-planet baseline. Here the effect of decreasing the mantle reference viscosity due to a higher water content dominates because the solidus temperature does not decrease significantly from that of dry peridotite until the water content approaches ≈0.05–0.1 wt%. When the water content increases beyond this point (χ_m ≳ 0.05–0.1 wt%), the decrease in solidus temperature begins to dominate over the decrease in viscosity, and the change in Age_max compared to a dry-planet baseline begins to increase. The higher the r, the more the viscosity effect dominates, as mantle viscosity is more sensitive to water content in this case. For a 1 M_⊕ planet, Age${}_{\mathrm{Max}}^{\mathrm{Avg}}$ could decrease by up to ≈1.4 Gyr when χ_m = 0.05 wt% and r = 1.2. Meanwhile, Age${}_{\mathrm{Max}}^{\mathrm{Avg}}$ could increase by up to ≈1.6 Gyr for χ_m = 0.3 wt% and r = 0.7. The effects are larger for a 6 M_⊕ planet, with the biggest decrease in Age${}_{\mathrm{Max}}^{\mathrm{Avg}}$ being ≈2.1 Gyr and the biggest increase being ≈2.3 Gyr. Note again that these estimates assume fixed mantle water concentrations, when in reality, the water concentration would evolve over time as a result of in- and outgassing. On a stagnant-lid planet, ingassing will be limited (e.g., Spaargaren et al. 2020), so the mantle water content is likely to decrease over time; this means the changes to Age${}_{\max }$ presented here are likely upper limits.

Mantle water content, then, is an important parameter that can have a large effect on ${{\rm{Age}}}_{\max }$. Should a planet's mantle form with an Earth-like concentration (10⁻² wt%), Age${}_{\mathrm{Max}}^{\mathrm{Avg}}$ is roughly the same as or lower than the dry mantle model (Figure A5). It is only those planets with mantles enriched in water relative to the Earth by a factor of 10–100 where Age ${}_{\mathrm{Max}}^{\mathrm{Avg}}$ increases, according to our simple estimates. Planets forming with this much water is a possibility (e.g., Raymond et al. 2004; Unterborn et al. 2018). However, despite a lower mantle solidus promoting melting on such water-rich planets, there are other factors that could hamper their habitability and the possibility of detecting biosignatures on such planets, should they be inhabited. Exposed land may be essential for life on Earth, as weathering of subaerial land supplies critical nutrients to the oceans (e.g., Maruyama et al. 2013; Dohm & Maruyama 2015), and some leading theories for the origin of life rely on wet–dry cycles and therefore require exposed land (e.g., Bada & Korenaga 2018). For an Earth-mass planet, complete flooding of the surface is expected to occur when the mass of surface water is roughly that of 8–40 oceans. This equates to ∼0.2% of the mass of the planet, with this mass-fraction threshold decreasing with increasing planet mass due to higher gravity limiting surface topography (Kite et al. 2009; Cowan & Abbot 2014).

Furthermore, as the surface water content increases, the pressure at the water–rock interface increases, forcing melting to occur at higher temperatures, limiting mantle degassing (Kite et al. 2009; Cowan & Abbot 2014). This pressure effect would then counteract the effect of mantle water content on the solidus temperature and therefore potentially limit the increase in Age${}_{\max }$ at higher mantle water concentrations than in Figure B5 (where ocean bottom pressure is assumed to be negligible for mantle melting). In fact, for a very thick ocean layer, pressures may be high enough to shut down degassing entirely (Kite et al. 2009; Krissansen-Totton et al. 2021) by preventing melting or forcing melting to occur at such high pressures that the melt would no longer be buoyant. Moreover, even if life were to develop on an ocean-covered planet with mantle degassing, the detectability of life on these planets is more difficult (Glaser et al. 2020). This is because, for planets with enough surface water to submerge the continental crust, the weathering rate of the biocritical element phosphorus is significantly reduced. With less bioavailable P, the biogenic flux of O₂ is reduced, becoming comparable to the abiotic O₂ flux (Glaser et al. 2020). Even for those water worlds where the pressure at the water–mantle boundary allows for mantle degassing, as these planets age, degassing will wane. In this case, the source of reductants from the mantle will also decrease, allowing abiotic O₂ to build up in the atmosphere, potentially masking any biotic source of O₂ (Krissansen-Totton et al. 2021). Each of these factors will limit our ability to distinguish between biotic or abiotic O₂ through atmospheric observations, even though the planet may be habitable and hosting life.

Therefore, it is unclear whether high water concentrations (≳0.2 wt%) would actually promote long-lived mantle degassing, temperate climates, and habitability. Instead, for water contents ≲0.2 wt%, where continents may still be exposed on small exoplanets, increasing water content generally leads to lower degassing lifetimes and a lower Age${}_{\max }$. While more work is needed to quantify these differences as a function of planet mass, our results shown in Figure A5 may represent the most optimistic case for Age${}_{\max }$ with regard to mantle water content.