EFFECT OF SURFACE-MANTLE WATER EXCHANGE PARAMETERIZATIONS ON EXOPLANET OCEAN DEPTHS

To date, the suite of observed exoplanets from Kepler has proven that Earth-sized planets are common in the universe ($\approx 0.16$ per star, Fressin et al. 2013; Morton & Swift 2014). Though we do not yet have a detailed understanding of the atmospheric composition of an extrasolar terrestrial planet, spectra of many extrasolar gas giants (Kreidberg et al. 2015; Sing et al. 2015) and a smaller Neptune-sized planet (Fraine et al. 2014) have shown that water is likely abundant in other solar systems. Calculations of volatile delivery rates to terrestrial planets via planetesimals (e.g., Raymond et al. 2004; Ciesla et al. 2015) have shown that planets can have a wide range of initial water fractions, with some planets being 1% water by mass or more. Both observations and simulations, hence, point toward the likelihood that terrestrial planets are also born with abundant water. However, the intertwined effects of climate (Kasting et al. 1993) and mantle-surface volatile interchange (Hirschmann 2006; Cowan 2015) determine whether there is abundant liquid water on the present-day surfaces of terrestrial exoplanets. Additionally, atmospheric escape (especially early in the atmospheric evolution) can cause loss of copious amounts of water (Ramirez & Kaltenegger 2014; Luger & Barnes 2015; Tian & Ida 2015; Schaefer et al. 2016), with $\gtrsim 10$ Earth oceans possibly lost from planets in the habitable zone of M-dwarfs.

The extent of the traditional habitable zone is determined by the continental silicate weathering thermostat (Kasting 1988), in which silicate minerals react with ${\mathrm{CO}}_{2}$ and rainwater to produce carbonates (Walker et al. 1981). Silicate weathering is extremely efficient at stabilizing the climate because the process runs faster with increasing temperature. This is due to faster reaction rates and increased rain in warmer climates. However, the silicate weathering thermostat itself depends on the surface water abundance.

If there is no surface water, the silicate weathering thermostat cannot operate due to the lack of reactants, and if the planet surface is completely water-covered the negative feedback does not operate unless seafloor weathering is also temperature-dependent (Abbot et al. 2012). Note that even if seafloor weathering is temperature-dependent, it might be insufficient to stabilize the climate (Foley 2015). A waterworld state is likely stable (Wordsworth & Pierrehumbert 2013) since water loss rates would be low because the atmosphere would be ${\mathrm{CO}}_{2}$-rich due to the lack of silicate weathering feedback. However, if water loss rates remain high due to a large incident stellar flux, it is possible that brief exposures of land can allow for a "waterworld self-arrest" process in which the planet adjusts out of the moist greenhouse state (Abbot et al. 2012). This can occur if the timescale for ${\mathrm{CO}}_{2}$ drawdown by the silicate weathering feedback is shorter than the timescale for water loss to space, which is probable for Earth parameters.

From the above discussion, we conclude that although waterworlds are by definition in the habitable zone (having liquid water on the surface), they may not actually be temperate and conducive to life. It is instead likely that waterworlds are less habitable than worlds with continents, and so determining whether or not waterworlds are common is important. To determine whether or not waterworlds should be common, we must look to the deep-water cycle, that is, the mantle-surface interchange of water over geologic time.

To understand the deep-water cycle on exoplanets, we look to Earth as an analogue because it is the only planet known with continuous (not episodic) mantle-surface water interchange due to plate tectonics. On present-day Earth, water is largely expelled from the mantle to the surface (degassed) through volcanism at mid-ocean ridges and volcanic arcs (Hirschmann 2006). Water is lost from the surface to the mantle (regassed) through subduction of hydrated basalt. The relative strength of regassing and degassing determines whether the surface water abundance increases or decreases with time.

It has long been suggested that Earth's surface water fraction is in an effective steady state (McGovern & Schubert 1989; Kasting & Holm 1992), due to the constancy of continental freeboard since the Archean (∼2.5 Gya). However, this may simply be due to isostasy, that is, the adjustment of the continental freeboard under varying surface loads (Rowley 2013; Cowan & Abbot 2014). A more convincing argument is that the degassing and regassing rates on Earth are high enough that if they did not nearly balance each other, the surface would have long ago become either completely dry or water-covered (Cowan & Abbot 2014). However, some studies of volatile cycling on Earth that utilized parameterized convection to determine the upper mantle temperature and hence the degassing and regassing rates have not found such a steady state (McGovern & Schubert 1989; Crowley et al. 2011; Sandu et al. 2011). If the Earth is indeed near steady state, this mismatch could be because there are many secondary processes, e.g., loss of water into the transition zone (Pearson et al. 2014) and early mantle degassing (Elkins-Tanton 2011), that are difficult to incorporate into a simplified volatile cycling model. Also, it is possible that our understanding of what processes control the release of water from the mantle and return of water to it via subduction is incomplete.

Using the maximum allowed fraction of water in mantle minerals (Hauri et al. 2006; Inoue et al. 2010), Cowan & Abbot (2014) estimate that Earth's mantle water capacity is $\approx 12$ times the current surface water mass. However, measurements of the electrical conductivity of Earth's mantle (Dai & Karato 2009) have found only $\sim 1\mbox{--}2$ ocean masses of water in the mantle, which is much less than the maximally allowed value. This measurement may vary spatially (Huang et al. 2005) and by method (Khan & Shankland 2012), but it is likely constrained to within a factor of a few. This implies that dynamic effects lead to a first-order balance between degassing and regassing on Earth, rather than the surface water complement being in steady state simply because the mantle is saturated.

Using a steady-state model wherein the degassing and regassing of water is regulated by seafloor pressure, Cowan & Abbot (2014) applied our knowledge of Earth's deep-water cycle to terrestrial exoplanets. They showed that terrestrial exoplanets require large amounts ($\sim 1 \%$ by mass) of delivered water to become waterworlds. Applying a time-dependent model and including the effects of mantle convection, Schaefer & Sasselov (2015) found that the amount of surface water is strongly dependent on the details of the convection parameterization. These works rely on other planets being in a plate-tectonic regime similar to Earth. However, it is important to note that there is debate about whether or not plate tectonics is a typical outcome of planetary evolution (e.g., O'Neill & Lenardic 2007; Valencia et al. 2007a; Valencia & O'Connell 2009; Korenaga 2010), potentially because plate tectonics is a history-dependent phenomenon (Lenardic & Crowley 2012). In this work, we also assume plate tectonics. We do so because our understanding of habitability is most informed by Earth and it enables us to examine how processes that are known to occur on Earth affect water cycling on exoplanets. As a result, we assume that continents are present, and that isostasy determines the depths of ocean basins. In the future, exploring other tectonic regimes (e.g., stagnant lid) may be of interest to exoplanet studies and potential investigations of Earth's future evolution (Sleep 2015).

The studies of volatile cycling on super-Earths discussed above used drastically different approaches, with Cowan & Abbot (2014) applying a two-box steady-state model of volatile cycling, and Schaefer & Sasselov (2015) extending the time-dependent coupled volatile cycling-mantle convection model of Sandu et al. (2011) to exoplanets. As a result, these works made different assumptions about which processes control water partitioning between ocean and mantle. The degassing parameterization of Cowan & Abbot (2014), based on the model of Kite et al. (2009), utilized the negative feedback between surface water inventory and volatile degassing rate that results from pressure reducing degassing. Their regassing rate was also related to the surface water inventory, using the prediction of Kasting & Holm (1992) that the hydration depth increases with increasing surface water abundance up to the limit where the hydration depth is equal to the crustal thickness. Meanwhile, the degassing and regassing parameterizations of Schaefer & Sasselov (2015) were both related directly to the mantle temperature, with the degassing rate determined by the abundance of water in melt and the regassing rate set by the depth of the hydrated basalt (serpentinized) layer, which is determined by the depth at which the temperature reaches the serpentinization temperature.

In this work, we seek to identify how different assumptions about regassing and degassing determine the surface water mass fraction. To do so, we utilize simplified models of convection and volatile cycling that separately incorporate the key features of both the Cowan & Abbot (2014) and Schaefer & Sasselov (2015) volatile cycling parameterizations. The latter model builds upon the analytic work of Crowley et al. (2011), who developed an analytic model that captures the key processes in the numerical models of Sandu et al. (2011) and Schaefer & Sasselov (2015). However, here we further simplify and also non-dimensionalize the Crowley et al. (2011) model, enabling us to elucidate the dependencies of water abundance on mantle temperature and planetary parameters. We then combine the models of Cowan & Abbot (2014) and Schaefer & Sasselov (2015), utilizing surface water budget-dependent degassing and temperature-dependent regassing. We do so because it is likely the most physically relevant choice, as temperature affects serpentinization depths (and resulting regassing rates) more directly than seafloor pressure. Additionally, temperature-dependent degassing would become small at late times while seafloor pressure-dependent degassing would not, and it has been shown by Kite et al. (2009) that degassing should be pressure-dependent. This is more in line with the approximate steady-state water cycling that Earth is currently in, as if both regassing and degassing are temperature-dependent, regassing will dominate at late times. We find that the choice of volatile cycling parameterization greatly impacts the end-state surface water mass reservoir. We also find that, regardless of volatile cycling parameterization, the water partitioning reaches a steady state after a few billion years of evolution due to the cooling of the mantle to below the melting temperature, which causes the effective end of temperature-dependent degassing and regassing.

This paper is organized as follows. In Section 2, we describe our parameterized convection model and the various volatile cycling parameterizations we explore, along with the consequences these have for the temporal evolution of mantle temperature and water mass fraction. Detailed derivations of the volatile cycling models can be found in Appendix A. In Section 4, we explore where, in water mass fraction–planet mass parameter space, each volatile cycling model predicts the waterworld boundary to lie. We discuss our results in Section 5, performing a sensitivity analysis of the waterworld boundary on key controlling parameters, comparing this work to previous works, and discussing our limitations and potential avenues for future work. Importantly, we also show how our model with pressure-dependent degassing and temperature-dependent regassing could, in principle, be observationally distinguished from the models of Cowan & Abbot (2014) and Schaefer & Sasselov (2015). Lastly, we express our conclusions in Section 6.

