What Factors Affect the Duration and Outgassing of the Terrestrial Magma Ocean?

We calculate the thermal state of the solidifying MO without examining the preceding stage of its formation. We build a model that simulates the coupled evolution of the interior (Sections 2.2–2.7) and the atmosphere (Sections 2.8–2.9). The COnvective Magma Radiative Atmosphere and Degassing (COMRAD) model resolves the mantle interior profiles of temperature and liquid and solid fraction, along with the degassing process, starting from a fully molten mantle up to the end of the MO phase (see Section 2.10). We employ the mantle surface temperature iteration method developed by Lebrun et al. (2013), with differences in the calculation of the mantle adiabat (Section 2.4), of the liquid viscosity (Section 2.6), and of the volatile mass balance (Section 2.7). The outgassing of H₂O and CO₂ is calculated according to a melt solubility curve for each volatile (Section 2.7). The atmosphere is treated in two alternative ways (Section 2.8): (i) a gray atmosphere accounting for two greenhouse gases, H₂O and CO₂ (Abe & Matsui 1985; Elkins-Tanton 2008; Section 2.8.1) and (ii) a pure H₂O atmosphere with a spectrally resolved Outgoing Longwave Radiation at the Top Of the Atmosphere (henceforth named "OLR at TOA," OLR_TOA) by Katyal et al. (2019) in a companion paper (hereafter "companion paper"; Section 2.8.2). In the following sections, we separately introduce each model component.

We consider a spherically symmetric Earth with outer radius R_p and core radius R_b. This yields a mantle of thickness R_p − R_b, the physical properties of which are defined by the melting curves (solidus and liquidus) of KLB-1 peridotite (Figure 1). By comparing the interior thermal profile to these curves, we identify the phase (liquid, partially molten, or solid) of each mantle layer. The mantle is initially assumed to be fully molten and convecting, with an adiabatic temperature profile. As it cools, the liquid adiabat (dotted line in Figure 1) intersects with the melting curves (solid lines in Figure 1). Due to the steeper slope of the adiabat compared to the melting curves, the adiabat and the liquidus intersect first atop the core–mantle boundary (CMB). The mantle is fully molten from the surface until the depth of intersection between the liquid adiabat and the liquidus, it is partially molten between the liquidus and solidus, and completely solid below the solidus. The melt fraction ϕ at any given depth is calculated as (e.g., Solomatov & Stevenson 1993a, 1993b; Abe 1997; Solomatov 2007; Lebrun et al. 2013)

$$\begin{eqnarray}&&\phi =\displaystyle \frac{T-{T}_{\mathrm{sol}}}{{T}_{\mathrm{liq}}-{T}_{\mathrm{sol}}},\end{eqnarray} \tag{ 1 }$$

where T is the temperature of the mantle at a given depth, and T_sol and T_liq the corresponding solidus and liquidus temperatures. The partially molten region is further divided by comparing the melt fraction ϕ with the critical melt fraction ϕ_C of 40% that separates liquid-like from solid-like behavior (Costa et al. 2009). For ϕ > ϕ_C, the region is considered liquid-like and belongs to the MO-convecting domain of depth D. Note that ϕ_C varies among 30% (e.g., Hier-Majumder & Hirschmann 2017; Maurice et al. 2017), 40% (Solomatov 2007; Bower et al. 2017), and 50% (Monteux et al. 2016; Ballmer et al. 2017) in the geodynamic literature.

Zoom In Zoom Out Reset image size

We use solidus and liquidus curves of KLB-1 peridotite obtained from experimental data. Depending on pressure, we adopt different parameterizations for different parts of the mantle. For the solidus, we use data from Hirschmann (2000) for P ∈ [0, 2.7) GPa, Herzberg et al. (2000) for P ∈ [2.7, 22.5) GPa, and Fiquet et al. (2010) for P ≥ 22.5 GPa, while for the liquidus, from Zhang & Herzberg (1994) for P ∈ [0, 22.5) GPa and from Fiquet et al. (2010) for P ≥ 22.5 GPa. Because we employ data from multiple studies, we refer to the resulting set of melting curves as "synthetic." Such curves are adopted in our experiments unless otherwise specified.

As shown in Figure 1, we also tested the linear melting curves adopted by Abe (1997) and later by Lebrun et al. (2013), as well as those introduced by Andrault et al. (2011) that are representative of a chondritic composition (for analytical expressions for all melting curves, see Appendix A). Apart from the linear curves of Abe (1997), the experimental data that we considered require higher order polynomial fittings because their slope is not constant with depth. As we discuss in Section 2.4, this imposes a limitation in calculating the two-phase adiabat.

The interior temperature profile is calculated using the expression of the adiabatic temperature gradient for a one-phase system:

$$\begin{eqnarray}&&\displaystyle \frac{{dT}}{{dP}}=\displaystyle \frac{{\alpha }_{T}T}{\rho {c}_{P}},\end{eqnarray} \tag{ 2 }$$

where P is the pressure in GPa, c_P the thermal capacity at constant pressure, and α_T the pressure-dependent thermal expansivity given by (Abe 1997)

$$\begin{eqnarray}&&{\alpha }_{T}(P)={\alpha }_{0}{\left[\displaystyle \frac{{PK}^{\prime} }{{K}_{0}}\right]}^{-(m-1+K^{\prime} )/K^{\prime} },\end{eqnarray} \tag{ 3 }$$

where α₀ is the surface expansivity, K₀ the surface bulk modulus and K' its pressure derivative (see Appendix D, Table 4 for their values), and m = 0. The pressure is simply calculated assuming a hydrostatic profile.

Solomatov & Stevenson (1993a, 1993b) derived the adiabat for a two-phase system by introducing a modified thermal expansivity and thermal capacity that depends on the melt fraction. However, the expressions they derived are explicitly valid for a system with a constant rate of temperature drop with depth, which is equivalent to a constant phase boundary slope and thus applies only to linear melting curves such as those of Abe (1997). The modified adiabat then tends to align with the slope of constant melt fraction. Because they do not cover the higher order parameterizations of the experimental data that we adopted, we employ instead the one-phase adiabat of Equation (3) with constant thermal capacity and pressure-dependent expansivity.

Assuming that the mantle temperature profile T(r) is adiabatic as described in Section 2.4, the time evolution of the MO is obtained by integrating the energy balance equation (Abe 1997) over the evolving MO volume:

$$\begin{eqnarray}&&\rho \left({c}_{P}+{\rm{\Delta }}H\displaystyle \frac{d\phi }{{dT}}\right)\displaystyle \frac{{dT}}{{dt}}=-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial ({r}^{2}{F}_{\mathrm{conv}})}{\partial r}+\rho {q}_{r},\end{eqnarray} \tag{ 4 }$$

where ρ is the density, ΔH the specific enthalpy difference due to phase change, F_conv the MO convective cooling flux, and q_r the internal heat released by the decay of the radioactive elements (see Appendix G).

Because of its large depth extent and liquid-like viscosity which approaches that of water (see Section 2.6), the MO is expected to undergo highly turbulent convection (Solomatov 2007) that is neither attainable in the laboratory (Shishkina 2016) nor is numerically resolvable (Maas & Hansen 2015). Key to the evolution of the interior temperature profile is the parameterization of the convective heat flux F_conv of Equation (4). This is calculated with the aid of the Rayleigh (Ra) and the Prandtl (Pr) numbers:

$$\begin{eqnarray}&&{Ra}=\displaystyle \frac{\rho {\alpha }_{T}g({T}_{p}-{T}_{\mathrm{surf}}){D}^{3}}{{\kappa }_{T}\eta },\quad {\Pr }=\displaystyle \frac{\eta }{\rho {\kappa }_{T}},\end{eqnarray} \tag{ 5 }$$

where g is the gravity acceleration, T_p the mantle potential temperature, T_surf the surface temperature, D the depth of the convective layer, κ_T the thermal diffusivity, and η the dynamic viscosity.

We consider two different parameterizations for "soft" and "hard" turbulence. In the first one, turbulent dissipation at the boundary layers affects the heat flux, while the flow is laminar in the bulk of the fluid (Solomatov 2007):

$$\begin{eqnarray}&&{F}_{\mathrm{soft}}=0.089\displaystyle \frac{{k}_{T}({T}_{p}-{T}_{\mathrm{surf}}){{Ra}}^{1/3}}{D},\end{eqnarray} \tag{ 6 }$$

where k_T = κ_Tρc_p is the thermal conductivity. With the above formulation, the heat flux becomes independent of the depth of the convective layer because the Rayleigh number is proportional to D³. The second parameterization depends additionally on the inverse of the Prandtl number. It represents a regime where the heat flux is assumed to be controlled not only by boundary layer friction but also has a contribution from turbulence generated in the bulk volume of the fluid (Solomatov 2007):

$$\begin{eqnarray}&&{F}_{\mathrm{hard}}=0.22\displaystyle \frac{{k}_{T}({T}_{p}-{T}_{\mathrm{surf}})}{D}{{Ra}}^{2/7}{{\Pr }}^{-1/7}{\lambda }^{-3/7},\end{eqnarray} \tag{ 7 }$$

where λ = basin length/depth is the aspect ratio of the mean flow. For very high Ra, the hard turbulence parameterization is suggested (Solomatov 2007). There, increasing values of Pr in a progressively more viscous fluid yields a lower heat flux if all other parameters are left unchanged. Both parameterizations were implemented and tested.

The melt viscosity prominently factors into the thermal evolution of the MO (Equations (5) and (6)). We separate the mantle into two regimes, namely liquid-like and solid-like. The transition between them upon cooling and solidification is a complex phenomenon that depends, among other factors, on the composition and cooling rate (Speedy 2003). Although the viscosity dependence on temperature follows an Arrhenius law below the solidus temperature (Kobayashi et al. 2000), nonlinear effects take place near and above it (Dingwell 1996; Kobayashi et al. 2000; Speedy 2003). In order to mitigate the transition to solid state, Salvador et al. (2017) proposed a "smoothening" of the sharp viscosity jump that occurs in earlier MO models (e.g., Lebrun et al. 2013; Schaefer et al. 2016). However, during cooling, the crystal content increases, and the melt rheology is expected to make a discontinuous jump from the liquid- to solid-like state over a short crystallinity range (Marsh 1981). Continuous variation in viscosity across five orders of magnitude is seen only in the case of glasses (Kobayashi et al. 2000), although such behavior is not consistent with common MO solidification assumptions.

The Vogel–Fulcher–Tammann, henceforth referred to as VFT, equation is employed in our work. By the definition of VFT, the obtained viscosity tends to an infinite value (hence not physically meaningful) at a threshold temperature T = C. We consider two such expressions for the liquid dynamic viscosity: one that depends only on the temperature η_l = f(T) (Karki & Stixrude 2010) and a second one that depends on both temperature and water content ${\eta }_{l}=f(T,{X}_{{{\rm{H}}}_{2}{\rm{O}}})$ (Giordano et al. 2008).

Karki & Stixrude (2010) found that for hydrous melts, the viscosity at a given potential temperature can be well fitted to the following VFT equation:

$$\begin{eqnarray}&&{\eta }_{l}(T)={A}_{K}\exp \left[\displaystyle \frac{{B}_{K}}{T-{C}_{K}}\right],\end{eqnarray} \tag{ 8 }$$

where A_K, B_K, and C_K (see Appendix D, Table 4 for their values) are calculated for a fixed water content of 10 wt%.

However, Equation (8) does not explicitly include the effect of water content, which is expected to vary during the simulation (see Section 4.2). The presence of water tends to lower the melt viscosity (Marsh 1981; Dingwell 1996; Giordano et al. 2008; Karki & Stixrude 2010), and it is important to include it as a time-dependent variable in our calculations. Hence, we implemented the empirical model of Giordano et al. (2008), which explicitly uses the water concentration, together with the concentration of 13 different oxides in the silicate melt. It calculates two of the three VFT parameters (B_G and C_G in Equation (9) below). The viscosity for a given temperature T and water concentration ${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ is then given by

$$\begin{eqnarray}&&{\eta }_{l}(T,{X}_{{{\rm{H}}}_{2}{\rm{O}}})={10}^{{A}_{G}+\left[\tfrac{{B}_{G}}{T-{C}_{G}}\right]},\end{eqnarray} \tag{ 9 }$$

where the parameter A_G = −4.55 ± 1 is a constant pre-exponential factor. The parameters ${B}_{G}({X}_{{{\rm{H}}}_{2}{\rm{O}}})$ and ${C}_{G}({X}_{{{\rm{H}}}_{2}{\rm{O}}})$ are calculated by the model at each time step according to the evolving concentration of water (see Section 2.7). Even though COMRAD is unable to resolve the evolution of the melt composition, it does resolve the melt water concentration with time, which allows us to evaluate the Giordano et al. (2008) model at each time step. In order to use this model, it is required that we choose a suitable composition as a constant, non-evolving basis. We found that the composition of basanite (Giordano & Dingwell 2003; Giordano et al. 2008) is able to reproduce the experimental values for the temperature-dependent viscosity for the anhydrous case (Urbain et al. 1982) because it is one of the least evolved in terms of silicate content. After calibration of the prefactor A_G (see Appendix F, Table 6), the error lies within ±10%. This fitting provides us with a parameterized description of the composition that allows us to treat the decrease of the melt viscosity with increasing water content (Figure 2(A)). For anhydrous melt, the viscosity calculated with Equation (9) and ${X}_{{{\rm{H}}}_{2}{\rm{O}}}=0$ yields similar values to those proposed by Karki & Stixrude (2010) at high temperatures, though the two tend to depart significantly for temperatures near the solidus (Figure 2(B)).

Zoom In Zoom Out Reset image size

The pressure dependence of the viscosity for the hydrous case is not explicitly provided by Karki & Stixrude (2010). The authors report that it varies by a relatively small factor between 2.5 and 10 over the pressure range [0, 140] GPa spanned by a global terrestrial MO. In the following, for simplicity, we will ignore such dependence and use the liquid viscosity evaluated at the potential temperature T_p of the MO as representative of the fully molten part.

The melt viscosity η_l is further corrected for the crystal fraction content in each layer. In the liquid-like partially molten region, the effect of crystals is taken into account with the following expression (Roscoe 1952):

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\eta }_{l}}{{\left(1-\tfrac{1-\phi }{1-{\phi }_{C}}\right)}^{2.5}}.\end{eqnarray} \tag{ 10 }$$

By combining Equation (8) or (9) with Equation (10) for each layer that belongs to the MO, we obtain the volumetric harmonic mean viscosity η that is then used in the calculation of the parameterized convective heat flux in Equation (6). The effect of the layer with the lowest viscosity value is prioritized in the calculation of the harmonic mean viscosity and sets the leading order of magnitude for the value used in the calculation of the Rayleigh number (5).

The viscosity employed during the MO lifetime is the liquid-like viscosity (η_l). Note that the solid-like viscosity (η_s) is employed only after the MO phase ends, and it is expressed following the Karato & Wu (1993) formulation. The viscosity of the partially molten solid-like region below the rheology front (i.e., for ϕ < ϕ_C) is modified by the presence of crystals according to Solomatov (2007).

Along with the thermal evolution, the concentration of volatiles in the mantle and their outgassing into the atmosphere are calculated. We use solubility curves to calculate the concentration and gas pressure in the melt for each volatile. For H₂O, we use (Caroll & Holloway 1994)

$$\begin{eqnarray}&&{P}_{\mathrm{sat},{{\rm{H}}}_{2}{\rm{O}}}={\left(\displaystyle \frac{{X}_{{{\rm{H}}}_{2}{\rm{O}}}}{6.8\cdot {10}^{-8}}\right)}^{(1/0.7)},\end{eqnarray} \tag{ 11 }$$

and for CO₂ (Pan et al. 1991),

$$\begin{eqnarray}&&{P}_{\mathrm{sat},{\mathrm{CO}}_{2}}=\displaystyle \frac{{X}_{{\mathrm{CO}}_{2}}}{4.4\,\cdot \,{10}^{-12}}.\end{eqnarray} \tag{ 12 }$$

In this way, we obtain the saturation vapor pressure over melt with a given volatile concentration. Due to efficient mixing, the volatile concentration is homogeneous throughout the MO.

Upon solidification, part of the volatile budget remains in the solid mantle according to the partition coefficients of lherzolite for the upper mantle (κ_vol,lhz) and of perovskite for the lower mantle (κ_vol,pv; see Appendix D, Table 4 for their values). By calculating the volatile content stored in the liquid and solid phases of the mantle, we estimate the mass balance for each volatile at each time iteration t as follows:

$$\begin{eqnarray}&&\begin{array}{l}{M}_{l,{t}_{0}}{X}_{\mathrm{vol},{t}_{0}}=P({X}_{\mathrm{vol},t})\displaystyle \frac{4\pi {R}_{p}^{2}}{g}\\ \quad +\ {M}_{s,{\rm{pv}}}{\kappa }_{\mathrm{vol},\mathrm{pv}}{X}_{\mathrm{vol},t}+{M}_{s,{\rm{lhz}}}{\kappa }_{\mathrm{vol},\mathrm{lhz}}{X}_{\mathrm{vol},t}\\ \quad +\ {M}_{l,z\lt {z}_{\mathrm{RF}}}{X}_{\mathrm{vol},t}+{M}_{l,z\gt {z}_{\mathrm{RF}}}{X}_{\mathrm{vol},t}{\rm{\Pi }}\\ \quad +\ {M}_{l,z\gt {z}_{\mathrm{RF}}}{X}_{\mathrm{vol},t-{dt}}(1-{\rm{\Pi }}),\end{array}\end{eqnarray} \tag{ 13 }$$

where ${M}_{l,{t}_{0}}$ is the initial (time = t₀) mass of the liquid mantle, ${X}_{\mathrm{vol},{t}_{0}}$ the initial volatile concentration in the melt, P(X_vol,t) the saturation pressure of the volatile for the respective concentration X_vol,t at time t, M_s,pv and M_s,lhz are the masses of solid mantle in perovskite and in lherzolite, respectively, M_l the mass of the melt at depth z, which is either shallower (z < z_RF) or deeper than the rheology front (z > z_RF), and X_vol,t−dt the concentration of the volatile in the previous time step. Comparing the melt percolation velocity (Solomatov 2007) to the rheology front velocity, we calculate a volumetric fraction Π of the total melt volume that upwells across the rheology front. Π takes values within 0 and 1. The last term on the right-hand side (RHS) of Equation (13) represents the volatile mass trapped in the liquid below the rheology front and is evaluated at the concentration of the previous time step.