Parameterized convection models utilize scalings from numerical calculations to relate the Nusselt number (the ratio of outgoing heat flux from the mantle to that which would be conducted if the entire mantle were not convecting) to the Rayleigh number of the mantle (Turcotte & Schubert 2002). We here consider a standard boundary-layer convection model, as in McGovern & Schubert (1989) and Sandu et al. (2011), with a top boundary layer of thickness δ and one characteristic mantle temperature T. Figure 1 shows a schematic of the temperature profile relevant for this convection parameterization. We can determine the boundary-layer thickness by the depth at which the boundary layer peels away and convects. To zeroth order, this peel off occurs where the timescale for advection of the boundary layer is shorter than the timescale for heat to diffuse out of the boundary layer. The time it takes for the boundary layer to overturn via advection is

$$\begin{eqnarray}&&{\tau }_{\mathrm{over}}\approx \displaystyle \frac{\eta (T,x)}{{\rm{\Delta }}\rho g\delta },\end{eqnarray} \tag{ 1 }$$

where $\eta (T,x)$ is the temperature and mantle water fraction-dependent viscosity of the boundary layer (the viscosity parameterization will be discussed further in Section 2.1.2), $g=g(M)$ is gravity (see Section 2.1.1 for how g and other variables scale with planet mass), and we take the density contrast ${\rm{\Delta }}\rho \approx \alpha {\rho }_{{\rm{m}}}{\rm{\Delta }}T$, where α is a characteristic thermal expansivity, ${\rm{\Delta }}T$ is the temperature contrast across the boundary layer, and ${\rho }_{{\rm{m}}}$ is the density of the upper mantle. The heat diffusion timescale is then

$$\begin{eqnarray}&&{\tau }_{\mathrm{diff}}\approx \displaystyle \frac{{\delta }^{2}}{\kappa },\end{eqnarray} \tag{ 2 }$$

where κ is the thermal diffusivity of the boundary layer (assumed equal to that of the upper mantle). Taking the ratio of Equations (1) and (2) defines the local boundary-layer Rayleigh number

$$\begin{eqnarray}&&{\mathrm{Ra}}_{\mathrm{loc}}=\displaystyle \frac{\alpha {\rho }_{{\rm{m}}}g{\rm{\Delta }}T{\delta }^{3}}{\eta \kappa }.\end{eqnarray} \tag{ 3 }$$

Then, boundary-layer peel-away occurs when the local Rayleigh number is greater than the critical Rayleigh number for convection (i.e., ${\mathrm{Ra}}_{\mathrm{loc}}\gt {\mathrm{Ra}}_{\mathrm{crit}}$, where ${\mathrm{Ra}}_{\mathrm{crit}}\sim 1100$). Setting ${\mathrm{Ra}}_{\mathrm{loc}}={\mathrm{Ra}}_{\mathrm{crit}}$, solving for the boundary-layer thickness, and substituting in the mantle Rayleigh number

$$\begin{eqnarray}&&\mathrm{Ra}=\displaystyle \frac{\alpha {\rho }_{{\rm{m}}}g{\rm{\Delta }}{{Th}}^{3}}{\eta \kappa },\end{eqnarray} \tag{ 4 }$$

where $h=h(M)$ is the mantle thickness, we find

$$\begin{eqnarray}&&\delta \sim h{\left(\displaystyle \frac{{\mathrm{Ra}}_{\mathrm{crit}}}{\mathrm{Ra}}\right)}^{1/3}.\end{eqnarray} \tag{ 5 }$$

Note that h itself cancels out in Equation (5) when inserting in Equation (4) because we have substituted in the mantle Rayleigh number in order to motivate the scaling relationship derived from numerical simulations (see Equation (7) below).

Zoom In Zoom Out Reset image size

Using Equation (5), we can find the conducted flux through the boundary layer

$$\begin{eqnarray}&&F=\displaystyle \frac{k{\rm{\Delta }}T}{\delta }\sim \displaystyle \frac{k{\rm{\Delta }}T}{h}{\left(\displaystyle \frac{\mathrm{Ra}}{{\mathrm{Ra}}_{\mathrm{crit}}}\right)}^{1/3},\end{eqnarray} \tag{ 6 }$$

where k is the thermal conductivity of the boundary layer. In this work, we take a more general power-law form for our Nusselt number scaling, which relates the outgoing flux from the mantle to the Rayleigh number:

$$\begin{eqnarray}&&\displaystyle \frac{F}{{F}_{\mathrm{cond}}}\equiv \mathrm{Nu}=a{\left(\displaystyle \frac{\mathrm{Ra}}{{\mathrm{Ra}}_{\mathrm{crit}}}\right)}^{\beta }.\end{eqnarray} \tag{ 7 }$$

In Equation (7), $\mathrm{Nu}$ is the Nusselt number, F is the convected flux, and ${F}_{\mathrm{cond}}=k{\rm{\Delta }}T/h$ is the flux that would be conducted if the mantle were not convecting. It is expected from numerical studies of convection (Schubert 1979) that $\beta \approx 0.3$, which we take as our nominal value. Note that if β is smaller than the value assumed here, planetary thermal evolution would be slower. As in Schaefer & Sasselov (2015), we set a = 1, because a is an order-unity parameter and using our theoretically derived outgoing flux from Equation (6) we expect that a should be equal to one. Note that this model only requires the characteristic temperature at the interface between the upper boundary layer and mantle. As a result, we do not consider the actual (nearly adiabatic) temperature profile of the mantle. Additionally, the temperature contrast across the boundary layer is much greater than that between the top and bottom of the mantle. Given that the argument for convection driven by boundary-layer peel off requires local quantities (e.g., κ, α relevant for the boundary layer itself) rather than global quantities, we consider the upper-mantle viscosity η in our model. This results in a pressure-independent viscosity, as will be discussed further in Section 2.1.2.

Given the flux conducted out of the mantle from Equation (7), we can write down a thermal evolution equation that allows us to solve for the mantle temperature as a function of time and mantle water mass fraction. This is

$$\begin{eqnarray}&&{\rho }_{{\rm{m}}}{c}_{p}\displaystyle \frac{{dT}}{{dt}}=Q-\displaystyle \frac{A(M)F(T,x)}{V(M)},\end{eqnarray} \tag{ 8 }$$

where c_p is the mantle heat capacity, $Q={Q}_{0}{e}^{-t/{\tau }_{\mathrm{decay}}}$ is the heating rate from radionuclides with ${\tau }_{\mathrm{decay}}=2\,\mathrm{Gyr}$, A(M) is the planet surface area, and V(M) is the mantle volume. We do not include the Kelvin–Helmholtz contraction term, which is small at late times. We will non-dimensionalize Equation (8) in Section 2.1.3 to elucidate its dependence on temperature and mantle water mass fraction.

2.1.1. Scaling with Planet Mass

To calculate mass-dependent planetary parameters ($h,A,V,g$), we use the scaling laws of Valencia et al. (2006) that take into account internal compression effects on the radius. These scaling relations utilize a constant core mass fraction to relate the planetary radius R and core radius R_c to planetary mass

$$\begin{eqnarray}R & = & {R}_{\oplus }{\left(\displaystyle \frac{M}{{M}_{\oplus }}\right)}^{p},\\ {R}_{{\rm{c}}} & = & {{cR}}_{\oplus }{\left(\displaystyle \frac{M}{{M}_{\oplus }}\right)}^{{p}_{{\rm{c}}}},\end{eqnarray} \tag{ 9 }$$

where p = 0.27, c = 0.547, ${p}_{{\rm{c}}}=0.25$. Using Equation (9), we can then calculate $h=R-{R}_{{\rm{c}}}$, $A=4\pi {R}^{2}$, $V=4\pi /3({R}^{3}-{R}_{{\rm{c}}}^{3})$, $g={GM}/{R}^{2}$.

2.1.2. Viscosity

The mantle viscosity depends both on temperature and the mantle water fraction. We use a similar parameterization as Sandu et al. (2011) and Schaefer & Sasselov (2015) for the mantle viscosity; however, we choose not to incorporate the pressure-dependence of viscosity. As discussed in Section 2.1, we do so because we are interested in convection driven by upper boundary peel off, which occurs in the upper mantle where pressures are relatively small. Additionally, the high-viscosity case of the Schaefer & Sasselov (2015) water cycling model did not reproduce Earth's near steady state or present ocean coverage. This is because the evolution timescales are too long to reach steady state in the high-viscosity case. However, when utilizing low viscosities, the system does converge to an approximate steady-state surface water mass fraction in all cases. We show in Section 3.1 that this choice of viscosity approximates Earth's mantle temperature well when we choose Earth-like parameters.

Our viscosity is hence parameterized as

$$\begin{eqnarray}&&\eta \approx {\eta }_{0}{f}_{{\rm{w}}}^{-r}\exp \left[\displaystyle \frac{{E}_{{\rm{a}}}}{{R}_{\mathrm{gas}}}\left(\displaystyle \frac{1}{T}-\displaystyle \frac{1}{{T}_{\mathrm{ref}}}\right)\right],\end{eqnarray} \tag{ 10 }$$

where ${\eta }_{0}$ gives the viscosity scale, E_a is activation energy, ${R}_{\mathrm{gas}}$ is the universal gas constant, ${T}_{\mathrm{ref}}$ is the reference mantle temperature, and f_w is the water fugacity. We assume throughout this work that r = 1, which is the nominal value used by Schaefer & Sasselov (2015) and that expected from experiments on wet diffusion in olivine (Hirth & Kohlstedt 2003). As in Schaefer & Sasselov (2015), we relate the water abundance to fugacity using experimental data on the concentrations of water in olivine from Li et al. (2008) as

$$\begin{eqnarray}\mathrm{ln}{f}_{{\rm{w}}} & = & {c}_{0}+{c}_{1}\mathrm{ln}\left(\displaystyle \frac{{Bx}{\mu }_{\mathrm{oliv}}/{\mu }_{{\rm{w}}}}{1-x{\mu }_{\mathrm{oliv}}/{\mu }_{{\rm{w}}}}\right)\\ & & +{c}_{2}{\mathrm{ln}}^{2}\left(\displaystyle \frac{{Bx}{\mu }_{\mathrm{oliv}}/{\mu }_{{\rm{w}}}}{1-x{\mu }_{\mathrm{oliv}}/{\mu }_{{\rm{w}}}}\right)+{c}_{3}{\mathrm{ln}}^{3}\left(\displaystyle \frac{{Bx}{\mu }_{\mathrm{oliv}}/{\mu }_{{\rm{w}}}}{1-x{\mu }_{\mathrm{oliv}}/{\mu }_{{\rm{w}}}}\right),\end{eqnarray} \tag{ 11 }$$

where ${c}_{0}=-7.9859$, ${c}_{1}=4.3559$, ${c}_{2}=-0.5742$, ${c}_{3}=0.0227$, and $B=2\times {10}^{6}$ is a conversion to number concentration (H atom/10⁶ Si atoms), ${\mu }_{\mathrm{oliv}}$ is the molecular weight of olivine and ${\mu }_{{\rm{w}}}$ is the molecular weight of water. As in Schaefer & Sasselov (2015), we choose ${\eta }_{0}$ such that $\eta (x={x}_{\oplus },T={T}_{\mathrm{ref}})={10}^{21}\,\mathrm{Pa}\ {\rm{s}}$ which yields mantle temperatures that approximately reproduce those on Earth.

2.1.3. Non-dimensional Thermal Evolution Equation

Throughout the remainder of this paper, we will work with non-dimensional versions of the thermal evolution and volatile cycling equations. We do so because it elucidates the essential physical processes and non-dimensional control variables. Substituting our scaling for mantle heat flux from Equation (7) into (8) and using our prescription for viscosity from Equation (10), we can non-dimensionalize the thermal evolution equation as

$$\begin{eqnarray}\displaystyle \frac{d\tilde{T}}{d{\tau }_{\mathrm{heat}}} & = & \ \tilde{Q}({\tau }_{\mathrm{heat}})\\ & & -\tilde{{F}_{0}}{\tilde{f}}_{{\rm{w}}}^{\beta }(\tilde{x}){(\tilde{T}-\tilde{{T}_{{\rm{s}}}})}^{\beta +1}\exp \left[-\displaystyle \frac{\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right],\end{eqnarray} \tag{ 12 }$$

where the non-dimensional temperature is $\tilde{T}=T/{T}_{\mathrm{ref}}$, the non-dimensional mantle water mass fraction is $\tilde{x}={{xf}}_{{\rm{m}}}/({\omega }_{0}{\tilde{f}}_{{\rm{b}}}),$ and $\tilde{{F}_{0}}={F}_{0}/{Q}_{0}$, where Q₀ is a constant and

$$\begin{eqnarray}&&{F}_{0}=\displaystyle \frac{{{kT}}_{\mathrm{ref}}^{1+\beta }A}{{hV}}{\left(\displaystyle \frac{\alpha {\rho }_{{\rm{m}}}{{gh}}^{3}{f}_{{\rm{w}}}(\tilde{x}=1)}{{\mathrm{Ra}}_{\mathrm{crit}}\kappa {\eta }_{0}}\right)}^{\beta },\end{eqnarray} \tag{ 13 }$$

where the non-dimensional fugacity is ${\tilde{f}}_{{\rm{w}}}={f}_{{\rm{w}}}/{f}_{{\rm{w}}}(\tilde{x}=1)$, the reference temperature is ${\tilde{T}}_{{\rm{m}}}={T}_{\mathrm{ref}}{R}_{\mathrm{gas}}/{E}_{{\rm{a}}}$, and the surface temperature is ${\tilde{T}}_{{\rm{s}}}={T}_{{\rm{s}}}/{T}_{\mathrm{ref}}$. Lastly, the non-dimensional heating timescale is ${\tau }_{\mathrm{heat}}={{tQ}}_{0}/({\rho }_{{\rm{m}}}{c}_{p}{T}_{\mathrm{ref}})$. The typical values of these non-dimensional parameters are shown in Table 1.

Table 1.  Non-dimensional Variables and Parameters Used in This Paper, Their Symbols, and Their Value for Earth-like Parameters

Quantity	Symbol	Fiducial value
Mantle water mass fraction${}^{\star }$	$\tilde{x}$	${\tilde{x}}_{\oplus }=1.32$
Mantle temperature${}^{\star }$	$\tilde{T}$	${\tilde{T}}_{\mathrm{ref}}=1$
Planet mass	$\tilde{M}$	1
Total water mass fraction	$\tilde{\omega }$	2.07
Heating timescale	${\tau }_{\mathrm{heat}}$	$t({\tau }_{\mathrm{heat}}=1)=4.02\,\mathrm{Gyr}$
Pressure-dependent volatile cycling timescale	τ	$t(\tau =1)=2.87\,\mathrm{Gyr}$
Temperature-dependent volatile cycling timescale	${\tau }_{\mathrm{SS}}$	$t({\tau }_{\mathrm{SS}}=1)=2.22\,\mathrm{Gyr}$
Hybrid volatile cycling timescale	${\tau }_{\mathrm{hyb}}$	$t({\tau }_{\mathrm{hyb}}=1)=2.22\,\mathrm{Gyr}$
Heat flux	${\tilde{F}}_{0}$	0.531
Heat flux scaling coefficient	β	0.3
Critical Rayleigh number	${\mathrm{Ra}}_{\mathrm{crit}}$	1100
Water fugacity	${\tilde{f}}_{{\rm{w}}}$	1
Surface temperature	${\tilde{T}}_{{\rm{s}}}$	0.175
Reference temperature	${\tilde{T}}_{{\rm{m}}}$	0.040
Mantle water mass fraction of Earth	${\tilde{X}}_{\oplus }$	1.32
Seafloor pressure degassing exponent	μ	1
Seafloor pressure regassing exponent	σ	1
Solidus temperature	${\tilde{T}}_{\mathrm{sol},\mathrm{dry}}$	0.780
Liquidus temperature	${\tilde{T}}_{\mathrm{liq},\mathrm{dry}}$	0.936
Temperature-dependent degassing coefficient	$\tilde{{\rm{\Pi }}}$	0.102
Solidus depression constant	$\tilde{\lambda }$	$8.16\times {10}^{-5}$
Solidus depression coefficient	γ	0.75
Melt fraction exponent	θ	1.5
Pressure-dependent degassing coefficient	$\tilde{E}$	0.473
Maximum mantle water mass fraction	${\tilde{x}}_{\max }$	15.9

Note. Stars denote model state variables.

Download table as:  ASCIITypeset image

We seek to explore a variety of different volatile cycling parameterizations, each of which relies on the following expression for the time rate of change of mantle water mass fraction (Cowan & Abbot 2014)

$$\begin{eqnarray}&&\displaystyle \frac{{dx}}{{dt}}=\displaystyle \frac{{L}_{\mathrm{MOR}}S(T)}{{f}_{{\rm{m}}}M}({w}_{\downarrow }-{w}_{\uparrow }),\end{eqnarray} \tag{ 14 }$$

where S(T) is the temperature-dependent spreading rate (discussed further in Section 2.2.2), ${f}_{{\rm{m}}}M$ is the mantle mass (where f_m is the mantle mass fraction), w_↓ is the regassing rate, and w_↑ is the degassing rate. Each of the volatile cycling parameterizations we consider utilizes different regassing and degassing rates, which we explore in the following Sections 2.2.1–2.2.3.

2.2.1. Seafloor Pressure-dependent Degassing and Regassing

In this section, we construct a non-dimensional version of Equation (14) corresponding to the volatile cycling model of Cowan & Abbot (2014). This model determines the water mass fraction of the mantle independent of the mantle temperature. We utilize their expressions for the regassing and degassing rates:

$$\begin{eqnarray}&&{w}_{\downarrow }={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}{d}_{{\rm{h}}}(P)\chi ,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{w}_{\uparrow }=x{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas}}(P),\end{eqnarray} \tag{ 16 }$$

where x_h is the mass fraction of water in the hydrated crust, ${\rho }_{{\rm{c}}}$ is the density of the oceanic crust, χ is the subduction efficiency, ${\rho }_{{\rm{m}}}$ is the density of the upper mantle, and ${d}_{\mathrm{melt}}$ is the depth of melting below mid-ocean ridges. We take ${d}_{\mathrm{melt}}$ and χ as constants, with their fiducial value equal to their fiducial value in Cowan & Abbot (2014). As in Cowan & Abbot (2014), we take the depth of the serpentinized layer ${d}_{{\rm{h}}}(P)$ and the fraction of the water in the melt that is degassed ${f}_{\mathrm{degas}}(P)$ to be power laws with seafloor pressure, with d_h increasing with increasing pressure and ${f}_{\mathrm{degas}}$ decreasing with increasing pressure. See Appendix A.1 for a thorough explanation of these parameters and the derivation that follows to give the mantle water mass fraction rate of change with time. Inserting Equations (15) and (16) into (14) and non-dimensionalizing, we find