Equation (13) combined with either Equation (11) or (12) forms a system of two equations in two unknowns (P, X) for each species. Equation (11) is nonlinear with respect to ${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ and is solved iteratively with 0.01 bar tolerance.

The thermal evolution of the system is coupled to the outgassing process. Upon cooling, the mantle volume is redistributed into a solid and a liquid phase (see Figure 1). Consequently, the masses of the solid and liquid reservoirs M_s and M_l where the volatiles are stored change continuously. For every new mantle layer that solidifies, volatile enrichment is ensured in the remaining melt. As the saturation pressure increases monotonically with concentration (Equations (11) and (12)), the equilibrium gas pressure also increases, resulting in a progressive build-up of atmospheric mass at the surface of the planet.

We adopt two alternative approaches to model the atmosphere generated upon MO outgassing: (i) a gray approximation after Abe & Matsui (1985) that treats two gas species H₂O and CO₂ and (ii) a line-by-line approach that calculates a spectrally resolved OLR (companion paper) that assumes a pure H₂O vapor composition.

2.8.1. Gray Atmospheric Model

The gray approximation that we use is derived in Abe & Matsui (1985). It considers the absorption of thermal radiation independently of the wavelength. Both outgassed species H₂O and CO₂ absorb significantly in the spectral region where thermal energy is emitted from the Earth surface and are therefore greenhouse contributors. By absorbing radiative energy, they exert direct control on the surface temperature. Water is the more potent greenhouse agent of the two under normal atmospheric conditions (P₀ = 101,325 Pa, T₀ = 293 K; see, e.g., Pierrehumbert 2010). For the H₂O absorption coefficient, the value ${k}_{0,{{\rm{H}}}_{2}{\rm{O}}}=0.01$ m² kg⁻¹ in the mid-infrared window region 1000 cm⁻¹ is adopted, after Abe & Matsui (1988). CO₂ is accounted for with absorption coefficient ${k}_{0,{\mathrm{CO}}_{2}}=0.001$ m² kg⁻¹ (Yamamoto 1952). A higher value ${k}_{0,{\mathrm{CO}}_{2}}=0.05$ m² kg⁻¹ has been employed by Elkins-Tanton (2008; along with lower ${k}_{0,{\mathrm{CO}}_{2}}$ values) and by Lebrun et al. (2013), which was calculated by Pujol & North (2003) in order to reproduce the present-day Earth climate sensitivity (ECS). ECS refers to the combined response of the climate system to the radiative forcing from doubling the atmospheric CO₂ abundance relative to its pre-industrial levels and corresponds to an increase of about 2 °C in surface mean temperature (Flato et al. 2013). Using it is a good practice for studying the role of CO₂ radiative forcing on today's temperate Earth climate, within which the water vapor is not saturated over the atmospheric column. The presence of a liquid ocean is a strong constraint on the ECS and affects the overlying atmospheric profile through the intercomponent exchange of vapor or "hydrological cycle" (Held & Soden 2006), provided that no RG regime (Pierrehumbert 2010) ensues. Extrapolating the climate sensitivity to MO mean surface temperature (>1000 K) well above today's (300 K) is unsuitable for our study. We therefore avoid using the ECS-based value for ${k}_{0,{\mathrm{CO}}_{2}}$ because it could overestimate CO₂'s radiative forcing on a planet with a qualitatively different surface and atmospheric dynamics.

In Abe & Matsui (1985), the downward radiation at the TOA is set to the incoming stellar flux F_Sun, which depends on the incident radiation S₀ at the assumed orbital distance. It relates to the blackbody equilibrium temperature T_eq of the planet through the Stefan–Boltzmann law:

$$\begin{eqnarray}&&{F}_{\mathrm{Sun}}=\left(1-\alpha \right)\displaystyle \frac{{S}_{0}}{4}=\sigma {T}_{\mathrm{eq}}^{4},\end{eqnarray} \tag{ 14 }$$

where α is the albedo and σ the Stefan–Boltzmann constant. The resulting net upward flux at the top of the atmosphere (F_gray) is given by (Abe & Matsui 1985)

$$\begin{eqnarray}&&{F}_{\mathrm{grey}}=\sigma \epsilon \left({T}_{\mathrm{surf}}^{4}-{T}_{\mathrm{eq}}^{4}\right)={F}_{\mathrm{conv}}.\end{eqnarray} \tag{ 15 }$$

According to Equation (15), the net radiative flux at the TOA is positive for T_surf > T_eq. We adopt the convention of positive flux to represent planetary cooling. In order to find a state of the system that satisfies the energy balance, assuming that the radiative atmospheric adjustment is instantaneous, we require that the convective heat flux F_conv at the top of the MO is equal to the flux at the TOA.

For a given potential temperature, we solve the system of Equations (6) and (15) using an iterative scheme built according to the method of Lebrun et al. (2013) with an accuracy of 10⁻² W m⁻².

2.8.2. Line-by-line Atmospheric Model

An alternative approach is to employ a line-by-line (lbl) code to calculate the atmospheric outgoing radiation. The respective model is described in the companion paper. The advantage of this approach is that it provides a detailed calculation of wavelength-dependent longwave emission. We assume a 100% water vapor atmosphere ("steam atmosphere") as commonly used by various authors when treating MO planets or exoplanets with water-dominated atmospheres (e.g., Hamano et al. 2013; Massol et al. 2016; Schaefer et al. 2016). A dry adiabatic temperature profile is used for the troposphere when the surface temperature is above the critical point of water (${T}_{{{\rm{H}}}_{2}{\rm{O}},{\rm{crit}}}$ = 647 K). When the dry adiabat intersects the saturation vapor pressure of water, a moist adiabatic temperature profile is assumed (Kasting 1988).

Using the temperature profiles of the atmosphere for a range of surface temperatures T_surf and surface pressures ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$, the emitted radiation is calculated for the spectral range from 20–29,995 cm⁻¹ using the lbl model GARLIC (Schreier et al. 2014) with HITRAN 2012 (Rothman et al. 2013). Integrating the emitted radiation at the TOA (corresponding to pressure 1 Pa), the outgoing radiation flux OLR_TOA is obtained.

Katyal et al. (2019) determines OLR_TOA for various H₂O surface pressures and temperatures on a (${P}_{{{\rm{H}}}_{2}{\rm{O}},0},{T}_{\mathrm{surf}}$) grid that covers pressures between 4 and 300 bar and temperatures between 650 and 4000 K (Appendix C.1). For the surface vapor pressure, we use the values 4, 25, 50, 100, 200, and 300 bar. For the surface temperature, we sample our calculations with a resolution ΔT = 100 K. We employ ΔT = 20 K in the region T ∈ [1400, 1800] K where the highest variability of outgassing with respect to surface temperature occurs for the melting curves used in this study. We obtain values of OLR_TOA that correspond to conditions intermediate to the grid points via bilinear interpolation. The relative interpolation error ranges from 10% for fluxes of the order of 10⁶ W m⁻² to about 1% for fluxes of the order of 10² W m⁻².

In order to couple the lbl model results with the MO thermal evolution, we impose a balance between the net energy flux at TOA and the MO cooling flux F_conv such that

$$\begin{eqnarray}&&{F}_{\mathrm{conv}}={{\rm{OLR}}}_{\mathrm{TOA}}-{F}_{\mathrm{Sun}}.\end{eqnarray} \tag{ 16 }$$

Equations (6) and (16) form a system of two unknowns T_surf and F_conv, which we solve iteratively with a tolerance of 10⁻¹ W m⁻². The resulting flux needed to balance the RHS of Equation (16) can be either positive or negative, corresponding to a cooling or warming case, respectively.

For the young Sun, we used a lower irradiation value following the expression for the time dependence of the solar constant of Gough (1981); otherwise, we used today's value (1361 Wm⁻²; see Table 1).

Table 1.  Parameter Values and Components of the Reference-A Model

Parameter	Description	Value/Type	Unit/Info
Atmosphere	Type of approximation	Gray	Equation (15); Abe & Matsui (1985)
${k}_{0,{{\rm{H}}}_{2}{\rm{O}}}$	Absorption coeff. at normal atmospheric conditions	0.01	m2 kg−1
${k}_{0,{\mathrm{CO}}_{2}}$	Absorption coeff. at normal atmospheric conditions	0.001	m2 kg−1
H2O content	Total water reservoir	300	bar
${X}_{{{\rm{H}}}_{2}{\rm{O}},0}$	Initial H2O mantle abundance	410	ppm
CO2 content	Total CO2 reservoir	100	bar
${{\rm{X}}}_{{\mathrm{CO}}_{2},0}$	Initial CO2 mantle abundance	130	ppm
S	Solar constant (S0)	1361	W m−2
α	Planetary albedo	0.30	#
Tp,0	Initial potential temperature	4000	K
D	Initial MO depth	2890	km
ηl	Melt viscosity parameterization	ηl = f(T)	Equation (8)
qr	Radioactive heating	0	Not included
tplanet	Planet accretion time	100	Myr (employed for qr only)
Tsol, Tliq	Melting curves	"Synthetic"	Herzberg et al. (2000); Hirschmann (2000);
			Fiquet et al. (2010) (Section 2.3)
TRF,0	Temperature of rheology front at z = 0	1645	K
Fconv	Convective heat flux parameterization	Soft turbulence	Equation (6); Solomatov (2007)

Note. Use this table to complement the results presented in Table 2.

Download table as:  ASCIITypeset image

Calculating the planetary albedo is outside the scope of this study. α is instead used as an input parameter. The suggested albedo for a cloudless steam atmosphere lies within the range [0.15, 0.40] (Kasting 1988; Goldblatt et al. 2013; Leconte et al. 2013; Pluriel et al. 2019). We employ the value 0.30 unless otherwise stated.

The end of the MO phase is defined as the point in time when the rheology front reaches the surface. At that stage, the mantle adiabat has potential temperature T_RF,0 such that all mantle layers have a melt fraction lower than ϕ_C, a condition termed "mush stage" by Lebrun et al. (2013). Although some melt still remains enclosed in the solid matrix, the mantle subsequently behaves as a solid. Moreover, the adiabatic profile used for the solid mantle implies that solid-state convection has fully developed. Although this is a likely scenario for slowly solidifying MOs, establishing whether and to what extent the solid mantle convects during its early stages is beyond the scope of this study as this would require the use of fully dynamic simulations (e.g., Ballmer et al. 2017; Maurice et al. 2017). We therefore consider the established convection considered in this work to be an end-member case within the geodynamic assumptions.

However, we stress that the thermal evolution model is designed to cover only the time until the MO end is reached via bottom-up crystallization.

In a similar fashion to prior studies of the MO solidification (Zahnle et al. 1988; Abe 1997; Elkins-Tanton 2008; Hamano et al. 2013, 2015; Lebrun et al. 2013; Monteux et al. 2016; Schaefer et al. 2016; Hier-Majumder & Hirschmann 2017), we present the thermal evolution using the state variables: surface temperature, potential temperature, heat flux, Ra number, and MO depth evolution. This enables both model validation and comparison. We adopt the multipanel approach of Lebrun et al. (2013), which is convenient for comparisons between varying modeling approaches. We performed four simulations for the following cases: (i) absence of atmosphere, referred to as the blackbody "bb" case, (ii) a gray atmosphere composed of both H₂O and CO₂, "gr-H₂O/CO₂," (iii) a gray atmosphere composed of only H₂O, "gr-H₂O," and (iv) a H₂O atmosphere treated with an lbl model (Figure 3). Apart from the representation of the atmosphere or absence thereof, all aspects of the model follow the Ref-A case (Table 1).

Zoom In Zoom Out Reset image size

Commonly in all simulations, T_p and T_surf coevolve until an abrupt difference between the two marks the end of the MO (Figure 3(B), dashed lines); the liquid-like behavior comes to an end and a layer with a melt fraction of 40% or lower remains. The average viscosity increases by more than eight orders of magnitude across the critical melt fraction, taking values from 10⁸ to 10¹⁸ Pa s (not shown). A smooth variation across this interval would be difficult to justify under the assumption of fractional crystallization (e.g., Marsh 1981). At this point, the whole domain switches to the low cooling flux that characterizes solid-like convection (Equation (6)) and the surface temperature drops abruptly, while the potential temperature remains unaffected.

The "bb" thermal evolution compares well with that presented by Lebrun et al. (2013). It demonstrates the highest T_p−T_s difference. The mantle consequently cools rapidly (0.002 Myr) with the highest convective flux (F = 5 · 10⁶–10⁴ W m⁻²), caused by this large temperature difference. Longer solidification times (0.15 Myr) are found by Monteux et al. (2016), who assume a bb radiative cooling (F = 10⁵–10² W m⁻²) but a different interior model with a heat contribution from the core. The bb case is only relevant for planets that lose their outgassed atmosphere instantaneously.

For a planet that retains its atmosphere, the gray approximations show that the presence of the additional greenhouse species CO₂ contributes only 0.05 Myr to the solidification time and is less significant in comparison to water (0.16 Myr versus 0.21 Myr MO duration). The longer solidification time (≈0.4 Myr) obtained by Lebrun et al. (2013) for their gray two-species case is consistent with absorption coefficient ${k}_{0,{\mathrm{CO}}_{2}}$ = 0.05, which is likely to be rather high for those climates (see Section 2.8.1). The gray approach employed in this study and theirs follows Abe & Matsui (1985) and should not be identified with the other gray models used in the literature: the study of Hamano et al. (2013) expands on the Nakajima et al. (1992) gray model and employs supercritical water thermal capacities. Hier-Majumder & Hirschmann (2017) formulated a hybrid energy balance for the atmosphere, employing elements from both Abe & Matsui (1985) and Hamano et al. (2013). For potential temperature equal to the equilibrium temperature, it results in net radiative warming of the planet. Moreover, the study of Hier-Majumder & Hirschmann (2017) preserves the mantle fully molten for the majority of the MO period, due to the lack of a convective cooling sink. The slow solid-matrix compaction process provided from their detailed melt/volatile percolation model further increases the solidification time (3 Myr) in comparison to our study. Lastly, we obtain a lower solidification time in comparison to Hamano et al. (2013), who define the MO end at the surface solidus and not at the higher temperature of the critical melt fraction.

The cooling path can be followed from the convective heat flux, the MO depth, and the Ra number (Figures 3(B)–(D)). For about 50% of its lifetime, the MO has a depth equal to or smaller than 50 km for the Ref-A case. The intersection of the adiabat with the rheology front at the two pressure depths, where switches in the parameterization of the melting curves occur (see Section 2.3 and Appendix A), results in equal characteristic jumps in MO evolution. The decrease in the cooling flux is independent of the decrease in depth D, since D is explicitly overwritten when using the soft-turbulence parameterization (see Equations (5) and (6)). Nevertheless, D defines the average viscosity within the convecting domain, which enters the Ra calculation (Equation (5)). Ra is ultimately responsible for the decrease in heat flux.

The role of radioactive decay as an energy source in the MO evolution of an Earth-sized planet is insignificant (Figure 3(C)), unless the planet is formed within a few Myr, which includes the contribution of the short-lived elements ²⁶Al and ⁶⁰Fe. Their contribution becomes comparable to the long-lived element contribution after 9 Myr and insignificant to the MO evolution by 7 Myr after calcium–aluminium inclusions (CAI) formation, in agreement with the Elkins-Tanton (2012) findings.

Comparing the two atmospheric approaches, we find that the pure water vapor gray approximation underestimates the thermal blanketing in comparison to the lbl model because of the low absorption coefficient used to represent the whole thermal radiation spectra (${k}_{0,{{\rm{H}}}_{2}{\rm{O}}}=0.01$). The lbl pure H₂O model resolves the steam IR absorption better, although it overlooks the role of CO₂.

The "bb–gray–lbl atmosphere" comparison captures the decreasing convective fluxes at the last MO time step as in the "bb–gray–spectrally resolved atmosphere" comparison of the Lebrun et al. (2013) model. In our approach, this is due to the decrease in temperature difference and in Ra toward the MO end. However, upon evaluation at the last time step before solidification, the Ra drop to 10¹⁰ is not seen in their work (where Ra = 10¹⁴ = const), likely due to differing average viscosity calculation and spatial resolution.

On a technical note, a Ra overshoot toward lower values is observed near the end of the MO phase (Figure 3(D)). Consequently, the switch to solid occurs abruptly from Ra = 10¹⁰ to 10¹², two orders of magnitude higher than the value obtained during convection of the last 1 km-deep liquid-like layer of MO. This is a numerical artifact that correlates with the high radial resolution of the model layers (≈1 km). Therefore, care should be taken when using convective heat flux parameterizations with high spatial resolution very close to the critical melt fraction, because the rheology becomes more complex at high crystal values.

The assumption of greenhouse gases H₂O and CO₂ as major species is in accordance with an oxidized MO surface (Hirschmann 2012; Zhang et al. 2017) and bulk silicate Earth (Lupu et al. 2014).

As far as volatile solubility is concerned, molten silicate is a poor CO₂ solvent. It thus operates as a "CO₂ pump" into the atmosphere. In contrast, H₂O is highly soluble in the silicate melt and does not leave the mantle until the latest stage of the MO, where the enrichment in the melt peaks. The evolving atmospheric composition reflects those features as it transitions from a CO₂-dominated to a H₂O-rich one (Figure 4). The major release (from 2.5 to 220 bar) of the water vapor occurs when the total melt fraction of the mantle reduces from 30% to 2% or as the potential temperature drops from ∼2200 to ∼1650 K (Figure 5). This effect is the basis for the so-called "catastrophic" outgassing of a steam atmosphere (e.g., Lammer et al. 2013). It reflects the progressive replacement of the melt volume with solid volume that has a small capacity for storing volatiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, the choice of melting curves defines the degree of melting throughout the MO lifetime and similarly affects the accompanying outgassing process. Comparing chondritic to peridotitic composition for the lower mantle, we find that the melt fraction differs by 10%–43% for the same potential temperature over the range T_p ∈ [3000, 2200] K (Figure 5). The choice of lower mantle melting curves does not affect the final outgassing but modifies the onset of catastrophic outgassing by a maximum 5% of the total volatile volume. Therefore, chondritic composition for the lower mantle disfavors early water release for a cooling MO for potential temperatures above 2200 K.

Simultaneous to the final outgassed quantity, we also calculate the relative volatile inventory extracted from the mantle assuming different initial concentrations (Figure 6). As expected, the higher initial concentration results in higher outgassing. However, the relative quantity varies as follows: we find that [45%, 10%] of the initial water reservoir remains in the mantle for the examined range X ∈ [10⁻⁵, 10⁻¹] respectively, while the rest [55%, 90%] is in the atmosphere. This suggests that the lower the initial mantle abundance, the larger the relative amount of water stored in the planet's interior after the MO ends. By contrast, only ≈6% of CO₂ remains in the mantle for an Earth-sized planet independently of the initial concentration assumed.

Zoom In Zoom Out Reset image size

The combined H₂O/CO₂ inventory was found to delay the MO termination in prior works (e.g., Zahnle et al. 1988; Abe 1997; Lebrun et al. 2013). The Elkins-Tanton (2008) work considers different MO depths (2000, 1000, 500 km) from our global MO for Earth (2890 km). Consequently, the volatile masses differ for the same assumed concentration, and a direct comparison is not possible. Recently, Salvador et al. (2017) studied the effect of water abundances on the global MO solidification time, yielding longer durations likely due to the use of a non-gray atmospheric model. In our study, we quantify the solidification time (t_s) by sampling a larger domain of initial abundances for the two species and assuming a gray atmosphere (Figure 7). The MO duration amounts to ≈0.21 Myr for conservative Earth volatile abundances, while it would reach 5–10 Myr for an (unlikely) Earth-sized planet made entirely out of carbonaceous chondritic (CC) material with 1 wt% of H₂O. Our results confirm that the atmosphere is the most important solidification delaying factor.

Zoom In Zoom Out Reset image size

The effect of each separate interior process on the duration of the MO stage remains difficult to disentangle; however, it would help clarify future modeling priorities. In Table 2, we present an overview of the effect of additional factors and parameters on the MO solidification time (t_s). Each t_s is obtained through varying parameters and/or including a different process (first column). The second column states the number of parameters (three at most) that have been modified in each experiment with respect to a reference case. The third column gives the details on the experiment changes with respect to the reference case. We calculate the solidification time (t_s, fourth column), as well as the absolute (Δt_s, fifth column) and relative difference (Δt_s/t_s,ref, sixth column) with respect to the solidification times (t_s,ref) obtained during two reference simulations: Ref-A (Table 1) and Ref-B. The latter uses the same parameters as the Ref-A settings but does not include CO₂. We thus obtain the tendency of each factor to increase or decrease the solidification time ("+" or "−" sign, respectively) as well as its magnitude. Below we discuss only the most crucial contributions.

Table 2.  Overview of the Effects of Various Parameters on the Solidification Time

Modified parameter	#	Value/Description	ts (yr)	Δts (yr)	Δts/ts,ref
Reference-A	0	Setting described in Table 1	208,600	0	0
1a: H2O content	1	${X}_{{{\rm{H}}}_{2}{\rm{O}},0}$ = 10 ppm	58,900	−149,700	−72%
1b:	1	${X}_{{{\rm{H}}}_{2}{\rm{O}},0}={10}^{5}\,\mathrm{ppm}$	69,699,000	+69,490,400	+33000%
2a: CO2 content	1	${X}_{{\mathrm{CO}}_{2},0}=10\,\mathrm{ppm}$	160,500	−48,100	−23%
2b:	1	${X}_{{\mathrm{CO}}_{2},0}={10}^{5}\,\mathrm{ppm}$	3,919,000	+3,710,400	+18000%
3a: Liquid viscosity	1	${\eta }_{l}=f\left(T,{X}_{{{\rm{H}}}_{2}{\rm{O}}}\right)$	213,400	+4,800	+2%
3b:	2	${\eta }_{l}=f\left(T,{X}_{{{\rm{H}}}_{2}{\rm{O}}}\right)$, ${X}_{{{\rm{H}}}_{2}{\rm{O}},0}={10}^{5}\,\mathrm{ppm}$	69,711,000	+69,502,400	+33000%
4a: Radioactive sources	1	tplanet = 100 Myr	208,600	0	+0%
4b:	2	tplanet = 2 Myr	6,036,780	+5,828,180	+2793%
5: Heat flux parameterization	1	${F}_{\mathrm{hard}};{\tfrac{L}{D}}_{\max }=1$	260,770	+52,170	+25%
6a: Upper-mantle solidus	1	Tsol − 20 K $\Rightarrow$ TRF,0 = 1625 K	216,400	+7,800	+4%
6b:	1	Tsol − 50 K $\Rightarrow$ TRF,0 = 1595 K	228,700	+20,100	+10%
6 c:	1	Tsol − 100 K $\Rightarrow$ TRF,0 = 1545 K	250,600	+42,000	+20%
6d:	1	Tsol − 400 K $\Rightarrow$ TRF,0 = 1245 K	434,600	+226,000	+108%
7: Lower-mantle melting curves	2	Tsol,liq; (Andrault et al. 2011)	207,100	−1,500	−1%
8: Alternative melting curves	2	Tsol,liq; Linear (A) $\Rightarrow$ TRF,0 = 1360 K	126,670	−81,930	−39%
9: Irradiation	1	72%S0	208,500	−100	+0%
10a: Irradiation & albedo	2	S0, α = 0.15	208,600	0	+0%
10b:	2	72%S0, α = 0.60	208,500	−100	+0%
11: No atmosphere and ${\eta }_{l}\left(T\right)$	1	No atmosphere	2000	−206,600	−99%
12a: No atmosphere and ${\eta }_{l}\left(T,{X}_{{{\rm{H}}}_{2}{\rm{O}}}\right)$	2	No atmosphere, ${X}_{{{\rm{H}}}_{2}{\rm{O}},0}$ = 410 ppm	2958	−205,642	−99%
12b: No atmosphere and ${\eta }_{l}\left(T,{X}_{{{\rm{H}}}_{2}{\rm{O}}}\right)$	3	No atmosphere, ${X}_{{{\rm{H}}}_{2}{\rm{O}},0}$ = 104 ppm	2,713	−205,887	−99%
Reference-B	0	Same as ref-A but with ${X}_{{\mathrm{CO}}_{2},0}=0$	156,700	0	0
13: Lbl atmosphere	1	Steam lbl	736,100	+579,400	+278%

Note. Different scenarios are compared to a reference case. The scenarios consist of varying or replacing a parameter or physical process as indicated in the first column. The total number of changed parameters (three at most) with respect to the reference scenario is indicated in the second column. The employed parameter values and/or the description of the process are in the third column. The fourth column shows the solidification time t_s, and the fifth and sixth columns the absolute and relative difference of t_s with respect to reference cases A or B. The Reference case A setting in bold is described in detail in Table 1. The Reference case B in bold is identical to case A except for X_CO2 = 0. The t_s values in bold represent the magma ocean duration (yr) for each case.

Download table as:  ASCIITypeset image

When accounting for the water dependence of the melt viscosity in experiment 3, we expect a shorter solidification time that reflects the more efficient convection due to lower viscosity. η_l decreases due to the progressive enrichment of water concentration in the melt during the MO evolution (from 410 to ≈10⁴ ppm (Ref-A)). The atmospheric radiative forcing remains identical to the Ref-A case. The expected cooling acceleration is counteracted by the delaying role of the outgassed vapor atmosphere (experiment 3), even so for particularly water-rich settings (as seen by the almost identical t_s of the water-rich experiments 1b and 3b). The effect of viscosity on t_s becomes evident in the blackbody cases (experiments 11 and 12). With respect to the blackbody case of experiment 11 (t_s = 2000 yr) that uses a constant 10 wt% water content (Karki & Stixrude 2010), we observe an increase in the solidification time (t_s = 2713 yr) in experiment 12b, which uses water-dependent viscosity. This is explained by the fact that in experiment 12b, the 10% water enrichment occurs only at the latest MO stage and not throughout the whole run. Our parameterizations show that one order of magnitude enrichment in H₂O in the melt causes a decrease of up to two orders of magnitude in the viscosity (Figure 2). This becomes important at lower melting temperatures T_RF,0 < 1400 K, which correspond to evolved silicate melts (Parfitt & Wilson 2008). Experiments 12a and 12b confirm the tendency we hypothesized for the role of viscosity in decreasing t_s with increasing water content (410 ppm and 10⁴ ppm accordingly). Therefore, the water-enriched melt accelerates the solidification process, and it should be taken into account for evolved surface compositions or planets around EUV and XUV active host stars that lose their atmospheres. According to Abe (1997), low viscosity enhances the differentiation of minerals. Therefore, such an η_l parameterization is also vital in better modeling the mineral solidification sequence.

Using the hard turbulence approximation for the convective flux rather than the soft approximation yields a slight increase in the solidification time (experiment 5). The abrupt decrease of ≈1000 K in the surface temperature at the MO termination is reduced by up to 300 K by employing the hard turbulence parameterization. During this, the Pr number is updated according to the evolution of the liquid viscosity, and the flow aspect ratio (λ) takes values between 1 and 2. Significant work that has been done in this direction shows numerical proof of the hard turbulence regime (Grossmann & Lohse 2011) and suggests that it could affect the thermal transport controlled by the boundary layers (Grossmann & Lohse 2003; Lohse & Toschi 2003).

In experiment 6, we examine the role of uncertainty in the upper mantle (0–22.5 GPa) solidus. The ±20 K error estimated in the solidus expression of Herzberg et al. (2000) has a measurable impact (+4%) on the solidification time. The mere uncertainty in the experimental data can thus affect the MO solidification time by a few thousands of years.

Further decreasing the upper mantle solidus by 50, 100, and 400 K causes the solidification time to decrease by 10%, 20%, and 108% respectively. Compositions more silicate-evolved compared to the KLB-1 peridotite have such lower melting temperatures. The −400 K value corresponds to rhyolite (Parfitt & Wilson 2008). Lebrun et al. (2013) and Salvador et al. (2017) previously acknowledged that the chemical composition of the MO at its latest stages would be a decisive factor in the evolution. Schaefer et al. (2016) and Wordsworth et al. (2018) further resolved the chemical evolution for specific compositions. Our result emphasizes the controlling role of the surface melting temperature in the solidification duration and reveals a linear dependence between them.

The solidification time is however insensitive to changes in the lower mantle melting curves (experiment 7) as long as bottom-up solidification is ensured. The reason is that they affect neither the amount of CO₂ in the atmosphere, the majority of which is degassed at the beginning of the MO phase, nor the water enrichment, which does not occur at high MO depths.

In experiment 8, we test the effect of linearizing the melting curves of Abe (1997), where the solidification time decreases significantly (−39%). The higher melt fraction preserved at the end of the MO is tied to lower final outgassing, which explains the difference from the Ref-A setting. Lebrun et al. (2013) previously discussed a similar effect of the curve linearization. A quantitative comparison is however inconclusive, due to the different atmospheres used.

We clarify a fundamental difference between the atmospheric approximations that were implemented in this work. We illustrate this by assuming a high (F_Sun(S = 1361 Wm⁻², α = 0.11) = 303 W m⁻²) and a low (F_Sun(S = 1361 Wm⁻², α = 0.30) = 238 W m⁻²) incoming solar radiation (Figure 8). The difference is only in the assumed albedo value, 0.11 or 0.30.

Zoom In Zoom Out Reset image size

In the lbl approximation (Figures 8(A), (B)), the color map combinations of ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ and T_surf lead to planetary cooling. In the high F_Sun case, for each value of the surface temperature T_surf between 700 and ≈1700 K, there exists a threshold value of outgassed water ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ across which the net radiation balance at TOA is negative, and the planet warms. This effect is absent in the low F_Sun case, which yields a cooling regime for all combinations of ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ and T_surf. On the contrary, the gray approximation shows a negligible difference of the MO cooling flux of the order of 10⁻¹ W m⁻², accounting for the T_eq of our solar system's inner planet orbits (Figure 8(C)). In fact, the gray atmosphere is insensitive to variations in the incoming stellar radiation.

The reason is that in the gray energy balance (Equation (15)), the incoming solar flux enters only in the calculation of the equilibrium temperature. The latter does not vary by more than a factor of 2 over the insolation range in our solar system history (T_eq = 144 K for the case of the young Sun and T_eq = 256 K for today's Sun at 1 au). The fourth power of T_eq has a minor contribution compared to the fourth power of the surface temperature of the MO, which is higher than T_RF,0 = 1645 K (Ref-A) throughout the evolution.

In the limit of our convecting MO model, we only explore cooling regimes and obtain the relevant solidification times. The convective cooling flux out of the MO, F_conv, requires T_surf < T_p to ensure the necessary gravitational instability for convection to occur (see Equation (6)). However, if the flux at the TOA becomes negative (RHS of Equation (16)), the system would warm, resulting in T_surf > T_p, a condition that describes a stably stratified system that will not overturn.

Sections 4.5–4.7 focus on the cooling/warming limit found with the lbl atmosphere.

The lifetime of an MO with a steam atmosphere is controlled by the longwave radiation through its steam layer, the energy received from the star, and the melting temperature of the mantle at its surface. All above factors combine into a comprehensive mechanism that distinguishes between a "transient" (or "short term," or "type I" after Hamano et al. 2013) and "continuous" (or "long term," or "type II" after Hamano et al. 2013) MO evolution path. Goldblatt (2015) and Ikoma et al. (2018) have discussed the warming/cooling distinction, always in relation to the constant radiation limit for RG ≈300 W m⁻². We exemplify this idea with an emphasis on the additional role of T_RF,0.

We use two simulations that are subject to different insolation conditions, namely F_Sun,low = 238 W m⁻² and F_Sun,high = 563 W m⁻² (Figure 9(a) black solid line and black dashed line, respectively), leaving all other parameters unchanged. F_Sun,high is obtained using S = 2648 W m⁻², which corresponds to the incident radiation at the orbital distance of Venus for today's Sun and α = 0.15, while F_Sun,low is equal to the incoming radiation at Earth's orbit today. Because F_Sun is independent of T_surf, it is plotted as a line parallel to the T_surf axis (Figure 9(a)). Both simulations have the same water reservoir (405 bar or 550 ppm initial concentration) to ensure outgassing of one Earth ocean (300 bar) at the end of the MO stage. OLR_TOA as a function of T_surf is plotted for three values of atmospheric water content (4, 100, and 300 bar), which we term "isovolatiles" (gray lines). F_Sun intersects with each isovolatile over a temperature value ${T}_{\mathrm{surf}}^{{\prime} }$. The cooling flux F_conv (read on the right axis) becomes zero for that specific water content, and the planet ceases to cool. If ${T}_{\mathrm{surf}}^{{\prime} }$ is higher than the mantle rheology front temperature at the surface (T_RF,0), the steam quantity indicated by the respective isovolatile balances the energy flux from the star, and the MO does not solidify.

Zoom In Zoom Out Reset image size

First, we examine the trajectory of the convective flux of the transient MO on the (T_surf, F_conv) plane as it cools from T_surf = 3000 K and with an insolation F_Sun,low (Figure 9(a), red solid line). F_conv progressively crosses isovolatiles of higher water content. As it approaches the highest outgassed quantity of 300 bar, the difference OLR_TOA − F_Sun,low = F_conv always remains positive since the 300 bar isovolatile allows the system to dispose of heat at a higher rate than it receives solar radiation. The high convective flux value ensures cooling until T_surf = T_RF,0, which marks the end of the MO. The abrupt cooling after the end of the MO stage and the final outgassing quantity are shown in the evolution of T_surf(t) and ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$(t) (Figures 9(b), (c)).

Second, we obtain a long-term MO (Figure 9(a), red dashed line) in a scenario that assumes F_Sun,high. Initially, for high values of T_surf, the same amount of water as before is outgassed, and its F_conv almost coincides with the one of the short-term case (the difference, ≈10³ W m⁻², is hardly noticeable on the logarithmic graph). During evolution, the outgassing proceeds, and the simulation trajectory crosses isovolatiles of higher water content. F_conv drops to very low values that tend numerically to zero for T_surf = ${T}_{\mathrm{surf}}^{{\prime} }$ ≈ 1915 K. The intersection of the incoming radiation F_Sun,high with the respective isovolatile over ${T}_{\mathrm{surf}}^{{\prime} }$ reflects the steam atmosphere already outgassed when the system ceased to cool. We obtain a point that falls between the isovolatiles of 100 and 300 bar (167 bar read in Figure 9(b)). Consequently, a continuous MO is maintained at potential temperature $\approx {T}_{\mathrm{surf}}^{{\prime} }$ (Figure 9(c)), due to a specific combination of incoming solar radiation, its intersection with the 167 bar isovolatile, and the solidification temperature (Figure 9(a)). Note that the long-term MO ocean is maintained with less water than one Earth ocean and at an insolation higher than the RG limit.

The prominent role of T_RF,0 on the MO type becomes evident when comparing the point (${T}_{\mathrm{surf}}^{{\prime} }$, ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$), where each isovolatile intersects F_Sun, with T_RF,0. For the short-term MO, the intersection point A occurs well below T_RF,0. That MO stage will be transient for every possible outgassing scenario within the [4, 300] bar range. In the case of higher solar irradiation, we have intersection points with each isovolatile (B, C, and D), which indicate different thermal evolution paths. On the one hand, points B and C are located at surface temperatures higher than T_RF,0, which means that if the MO has outgassed the respective quantities of 300 and 100 bar by the time ${T}_{\mathrm{surf}}^{{\prime} }$ is reached, it will cease cooling. On the other hand, point D corresponds to a much lower temperature than T_RF,0, which means that a steam atmosphere of 4 bar under those insolation conditions can counteract the cooling process only if T_surf decreases to 900 K. The respective MO stage is transient, because it solidifies at a much higher temperature (i.e., 1645 K). The variation of the OLR as a function of P, T is explored in detail in the companion work.

Clearly, the essential quantity regarding the planetary heat budget is F_Sun, as it determines the fate of the MO between transient and continuous. Below we refer to this limiting incoming flux as F_lim (where F_conv = 0), and we specify the incident solar radiation and albedo combinations that satisfy it.

On combining the stellar luminosity (L_Star),

$$\begin{eqnarray}&&{L}_{\mathrm{Star}}={\left(4\pi {R}_{\mathrm{Star}}\right)}^{2}\sigma {T}_{\mathrm{eff},\mathrm{Star}}^{4},\end{eqnarray} \tag{ 17 }$$

where R_Star is the stellar radius and T_eff,Star the effective temperature at the star's photosphere; with the expression of Gough (1981) for the evolution of solar luminosity, we get T_eff,Star(t_star). Combining with the blackbody radiation law for the equilibrium temperature of a planet, we obtain the following equation:

$$\begin{eqnarray}&&R=\displaystyle \frac{{R}_{\mathrm{Star}}\cdot {T}_{\mathrm{eff},\mathrm{Star}}^{2}({t}_{{\rm{star}}})\cdot \sqrt{1-{\alpha }_{\max }}}{2\sqrt{\tfrac{{F}_{\mathrm{lim}}({P}_{{{\rm{H}}}_{2}{\rm{O}}},{T}_{\mathrm{RF},0})}{\sigma }}}.\end{eqnarray} \tag{ 18 }$$

Equation (18) relates to F_lim the maximum albedo α_max that an Earth-sized planet at orbital distance R from a star of effective temperature T_eff,Star can possess in order to maintain a continuous MO stage. The limiting flux included in the denominator of Equation (18) is not constant but equal to

$$\begin{eqnarray}&&{F}_{{\rm{lim}}}({P}_{{{\rm{H}}}_{2}{\rm{O}}},{T}_{{\rm{surf}}}^{{\rm{{\prime} }}})=\displaystyle \frac{\left(1-{\alpha }_{c}\right)S{\left({t}_{{star}}\right)}_{1{\rm{au}}}}{4},\end{eqnarray} \tag{ 19 }$$

where α_c is the critical albedo found with a sensitivity experiment for a given planetary volatile inventory. There, S(t_star)_1au is the solar constant at an orbital distance of 1 au and stellar age t_star, ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ (in bar) the mass of the water vapor outgassed, and ${T}_{\mathrm{surf}}^{{\prime} }$ the temperature over which the stellar insolation crosses the ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ isovolatile (Section 4.5). The obtained limiting flux F_lim maps to the data product ${{\rm{OLR}}}_{\mathrm{TOA}}({T}_{\mathrm{surf}}^{{\prime} },{P}_{{{\rm{H}}}_{2}{\rm{O}}})$. Using the Katyal et al. (2019) values of the limiting flux, we compared our Equation (18) with an equivalent expression calculated by Hamano et al. (2013). The solution is similar with minor differences due to the astrophysical properties assumed. We generalize the formulation to cover our Sun or any other host star with a known photospheric temperature T_eff and radius.

The irradiation conditions that can pinpoint to an Earth-sized planet stalling in a MO stage just above the 40% melting temperature are extended here to include a range of steam atmosphere masses that span [4, 300] bar for two different T_RF,0 values (the lower melting temperature is representative of a more evolved composition than the KLB-1 peridotite of Ref-A).

Different radiative flux limits F_lim are obtained for the two values T_RF,0 = 1370 and 1645 K, depending on the vapor amount (Figure 10). From the superposition of points that correspond to 200 and 300 bar at T_RF,0 = 1645 K in Figure 10, we note the tendency of steam atmospheres exceeding 200 bar to converge to the constant RG limit (RG = 282 W m⁻²; Katyal et al. 2019). For T_RF,0 = 1370 K, atmospheres already equal to or higher than 100 bar suffice to reach the RG limit. A similar tendency is shown in Hamano et al. (2013).

Zoom In Zoom Out Reset image size

We also find that at lower steam contents, F_lim is greater than the RG limit. All F_lim values can be found in Table 3. Using Equation (18), we calculate the orbital distance–albedo combinations for which the radiation limits (Table 3) of different isovolatiles are attained (Figure 11). We assume solar luminosity at the beginning of its main-sequence evolution at τ = 100 Myr (72% of today's value; Gough 1981). Not all of the calculated albedo values are realistic. The albedo for a cloudless steam atmosphere, based on 1D models and 3D global circulation model calculations, lies between [0.15, 0.40] (Kasting 1988; Goldblatt et al. 2013; Leconte et al. 2013; Pluriel et al. 2019).

Zoom In Zoom Out Reset image size

Table 3.  ${F}_{\mathrm{lim}}({P}_{{{\rm{H}}}_{2}{\rm{O}}},{T}_{\mathrm{RF},0})$ for Indicative T_RF,0 Cases Calculated with the Katyal et al. (2019) Data

${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ (bar)	Flim(TRF,0 = 1645 K) (W m−2)	Flim(TRF,0 = 1370 K) (W m−2)
4	16724.2	5885.6
25	1362.8	519.4
50	677.3	332.8
100	398.5	286.6
200	298.1	283.6
300	285.8	283.6

Download table as:  ASCIITypeset image

Apart from the new radiation limits found, our results are in line with those of previous studies, as far as the role of insolation is concerned. In Hamano et al. (2013), the threshold distance between continuous and transient MO types for albedo 0.3 and solar constant 0.72S₀ is 0.77 au, whereas under the same conditions, our calculations show 0.79 au. The difference is due to the lower absolute OLR steam atmosphere limit of 282 ± 1 W m⁻² that we obtain compared to the 294 W m⁻² limit employed in that study.

By raising the albedo to the critical value α_c = 0.146 found in our simulations for an Earth-sized planet at 1 au that outgasses one Earth ocean at today's Sun, we obtain ≈12 Myr MO duration. The longer solidification times of 10–30 Myr reported in Hamano et al. (2013) at the same F_Sun = 285.5 W m⁻² for the same steam atmosphere are due to the lower surface melting temperature used (T_RF,0 = 1370 K).

With an albedo of 0.63, a solar constant of 0.7S₀, and total planetary water content ${{\rm{X}}}_{{{\rm{H}}}_{2}{\rm{O}},0}$ = 5.53 · 10⁻² wt% (equivalent to 405 bar) and ${{\rm{X}}}_{{\mathrm{CO}}_{2},0}$ = 1.4 · 10⁻² wt% (equivalent to 100 bar), Lebrun et al. (2013) found that the distance at which the outgassed water vapor could no longer condense is 0.67 au. Under the same conditions, excluding the influence of CO₂, we find in our model that the atmosphere would exist in an RG state trapped in a continuous MO at a critical distance of 0.59 au. The reason for this discrepancy is twofold. First, our lbl model approach does not include CO₂, which also contributes to the greenhouse effect. Second, the absolute OLR limit used by Lebrun et al. (2013) is ≈200 W m⁻² (Marcq 2012; see Marcq et al. 2017 for an updated limit). This is substantially lower than the limit of 282 W m⁻² used in our study. Therefore, the shift of our limit inward toward the star corresponds to the higher critical flux that needs to be achieved in order to trigger the qualitative shift from a transient to a continuous MO regime. Considering the orbital distances of the inner terrestrial planets of the solar system (Mercury—0.38 au, Venus—0.72 au, Earth—1 au, and Mars—1.52 au), we find that the planets inwards of Earth could sustain a continuous MO within the range of albedos expected for a cloudless steam atmosphere (Figure 11). Moreover, a 100 Myr old Earth at 1 au around the Sun cannot exist in a continuous MO state under any albedo for a steam atmosphere of up to 300 bar (Figure 11) or of up to 1000 bar according to the recent study of Ikoma et al. (2018).

Note that in this work "continuous" MO refers to planets that would cease cooling if the amount of steam in the atmosphere was conserved. This cannot be ensured under atmospheric escape processes, which have not been accounted for, and as such the limits calculated here yield the farthest possible distance from the Sun for achieving a continuous MO with a constant atmospheric steam content. Using the database of F_lim that depends on both the atmospheric water content and T_RF,0 that we provide in this (Table 3) and in the companion work, Equation (18) is qualitatively extended. It covers MO-type transitions for intermediate levels of outgassing below the 300 bar reference value, and hence has higher F_lim. This database is backwards compatible and can also be used in the Hamano et al. (2013) equation.

We draw specific attention to the difference between "evolutionary" and "permanent" MO, which is studied in other works (e.g., Hammond & Pierrehumbert 2017). Both the transient and continuous MO that we study belong to the so-called "evolutionary" MOs generated during the accretional process. However, within a certain orbital distance, the energy for melting the mantle is already provided by the solar irradiation alone, and the atmosphere blanketing effect becomes irrelevant. This is the "permanent" MO caused by the star. Radioactivity and delivery of kinetic energy cause the evolutionary MO. To help distinguish these, the minimum distance from the star for which the equilibrium temperature T_eq equals the 40% melt fraction surface temperature (where permanent MO stage ensues) is indicated by the black lines in Figure 11. Further research, which is beyond the scope of this study, could lead to an expansion of the orbital distance defining the permanent MO, upon accounting for climate feedbacks that raise the surface temperature above the T_RF,0 melting point.

F_lim is not significantly affected by the gravitational acceleration of the planet as long as this has between 0.1 and 2 Earth masses (Goldblatt et al. 2013). For greater planetary masses, the pressure levels in the atmosphere change height, as does the level of opacity depth, which is crucial for the calculations of the outgoing radiation. Goldblatt et al. (2013) also calculated that for a planet of half the Earth mass, the OLR-limiting flux is lowered by only 5 W m⁻². In comparison, F_lim ≈ 282 W m⁻² ± 1 W m⁻² as calculated for Earth by Katyal et al. (2019) has a lower uncertainty. Therefore, Equation (18) can be applied without loss of generality to planets between 0.1 and 2 Earth masses, using the ${F}_{\mathrm{lim}}({P}_{{{\rm{H}}}_{2}{\rm{O}}},{T}_{\mathrm{RF},0})$ (Table 3) calculations by Katyal et al. (2019).

Given their similarity in mass and radius, the criteria for a continuous MO applied for Earth can be extended to Venus. A continuous MO could not have been possible for Earth during the young Sun period for any bulk water abundance. We find however, that the Venus orbit qualifies for a long-term MO within a wide range of planetary albedos [0.15, 0.40] proposed for cloud-free steam atmospheres, as long as its outgassed steam atmosphere amounts to 200 bar or more for a surface solidification temperature of 1645 K (Figure 11). In the case of the lowest solidification temperature (T_RF,0 = 1370 K), the minimum atmosphere required for a continuous MO at Venus orbit is 50 bar (Figure 11). This highlights the fact that the melt composition alone could dictate a different MO evolution path for two hypothetical planets with equal water vapor atmosphere masses.

It is additionally important to consider whether a planetary body has had a long impact history or has chemically evolved before impacts remelt it into an MO (Lammer et al. 2018). Such bodies could more easily maintain a secondary continuous MO. Due to their lower T_RF,0, they would require smaller steam atmospheric mass, instead of the reference one Earth ocean (300 bar) usually assumed in RG studies. On the contrary, a chemically unevolved silicate primitive composition that melts at high temperatures would require a massive steam atmosphere >100 bar in order to maintain a continuous MO. We conclude that past events of chemical alteration may influence the fate of the MO under the same orbital configuration. Therefore, the age of the star and of its planetary system matters. Evolution of the mantle composition during the MO solidification (Elkins-Tanton 2008; Schaefer & Fegley 2010) will be an additional factor that prolongs the MO lifetime if it results in decreased T_RF,0.

Interestingly, the drop in surface temperature during cooling combined with the tendency of stellar environments to gradually strip planets of their atmospheres (Johnstone et al. 2015; Odert et al. 2017; Lammer et al. 2018; therefore lowering the surface pressure ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$) could result in the same outgoing radiation limit during planetary evolution. We see this in Figure 10 taking any constant OLR value that crosses multiple isovolatiles. In that case, the stellar evolution plays a primary role in the fate of the continuous MO. A G-star with luminosity increasing with time (Gough 1981) favors the maintenance of an existing MO because it contributes warming at the critical distance. In contrast, continuous MOs will be more elusive around M-stars, the luminosity of which decreases with time (Baraffe et al. 2015). This is because a continuous MO close to its critical distance will receive less and less stellar radiation, eventually creating a window of cooling. A buffer against this effect is the additional vapor outgassing that increases the opacity and lowers the required F_lim. However, during progressive cooling, the interior will exhaust its water supply into the atmosphere. Under these conditions, only water-rich planets can sustain a continuous MO. This shows that there are numerous processes that affect MO feasibility. Consider also the possible Trappist-1 exoplanet migration scenarios (Unterborn et al. 2018; suggested in order to justify water-rich composition).

We explore our findings in view of potentially rocky exoplanets having radius and/or mass within a few Earth units (parameters in Appendix E). Results suggest that there are orbital regions where the MO can be transient or permanent, and an intermediate region where it is "conditionally" continuous (Figure 12). "Conditionally" here refers to the dependence on water content and rheology front temperature. We observe the overlap in the regions of continuous MO for different T_RF,0. Considering the interior composition adds a measurable level of uncertainty since different planets with different atmospheric water content and solidification temperatures can be characterized by the same outgoing OLR. Note that unless we are able to constrain the surface pressure of water vapor on an exoplanet, which is not feasible with the current observational capabilities (Madhusudhan et al. 2014), we will not be able to constrain the type of evolutionary (continuous, transient) MO. However, hypotheses and proxies concerning planetary water abundance could break this OLR degeneracy that disappears at low vapor pressures close to 4 bar (see Figure 10). A water-poor planet with a thin atmosphere of 4 bar water would be sensitive to the T_RF,0 value for developing a continuous/transient MO close to their separation limit. Such could be the case for distinguishing the compositions of HD 219134 b and c if one is found in MO state and the other is not (Figure 12). Kepler 36 b' s orbit is farther than this distinction possibility and receives enough energy from the star to be in continuous MO as long as it has at least 4 bar water. As soon as its atmosphere is lost, it would resemble a blackbody on which the liquid viscosity effect of any water present would ensure rapid MO solidification, as we showed earlier. Planets Kepler 236c, Ross 128 b, and LHS 1140c, on the contrary, are located in the F_lim region for relatively high vapor pressure. Assuming ≥200 bar, the system converges to the minimum OLR solution of 282 W m⁻² (see Figure 10), which is maintained for up to 1000 bar (Ikoma et al. 2018). Detecting any MO state on those planets would be difficult because of the opaque atmosphere. However, if detected, it would mean that the planet formed within a water-rich environment that ensured the minimum atmospheric 200 bar required for the continuous MO. Especially for LHS 1140c, the planet LHS 1140 b located in the transient MO region of the same system could provide complementary information for the likelihood of high water content. GJ 1132 b is located at the compositional distinction limit. Its potential MO has been studied before by Schaefer et al. (2016). A low atmospheric water content in its MO state would be a proxy of primitive silicate composition. Any of the continuous MOs on those planets would eventually solidify if their atmospheric water were lost and were not replenished by the interior.

Zoom In Zoom Out Reset image size

The possibility of observing a transient MO system is insignificant, due to the order of million years duration that we find for them, which is very short compared to the ages of observable systems. Detection of continuous MOs on candidate planets (at orbits receiving 282 W m⁻² or more; see Figure 12) is challenging but is aided by the fact that the planet's MO brightness temperature would be much higher than that corresponding to its equilibrium temperature, yielding OLR of up to 16,000 W m⁻² (see Table 3). Such measurements require secondary transit observations as carried out for 55 Cnc e with the Spitzer telescope (Demory et al. 2016) aided by the longer wavelength coverage of JWST. A low brightness temperature, in agreement with a low OLR of 282 W m⁻², would be an indication toward high steam pressures (see companion paper for possible emission spectra). The surface pressure is not retrievable with the current capabilities, but promising methods are being developed for low pressure atmospheres (10 bar) that demonstrate pressure broadening of absorbers such as CO₂ and O₂ (Misra et al. 2014). Transmission methods could not probe high surface pressure atmospheres but the OLR would be already near the RG limit in those cases so one should focus on retrieving the latter. A measured OLR = 282 W m⁻² would be indicative of MOs with high steam pressures. We suggest the auxiliary/complementary use of observations obtained from the permanent MO type, such as potentially on 55 Cancri e (Demory et al. 2016; Angelo & Hu 2017) and Kepler 78 b. From there one could isolate characteristic atmospheric signatures such as the atmospheric effects of evaporated silicate species that develop over the molten rocky surface (Fegley et al. 2016; Kite et al. 2016; Hammond & Pierrehumbert 2017) and the oxides in the presence of a steam atmosphere (Fegley et al. 2016). Detecting similar silicate cloud signatures on planets close to the continuous MO compositional distinction that is observed at low vapor pressures (4 bar) would serve as a proxy of their composition (T_RF,0) and of their water content.

Detection of evolutionary MOs additionally requires stellar ages in order to focus on systems with ongoing planetary formation, preferably after recently completed accretion. Constraining the albedo from observations is a possibility given favorable orbital configurations (Madhusudhan et al. 2014; Kite et al. 2016) and would help define the range of orbital distances for a conditionally continuous MO.

For the solidus temperature (T_sol) of nominally anhydrous peridotite and pressures 0 ≤ P ≤ 2.7 GPa, we use (Hirschmann 2000)

$$\begin{eqnarray}&&{T}_{\mathrm{sol}}=1120.661+273.15+132.899P-5.904{P}^{2},\end{eqnarray} \tag{ 20 }$$

with reported error by the authors of ±20 K. For 2.7 < P ≤ 22.5 (Herzberg et al. 2000),

$$\begin{eqnarray}&&{T}_{\mathrm{sol}}=1086+273.15-5.7P+390\mathrm{log}(P),\end{eqnarray} \tag{ 21 }$$

with reported error by the authors of ±68 K. At lower mantle pressures, for P > 22.5 GPa, we use a quadratic fit to the data of Fiquet et al. (2010) for fertile peridotite:

$$\begin{eqnarray}&&{T}_{\mathrm{sol}}=1762.722+31.595P-0.102{P}^{2}.\end{eqnarray} \tag{ 22 }$$

For the liquidus of fertile peridotite, we use a fit to the data of Zhang & Herzberg (1994) for 0 ≤ P ≤ 22.5 GPa,

$$\begin{eqnarray}&&{T}_{\mathrm{liq}}=2014.497+37.743P-0.472{P}^{2},\end{eqnarray} \tag{ 23 }$$

and for 22.5 > P GPa, again a quadratic fitting to the data of Fiquet et al. (2010),

$$\begin{eqnarray}&&{T}_{\mathrm{liq}}=1803.547+50.810P-0.185{P}^{2}.\end{eqnarray} \tag{ 24 }$$

A linear approximation for the melting curves presented in Abe (1997) is used for both solidus and liquidus:

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{sol}} & = & 1500.0+\displaystyle \frac{4450.0}{3000.0}z,\\ {T}_{\mathrm{liq}} & = & 2000.0+\displaystyle \frac{4560.0}{3000.0}z,\end{array}\end{eqnarray} \tag{ 25 }$$

where z is the depth from the surface in kilometers.

The following quadratic fitting to the data of Andrault et al. (2011) for a chondritic composition is employed for the lower mantle only for P > 22.5 GPa:

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{sol}} & = & 2056.489+15.801P-0.003{P}^{2},\\ {T}_{\mathrm{liq}} & = & 2049.555+24.671P-0.035{P}^{2}.\end{array}\end{eqnarray} \tag{ 26 }$$

Melting curves for the upper mantle are identical to those described above as "synthetic," unless otherwise specified.

In Figure 14, we show the OLR on the 50 × 8 grid of $\left({T}_{\mathrm{surf}},{P}_{{{\rm{H}}}_{2}{\rm{O}}}\right)$ points that we used as input for our simulations. The OLR data at each grid point have been obtained with the method described in the companion paper using an lbl code (GARLIC) of Schreier et al. (2014).

Zoom In Zoom Out Reset image size

The grid spans surface temperatures from 650 to 4000 K and water vapor surface pressures from 4 to 300 bar. It is irregularly spaced and denser over the temperature range where the highest rate of enrichment and outgassing takes place. This range is obtained by performing simulations using synthetic melting curves with the interior model coupled to the gray atmosphere model (see Figure 5). For ${T}_{s}\in \left[1400,1800\right]$ K, the OLR was sampled with a resolution of 20 K, while a 100 K resolution was employed outside this range. The sampling is sparser (eight values) on the pressure axis for $P\in \left[\mathrm{4,300}\right]$ bar. In order to obtain the OLR values at intermediate (P, T) points of the above data set, a bilinear interpolation method was used (van Rossum & Drake 2001). In order to estimate the interpolation error, we compared the interpolated field with an independent set of intermediate data points obtained from the atmospheric model. The relative interpolation error amounts to between 1% and 10%. The maximum of 10% occurs at pressures lower than 10 bar and high temperatures. The minimum occurs for high pressures and temperatures at the lower end of the data set. Therefore, the quality of the result is acceptable for this study that focuses on the coolest end of the MO phase where the outgassed atmosphere has high pressure, and the errors are minimal (1–2 W m⁻²).

The data of Figure 14 represent the OLR at the TOA with a viewing angle of 38° and thus differ from the field shown in Figure 8, which represents the net planetary flux at the TOA. More details on the OLR value calculation can be found in the companion paper.

In order to satisfy the requirement of our iteration algorithm for surface temperatures lower than ${T}_{{{\rm{H}}}^{2}{\rm{O}},\mathrm{crit}}$ = 647 K, which are not covered by our grid, we use a fit to the OLR data of (Nakajima et al. 1992). This aspect does not affect our results for the solidification process, which occurs for T_surf ≈ T_RF,0 ≫ ${T}_{{\rm{H}}2{\rm{O}},\mathrm{crit}}$, but ensures that the iteration algorithm runs unhindered until convergence to the solution.

We use Equation (18) to estimate the orbital distance for which a planet of a given albedo is located at the boundary that separates a long-term and short-term MO. To this end, the value of the limiting radiation F_lim corresponding to a specific water vapor pressure ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ and rheology front temperature at the surface T_RF,0 is needed. In Table 3, we report ${F}_{\mathrm{lim}}({P}_{{{\rm{H}}}_{2}{\rm{O}}},{T}_{\mathrm{surf}})$ for two different rheology front surface temperatures as obtained by interpolating the OLR data points of the companion paper (Figure 14). The same values are plotted in Figure 10 and used to calculate the critical distances for the young Sun in Figure 11.