$$\begin{eqnarray}\displaystyle \frac{d\tilde{x}}{d\tau } & = & {[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{\sigma }-{\tilde{X}}_{\oplus }^{-1}\tilde{x}{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{-\mu }\\ & = & {\tilde{F}}_{\downarrow ,\mathrm{CA}}-{\tilde{F}}_{\uparrow ,\mathrm{CA}}.\end{eqnarray} \tag{ 17 }$$

In Equation (17),

$$\begin{eqnarray}&&{\tilde{X}}_{\oplus }=\displaystyle \frac{{x}_{h}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }{f}_{{\rm{M}}}}{{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}\end{eqnarray} \tag{ 18 }$$

is a degassing coefficient identified by Cowan & Abbot (2014) as the mantle water mass fraction of Earth, the non-dimensional mantle water mass fraction is (as before) $\tilde{x}={{xf}}_{{\rm{m}}}/({\omega }_{0}{\tilde{f}}_{{\rm{b}}})$, $\tilde{\omega }=\omega /({\omega }_{0}{\tilde{f}}_{{\rm{b}}})$ is the normalized total water mass fraction, $\tilde{g}=g/{g}_{\oplus }$, and

$$\begin{eqnarray}&&{\tau }_{\mathrm{CA}}=\tau =t\displaystyle \frac{{L}_{\mathrm{MOR}}{{Sx}}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }}{M{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}\end{eqnarray} \tag{ 19 }$$

is the non-dimensional time, which is inversely related to the seafloor overturning timescale $A/({L}_{\mathrm{MOR}}S)$. The first term on the right-hand side of Equation (17) is the regassing flux ${\tilde{F}}_{\downarrow ,\mathrm{CA}}$ and the second term is the degassing flux ${\tilde{F}}_{\uparrow ,\mathrm{CA}}$. In this model, the spreading rate S does not depend on mantle temperature, but it will in Sections 2.2.2 and 2.2.3. We write the non-dimensional timescale here as τ because it will be the timescale that all of our solutions are converted to for inter-comparison.

2.2.2. Temperature-dependent Degassing and Regassing

In this section, we write down a simplified, non-dimensional form of Section 2.3 in Schaefer & Sasselov (2015). Their regassing and degassing rates are

$$\begin{eqnarray}&&{w}_{\downarrow }={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}}}(T),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{w}_{\uparrow }={\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{f}_{\mathrm{melt}}(T)x.\end{eqnarray} \tag{ 21 }$$

Equation (20) is identical to Equation (15) except now the hydrated layer depth is a function of temperature (see Appendix A.2 for details), and Equation (21) is similar to Equation (16), except ${f}_{\mathrm{degas}}(P)$ has been replaced by ${f}_{\mathrm{degas},\oplus }{f}_{\mathrm{melt}}(T)$ with ${f}_{\mathrm{melt}}(T)$ the temperature-dependent mass fraction. In Equation (21), we have assumed that the mass fraction of water in melt is the same as the mass fraction of water in the mantle due to the extremely low ($\approx 1 \%$) difference in water partitioning between melt and mantle rock.

Inserting our expressions (20) and (21) for regassing and degassing rates into Equation (14) and non-dimensionalizing, we find (see Appendix A.2 for the steps and parameterizations of $S(T),{d}_{{\rm{h}}}(T),{f}_{\mathrm{melt}}(T)$)

$$\begin{eqnarray}\displaystyle \frac{d\tilde{x}}{d{\tau }_{\mathrm{SS}}} & = & {\tilde{f}}_{{\rm{w}}}^{\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{\beta -1}\exp \left[\displaystyle \frac{-\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & -\tilde{{\rm{\Pi }}}{\tilde{f}}_{{\rm{w}}}^{2\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{2\beta }\exp \left[\displaystyle \frac{-2\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & \times \tilde{x}{(\tilde{T}-{\tilde{T}}_{\mathrm{sol},\mathrm{dry}}+\tilde{\lambda }{\tilde{x}}^{\gamma })}^{\theta }\\ & = & {\tilde{F}}_{\downarrow ,\mathrm{SS}}-{\tilde{F}}_{\uparrow ,\mathrm{SS}}.\end{eqnarray} \tag{ 22 }$$

In Equation (22), the non-dimensional solidus depression coefficient is $\tilde{\lambda }=\tilde{K}{({\omega }_{0}{\tilde{f}}_{{\rm{b}}}/{f}_{{\rm{m}}})}^{\gamma }$, the degassing coefficient is $\tilde{{\rm{\Pi }}}={\rm{\Pi }}/D$, where

$$\begin{eqnarray}&&{\rm{\Pi }}={\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }\displaystyle \frac{{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}{{f}_{{\rm{m}}}}{({\tilde{T}}_{\mathrm{liq},\mathrm{dry}}-{\tilde{T}}_{\mathrm{sol},\mathrm{dry}})}^{-\theta },\end{eqnarray} \tag{ 23 }$$

and the regassing coefficient (related to the hydrated layer depth) is

$$\begin{eqnarray}&&D={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}}{h}^{(1-3\beta )}\displaystyle \frac{({T}_{\mathrm{serp}}-{T}_{{\rm{s}}})}{{T}_{\mathrm{ref}}^{1+\beta }}{\left(\displaystyle \frac{\kappa {\mathrm{Ra}}_{\mathrm{crit}}{\eta }_{0}}{\alpha {\rho }_{{\rm{m}}}{{gf}}_{{\rm{w}}}(\tilde{x}=1)}\right)}^{\beta },\end{eqnarray} \tag{ 24 }$$

and ${\tau }_{\mathrm{SS}}={tD}/{\rm{\Sigma }}$, where

$$\begin{eqnarray}&&{\rm{\Sigma }}=M{\omega }_{0}{\tilde{f}}_{{\rm{b}}}\displaystyle \frac{{h}^{(1-6\beta )}{\kappa }^{(2\beta -1)}}{10.76{L}_{\mathrm{MOR}}}{\left(\displaystyle \frac{{\mathrm{Ra}}_{\mathrm{crit}}{\eta }_{0}}{\alpha {\rho }_{{\rm{m}}}{{gT}}_{\mathrm{ref}}{f}_{{\rm{w}}}(\tilde{x}=1)}\right)}^{2\beta }\end{eqnarray} \tag{ 25 }$$

is related to the spreading rate. Additonally, we have re-expressed the first term on the right-hand side of Equation (22) as the regassing flux ${\tilde{F}}_{\downarrow ,\mathrm{SS}}$ and the second term on the right-hand side as the degassing flux ${\tilde{F}}_{\uparrow ,\mathrm{SS}}$.

In our coupled integrations of Equations (12) and (22), we ensure that the hydrated layer does not contain more water than the surface in order to maintain water mass balance (Schaefer & Sasselov 2015). In terms of our analytic model, this is equivalent to ensuring that the regassing coefficient D (which is related inversely to the non-dimensional degassing coefficient $\tilde{{\rm{\Pi }}}$) never exceeds a critical value, which is written in Equation (58). As a result, this is a constraint on the rate of subduction of water that ensures that the amount of water in the mantle never exceeds the total amount of water in the planet.

2.2.3. Seafloor Pressure-dependent Degassing and Temperature-dependent Regassing

In this section, we construct a model where the degassing rate is determined by seafloor pressure (because volcanism rates will be lower if overburden pressure is higher) and the regassing rate is determined by the mantle temperature (as the depth of serpentinization will be lower if temperature is higher). We construct such a model because serpentinization can only happen below a critical temperature, whereas it has not been conclusively shown to depend on overburden pressure. Meanwhile, it has been shown that volcanism rates on exoplanets should be inversely related to the overburden pressure (Kite et al. 2009). In this model, the degassing rate is taken from Equation (16) with $\mu =1$ (the value expected from Kite et al. 2009) and the regassing rate from Equation (20). Using the same method as in Sections 2.2.1 and 2.2.2, we substitute these into Equation (14) and non-dimensionalize (see Appendix A.3 for more details). Doing so, we find

$$\begin{eqnarray}\displaystyle \frac{d\tilde{x}}{d{\tau }_{\mathrm{hyb}}} & = & {\tilde{f}}_{{\rm{w}}}^{\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{(\beta -1)}\exp \left[\displaystyle \frac{-\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & -\tilde{E}{\tilde{f}}_{{\rm{w}}}^{2\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{(2\beta )}\exp \left[\displaystyle \frac{-2\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & \times \tilde{x}{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{-1}\\ & = & {\tilde{F}}_{\downarrow ,\mathrm{hyb}}-{\tilde{F}}_{\uparrow ,\mathrm{hyb}},\end{eqnarray} \tag{ 26 }$$

where $\tilde{E}=E/D$, $E={\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{\omega }_{0}{\tilde{f}}_{{\rm{b}}}/{f}_{{\rm{m}}}$, and ${\tau }_{\mathrm{hyb}}={\tau }_{\mathrm{SS}}={tD}/{\rm{\Sigma }}$. As before, we have re-written the first term on the right-hand side as the regassing flux ${\tilde{F}}_{\downarrow ,\mathrm{hyb}}$ and the second term on the right-hand side as the degassing flux ${\tilde{F}}_{\uparrow ,\mathrm{hyb}}$.

Before turning to the steady-state solutions, we compare directly the time-dependent evolution of the three models. Figure 2 shows such a comparison for Earth-like parameters. The hybrid and pressure-dependent models reach a steady state after a time $\tau \approx 1$, with a value that is independent of initial conditions (not shown). Meanwhile, the solely temperature-dependent model reaches a near steady state where a tiny amount of net regassing still occurs. Figure 2 shows that although the mantle temperature evolution does not vary by more than $\sim 15 \%$ among models, the evolution and steady-state value of mantle water mass fraction varies greatly. Notably, the degassing parameterizations lead to different values of the end-state mantle water mass fraction $\tilde{x}$ even though both the hybrid and solely temperature-dependent models have their late volatile evolution determined by water mass balance between the hydrated layer and surface. We will explain this in detail in the steady-state solutions of Section 3.3.

Zoom In Zoom Out Reset image size

Figure 2 also shows the evolution of the non-dimensional ocean depth for each of the models considered. To determine the ocean depth for a given $\tilde{x}$, we utilize Equation (15) of Cowan & Abbot (2014). This relates ocean depth to seafloor pressure by ${d}_{{\rm{w}}}=P/(g{\rho }_{{\rm{w}}}$), where d_w is the ocean depth and ${\rho }_{{\rm{w}}}$ the density of water. Non-dimensionalizing, we find

$$\begin{eqnarray}&&{\tilde{d}}_{{\rm{w}}}=\tilde{g}(\tilde{\omega }-\tilde{x}),\end{eqnarray} \tag{ 27 }$$

where ${\tilde{d}}_{{\rm{w}}}={d}_{{\rm{w}}}/{d}_{{\rm{w}},\oplus }$ is the ocean depth normalized to that of Earth, ${d}_{{\rm{w}},\oplus }=4$ km. Note that gravity comes into Equation (27) from the ratio of planetary mass to area, which is proportional to gravity. Because we use an Earth-like $\tilde{g}=1$ in these calculations, the ocean depth is equivalent to the surface water mass fraction $\tilde{\omega }-\tilde{x}$. The evolution of ocean depth with time hence has an opposite sign to that for the mantle water mass fraction, with the hybrid model having the deepest oceans and the Schaefer & Sasselov (2015) model the least end-state surface water.

Alongside the ocean depth (equivalently surface water mass fraction) we show the evolution of the mass fraction of water in the hydrated layer ${\tilde{\omega }}_{\mathrm{hyd}}={M}_{\mathrm{hyd}}/({\omega }_{0}{\tilde{f}}_{{\rm{b}}}M)$, where ${M}_{\mathrm{hyd}}$ is the mass of water in the hydrated layer. We can relate the mass fraction of water in the hydrated layer and the hydrated layer depth as

$$\begin{eqnarray}&&{\tilde{\omega }}_{\mathrm{hyd}}=\displaystyle \frac{4\pi }{3{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}[{R}^{3}-{(R-{d}_{{\rm{h}}})}^{3}]\displaystyle \frac{{x}_{{\rm{h}}}{\rho }_{{\rm{m}}}}{M},\end{eqnarray} \tag{ 28 }$$

where as in Cowan & Abbot (2014) we use ${x}_{{\rm{h}}}=0.05$ and ${\rho }_{{\rm{m}}}=3.3\times {10}^{3}\ \mathrm{kg}\ {{\rm{m}}}^{-3}$. For the pressure-dependent model, we calculate the hydrated layer thickness from Equation (35) using ${d}_{{\rm{h}},\oplus }=3\ \mathrm{km}$. For the temperature-dependent and hybrid models, we calculate the hydrated layer thickness using Equation (50) unless the hydrated layer thickness limit is violated, in which case we calculate it from Equations (57) and (58). As shown in Figure 2, the amount of water in the hydrated layer increases drastically in the temperature-dependent and hybrid models when regassing begins to dominate over degassing. At late times in the temperature-dependent and hybrid models, the hydrated layer water mass fraction is equal to the surface mass fraction. As discussed further in Section 3.2, this is not representative of present-day Earth. However, the hydrated layer water mass fraction in the pressure-dependent model stays small at all times, as in this model the hydrated layer thickness is a simple power-law with pressure.

Figure 2 also shows the individual regassing and degassing rates (the first and second terms on the right-hand side of Equations (17), (22), and (26)). The models in Section 2.2.2 and Section 2.2.3 have an initial phase of degassing from the mantle followed by strong regassing of water back to the mantle. Note that this initial phase of degassing occurs over a shorter timescale than the low-viscosity ("boundary-layer") model of Schaefer & Sasselov (2015), as we assume that all of the water is in melt (thereby increasing the amount able to be outgassed) and use a higher initial mantle temperature. We can qualitatively understand the varied evolution in different models by examining how the degassing and regassing rates vary with temperature and/or pressure. Note that the initial phase of degassing from the mantle occurs whether the degassing is temperature or seafloor pressure-dependent, as long as regassing of water back into the mantle is temperature-dependent. This is because initially the regassing rate is very low due to the smaller hydrated layer thickness when the mantle is hot. Regassing then becomes more efficient as the mantle temperature drops and the hydrated layer thickness grows. Similarly, degassing becomes less efficient at later times. This is because the mantle is cooler and hence has a lower melt fraction and the seafloor pressure is greater, both decreasing the rate of volcanism.

However, the end-state evolution for these models is slightly different. For the Schaefer & Sasselov (2015) model, the degassing rate drops to zero at late times because of the lack of melt available for degassing, while the regassing rate (though small) is non-negligible. As a result, there is net regassing at late times in our temperature-dependent model. The hybrid model, meanwhile, does reach a true steady state. This is because the degassing rate stays large at late times, as it is not dependent on the mantle convection itself, and because the hydrated layer thickness limit is reached, which sharply decreases the regassing rate. Meanwhile, the regassing rate is limited by the hydrated layer thickness, which reaches the limit given by Equation (58) after ∼1 Gyr of evolution. This regassing rate then decreases more strongly with time due to the constancy of the maximum hydrated layer thickness and the decreasing spreading rate with decreasing mantle temperature, leading to convergence of the degassing and regassing rates and entrance into steady state. Such a steady state was not found in the volatile cycling models of Sandu et al. (2011) and Schaefer & Sasselov (2015), which utilized parameterized convection. This is because the evolution was either too slow to reach steady state over the age of the observable universe or because slow ingassing continued to deplete the surface water reservoir. The former occurs in models that consider the viscosity as an average mantle viscosity rather than that relevant for the interface between the boundary layer and mantle interior. We do the latter in this work. As discussed in Section 1.2, Earth is likely currently at or near a steady state in surface water mass fraction. We can estimate the steady-state ocean depths for our models as a function of planetary parameters, and will do so in Section 3.3.

Given that each of the models considered produces different mantle and volatile evolution, here we compare their results to that of Earth. We do so in order to understand which model(s) may be most physical for application to exoplanets. We compare the key variables (mantle water mass fraction, surface water mass fraction, upper mantle temperature, ocean depth, hydrated layer mass fraction, degassing and regassing rates) derived from our Earth-like model after 4.5 Gyr of evolution to those of Earth itself in Table 2. Note that we are using an ocean basin covering fraction that is 1.3 times that of Earth, so we calculate equivalent ocean depths for this increased ocean basin area.

Table 2.  Values of Key Model Variables (Mantle Water Mass Fraction, Surface Water Mass Fraction, Upper Mantle Temperature, Ocean Depth, Mass Fraction of Water in the Hydrated Layer, Degassing and Regassing Fluxes) after $4.5\ \mathrm{Gyr}$ of Evolution for the Earth-like Case Shown in Figure 2

	Pressure-dependent	Temperature-dependent	Hybrid	Earth
Mantle water mass fraction $\tilde{x}$	1.14 (0.050%)	1.65 (0.072%)	1.06 (0.047%)	1.3 (0.057%)
Surface water mass fraction $\tilde{\omega }-\tilde{x}$	0.931 (0.028%)	0.424 (0.013%)	1.01 (0.030%)	0.75 (0.022%)
Upper mantle temperature $\tilde{T}$	0.866 (1386 K)	0.844 (1350 K)	0.875 (1400 K)	$0.8\mbox{--}1.2$ (1280–1920 K)
Ocean depth ${\tilde{d}}_{{\rm{w}}}$	0.931 (3.72 km)	0.424 (1.70 km)	1.01 (4.04 km)	0.75 (3.0 km)
Hydrated layer mass fraction ${\tilde{\omega }}_{\mathrm{hyd}}$	0.132 ($3.9\times {10}^{-3}$%)	0.423 (0.013%)	1.01 (0.030%)	0.14 ($4.23\times {10}^{-3}$%)
Degassing flux ${\tilde{F}}_{\uparrow }$	0.928	$2.55\times {10}^{-3}$	0.105	0.03
	($6.06\times {10}^{12}$ kg yr−1)	($1.67\times {10}^{10}$ kg yr−1)	($6.86\times {10}^{11}$ kg yr−1)	($2\times {10}^{11}$ kg yr−1)
Regassing flux ${\tilde{F}}_{\downarrow }$	0.932	0.0973	0.119	0.1–0.4
	($6.09\times {10}^{12}$ kg yr−1)	($6.36\times {10}^{11}$ kg yr−1)	($7.78\times {10}^{11}$ kg yr−1)	($[0.7\mbox{--}2.9]\times {10}^{12}$ kg yr−1)

Note. Approximate present-day Earth values from Hirschmann (2006) and Cowan & Abbot (2014) are shown for reference. Dimensionful values are shown in parentheses. The value of Earth's ocean depth is computed given an ocean basin covering fraction (${\tilde{f}}_{{\rm{b}}}$) of 1.3, which is larger than for current Earth and leads to a smaller ocean depth than seen. Note that Earth's mantle water mass fraction and regassing and degassing rates are very approximate, accurate at best to within a factor of ∼2 (Cowan & Abbot 2014). The uncertainty in characteristic Earth upper mantle temperatures stems from the range of possible relevant depths.

Download table as:  ASCIITypeset image

No model matches well Earth's present-day water partitioning between mantle and surface, with the estimated Earth mantle water mass fraction lying between those of the temperature-dependent and pressure-dependent models. As mentioned previously, the temperature-dependent and hybrid models well over-predict the amount of water in the hydrated layer. However, the pressure-dependent model matches well the estimated mass fraction of water in the hydrated layer. Though there is considerable uncertainty in Earth's mantle water mass fraction and regassing and degassing rates, the hybrid model has reasonably similar values to both of these. The pressure-dependent model over-predicts both the degassing and regassing fluxes, as these do not decrease strongly with time in this model. However, note that the degassing and regassing fluxes estimated from observations for Earth do not match, so if Earth water cycling is in steady state one of these must be erroneous by approximately an order of magnitude. The temperature-dependent model reasonably matches Earth's present-day regassing rate, but well over-predicts the mantle water mass fraction and under-predicts the degassing rate, due to the lack of degassing at late times.

In general, none of the models alone match all of the Earth constraints, though each model does approximate at least one constraint. Given that the temperature-dependent and hybrid models have almost all of their surface water in the hydrated layer, the pressure-dependent model is most representative of present-day Earth. The pressure-dependent model is also closest to the present-day surface water mass fraction of Earth. Though it matches Earth's regassing flux within a factor of two, it over-predicts the degassing flux. However, if the degassing and regassing fluxes of present-day Earth are in steady state, including a necessary increase in degassing flux such that it matches the regassing flux would allow the pressure-dependent model to match all available constraints. If not, future work is needed to develop a model that matches well all of the available constraints from Earth.

Given that all of our models reach a steady state in water partitioning on the timescale of a few billion years, we examine steady-state solutions to the models in Sections 2.2.1–2.2.3. We do so because these steady states are the most observationally relevant, as most planets in the habitable zone will lie around ∼Gyr-age main-sequence stars. We note that due to continuous regassing, the steady state for the Schaefer & Sasselov (2015) model is one where the amount of surface water is simply determined by the amount of water that can be incorporated into the mantle of a planet. This is the "petrological limit" of the mantle, and will be discussed in detail in Section 4 below. Note that if the total water mass fraction is less than the petrological limit, the mantle holds all of the water except that which remains on the surface due to mass balance with the hydrated layer. Solving for the steady state of the pressure-dependent and hybrid models using Equations (17) and (26) gives

$$\begin{eqnarray}&&\tilde{x}={\tilde{X}}_{\oplus }{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{\mu +\sigma }\end{eqnarray} \tag{ 29 }$$

for the Cowan & Abbot (2014) model in Section 2.2.1 and

$$\begin{eqnarray}&&\tilde{x}=\tilde{\omega }{\left(1+\tilde{E}{\tilde{f}}_{{\rm{w}}}^{\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{1+\beta }\exp \left[\displaystyle \frac{-\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]{\tilde{g}}^{-2}\right)}^{-1}\end{eqnarray} \tag{ 30 }$$

for the hybrid model in Section 2.2.3. We choose $\tilde{T}={\tilde{T}}_{\mathrm{sol},\mathrm{dry}}$ to calculate the steady-state mantle water mass fraction for the hybrid model. This is because the steady state is nearly independent of temperature for $\tilde{T}\lesssim {\tilde{T}}_{{\rm{m}}}$ (see Figure 2). Equation (29) reproduces Equation (20) of Cowan & Abbot (2014). Note that the steady-state value of $\tilde{x}$ for the model in Section 2.2.2 is independent of $\tilde{\omega }$. This greatly limits the relative amount of water that can be put into the mantle of planets with large total water fractions.

Equations (29) and (30) give us transcendental expressions for the steady-state mantle water mass fraction as a function of mantle temperature and planet mass for each model. In Section 4, we solve Equations (29) and (30) and relate the mantle water mass fraction to the surface ocean depth to determine the waterworld limit for various assumptions about the processes that control volatile cycling on exoplanets.

In this section, we perform a sensitivity analysis to determine how the non-dimensional parameters affect our steady-state solutions from Section 3.3. The key unknown parameters that affect our solutions are the maximum mantle water mass fraction ${\tilde{x}}_{\max }$ (which affects all models), Earth mantle water mass fraction ${\tilde{X}}_{\oplus }$ and seafloor pressure power-law exponents $\mu +\sigma$ for the model of Cowan & Abbot (2014), degassing coefficient $\tilde{{\rm{\Pi }}}$ for the Schaefer & Sasselov (2015) model, degassing coefficient $\tilde{E}$ in the hybrid model, and surface temperature ${\tilde{T}}_{{\rm{s}}}$ for both the Schaefer & Sasselov (2015) and hybrid models. Importantly, our steady states are independent of the abundance of radiogenic elements, eliminating some of the natural variation between planetary systems. Though the abundance of radiogenic elements affects the time it takes to reach steady state, the steady-state volatile cycling is independent of the decreasing mantle temperature at late times. Figure 4 shows how varying these parameters in each model affects our derived waterworld boundary.

Zoom In Zoom Out Reset image size

Though the changes in the waterworld boundary with changing $\mu +\sigma$ have been explored in Cowan & Abbot (2014), we reproduce them here for comparison with the other models. Decreasing the dependencies of degassing and regassing on seafloor pressure reduces the water mass fraction at which the surface is completely water-covered, with a maximum decrease of a factor of two between the $\mu +\sigma =2$ and $\mu +\sigma =1$ cases. Similarly, increasing the dependence to $\mu +\sigma =3$ increases the limiting water mass fraction to become a waterworld, but the model reaches the maximum mantle water mass fraction. If ${\tilde{X}}_{\oplus }$ is a factor of 10 lower than used here, the waterworld boundary decreases by a comparably large fraction, especially for super-Earths. If ${\tilde{X}}_{\oplus }$ is much larger than assumed here, the mantle will be at its petrological limit of water intake and the waterworld boundary will be determined by the maximum mantle water mass fraction.

For the model of Schaefer & Sasselov (2015), which is at the petrological limit of maximum mantle water mass fraction, varying $\tilde{{\rm{\Pi }}}$ by an order of magnitude in either direction does not change the waterworld boundary. However, increasing the maximum mantle water mass fraction by a given value increases the total water mass fraction needed to become a waterworld by a comparable amount. For both the hybrid model and the Schaefer & Sasselov (2015) model, changing the surface temperature only leads to minute changes in the waterworld boundary. This is because the surface temperature cannot vary by more than a factor of a few or else liquid water would not be stable on the surface. Increasing $\tilde{E}$ in the hybrid model decreases the total mass fraction needed to become a waterworld, but by less than a factor of two for all masses.

From Figure 4, we can identify four key non-dimensional parameters that change the waterworld boundary by a sizeable amount: ${\tilde{x}}_{\max }$ (which mainly affects the Schaefer & Sasselov 2015 model), $\tilde{E}$ for the hybrid model, and $\sigma +\mu$ and ${\tilde{X}}_{\oplus }$ for the Cowan & Abbot (2014) model. Figure 5 shows how continuously varying these parameters by one order of magnitude around their fiducial value with a planet mass fixed equal to that of Earth affects the water mass fraction at which planets become waterworlds. We also consider varying β, which could affect the solution since a reduced outgoing flux would lead to larger mantle temperatures and hence larger steady-state mantle water mass fractions. As mentioned above, our results are very sensitive to the petrological limit of the mantle water mass fraction, but an order of magnitude increase in $\tilde{E}$, $\sigma +\mu$, and ${\tilde{X}}_{\oplus }$ leads to only a factor of approximately two or less increase in the waterworld limit. The waterworld boundary is also largely insensitive to β, which should not vary by more than a factor of two from its nominal value of 0.3. Note that increasing both $\sigma +\mu$ and ${\tilde{X}}_{\oplus }$ cannot lead to continuous increases in the waterworld limit, as the petrological limit of water mantle mass fraction is reached just above our fiducial values for these parameters. As a result, the Cowan & Abbot (2014) model is, like the Schaefer & Sasselov (2015) model, sensitive to the maximum mantle water mass fraction ${\tilde{x}}_{\max }\,.$

Zoom In Zoom Out Reset image size

Our results are much less sensitive to $\tilde{E}$, for which an order-of-magnitude increase only decreases the waterworld boundary by $\sim 10 \%$. Note that the non-dimensional degassing rate $\tilde{E}$ scales with the normalization of mantle viscosity as $\tilde{E}\propto {\eta }_{0}^{-\beta }$. As a result, if the viscosities were increased (for instance, in the case that the middle-mantle viscosity is more relevant for the convection parameterization), the degassing rate would decrease as a power law with increasing viscosity. However, because our solutions are only weakly dependent on $\tilde{E}$, the choice of mantle viscosity does not greatly affect the waterworld boundary itself. We explore the effects of a larger viscosity further in Section 5.2, as it will substantially affect the evolutionary timescales for water cycling. In general, the conclusion that super-Earths are more likely to be waterworlds if degassing is temperature-independent is robust to order-of-magnitude uncertainties in our non-dimensional parameters.

Though the waterworld boundary itself is largely robust to varying our non-dimensional parameters, the timescale to reach steady state depends on the planet mass and mantle viscosity. In this section, we derive how the evolution timescale varies with these parameters in order to determine the mass regime at which planets may not reach steady state. The evolution timescale scales with planet mass as $t(\tau =1)\propto {M}^{1-p}$, if S is independent of mass (as it is in the pressure-dependent model). Similarly, $t(\tau =1)\propto {M}^{1-2\beta p-2\beta }$ if S is dependent on mass (as it is in the temperature-dependent and hybrid models, see Equation (48)) where p = 0.27 and $\beta \approx 0.3$. As a result, the evolution of the mantle water fraction is slower for larger planets and is only weakly sensitive to mass, with $t(\tau =1)\propto {M}^{0.24}$, in the more realistic case where S depends on planet mass.

The evolution timescale scales with the viscosity as $t({\tau }_{\mathrm{SS}}=1)\propto {\eta }_{0}^{\beta }$. As a result, the evolution timescale is also a power-law in viscosity, increasing with increasing viscosity. The choice of a characteristic mantle viscosity hence may affect the resulting mantle evolution, with an order of magnitude increase in viscosity leading to a factor of approximately two increase in the timescale to reach steady state.

Based on the scaling of the evolutionary timescales with mass alone, our conclusion that water cycling reaches steady state is unaffected. This is because a $5{M}_{\oplus }$ planet would only take $\approx 1.5$ times longer to reach steady state than Earth. However, if the viscosity normalization is more than a factor of $\approx 5$ larger than assumed here, the evolution would take longer than the age of the solar system to reach steady state. This could occur if the viscosity of the deep mantle is relevant for our parameterized convection scheme, and is similar to the conclusion from the high-viscosity models of Schaefer & Sasselov (2015). However, these models were not shown to produce water cycling evolution similar to Earth, while we showed in Section 3.2 that a boundary-layer viscosity can match some of the constraints from Earth. In general, it is clear that understanding in detail which characteristic viscosity is relevant for parameterized convection is necessary to make more robust predictions of water cycling on exoplanets.

In this work, we developed simplified models for water cycling between the mantle and the surface based on previous models in the literature. We did so in order to compare their predictions for whether or not terrestrial exoplanets will be waterworlds. This is the first such test of physical assumptions that has been performed for volatile cycling on planets with varying mass, though Sandu et al. (2011) explored how models of varying complexity affect volatile cycling on Earth. Our models find that the surface water mass fraction reaches a steady state after ∼2 Gyr of evolution. Though this has been found when keeping the ratio of the degassing and regassing rates fixed in time (McGovern & Schubert 1989), no steady state has previously been found when these rates are dependent on the mantle temperature and allowed to separately evolve.

We find a steady state in our models for two reasons: the degassing rate is initially larger than the regassing rate (leading to convergent evolution of the two rates), and the volatile evolution is relatively quick. In our temperature-dependent and hybrid models, the volatile evolution is quick because we use a viscosity relevant for the upper mantle, leading to faster evolution than for a viscosity relevant in the deep mantle (Schaefer & Sasselov 2015). The use of an upper mantle viscosity here is reasonable since it is physically motivated from boundary-layer theory (see Section 2.1) and better matches Earth's current near steady state. Additionally, the regassing rate is limited because the mass of water in the hydrated layer cannot be greater than that on the surface, which as in Schaefer & Sasselov (2015) leads to a sharp decrease in the regassing rate at late times. Given that the degassing rate is also small due to the low temperatures, these rates balance to determine our steady states and hence waterworld boundary limits from Section 4.

We agree with the conclusion of Cowan & Abbot (2014) that super-Earths are unlikely to be waterworlds for both of the models with solely seafloor pressure-dependent and temperature-dependent water cycling. We also similarly find that conclusion is likely unaffected by parameter uncertainties. However, if these models themselves are less physical than a model with seafloor pressure dependencies dominating the degassing rate and temperature dependencies controlling regassing, that conclusion may change. In this hybrid model, the total mass fraction needed to become a waterworld is much smaller than that in both the Cowan & Abbot (2014) and Schaefer & Sasselov (2015) models. Ideally, future work will help distinguish between the viability of the three models considered here.

In this work, we considered three separate parameterizations for water cycling between ocean and mantle in order to make predictions for how they might affect exoplanet surface water abundances. We did so because the processes that control volatile cycling on Earth are not understood to the level of detail needed to make predictions for exoplanets with varying masses, total water mass fractions, compositions, and climates. Due to this, we utilized a simplified semi-analytic model and parameterized volatile cycling rates as either a power law in pressure or a function of temperature. In general, though this simplified model is powerful for understanding how a given process changes the surface water budget of the suite of exoplanets, studying surface water evolution on a given planet enables the use of more detailed coupling of parameterized convection and volatile cycling as in Sandu et al. (2011) and Schaefer & Sasselov (2015).

There also remain important parameters that do not have well-characterized dependencies with planet mass. Similarly to Schaefer & Sasselov (2015), we identified that the characteristic mantle viscosity is an important unknown in the problem because it can affect the evolutionary timescales. Additionally, the maximum mantle water mass fraction alone determines the waterworld boundary for the temperature-dependent model, and it is not known exactly how this should depend on planet mass. Understanding these parameters in detail will be necessary to make more detailed predictions of volatile cycling on terrestrial exoplanets.

It is clear that there is a dichotomy in the waterworld boundary based on whether or not one assumes that volatile cycling is temperature-dependent and/or pressure-dependent. As a result, understanding better which processes control volatile cycling on Earth is important to make more stringent predictions of whether terrestrial exoplanets should be waterworlds. Alternatively, observations with post-JWST-era instruments may be able to determine whether or not there is exposed land through either infrared spectra (if the atmosphere is not too optically thick, Abbot et al. 2012) or photometric observations over an entire planetary orbit in many wavelengths (Cowan et al. 2009; Kawahara & Fujii 2010; Cowan & Strait 2013; Cowan & Abbot 2014). This would serve as a test of the different volatile cycling parameterizations. If some super-Earths are found to have non-zero land fractions, the hybrid model considered here is not important or volatile delivery is inefficient for these objects. If, on the other hand, super-Earths are found to all be waterworlds, considering the combined effects of seafloor pressure-limited degassing and mantle convection may be necessary to explain volatile cycling on terrestrial planets.

In the future, one could use sophisticated multi-dimensional calculations of mantle convection including degassing through mid-ocean ridge volcanism and regassing through subduction of hydrated basalt, but this would be computationally expensive. However, these sophisticated calculations will not be worthwhile until the specific processes that govern volatile cycling on Earth and terrestrial exoplanets are understood in detail. We propose, then, that future observations of terrestrial exoplanets will be able to distinguish between the various water cycling models considered in this work. This could help constrain theories for water cycling on Earth and enable more sophisticated models to make predictions for the surface water inventory of individual planets. However, we must first understand in detail the effects of early water delivery and loss, and the effects of various tectonic regimes on water cycling itself.

In this section, we write down a time-dependent version of the model from Cowan & Abbot (2014), where degassing and regassing are regulated by seafloor pressure. The regassing and degassing rates in this case are

$$\begin{eqnarray}&&{w}_{\downarrow }={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}{d}_{{\rm{h}}}(P)\chi ,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{w}_{\uparrow }=x{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas}}(P),\end{eqnarray} \tag{ 33 }$$

equivalent to Equations (15) and (16). Here x_h is the mass fraction of water in the hydrated crust, ${\rho }_{{\rm{c}}}$ is the density of the oceanic crust, χ is the subduction efficiency, ${\rho }_{{\rm{m}}}$ is the density of the upper mantle, ${d}_{\mathrm{melt}}$ is the depth of melting below mid-ocean ridges, ${d}_{{\rm{h}}}$ is the hydrated layer depth, and ${f}_{\mathrm{degas}}$ is the degassing efficiency. To derive Equation (17), we start with Equation (14) and substitute Equations (32) and (33)

$$\begin{eqnarray}&&\displaystyle \frac{{dx}}{{dt}}=\displaystyle \frac{{L}_{\mathrm{MOR}}S}{{f}_{{\rm{M}}}M}\left[{x}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }{\left(\displaystyle \frac{P}{{P}_{\oplus }}\right)}^{\sigma }\right.\\ &&\,\left.-x{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{\left(\displaystyle \frac{P}{{P}_{\oplus }}\right)}^{-\mu }\right],\end{eqnarray} \tag{ 34 }$$

where ${L}_{\mathrm{MOR}}=3\pi {R}_{{\rm{p}}}$ is the mid-ocean ridge length, S is the average spreading rate of Earth ($\approx 10\ \mathrm{cm}\ {\mathrm{year}}^{-1}$), and we have used the power laws

$$\begin{eqnarray}&&{d}_{{\rm{h}}}(P)={d}_{{\rm{h}},\oplus }{\left(\displaystyle \frac{P}{{P}_{\oplus }}\right)}^{\sigma },\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{degas}}(P)={f}_{\mathrm{degas},\oplus }{\left(\displaystyle \frac{P}{{P}_{\oplus }}\right)}^{-\mu }.\end{eqnarray} \tag{ 36 }$$

Cowan & Abbot (2014) chose power laws to illustrate how different strengths of seafloor pressure-dependence would operate. In Equation (35), ${d}_{{\rm{h}},\oplus }$ is the hydration depth on Earth, P_⊕ is Earth's seafloor pressure, and in Equation (36) ${f}_{\mathrm{degas},\oplus }$ is the melt degassing fraction on modern Earth. Note that seafloor pressure $P=g{\rho }_{{\rm{w}}}{d}_{{\rm{w}}}$, where d_w is the ocean depth and ${\rho }_{{\rm{w}}}$ the density of water.

Cowan & Abbot (2014) relate seafloor pressure to mantle water mass fraction by

$$\begin{eqnarray}&&P={P}_{\oplus }{\tilde{g}}^{2}\displaystyle \frac{(\omega -{{xf}}_{{\rm{M}}})}{{\omega }_{0}\tilde{{f}_{{\rm{b}}}}},\end{eqnarray} \tag{ 37 }$$

where ${\omega }_{0}=2.3\times {10}^{-4}$ is the fractional mass of Earth's surface water and $\tilde{{f}_{{\rm{b}}}}={f}_{b}/{f}_{b,\oplus }=1.3$ is the ocean basin covering fraction normalized to that of Earth. Plugging this expression for P into Equation (34), we find

$$\begin{eqnarray}\displaystyle \frac{{dx}}{{dt}} & = & \displaystyle \frac{{L}_{\mathrm{MOR}}S}{{f}_{{\rm{M}}}M}\left[{x}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }{\left({\tilde{g}}^{2}\displaystyle \frac{(\omega -{{xf}}_{{\rm{M}}})}{{\omega }_{0}\tilde{{f}_{{\rm{b}}}}}\right)}^{\sigma }\right.\\ & & -\left.{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{{xf}}_{\mathrm{degas},\oplus }{\left({\tilde{g}}^{2}\displaystyle \frac{(\omega -{{xf}}_{{\rm{M}}})}{{\omega }_{0}\tilde{{f}_{{\rm{b}}}}}\right)}^{-\mu }\right].\end{eqnarray} \tag{ 38 }$$

Note that we can write Equation (38) using $\tilde{\omega }=\omega /({\omega }_{0}{\tilde{f}}_{{\rm{b}}})$ and $\tilde{x}={{xf}}_{{\rm{m}}}/({\omega }_{0}{\tilde{f}}_{{\rm{b}}})$ as

$$\begin{eqnarray}\displaystyle \frac{{dx}}{{dt}} & = & \displaystyle \frac{{L}_{\mathrm{MOR}}S}{{f}_{{\rm{M}}}M}[{x}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{\sigma }\\ & & -{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{{xf}}_{\mathrm{degas},\oplus }{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{-\mu }].\end{eqnarray} \tag{ 39 }$$

Non-dimensionalization of Equation (39) then gives

$$\begin{eqnarray}&&\displaystyle \frac{d\tilde{x}}{d\tau }={[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{\sigma }-{\tilde{X}}_{\oplus }^{-1}\tilde{x}{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{-\mu },\end{eqnarray} \tag{ 40 }$$

equivalent to Equation (17). In Equation (40),

$$\begin{eqnarray}&&{\tilde{X}}_{\oplus }=\displaystyle \frac{{x}_{h}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }{f}_{{\rm{M}}}}{{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}\end{eqnarray} \tag{ 41 }$$

is the non-dimensionalized mantle water mass fraction of Earth, $\tilde{\omega }=\omega /({\omega }_{0}{\tilde{f}}_{{\rm{b}}})$ is the non-dimensionalized total water mass fraction, $\tilde{g}=g/{g}_{\oplus }$, and

$$\begin{eqnarray}&&{\tau }_{\mathrm{CA}}=\tau =t\displaystyle \frac{{L}_{\mathrm{MOR}}{{Sx}}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}},\oplus }}{M{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}\end{eqnarray} \tag{ 42 }$$

is the non-dimensional time, which is inversely related to the seafloor overturning timescale $A/({L}_{\mathrm{MOR}}S)$.

In this section, we derive a simplified version of the Schaefer & Sasselov (2015) model, where volatile cycling rates are determined by the mantle temperature. The regassing and degassing rates in this case are

$$\begin{eqnarray}&&{w}_{\downarrow }={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}}}(T),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{w}_{\uparrow }={\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{f}_{\mathrm{melt}}(T)x,\end{eqnarray} \tag{ 44 }$$

equivalent to Equations (20) and (21). Here we have written the hydrated layer depth as a function of temperature. We have written ${f}_{\mathrm{degas}}$ as ${f}_{\mathrm{degas},\oplus }{f}_{\mathrm{melt}}(T)$ where ${f}_{\mathrm{melt}}(T)$ is the temperature-dependent melt fraction. Inserting Equations (43) and (44) into Equation (14), the dimensionful time-derivative of mantle water mass fraction is

$$\begin{eqnarray}&&\displaystyle \frac{{dx}}{{dt}}=\displaystyle \frac{{L}_{\mathrm{MOR}}S(T)}{{f}_{{\rm{m}}}M}[{x}_{{\rm{h}}}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}}{d}_{{\rm{h}}}(T)-{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{f}_{\mathrm{melt}}(T)x].\end{eqnarray} \tag{ 45 }$$

The functional forms of $S,{d}_{{\rm{h}}},{f}_{\mathrm{melt}}$ are developed in Section 2.3 of Schaefer & Sasselov (2015). Here we simplify them in order to obtain an analytically tractable version of Equation (45). First, the spreading rate is defined as

$$\begin{eqnarray}&&S=2{u}_{\mathrm{conv}}=2\displaystyle \frac{5.38\kappa h}{{\delta }^{2}},\end{eqnarray} \tag{ 46 }$$

The boundary-layer thickness δ is

$$\begin{eqnarray}&&\delta =h{\left(\displaystyle \frac{{\mathrm{Ra}}_{\mathrm{crit}}}{\mathrm{Ra}}\right)}^{\beta }.\end{eqnarray} \tag{ 47 }$$

Substituting δ from Equation (47),

$$\begin{eqnarray}S & = & \displaystyle \frac{10.76\kappa }{h}{\left(\displaystyle \frac{\mathrm{Ra}}{{\mathrm{Ra}}_{\mathrm{crit}}}\right)}^{(2\beta )}\\ & = & 10.76{\kappa }^{(1-2\beta )}{h}^{(6\beta -1)}{\left(\displaystyle \frac{\alpha {\rho }_{{\rm{m}}}g(T-{T}_{{\rm{s}}})}{\eta (T,x){\mathrm{Ra}}_{\mathrm{crit}}}\right)}^{2\beta }.\end{eqnarray} \tag{ 48 }$$

The hydration depth (the depth to which rock can be serpentinized) is defined as

$$\begin{eqnarray}&&{d}_{{\rm{h}}}=k\displaystyle \frac{({T}_{\mathrm{serp}}-{T}_{{\rm{s}}})}{{F}_{{\rm{m}}}}.\end{eqnarray} \tag{ 49 }$$

Using the mantle heat flux from Equation (7), we find

$$\begin{eqnarray}{d}_{{\rm{h}}} & = & h\displaystyle \frac{({T}_{\mathrm{serp}}-{T}_{{\rm{s}}})}{(T-{T}_{{\rm{s}}})}{\left(\displaystyle \frac{{\mathrm{Ra}}_{\mathrm{crit}}}{\mathrm{Ra}}\right)}^{\beta }\\ & = & {h}^{(1-3\beta )}{(T-{T}_{{\rm{s}}})}^{-(1+\beta )}({T}_{\mathrm{serp}}-{T}_{{\rm{s}}}){\left(\displaystyle \frac{\eta (T,x)\kappa {\mathrm{Ra}}_{\mathrm{crit}}}{\alpha {\rho }_{{\rm{m}}}g}\right)}^{\beta }.\end{eqnarray} \tag{ 50 }$$

Lastly, we use the same expression for the melt fraction as Schaefer & Sasselov (2015), which relates the melt fraction to mantle temperature through a power law, taking into account the solidus depression of the wet mantle

$$\begin{eqnarray}&&{f}_{\mathrm{melt}}={\left(\displaystyle \frac{T-{T}_{\mathrm{sol},\mathrm{wet}}(x)}{{T}_{\mathrm{liq},\mathrm{dry}}-{T}_{\mathrm{sol},\mathrm{dry}}}\right)}^{\theta }.\end{eqnarray} \tag{ 51 }$$

Here, we take ${T}_{\mathrm{liq},\mathrm{dry}}\approx 1498\ {\rm{K}}$, ${T}_{\mathrm{sol},\mathrm{dry}}\approx 1248\ {\rm{K}}$ as constants, and ${T}_{\mathrm{sol},\mathrm{wet}}={T}_{\mathrm{sol},\mathrm{dry}}-{{Kx}}^{\gamma }$, assuming that the mass fraction of water in melt is the same as the mass fraction of water in the mantle. We assume so because the partitioning coefficient of water in the mantle is thought to be extremely small ($\approx 1 \%$). Plugging in Equations (48), (50), and (51) into (45) and non-dimensionalizing gives

$$\begin{eqnarray}\displaystyle \frac{d\tilde{x}}{d{\tau }_{\mathrm{SS}}} & = & {\tilde{f}}_{{\rm{w}}}^{\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{\beta -1}\exp \left[\displaystyle \frac{-\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & -\tilde{{\rm{\Pi }}}{\tilde{f}}_{{\rm{w}}}^{2\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{2\beta }\exp \left[\displaystyle \frac{-2\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & \times \tilde{x}{(\tilde{T}-{\tilde{T}}_{\mathrm{sol},\mathrm{dry}}+\tilde{\lambda }{\tilde{x}}^{\gamma })}^{\theta },\end{eqnarray} \tag{ 52 }$$

equivalent to Equation (22). Here the non-dimensional solidus depression coefficient is $\tilde{\lambda }=\tilde{K}{({\omega }_{0}{\tilde{f}}_{{\rm{b}}}/{f}_{{\rm{m}}})}^{\gamma }$ and the degassing coefficient is $\tilde{{\rm{\Pi }}}={\rm{\Pi }}/D$, where

$$\begin{eqnarray}&&{\rm{\Pi }}={\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }\displaystyle \frac{{\omega }_{0}{\tilde{f}}_{{\rm{b}}}}{{f}_{{\rm{m}}}}{({\tilde{T}}_{\mathrm{liq},\mathrm{dry}}-{\tilde{T}}_{\mathrm{sol},\mathrm{dry}})}^{-\theta }.\end{eqnarray} \tag{ 53 }$$

The regassing coefficient (related to the hydrated layer depth) is

$$\begin{eqnarray}&&D={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}}{h}^{(1-3\beta )}\displaystyle \frac{({T}_{\mathrm{serp}}-{T}_{{\rm{s}}})}{{T}_{\mathrm{ref}}^{1+\beta }}{\left(\displaystyle \frac{\kappa {\mathrm{Ra}}_{\mathrm{crit}}{\eta }_{0}}{\alpha {\rho }_{{\rm{m}}}{{gf}}_{{\rm{w}}}(\tilde{x}=1)}\right)}^{\beta },\end{eqnarray} \tag{ 54 }$$

and ${\tau }_{\mathrm{SS}}={tD}/{\rm{\Sigma }}$, where

$$\begin{eqnarray}&&{\rm{\Sigma }}=M{\omega }_{0}{\tilde{f}}_{{\rm{b}}}\displaystyle \frac{{h}^{(1-6\beta )}{\kappa }^{(2\beta -1)}}{10.76{L}_{\mathrm{MOR}}}{\left(\displaystyle \frac{{\mathrm{Ra}}_{\mathrm{crit}}{\eta }_{0}}{\alpha {\rho }_{{\rm{m}}}{{gT}}_{\mathrm{ref}}{f}_{{\rm{w}}}(\tilde{x}=1)}\right)}^{2\beta }.\end{eqnarray} \tag{ 55 }$$

To ensure water mass balance in their time-dependent solutions, Schaefer & Sasselov (2015) force the hydrated layer to hold no more water than the surface itself. Formally, this ensures that

$$\begin{eqnarray}&&{x}_{h}{\rho }_{{\rm{m}}}\displaystyle \frac{4\pi }{3}({R}^{3}-{(R-{d}_{{\rm{h}}})}^{3})\leqslant M{\omega }_{0}{\tilde{f}}_{{\rm{b}}}(\tilde{\omega }-\tilde{x}).\end{eqnarray} \tag{ 56 }$$

Noting that we can re-write the hydrated layer depth from Equation (50) as

$$\begin{eqnarray}&&{d}_{{\rm{h}}}={D}_{2}\exp \left[\displaystyle \frac{\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\,{\tilde{f}}_{{\rm{w}}}^{-\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{-(1+\beta )},\end{eqnarray} \tag{ 57 }$$

where ${D}_{2}=D/({x}_{h}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}})$, we find a constraint for D₂ to ensure that the hydrated layer water mass is less than or equal to that on the surface:

$$\begin{eqnarray}&&{D}_{2}\leqslant \left[R-{\left({R}^{3}-\displaystyle \frac{3{\omega }_{0}{\tilde{f}}_{{\rm{b}}}M(\tilde{\omega }-\tilde{x})}{4\pi {x}_{{\rm{h}}}{\rho }_{{\rm{m}}}}\right)}^{1/3}\right]\\ &&\quad \times \exp \left[\displaystyle \frac{-\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\,{\tilde{f}}_{{\rm{w}}}^{\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{(1+\beta )}.\end{eqnarray} \tag{ 58 }$$

We force the constraint from Equation (58) in each timestep to ensure stability.³ Using the maximum value of D₂, we can find the maximum value of $\tilde{{\rm{\Pi }}}$ for use in Equation (22)

$$\begin{eqnarray}&&{\tilde{{\rm{\Pi }}}}_{\max }=\displaystyle \frac{{\rm{\Pi }}}{{D}_{2,\max }{x}_{h}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}}}.\end{eqnarray} \tag{ 59 }$$

Given the above models with either temperature or seafloor pressure-dependent volatile cycling rates, one can envision a model where surface water abundance is regulated by both seafloor pressure and mantle temperature. Here we consider a hybrid model where seafloor pressure regulates the degassing rate (because volcanism is less efficient with greater overburden pressure) and mantle temperature regulates the regassing rate (because serpentinization cannot occur if temperatures are too high). This hybrid model follows similarly from our derivations in Appendices A.1 and A.2. The regassing and degassing rates in this case are

$$\begin{eqnarray}&&{w}_{\downarrow }={x}_{{\rm{h}}}{\rho }_{{\rm{c}}}\chi {d}_{{\rm{h}}}(T),\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{w}_{\uparrow }=x{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas}}(P),\end{eqnarray} \tag{ 61 }$$

equivalent to Equations (32) and (44). Inserting these into Equation (14), we find the dimensional form of the time-derivative of the water mass fraction

$$\begin{eqnarray}&&\displaystyle \frac{{dx}}{{dt}}=\displaystyle \frac{S(T)}{{f}_{{\rm{m}}}M}[{x}_{{\rm{h}}}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}}{d}_{{\rm{h}}}(T)-{\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{{xf}}_{\mathrm{degas}}(P)].\end{eqnarray} \tag{ 62 }$$

We insert our prescriptions for S and ${d}_{{\rm{h}}}$ from Equations (48) and (50), respectively, and the seafloor pressure-dependence of ${f}_{\mathrm{degas}}$ from Equation (36) into Equation (62). Non-dimensionalizing, we find

$$\begin{eqnarray}\displaystyle \frac{d\tilde{x}}{d{\tau }_{\mathrm{hyb}}} & = & {\tilde{f}}_{{\rm{w}}}^{\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{(\beta -1)}\exp \left[\displaystyle \frac{-\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & -\tilde{E}{\tilde{f}}_{{\rm{w}}}^{2\beta }{(\tilde{T}-{\tilde{T}}_{{\rm{s}}})}^{(2\beta )}\exp \left[\displaystyle \frac{-2\beta }{{\tilde{T}}_{{\rm{m}}}}\left(\displaystyle \frac{1}{\tilde{T}}-1\right)\right]\\ & & \times \tilde{x}{[{\tilde{g}}^{2}(\tilde{\omega }-\tilde{x})]}^{-1},\end{eqnarray} \tag{ 63 }$$

equivalent to Equation (26). Here $\tilde{E}=E/D$, $E={\rho }_{{\rm{m}}}{d}_{\mathrm{melt}}{f}_{\mathrm{degas},\oplus }{\omega }_{0}{\tilde{f}}_{{\rm{b}}}/{f}_{{\rm{m}}}$, and ${\tau }_{\mathrm{hyb}}={\tau }_{\mathrm{SS}}={tD}/{\rm{\Sigma }}$.

As in to the solely temperature-dependent model, we restrict the hydrated layer depth using Equation (58). If ${D}_{2}={D}_{2,\max }$ the corresponding constraint on $\tilde{E}$ is

$$\begin{eqnarray}&&{\tilde{E}}_{\max }=\displaystyle \frac{E}{{D}_{2,\max }{x}_{h}{\rho }_{{\rm{c}}}{\chi }_{{\rm{r}}}}.\end{eqnarray} \tag{ 64 }$$