# Retention of Water in Terrestrial Magma Oceans and Carbon-rich Early
Atmospheres

Dan J. Bower Center for Space and Habitability, University of Bern
Gesellschaftsstrasse 6
3012 Bern, Switzerland  Kaustubh Hakim Center for Space and Habitability,
University of Bern
Gesellschaftsstrasse 6
3012 Bern, Switzerland  Paolo A. Sossi Institute of Geochemistry and
Petrology, Department of Earth Sciences, ETH Zurich
Clausiusstrasse 25
8092 Zurich, Switzerland  Patrick Sanan Institute of Geophysics, ETH Zurich
Sonneggstrasse 5
8092 Zurich, Switzerland

(Received July 14, 2021; Revised October 12, 2021; Accepted March 21, 2022)

###### Abstract

Massive steam and CO2 atmospheres have been proposed for magma ocean
outgassing of Earth and terrestrial planets. Yet formation of such atmospheres
depends on volatile exchange with the molten interior, governed by volatile
solubilities and redox reactions. We determine the evolution of magma
ocean–atmosphere systems for a range of oxygen fugacities, C/H ratios and
hydrogen budgets that include redox reactions for hydrogen (H2–H2O), carbon
(CO–CO2), methane (CH4), and solubility laws for H2O and CO2. We find that
small initial budgets of hydrogen, high C/H ratios, and oxidizing conditions,
suppress outgassing of hydrogen until the late stage of magma ocean
crystallization. Hence early atmospheres in equilibrium with magma oceans are
dominantly carbon-rich, and specifically CO-rich except at the most oxidizing
conditions. The high solubility of H2O limits its outgassing to melt fractions
below $\sim$30%, the fraction at which the mantle transitions from vigorous to
sluggish convection with melt percolation. Sluggish melt percolation could
enable a surface lid to form, trapping water in the interior and thereby
maintaining a carbon-rich atmosphere (equilibrium crystallization).
Alternatively, efficient crystal settling could maintain a molten surface,
promoting a transition to a water-rich atmosphere (fractional
crystallization). However, additional processes, including melt trapping and H
dissolution in crystallizing minerals, further conspire to limit the extent of
H outgassing, even for fractional crystallization. Hence, much of the water
delivered to planets during their accretion can be safely harbored in their
interiors during the magma ocean stage, particularly at oxidizing conditions.

Planetary interior — Planetary atmospheres — Planetary structure — Super
Earths — Extrasolar rocky planets — Exoplanet atmospheric composition —
Exoplanet evolution

††journal: The Planetary Science Journal††software: Simulating Planetary
Interior Dynamics with Extreme Rheology (SPIDER) (Bower et al., 2018, 2019) is
open source and hosted at https://github.com/djbower/spider. SPIDER Version
0.2.1 was used in this study and is also available on Zenodo (Bower et al.,
2021).

##  1 Introduction

The origin and distribution of light volatile elements such as C, H, and O, in
the major planetary reservoirs (core, mantle, and atmosphere), regulates how
terrestrial planet atmospheres form and evolve (Gaillard et al., 2021, and
references therein). During planetary accretion and contemporaneous iron-core
formation—both of which promote the formation of a global magma ocean—the
siderophile nature of both C (Corgne et al., 2008) and H (Tagawa et al., 2021)
strongly partitions them into the core. Nevertheless, the mantle retains some
fraction of C at core-forming conditions, in addition to replenishment of both
C and H after core formation (Hirschmann, 2016); hence, these volatiles
subsequently participate in short- and long-term geochemical cycling (e.g.,
Dasgupta, 2013). Here we examine how the inventories of C and H exchange
between the interior and atmosphere as the magma ocean evolves after core
formation.

The magma ocean outgasses volatiles as it cools from molten to mostly solid
silicate rock to form a secondary atmosphere. Therefore, chemical and physical
processes operating in the magma ocean, both at its surface and at depth,
exert strong controls on the size and composition of such early atmospheres.
Early atmospheres establish the initial conditions for the subsequent long-
term evolution of planetary interiors and their atmospheres through
interaction with the surface environment via geochemical cycles. Therefore,
magma ocean atmospheres provide the connection between the volatile budgets
emplaced during planet formation and the establishment of long-term climate
states.

The shock degassing of substantial H2O from hydrated minerals (Lange & Ahrens,
1982) to simulate impacts during planetary accretion motivated investigation
of the blanketing effect of a steam atmosphere above the molten early Earth
(Abe & Matsui, 1985). Additional justification for the widespread presence of
H2O and CO2 is provided by experiments and accompanying equilibrium models for
the devolatilization of ordinary and carbonaceous chondrites, which are shown
to produce oxidized species in abundance (Gooding & Muenow, 1977; Schaefer &
Fegley, 2010; Thompson et al., 2021). Consequently, coupled magma
ocean–atmosphere models have focused almost exclusively on the role of
oxidized species (H2O and CO2) in controlling the lifetime of magma oceans and
growth of early atmospheres (e.g., Lebrun et al., 2013; Hier-Majumder &
Hirschmann, 2017; Salvador et al., 2017; Nikolaou et al., 2019; Barth et al.,
2021).

More reduced materials (such as enstatite chondrites that are thought to have
contributed to Earth’s accretion on isotopic grounds) may instead produce
H2-rich or CO-rich gases (Schaefer & Fegley, 2010). This has partly inspired
calculations of the chemistry and cooling timescales of reducing atmospheres
(Zahnle et al., 2020; Lichtenberg et al., 2021). However, chondrites may not
be appropriate analogs for the materials that degas to produce early
atmospheres around terrestrial planets (Sossi, 2021). This is because mass
transfer is much more efficient when materials are liquid rather than solid,
allowing for the more rapid replenishment of surface material via convection
and thus facilitating atmosphere formation. This has further led coupled
models to understand the astrophysical observables of magma oceans to
interpret the spectroscopic data from next-generation telescopes (Hamano et
al., 2015; Bonati et al., 2019; Bower et al., 2019; Katyal et al., 2020).

The redox state (more precisely, the oxygen fugacity, $f{\rm O_{2}}$) of the
mantle depends on pressure, temperature, and composition and dictates the
speciation of outgassing products. It controls the relative fugacities of
reduced and oxidized species, which is equivalent to partial pressure at the
low total pressures of terrestrial atmospheres. Core formation on Earth
established $f{\rm O_{2}}$ around two log10 units below the iron-wüstite (IW)
buffer at equilibrium at depth, based on the Fe content of the mantle and core
(e.g., Rubie et al., 2015). However, the volume change of the reaction between
ferric and ferrous iron means that silicate melt transported to the surface
defines higher $f{\rm O_{2}}$ than that set by core formation (Armstrong et
al., 2019; Hirschmann, 2012). The surface $f{\rm O_{2}}$, rather than the
$f{\rm O_{2}}$ at depth, imposes the redox state at the magma ocean–atmosphere
interface and thereby controls outgassing chemistry (Sossi et al., 2020).

Evolutionary models of magma oceans are usually derived from geochemical or
geodynamic considerations. Chemical models focus on tracking the compositional
evolution of cumulate mineral assemblages that form as the magma ocean cools
(e.g., Elkins-Tanton, 2008). However, they ignore dynamics during the
crystallization process and often assume fractional crystallization. In
(stepwise) fractional crystallization, cumulates that form owing to cooling
are isolated from the molten magma ocean, and the evolving composition of the
melt is determined from mass balance. By contrast, dynamic models track
thermal energy transport but at the expense of reduced compositional
complexity. However, dynamic formulations that include a local representation
of melt–solid separation (Abe, 1993; Bower et al., 2018) can replicate
fractional crystallization. They can also replicate equilibrium
crystallization, in which melt and solid freeze together before any
significant relative motion has occurred.

Equilibrium chemistry calculations propose a $\sim$100 bar CO-rich atmosphere
for an Earth-like magma ocean around 2200 K (Sossi et al., 2020), although it
remains unclear how the atmospheric size and speciation change during the
cooling and crystallization of the magma ocean. In particular, it is unknown
whether all early atmospheres are expected to be CO-rich and whether they can
transition to atmospheres that are instead dominated by hydrogen species, such
as H2, H2O, and CH4—with implications for habitability and the formation of
surface water oceans. Furthermore, it remains to be determined how the
transition is influenced by planetary conditions (e.g., $f{\rm O_{2}}$, C/H
ratio, H budget), which is pertinent given the diversity of formation
environments for terrestrial planets outside the solar system.

To this end, we combine mass balance and equilibrium chemistry in a self-
consistent and time-evolved model to probe magma ocean–atmosphere exchange.
Unlike previous models, our model accounts for the CO2/CO, H2O/H2, and CH4/CO2
ratio set by the oxygen fugacity, C/H ratio and temperature of the magma
ocean, as well as their relative solubilities. We explore a range of initial
endowments of the hydrogen and carbon volatile inventory and quantify the
relationship between these variables and the partial pressures of the degassed
species. Finally, atmospheric escape is also modeled to determine the extent
to which the preferential loss of hydrogen can influence the subsequent
outgassing of dissolved volatiles.

##  2 Interior–Atmosphere Coupling

###  2.1 Overview

The SPIDER code (Bower et al., 2018, 2019) solves for the coupled evolution of
the silicate mantle and atmosphere, where the mantle can be molten, solid, or
a mixture (partial melt). It considers interior energy transport by
convection, conduction, and relative motion of melt and solid through mixing
and separation and is a true 1D model in the sense that energy fluxes are
determined locally (Abe, 1993; Bower et al., 2018). Cooling of the interior is
regulated by radiative transfer in the atmosphere, and the atmosphere itself
grows through outgassing of the interior as the mantle solidifies owing to
cooling. We consider a two-phase system of single composition (MgSiO3) for
simplicity, where the thermophysical properties of solid and melt are
determined from Mosenfelder et al. (2009) and Wolf & Bower (2018),
respectively. The melting curves adhere to peridotite melting data in the
upper mantle and measurements on chondritic material in the lower mantle
(Andrault et al., 2011); they are not perturbed according to the mantle
volatile content. The atmosphere is treated as gray with no scattering, and
the two-stream approximation is applied to solve for radiative transfer (Abe &
Matsui, 1985); Radiation limits are not considered. Opacities of the gas
species are provided in Appendix A.

We modified the SPIDER code to additionally account for mass exchange between
volatiles according to equilibrium chemistry (Section 2.3), as well as
atmospheric escape (Section 2.4). The mathematical description of the volatile
mass balance and its adherence to equilibrium chemistry is given in detail in
Appendix C. For each chemical reaction, we introduce a term into the volatile
mass balance to account for the exchange of mass necessary to retain chemical
equilibrium between participating gas species as dictated by oxygen fugacity
and temperature. Similarly, another term in the mass balance equation tracks
volatile loss due to atmospheric escape. Hence, chemical equilibrium and
escape are self-consistently determined as part of the same system of
equations that are integrated to describe the coupled interior–atmosphere
evolution.

Our parameter choices are similar to those in Bower et al. (2018, 2019), so we
only present pertinent details here. We consider a planet with Earth
dimensions and a magma ocean that crystallizes from the bottom up, an approach
that is justified owing to the steeper melting curves compared to the mantle
adiabat. Convection and melt–solid mixing are determined by an eddy
diffusivity that is based on a constant mixing length, and internal heat
sources are not considered because their influence is negligible over the
lifetime of a magma ocean. Viscosity controls the dynamic regime of the mantle
and varies from $10^{2}$ Pa s for pure silicate melt to $10^{22}$ Pa s for
pure solid, increasingly abruptly at a melt fraction of 40% (Costa et al.,
2009; Bower et al., 2018). Other controls on viscosity, such as temperature
and composition, are second order and therefore neglected. The crystal size,
which controls the efficiency of melt–solid separation, is set at either 1 mm
or 5 cm (Appendix B) to simulate equilibrium or fractional crystallization,
respectively.

We track the reservoir evolution of five volatile species during the magma
ocean stage: H2, H2O, CO, CO2, and CH4. For convenience, we refer to H2O and
H2 as hydrogen species and CO, CO2, and CH4 as carbon species. The dissolved
content of H2O and CO2 in the magma ocean is determined by their solubility,
and reduced species are set to zero solubility (Section 2.2). Redox reactions
impose H2–H2O, CO$-$CO2, and CO2–H2–CH4 equilibrium (Section 2.3).
Partitioning of volatiles into solids (crystals) is not considered because its
effect is small compared to the partitioning of volatiles between the melt and
atmosphere. Therefore, the amount of volatiles stored in the mantle during
solidification is a minimum estimate. We simulate cooling of a magma ocean
from an initial surface temperature of 2700 K until it reaches 1650 K. For
surface temperatures cooler than 1650 K the assumption of equilibrium
dissolution begins to break down. This is because a rheological transition
caused by increasing crystal fraction results in mass transfer that is too
sluggish to maintain equilibrium with the atmosphere; for conciseness we
subsequently refer to this event as ”surface lid formation.”

###  2.2 Solubility

A solubility law relates the dissolved volatile content of a particular
species in melt $X$ to its fugacity $f$. Pressures at the magma
ocean–atmosphere interface of an Earth-sized planet are expected to be around
several hundred bars, which motivates our subsequent choice of solubility
laws. In the following, the dissolved content of a given volatile in the magma
ocean is set solely by its fugacity at the interface and the relevant
solubility law. The volatile species are considered to be homogeneously
distributed in the magma ocean.

####  2.2.1 Solubility of H2O

We employ a general expression for solubility with a power-law dependence:

$X=\alpha f^{1/\beta},$ (1)

where $\alpha$ and $\beta$ are constants specific to a particular equilibrium
between a gas species and its dissolved component. In Equation 1, $\alpha$ is
an empirically determined coefficient that encompasses the temperature,
pressure, and composition of the liquid in which the solubility is being
determined, while $\beta$ reflects the stoichiometric relationship between the
mole fraction of the dissolved species and the fugacity of the gas species.
Henry’s law is recovered when $\beta=1$, which indicates no chemical reaction
(hence no change in speciation) when the gas dissolves in the solvent.
Throughout this work we assume gas ideality such that fugacity and partial
pressure are equivalent, and hence symbols for fugacity $f$ and partial
pressure $p$ are used interchangeably. In Equation 1, the relevant partial
pressure (fugacity) is determined at the interface of the magma ocean and
atmosphere.

At low H2O fugacities, water dissolves in silicate melts as hydroxyl (OH-)
group, so $\beta_{\rm H_{2}O}=2$ describes the solubility of H2O in melt
(e.g., McMillan, 1994). However, this relationship breaks down at high H2O
fugacities because H2O replaces OH- as the predominant melt species (above
$\sim 2000$ bars, Berndt et al., 2002; Stolper, 1982). Therefore, a single
$\beta_{\rm H_{2}O}$ determined by simply best-fitting data across large
$f$H2O ranges yields spurious values of $\beta_{\rm H_{2}O}$ because different
reactions govern solubility at different $f$H2O.

Table 1 summarizes recent experimental constraints on H2O solubility for
compositions, pressures, and temperatures relevant to the magma
ocean–atmosphere interface of a terrestrial planet, i.e. a few hundred bars
and greater than about 1500 K. The table is ordered based on the silica
content of the sample, increasing from less than 45 wt% SiO2 to more than 50
wt% SiO2. Experiments on evolved compositions (e.g. rhyolitic/granitic melts)
are not included because they are not representative of compositions during
the magma ocean stage: these are instead summarized in Iacono-Marziano et al.
(2012).

Peridotite is a rock composed predominantly of olivine, which is the major
constituent of Earth’s mantle and likely exoplanetary mantles (Putirka &
Rarick, 2019). The solubility of water in liquid peridotite (red line, Figure
1) is determined at high temperature and low $f$H2O, so its solubility law is
most appropriate for a terrestrial magma ocean (Sossi et al., 2022). For the
relevant pressure–temperature conditions of the magma ocean surface, the mole
fraction of H2O dissolved depends largely on variations in $f$H2O, while the
effect of temperature is a subordinate (Hamilton & Oxtoby, 1986), if not
poorly constrained, variable. As such, temperature does not explicitly appear
as a functional dependence in Equation 1, and $\alpha$ is taken to be
constant.

####  2.2.2 Solubility of CO2

There are no experimental constraints on the solubility of CO2 in peridotite,
so instead we use a solubility law formulated for basalt as a proxy. At low
CO2 fugacities, carbon dioxide dissolves in silicate melts as either molecular
CO2 or the CO${}_{3}^{2-}$ group (Stolper & Holloway, 1988). Because the
stability of the carbonate ion increases according to melt basicity (Mg2+,
Ca2+ cations, Holloway & Blank (1994)), the Mg-rich nature of planetary
mantles ensures that it occurs exclusively as CO${}_{3}^{2-}$ in such liquids.
Determining the mole fraction of dissolved carbonate for pressures up to 815
bars $f$CO2, Dixon et al. (1995) constrained the solubility of CO2 in mid-
ocean ridge basalt (MORB) melt:

$X^{m}_{\rm CO_{2}}=(3.8\times 10^{-7})f{\rm CO_{2}}\exp\left(-23\ \frac{f{\rm
CO_{2}}-1}{10RT}\right),$ (2)

where $X^{m}_{\rm CO_{2}}$ is the mole fraction of CO2 in the melt, $f{\rm
CO_{2}}$ is fugacity in bars, and $T$ temperature and $R$ gas constant in SI
units. The exponential function on the right-hand side is known as the
Poynting factor, which captures the influence of pressure and temperature on
CO2 fugacity. The abundance of CO2 (ppmw) in melt is then

$X_{\rm CO_{2}}=10^{4}\left(\frac{4400X^{m}_{\rm CO_{2}}}{36.594-44X^{m}_{\rm
CO_{2}}}\right).$ (3)

This solubility law compares favorably to experimental results at 10 kbar (Pan
et al., 1991). At water concentrations less than about 3 wt% there is a
negligible influence of the presence of dissolved H2O on CO2 solubility (Dixon
et al., 1995). Above about 6 wt%, water can enhance the solubility of CO2 by
around a factor of two (Iacono-Marziano et al., 2012), although CO2 remains
insoluble compared to H2O. In our model, the water concentration only
increases beyond 3 wt% at the end of the magma ocean stage for the most
oxidizing condition and largest initial hydrogen budget; otherwise, the water
concentration is significantly less. Hence, we ignore the minor influence of
both carbon and hydrogen on each other’s solubility.

####  2.2.3 Solubility of H2, CO, and CH4

We assume that the solubilities of H2, CO, and CH4 are negligible compared to
H2O and CO2, respectively, and therefore set $\alpha_{\rm H_{2}}=\alpha_{\rm
CO}=\alpha_{\rm CH_{4}}=0$ (Equation 1). This is motivated by high-pressure
experiments that find that the equilibrium constants for the dissolution of H2
gas in basaltic melts result in solubilities more than two orders of magnitude
lower than that for H2O (Li et al., 2015b; Hirschmann et al., 2012).
Similarly, the mole fraction of reduced carbon dissolved in silicate melts is
typically one to two orders of magnitude lower than for CO2 (Yoshioka et al.,
2019).

Methane dissolved in silicate melts has been identified by spectroscopy in
quenched glasses (Mysen et al., 2011; Ardia et al., 2013; Wetzel et al., 2013;
Armstrong et al., 2015; Dalou et al., 2019). These studies together
demonstrate that the fraction of dissolved methane increases in melts
equilibrated under increasingly reducing, high-pressure and H-rich conditions.
To quantify the solubility of methane in the melt at a given methane fugacity,
Ardia et al. (2013) performed experiments with Fe-free basaltic melts at
pressures between 0.7 and 3 GPa. The equilibrium constant of the reaction
CH4(l) = CH4(g) results in an order of magnitude less dissolved C for a given
$f$CH4 compared to an equivalent $f$CO2. Because estimated CH4 contents are
lower than several hundred ppmw even at $f$CH4 of 1000 bars (Ardia et al.,
2013), its solubility is neglected. It is noteworthy that solubility data for
CH4 in silicate melts at low total pressures relevant to magma oceans ($<$0.7
GPa) are lacking, preventing an accurate assessment of its solubility
behavior.

Table 1: H2O Solubilities Constrained by Experimental Studies

${\alpha}^{\ddagger}$ | $\beta$ | Composition | Temp. (K) | Pres. (bars) |  $f{\rm O_{2}}$ ($\Delta$IW) | Reference
---|---|---|---|---|---|---
534 | 2.0 | Peridotite | 2173 | 1.013 | –1.4 to +6.5 | Sossi et al. (2022)
683 | 2.0 | Lunar Glass | 1623 | 1.013 | –3.0 to +4.8 | Newcombe et al. (2017)
727 | 2.0 | Anorthite-Diopside | 1623 | 1.013 | –3.0 to +4.8 | Newcombe et al. (2017)
965 | 2.0 | Basalt (MORB) | 1473 | 176-717 | +3.9 to +13.2 | Dixon et al. (1995)
1007 | 2.0 | Basalt (MORB) | 1473 | 503-2021 | +3.5 and +7.9 | Berndt et al. (2002)
215 | 1/0.7 | Basalt | 1373 | 1034-6067 | unbuffered |  Wilson & Head III (1981)†

Note. — † Used in several modeling studies, e.g. Lebrun et al. (2013),
Salvador et al. (2017), Nikolaou et al. (2019), and Bower et al. (2019). ‡
$\alpha$ has units of ppmw/bar1/β and where relevant is determined by refits
to existing data by constraining $\beta=2$ and no solubility at zero fugacity.
We only consider experimental studies of basaltic or more mafic melts (i.e.,
Mg and Fe rich) at relatively low $f{\rm H_{2}O}$ ($\lesssim$2000 bars).

Figure 1: Experimentally determined solubility data of H2O for several melt
compositions and temperatures, for fugacities from 0 to 200 bars (Table 1).
Dotted lines denote $\beta=2.0$, and $\alpha$ is reported at the end of the
lines (Equation 1).

###  2.3 Redox Reactions

When at equilibrium, the magma ocean and atmosphere must be considered as a
single thermodynamic unit. This is because, based on the scaling of convective
velocity (Solomatov, 2000), the magma ocean mixing timescale is a few weeks,
which is considerably shorter than the cooling timescale of around 1 Myr.
Accordingly, $f{\rm O_{2}}$ defined by the magma ocean at its surface is
equivalent to the oxygen partial pressure (or fugacity) in the atmosphere.
Therefore, the fugacity ratio of gaseous redox couples (H2/H2O and CO/CO2), as
well as CO2–H2–CH4, satisfy chemical equilibrium at the $T$ and $f{\rm O_{2}}$
imposed by the magma ocean at its surface. We consider three reactions, where
all species exist in the gas phase, albeit partially dissolved in the melt
according to their respective solubilities (Sect. 2.2). For the hydrogen redox
couple,

$\displaystyle\rm{H}_{2}+0.5\rm{O}_{2}$ $\displaystyle=\rm{H}_{2}O,$ (4)
$\displaystyle\log_{10}{\rm K_{eq}}$ $\displaystyle=\frac{13152}{T}-3.039,$

for the carbon redox couple,

$\displaystyle\rm{CO}+0.5\rm{O}_{2}$ $\displaystyle=\rm{CO}_{2},$ (5)
$\displaystyle\log_{10}{\rm K_{eq}}$ $\displaystyle=\frac{14468}{T}-4.348,$

and for CO2–H2–CH4,

$\displaystyle\rm{CO}_{2}+2\rm{H}_{2}$ $\displaystyle=\rm{CH}_{4}+\rm{O}_{2},$
(6) $\displaystyle\log_{10}{\rm K_{eq}}$
$\displaystyle=\frac{-16276}{T}-5.4738.$

The equilibrium constants $\rm{K_{eq}}$ for the redox couples are determined
from the Gibbs free energy of reaction using data from the JANAF database fit
for temperatures between 1500 and 3000 K (Chase, 1998). The equilibrium
constant for the CO2–H2–CH4 reaction is taken from the IVTANTHERMO database
fit for temperatures between 300 and 2000 K (Schaefer & Fegley, 2017). The
most abundant redox-sensitive element in Earth’s mantle is iron, so the amount
of oxygen that is free to participate in reactions will be approximately
regulated by Gibbs free energy changes along the IW buffer during the magma
ocean stage. Due to the presence of melt and the difficulty of determining the
thermodynamic properties of a nonstoichiometric phase (wüstite), it is
preferable to consider equilibrium of solid metallic Fe with liquid FeO (the
IW buffer), where the Gibbs free energy of reaction constrains $f{\rm{O_{2}}}$
(O’Neill & Eggins, 2002):

$\displaystyle\rm{Fe}+0.5\rm{O}_{2}=\rm{FeO},$ (7)
$\displaystyle\log_{10}{f{\rm O_{2},IW}}=\frac{-244118+115.559T-8.474T\ln
T}{0.5\ln(10)RT}.$

It is convenient to define oxygen fugacity in log10 units relative to the IW
buffer ($\Delta{\rm IW}$):

$\Delta{\rm IW}=\log_{10}f{\rm O_{2}}-\log_{10}f{\rm O_{2},IW}.$ (8)

In our nominal models we set $\Delta{\rm IW}=0.5$ based on a recent
determination of the inferred oxygen fugacity of Earth’s magma ocean at its
surface using the Fe3+/Fe2+ ratio from a global compilation of peridotites
(Sossi et al., 2020). Hence, by fixing $f{\rm O_{2}}$ relative to the IW
buffer, H2O/H2 (Equation 4) and CO2/CO (Equation 5) are constrained as a
function of temperature, and by extension CH4/CO2 (Equation 6). Therefore,
while solubility controls the abundance of species in the atmosphere versus
the interior, redox reactions dictate the relative abundance of oxidized to
reduced species in the atmosphere. The thermodynamic coupling between the
atmosphere and the interior is set by the surface temperature at the magma
ocean–atmosphere interface. Appendices C.1 and C.3 describe how reactions are
self-consistently incorporated into the time stepper that evolves the magma
ocean–atmosphere system.

###  2.4 Atmospheric Escape

Atmospheric escape from terrestrial planets is challenging to parameterize
given the complex interaction of thermal and nonthermal processes, as well as
uncertainty regarding the history and details of the stellar environment and
planetary magnetic field (e.g., Gronoff et al., 2020). In the context of this
study, we concern ourselves with the impact of H2 loss from the atmosphere on
the evolution of the magma ocean–atmosphere system. We focus on hydrogen since
it is the lightest atmospheric component and therefore the most prone to
escape. Escape of H2 may be buffered by outgassing from the interior,
resulting in efficient depletion of hydrogen from the magma ocean. We apply a
model that obeys an upper limit for H2 loss based on the diffusion rate of H2
through a hydrostatic carbon-dominated atmosphere (Zahnle et al., 2019). Below
this limit, the escape flux follows the energy-limit that considers the
stellar energy input necessary for a volatile to escape Earth’s gravitational
field. Therefore, the H2 escape flux is (Zahnle et al., 2019)

$\phi_{\rm H_{2}}\approx\frac{\gamma 10^{16}\rm{[VMR]}_{\rm
H_{2}}S}{N_{A}\sqrt{1+0.006S^{2}}}\qquad\text{moles}\ {\rm m}^{-2}\ {\rm
s}^{-1},$ (9)

where $N_{A}$ is the Avogadro constant and ${\rm[VMR]}_{\rm H_{2}}$ is the
volume mixing ratio of H2 that provides coupling between the evolving
atmospheric speciation and H2 escape. Here $\gamma$ is a scaling factor (unity
by default), and $S$ is the normalized solar irradiation according to Earth’s
present-day value:

$S(t)=F_{\rm XUV}/F_{\rm XUV\odot},$ (10)

where XUV combines the influence of X-ray, extreme ultraviolet, and far
ultraviolet radiation. We adopt an upper estimate of $S=40$ at 4.0 Ga (Tu et
al., 2015). For ${\rm[VMR]}_{\rm H_{2}}$ of 30% and $\gamma=1$, the H2 mass-
loss rate is 2.1$\times 10^{5}$ kg s-1 (Equation C4). Both the escape
prefactor $\gamma$ (constant for a given case) and ${\rm[VMR]}_{\rm H_{2}}$
(time dependent) scale the mass-loss rate linearly.

##  3 Results

###  3.1 Outgassing of Hydrogen (One Earth Ocean)

Figure 2: Coupled interior and atmosphere evolution during magma ocean
outgassing at $\Delta$IW=0.5 with a hydrogen budget of one Earth ocean (no
carbon), where E1 shows equilibrium crystallization and F1 shows fractional
crystallization. For E1 and F1: (a, d) Interior temperature with gray showing
the mixed phase region where melt and solid coexist; (b, e) interior melt
fraction. (c) Atmospheric partial pressure of H2 and H2O for E1 and F1, where
the pressures of F1 (black lines) initially follow E1 (gray lines) before
departing at late time. Colored vertical ticks and lines show times in panels
(a), (b), (d), and (e), and labels show when a surface lid forms for E1 and
F1. (f) Interior depletion of melt (‘Melt’) and hydrogen (‘Hydr.’) for E1 and
F1.

We first present magma ocean outgassing models (Cases E1 and F1) where the
total hydrogen (H) budget is equivalent to the mass of H in Earth’s present-
day water ocean (hereinafter one ocean) (Figure 2). This budget is equivalent
to a fully outgassed atmosphere of 270 bars for a single species atmosphere
consisting of H2O. These models allow us to describe H2 and H2O outgassing
before subsequently introducing further complexity in the form of carbon
species, redox variation, and atmospheric escape. The total budget of H is
partitioned between a reduced (H2) and an oxidized (H2O) species according to
the $f{\rm O_{2}}$ buffer set at $\Delta$IW=0.5, with an initial surface
temperature of 2700 K. We recall that in our model H2 has zero solubility in
the melt and so can only reside in the atmosphere. The only parameter
difference between Case E1 (hereinafter E1) and Case F1 (hereinafter F1) is
the crystal size, where $a=$ 1 mm for E1 and $a=$ 5 cm for F1. Both of these
crystal sizes are reasonable for a crystallizing magma ocean and demonstrate
the end-member scenarios of equilibrium and fractional crystallization that
arise depending on the efficiency of melt–solid separation, which scales
according to the crystal size squared (Appendix B).

Both E1 and F1 chart a similar thermal evolution when the melt fraction is
high. Cooling proceeds quickly since most H remains in the magma ocean in the
form of H2O; this keeps the atmospheric pressure low ($<$ 1 bar). The
atmospheric composition (mole fraction of H2 compared to H2O) is determined by
chemical equilibrium, where $f{\rm O_{2}}$ decreases with decreasing surface
temperature (Equation 4) but parallels the IW buffer. Vigorous convection
keeps crystals mostly in suspension assuming efficient re-entrainment of
crystals (Solomatov et al., 1993; Solomatov & Stevenson, 1993a) and the
interior is approximately adiabatic. Due to the curvature of the melting
curves and the two-phase adiabat, the deepest mantle has the lowest melt
fraction and hence reaches the ‘rheological transition’ first. The rheological
transition is defined by an abrupt increase in viscosity around 30–40% melt
fraction, where the formation of an interconnected solid matrix begins to
dictate the convective timescale. The ”rheological front” is the interface
between melt-dominated dynamics above (lower pressure) and solid creep below
(higher pressure), and it moves upward through the mantle as cooling proceeds.
Hence, convection in the deep mantle becomes sluggish, akin to the present-day
solid mantle, acting as a brake on deep mantle cooling (Andrault et al.,
2016). Gravitational separation of melt and solid becomes a dominant driver of
energy transport at the rheological transition (Abe, 1993). Hence, the E1 and
F1 models exhibit different interior and outgassing evolution because of the
timescale of melt–solid separation (percolation) compared to net cooling rate.

E1 represents equilibrium crystallization, where the timescale for melt–solid
separation is longer than the timescale for the advancement of the rheological
front (Tonks & Melosh, 1990; Solomatov & Stevenson, 1993b). The rheological
front rapidly advances through the mantle, causing the thermal profile to
collapse on top of the rheological transition (Figure 2b). When it reaches the
surface, it initiates lid formation since convection has become so sluggish
that heat transport by convection to the surface cannot prevent top-down
cooling. Lid formation brings about the end of interior–atmosphere dissolution
equilibrium (‘E1 lid’ at 4.3 kyr; Figure 2). At this time, the mantle is
around 30% molten since it remains pinned at the rheological transition
because it has not had time to undergo any significant differentiation through
melt–solid separation prior to the formation of a surface lid. Therefore, the
mantle contains sufficient melt to trap the majority of H as dissolved H2O
beneath the surface lid (Figure 2f). Hence, the subsequent interaction of
these sequestered volatiles with the surface and atmosphere will be regulated
by geological processes operating over long timescales (millions to billions
of years), rather than by dissolution equilibrium with a comparatively short-
lived magma ocean.

By contrast, for F1 the melt percolation velocity keeps pace with the upward
velocity of the rheological front. This enables the complete solidification of
deep mantle layers below the rheological front by crystal settling, as well as
stalling the upward progression of the rheological front. This occurs because
efficient upward draining of melt keeps the upper regions molten at the
expense of cooling and fully crystallizing the lowermost mantle (Figure 2e).
Therefore, dissolution equilibrium with the atmosphere is maintained while
crystals form and are displaced to depth. Hence, F1 represents fractional
crystallization, where formed solids deep in the mantle can be considered
chemically isolated from the molten reservoir above. Below the rheological
transition, efficient melt–solid separation causes the thermal profile to
abruptly drop to the solidus, where now subsolidus cooling is greatly
restricted by the viscous transport timescale and partly by the core that
buffers the cooling of the mantle (Figure 2d). This ultimately enables a
substantial reduction of the melt reservoir ($98\%$ crystallized) at a surface
temperature of 1650 K compared to E1, where the crystal fraction only reaches
70% (Figure 2f). Due to the high solubility of H2O in melt, for F1 the
additional $\approx 30\%$ crystallization compared to E1 further depletes H
from the interior by $65\%$ before a surface lid forms. We note that
fractional crystallization could also occur for E1 after the surface lid has
formed, but in this case the melt reservoir is separated from the atmosphere
by a surface lid, and therefore outgassing is not driven by dissolution
equilibrium.

The atmospheric pressures of H2O and H2 increase in unison during cooling
(Figure 2c) because the relative abundance of these species is controlled by
chemical equilibrium at $\Delta$IW=0.5 (Equation 4). Since H2O is highly
soluble in silicate liquid, early in the lifetime of the magma ocean (prior to
1 kyr while melt fractions are high) the atmospheric mass remains mostly
constant and the partial pressures are controlled by the temperature
dependence of the $f{\rm O_{2}}$ buffer and equilibrium constant (Figure 2c,
Equation 4). After 1 kyr, the extent to which H outgasses is proportional to
the degree of crystallinity (i.e., continued cooling and crystallization)
according to the lever rule. Despite the mantle reaching 70% crystal fraction
for E1 before a surface lid forms, the atmospheric pressures of H2O and H2
only reach 4 and 2.4 bars, respectively (E1 lid, Figure 2c). For F1, the
atmospheric pressures of H2 and H2O (and hence their depletion in the magma
ocean) track E1 until around 2 kyr (gray lines, Figure 2c,f). At this time,
efficient melt–solid separation enables growth of a more substantial H2O (80
bars) and H2 (50 bars) atmosphere prior to surface lid formation (”F1 lid”,
Figure 2c). It takes around 4.3 and 26 kyr for a surface lid to form for E1
and F1, respectively, which compares favorably to a timescale of 10 kyr
determined by a model with a more advanced treatment of radiative transfer
(Lichtenberg et al., 2021). When a surface lid forms, the atmospheric pressure
of E1 and F1 differs by more than an order of magnitude.

###  3.2 Outgassing of Hydrogen (1–10 Earth Oceans)

Planet formation models predict substantial delivery of water to the inner
solar system (Raymond & Izidoro, 2017; Raymond et al., 2007), which can
account for the 1–10 oceans of water currently stored in Earth (Lécuyer et
al., 1998). Recent geochemical estimates propose that around 2.5 oceans are
currently stored in the mantle (Marty, 2012; Hirschmann, 2018). Therefore, we
now present cases the same as E1 and F1 but with an initial H budget up to 10
Earth oceans. Since the main features of the coupled evolution of the
interior–atmosphere are described for E1 and F1 (Section 3.1), only
differences relative to these fiducial cases are now presented. Figure 3
reveals the interior depletion of H as a percentage of the H inventory, where
the circles with a solid line and containing a 1 are E1 (Figure 3a) and F1
(Figure 3b).

For both equilibrium and fractional crystallization, as the H inventory is
increased from 1 to 10 oceans, the magma ocean lifetime increases by around
two orders of magnitude, reaching a maximum duration slightly less than 1 Myr.
This is because the optical depth of the atmosphere scales with pressure, and
a larger inventory results in a larger outgassed atmosphere and hence a higher
atmospheric pressure. Therefore, due to the thicker atmosphere, the magma
ocean takes longer to cool before the rheological front reaches the surface
and initiates lid formation. In addition, the interior depletion of H
increases for larger total inventories as a direct consequence of power-law
solubility where the exponent $\beta=2$. This is because, during cooling, the
atmospheric reservoir of H grows with the water mole fraction to the second
power (see Appendix C.6 for analysis). Hence, for H outgassing, both the
pressure of the atmosphere and relative depletion of the interior increase
with the total inventory at a given melt fraction. For a volatile species that
obeys a linear solubility relation ($\beta=1$), its relative depletion at a
given melt fraction is independent of its total inventory (Appendix C.6).
Water is expected to transition to such a linear solubility above around 10
ocean masses on Earth owing to the prevalence of dissolved H2O.

Figure 3: Depletion of hydrogen (H) from the interior during magma ocean
outgassing at $\Delta$IW=0.5 for H budgets from 1 to 10 Earth oceans: (a)
equilibrium crystallization, and (b) fractional crystallization. Lines show
constant melt fraction, and circles show the H budget. Magma ocean stage ends
at around 30% melt for equilibrium crystallization and 2% melt for fractional
crystallization. Cases E1 and F1 (Figure 2) correspond to 1 H ocean for
equilibrium and fractional crystallization, respectively. Dashed lines in
panel (a) with a ‘W’ suffix denote models that only include H as H2O, i.e. a
completely oxidized mantle.

For equilibrium crystallization, lid formation occurs soon after 30% melt is
reached (Figure 3a). For all equilibrium cases, melt–solid separation is slow
compared to the cooling timescale, so when a lid forms, the melt fraction is
only marginally reduced for cases with a long compared to a short cooling
timescale. At 30% melt, a maximum of around 23% of the total H budget is
outgassed for 10 oceans, although this decreases to about 4% for 1 ocean
(Figure 3a). For fractional crystallization, the melt fraction decreases below
30% until lid formation occurs around 2% melt. Compared to equilibrium
crystallization, this extra depletion of melt drives outgassing of hydrogen to
around 80% of its total budget (Figure 3b). Regardless of total inventory,
depletion must tend to 100% when the melt fraction reaches zero because we
assume that no volatiles are retained in solids. This explains the flattening
of the depletion curve for 2% melt fraction compared to, for example, 10% melt
fraction.

The magma ocean lifetimes for $\Delta{\rm IW}=0.5$ cases (solid lines, Figure
3a) are comparable to cases that are fully oxidized where H can only exist as
water (W cases, dashed lines, Figure 3a). This is because for all cases H2O
dominates the opacity of the atmosphere. However, H depletion of the
$\Delta$IW=0.5 cases and that of the fully-oxidized (W) cases differ by up to
about 11%. This is because at $\Delta$IW=0.5 outgassing of H2O must be
accompanied by commensurate growth of H2 in the atmosphere to maintain
chemical equilibrium. Since the H budget is conserved, increasing H2 in the
atmosphere necessitates a decrease in the H2O reservoir dissolved in the magma
ocean. In short, H2O solubility determines H2O outgassing, but oxygen fugacity
dictates the additional depletion of interior hydrogen to maintain equilibrium
between H2 and H2O in the atmosphere.

###  3.3 Outgassing of Hydrogen and Carbon

Figure 4: Atmospheric composition during magma ocean outgassing at
$\Delta$IW=0.5. Columns from left to right show hydrogen (H) budgets of 1, 3,
and 10 oceans, respectively. Rows from top to bottom show C/H by weight of 0,
0.1, 1, and 5. Magma ocean lifetime corresponds to the termination of the
colored region (at $\approx$ 2% melt for fractional crystallization, between
$10^{4}$ and $10^{7}$ yr). White dotted lines correspond to 30% melt, which
indicates both the duration and atmospheric composition if instead the mantle
underwent equilibrium crystallization. Initial partial pressures of volatile
species ($p_{\rm init}$) are for a completely molten magma ocean and $T_{\rm
surf}=2700$ K. Final partial pressures ($p_{\rm final}$) are for complete
solidification and outgassing of all volatile species and $T_{\rm surf}=1400$
K.

In contrast to the high solubility of water in silicate melt, carbon dioxide
is relatively insoluble (e.g., Dixon et al., 1995). Hence, due to the high
abundance of C in rocky planets, the presence of gaseous CO, CO2, and CH4 may
regulate the lifetime of a magma ocean and thereby the timescale over which H
outgasses. Carbon also suppresses the outgassing of H through its influence on
the mean molar mass of the atmosphere (Bower et al., 2019). Therefore, we now
present fractional crystallization cases that additionally include C species
(CO2, CO, and CH4), with their relative abundances in the gas phase determined
by redox reactions. The initial carbon-to-hydrogen ratio (C/H by weight) is
varied from 0.1 to 5, which is a range compatible with C/H ratios in
chondritic meteorites and in the bulk silicate Earth (1.1 and 1.4 by mass,
respectively (Hirschmann, 2016)).

We independently verified the results of our model by tabulating the final
total pressure and composition (in terms of moles of H, C, and O) of select
outgassed atmospheres calculated using our model (Appendix E). These were used
as input parameters for the Equilib module of FactSage 8.0 (Bale et al.,
2009), which calculates the equilibrium partial pressures of gases and stable
condensed phases using a Gibbs free energy minimizer with a database of more
than 40 gas species in the H-C-O system. For the sake of comparison, we
assumed ideal gas behavior. Provided that graphite is not predicted to
precipitate (Section 4.4) the agreement with our models, which utilize a
comparatively simple chemical network and five gas species, is commendable;
the partial pressures differ by at most a few percent (Table E).

Figure 4 compares the influence of the C/H ratio on magma ocean outgassing for
an H budget of 1, 3, and 10 oceans. Cases F1 (Figure 4a), F3 (Figure 4b), and
F10 (Figure 4c) serve as carbon-free reference cases. We first focus attention
on increasing C/H for an H budget of one ocean (left column, Figure 4). Even
for an addition of only 10% C by weight (i.e., C/H=0.1), the atmosphere is
dominated by C species (almost 90% volume mixing ratio) for a fully molten
mantle. This is due to the low solubility of CO2 compared to H2O, coupled with
equilibrium chemistry that establishes the partial pressure of CO around a
factor of seven greater than the partial pressure of CO2. During magma ocean
cooling, the CO pressure in the atmosphere steadily decreases. This occurs as
a result of the compounding effects of the equilibrium constant of Equation 5,
which results in a decrease of $f{\rm CO}/f{\rm CO_{2}}$ as the surface
temperature decreases, and the continued outgassing of CO2 as the melt
fraction decreases. The latter also drives an initial increase in CO2
pressure, although the CO2 pressure decreases later once H begins to outgas in
earnest. Hence, even though the ever-decreasing melt fraction drives
outgassing of both H and C throughout the magma ocean evolution, the partial
pressure of volatile species can either mostly decrease (CO), mostly increase
(H2, H2O), or increase and then decrease (CO2) (Figure 5).

Figure 5: Atmospheric pressures during magma ocean outgassing for different
C/H and $f{\rm O_{2}}$. Columns from left to right show hydrogen (H) budgets
of 1, 3, and 10 oceans, respectively. Rows from top to bottom show $f{\rm
O_{2}}$ of $\Delta$IW=$-2$, $\Delta$IW=0.5, and $\Delta$IW=4. Colors denote
volatiles (red: CO; green: CO2; orange: H2; blue: H2O; purple: CH4, where
solid lines show C/H=1 and dashed lines show C/H=5.

For an H budget of one ocean, the initial pressure of H2 and H2O is 0.5 and
0.4 bars, respectively, regardless of C/H. However, at higher initial C/H
ratios, the H species constitute a smaller fraction (volume mixing ratio) of
the atmosphere owing to the relative insolubility of C-bearing species
compared to H2O. For C/H=5, hydrogen species constitute less than 1% of the
atmosphere during the early evolution of the magma ocean ($>$30% melt). For
this case, the mixing ratio of H species increases by a factor of 230 as the
magma ocean cools from fully molten to fully solid (Figure 4j). In contrast,
for C/H=0.1 the mixing ratio of H species only increases by a factor of about
nine (Figure 4d). Even though the total inventory of H is fixed, the final
outgassed pressures of H2 and H2O depend on C/H, increasing from 62 and 133
bars (C/H=0) to 102 and 221 bars (C/H=5.0), respectively. White dotted lines
in Figure 4 denote a melt fraction of 30%, which corresponds to when a surface
lid forms for equilibrium crystallization. Hence, for equilibrium
crystallization, outgassed atmospheres are usually more dominated by C species
versus H species. For C/H=0.1, the mixing ratio of H species increases by
about a factor of two from 30% melt to 2% melt (Figure 4d). In contrast, for
C/H=5, H species increase from about 1% to more than 50% over the same range
of melt fraction (Figure 4j).

Figure 6: Depletion of hydrogen (H) versus melt fraction at $\Delta$IW=0.5 for
H budgets of 1, 3, and 10 oceans and different C/H by weight. For each shaded
region, the upper bound is C/H=0, lower bound is C/H=5, and black line is
C/H=1. Symbols indicate times of C/H=1 cases near 20% melt and 5% melt.

For an inventory of three oceans (middle column, Figure 4), a larger fraction
of the atmosphere consists of H species during the evolution compared to an
inventory of one ocean. This is a direct consequence of the power-law
solubility of H2O and the higher amount of H (Appendix C.6). When no melt
remains—as for the end of fractional crystallization—the entire initial
inventory of H and C sets the atmospheric composition where equilibrium
chemistry continues to dictate the partitioning between reducing and oxidizing
species. Therefore, for fixed C/H at the end of fractional crystallisation,
the volume mixing ratios are also set. If CH4 is insignificant, the final
partial pressures ($p_{\rm final}$) scale according to the initial inventory;
for example, the final pressures for C/H=0.1 are a factor of three larger for
a three-ocean inventory compared to a one-ocean inventory (compare $p_{\rm
final}$ in Figure 4d,e). For large H inventories, the volume mixing ratio of H
species changes less during outgassing because initially more H (compared to
C) resides in the atmosphere owing to the power-law solubility of H2O. The
case with a three-ocean budget and C/H=1 (Figure 4h, 5e) is similar to the
bulk silicate Earth calculated in Sossi et al. (2020). At high temperatures, a
CO-dominated atmosphere forms with $\sim$200-bar C-bearing species and
$\sim$10-bar H-bearing species in which fH2 and fH2O are subequal, while fCO
$>$ fCO2.

Figure 7: Depletion of hydrogen (H) versus time at $\Delta$IW=0.5 for
different C/H, where colored lines show constant melt fraction. Dashed lines
show C/H=0 for H budgets from 1 to 10 Earth oceans; these are the solid lines
in Figure 3b. Departing from the dashed lines, circles show 3 oceans and
C/H=X, where X is 0, 1, or 5. Similarly, rectangles show 10 oceans and C/H=X,
where X is 0, 1, or 5.

Figure 6 summarizes the interior depletion of H for 1–10 oceans and C/H by
weight varying from zero to five. To first order, depletion at a given melt
fraction is dictated by the total inventory of H due to the power-law
solubility of H2O. For 10% melt, this gives rise to around a 40% increase in
depletion as the H inventory increases from 1 to 10 oceans. A secondary
influence on depletion is C/H, where higher C/H suppresses the outgassing and
hence depletion of H. This effect is not connected to the solubility of H2O,
but rather is due to the dominant presence of C species in the atmosphere that
increases the mean molar mass of the atmosphere relative to H species
(Appendix C.1). Figure 7 reveals the influence of C/H on the depletion during
magma ocean outgassing and can be compared with Figure 3. For a fixed H
inventory, increasing C/H prolongs the time taken to reach a given melt
fraction since the optical depth of the atmosphere increases. Increasing C/H
increases the volume mixing ratio of C species in the atmosphere, and this
alone would actually decrease the opacity of the atmosphere since CO has the
lowest opacity of all the considered volatiles (Appendix A). However,
increasing C/H also increases the total surface pressure substantially because
CO2 is relatively insoluble compared to H2O. Hence, the optical depth
increases in response to the larger surface pressure, and this increase
predominates over the change in speciation; the net outcome is prolonged
cooling for larger C/H. Figure 7 further demonstrates the trend of Figure 6,
in which increasing C/H suppresses the outgassing and hence interior depletion
of H at a given melt fraction.

###  3.4 Oxygen Fugacity

Figure 8: Atmospheric composition during magma ocean outgassing for C/H=1 and
different $f{\rm O_{2}}$. Columns from left to right show hydrogen (H) budgets
of 1, 3, and 10 oceans, respectively. Rows from top to bottom show $f{\rm
O_{2}}$ of $\Delta$IW=$-2$, $\Delta$IW=0.5, $\Delta$IW=2, and $\Delta$IW=4.
Magma ocean lifetime corresponds to the termination of the colored region (at
$\approx$ 2% melt for fractional crystallization, between $10^{4}$ and
$10^{7}$ yr). White dotted lines correspond to 30% melt, which indicates both
the duration and atmospheric composition if instead the mantle underwent
equilibrium crystallization. Initial volatile partial pressures ($p_{\rm
init}$) are for a completely molten magma ocean and $T_{\rm surf}=2700$ K.
Final partial pressures ($p_{\rm final}$) are for complete solidification and
outgassing of all volatiles and $T_{\rm surf}=1400$ K. Figure 9: Depletion of
hydrogen (H) vs. melt fraction for H budgets of 1 and 10 oceans and different
$f{\rm O_{2}}$. For each shaded region, the upper bound is $\Delta$IW=$-2$,
lower bound is $\Delta$IW=4, and black line (same as Figure 6) is
$\Delta$IW=0.5. Symbols indicate times of $\Delta$IW=0.5 cases near 20% and 5%
melt. Figure 10: Depletion of hydrogen (H) vs. time for different $f{\rm
O_{2}}$, where colored lines show constant melt fraction. Dashed lines show
$\Delta$IW=0.5 and C/H=0 for H budgets from 1 to 10 Earth oceans; these are
the solid lines in Figure 3b. Departing from the dashed lines, circles show 3
oceans and C/H=1 for $\Delta$IW=X, where X is $-2$, 0.5, or 4. Similarly,
rectangles show 10 oceans and C/H=1 for $\Delta$IW=X, where X is $-2$, 0.5, or
4.

Iron-core formation established the deep Earth’s redox state at
$\Delta$IW=$-2$ whereas the present-day upper mantle redox is $\Delta$IW=4.
This indicates that mantle redox can vary with both time (e.g., Scaillet &
Gaillard, 2011) and pressure (Armstrong et al., 2019), potentially giving rise
to atmospheres that are more reduced or oxidized than our nominal cases at
$\Delta$IW=0.5 (e.g., Hirschmann, 2012). Hence, we supplement our calculations
at $\Delta$IW=0.5 by additionally considering outgassing scenarios at
$\Delta$IW=$-2$, $\Delta$IW=2, and $\Delta$IW=4. Figure 8 summarizes the
results for C/H=1, and detailed figures for all C/H cases at each $f{\rm
O_{2}}$ are presented in Appendix D. Increasing the oxygen fugacity from
$\Delta$IW=0.5 to $\Delta$IW=4 progressively increases the ratio of oxidized
to reduced species in the atmosphere at a given temperature. Since oxidized
species of C and H have a higher molar mass than their reduced counterparts,
this also increases the mean molar mass of the atmosphere which influences
partial pressures through mass balance (e.g., Equation C1). Cooling times are
generally extended for oxidized versus reduced atmospheres, due to the higher
surface pressure and the intrinsic higher opacity of oxidized species
(Appendix A). However, for $\Delta$IW=$-2$ and C/H=5, production of CH4
increases the cooling time to be comparable to an oxidized atmosphere
dominated by CO2 and H2O.

Oxidized interiors mitigate the interior depletion of H relative to reducing
conditions (e.g., compare $\Delta$IW=$-2$ and $\Delta$IW=4 in Figure 9 and
10). At 30% melt, H depletion is around a factor of 10 larger for
$\Delta$IW=$-2$ compared to $\Delta$IW=4. At 2% melt, H depletion at
$\Delta$IW=$-2$ is approximately 20% larger than $\Delta$IW=4; we recall that
depletion for all cases must reach 100% when no melt remains. Depletion of H
is largest for $\Delta$IW=$-2$ because it exists mainly as the less soluble H2
and CH4 in the gas, so only a small quantity of H is dissolved in the melt
reservoir as H2O. Whereas increasing C/H mitigates H depletion and extends
cooling times, decreasing $f{\rm O_{2}}$ greatly enhances H depletion (i.e.,
outgassing) and generally reduces cooling times. For C/H=1, early atmospheres
are almost always dominated by carbon species, as either CO (reduced cases) or
CO2 (oxidized cases). For 10 H oceans and $\Delta$IW=$-2$ the early atmosphere
is instead dominated by H2. Even for the most oxidized scenario when
$\Delta$IW=4, H2O is never the dominant species in the atmosphere until the
melt fraction drops below about 20% or even less. Hence, magma oceans that
undergo fractional crystallization are more likely to produce oxidized and
H2O-rich atmospheres versus equilibrium crystallization which produces reduced
and often CO-rich atmospheres.

###  3.5 Atmospheric Escape of H During Outgassing

Figure 11: Comparison of atmospheric composition and interior depletion at
$\Delta$IW=0.5 for (1) C/H=0.1, three H oceans, $\gamma=1$ (2) C/H=1, three H
oceans, $\gamma=1$ (3) C/H=0.1, one H ocean, $\gamma=1$ and (4) ”Large escape”
with C/H=0.1, three H oceans, and an escape prefactor $\gamma=1000$ (Equation
9). (a, b) H2 partial pressure and volume mixing ratio; (d, e) CO partial
pressure and volume mixing ratio; (c, f) interior depletion of H and C,
respectively.

During the magma ocean stage, escape of H2—the lightest atmospheric
component—could impact the reservoir evolution of volatiles, so we explore
this possibility in Figure 11. The reference case has C/H=0.1, three H oceans,
and unity escape prefactor ($\gamma=1$, Equation 9). However, an equivalent
case with no escape is visually indistinguishable from $\gamma=1$. This
immediately demonstrates that H2 escape due to irradiation of a
carbon–hydrogen atmosphere does not appreciably alter the volatile reservoirs
during the magma ocean stage for an Earth-sized body at 1 AU, largely because
of the magma ocean’s short lifetime ($\sim$105 yr). Even for a larger H budget
of 10 oceans where cooling is relatively prolonged, $\gamma\lesssim 100$ has
no appreciable impact on the evolution of the volatile reservoirs. Therefore,
we show a case in which escape is increased by setting $\gamma=1000$ (”large
escape”) to investigate an extreme end-member escape scenario, which may
arise, for example, as a result of tidal effects. For comparison,
$\gamma=1000$ gives a mass-loss rate around $2.1\times 10^{8}$ kg s-1, which
is within the range of rates determined for volatile stripping of the proto-
Moon by Earth (Charnoz et al., 2021) or for a Mars-sized body near 1 AU
undergoing EUV-driven thermal escape (Benedikt et al., 2020). For large escape
the inventory of hydrogen decreases by almost 50% over the lifetime of the
magma ocean of around $10^{5}$ yr. By comparison, the magma ocean lifetime for
an equivalent case without escape is around $1.5\times 10^{5}$ yr.

Early on (equivalently times at which melt $>$40%), the H2 partial pressure of
all cases with three H oceans follows a comparable trajectory. At later times
(equivalently $<$40% melt), the cases with three oceans and $\gamma=1$
continue to track a similar trend, although the final H2 partial pressure
depends on the initial C/H (e.g., Figure 4e, h). For high escape rate,
however, the growth of H2 in the atmosphere is mitigated by escape such that
the H2 partial pressure is reduced compared to the cases with an unchanging H
inventory of three oceans. Nevertheless, equilibrium chemistry continues to
dictate $f{\rm H_{2}}/f{\rm H_{2}O}$ so only the magnitudes of the partial
pressures of H species are modulated by escape. The buffering of the
atmosphere by the magma ocean presupposes that the loss of H2 is not
sufficient to influence the redox state of the magma ocean itself. In
practice, the preferential loss of H2 promotes oxidation of the residual
mantle (e.g., Olson & Sharp, 2019). The precise amount depends on the
abundance of H, but the net effect is to self-arrest loss of H as H2. This is
because as escape proceeds, necessarily more H2 is converted to H2O, given
that O-bearing species are heavier and thus escape less readily than H2. Since
for large escape almost 50% of the initial three-ocean inventory is lost,
$f{\rm H_{2}}$ is steered toward the final pressure (at 0% melt) of the one-
ocean case, which has the same C/H=0.1. For fixed C/H=0.1, the volume mixing
ratio of H2 is the same at 0% melt regardless of the final H inventory.
However, for large escape, H is lost relative to C, and hence the H2 volume
mixing ratio at the end of the magma ocean stage is slightly reduced compared
to cases with unchanging C/H=0.1.

Atmospheric escape modulates two controls on outgassing during magma ocean
cooling that we have previously investigated: (1) escape decreases the
inventory of H and hence modulates the partitioning of H between the interior
and atmosphere when melt is present according to the power-law solubility of
H2O, and (2) escape increases C/H, resulting in a final atmosphere richer in C
compared to H. For example, the loss of H2 for large escape increases the
volume mixing ratio of CO at later time compared to C/H=0.1. This abides by
our previous models that show a larger contribution of carbon species to the
atmosphere for a larger C/H. For large escape, C/H by weight increases from
0.1 initially to 0.2 at the end of the magma ocean stage. For minimal escape,
both decreasing the H inventory and/or increasing C/H suppresses the interior
depletion of H at a given melt fraction, relative to the total inventory
(Figure 6). However, for large escape, the depletion of H is enhanced relative
to other cases (Figure 11). This demonstrates that the direct loss of H due to
escape is a stronger control on its depletion in the magma ocean than the
secondary effects on its solubility (through H inventory) and the suppression
of H outgassing due to a heavier, carbon-rich atmosphere. In short, H2 escape
driven by irradiation given best estimates for the XUV flux of the young Sun
has an insignificant impact on volatile reservoir evolution of an Earth-sized
body at 1 AU during the magma ocean stage. However, more extreme escape
scenarios can reduce the partial pressure of H species in the atmosphere and
drive faster depletion of the interior H reservoir.

##  4 Discussion

###  4.1 Atmospheric Composition and Evolution

Figure 12: Summary of atmospheric composition and hydrogen depletion of the
interior during magma ocean cooling at 100%, 30%, and 0% melt, where the
surface temperature is given below the melt fraction. Transition from carbon-
dominated to hydrogen-dominated atmospheres depends on surface lid formation
and the style of magma ocean crystallization, which may prevent or hinder
interior outgassing below 30% melt fraction.

For a redox state inferred for the early Earth, the earliest outgassed
atmospheres of terrestrial planets are dominantly CO-rich by volume mixing
ratio when the mantle is mostly molten (Figure 12). This is borne out of the
equilibrium chemistry of CO–CO2 and the low solubility of C species compared
to H species. With C present, H species can in principle constitute more than
50% of the atmosphere if the H budget is greater than about 6 Earth oceans
owing to the power-law solubility of H2O. However, this is only possible if
C/H by weight is low—around an order of magnitude less than for the bulk
silicate Earth. As C/H increases to unity (Earth-like) and beyond, the
atmosphere is CO-rich regardless of the H budget. Hence, only for C/H$\lesssim
0.1$ is an H-rich atmosphere dominant during the early cooling phase, when
most of the mantle is molten. Throughout the magma ocean stage, H2 and H2O
have roughly constant mixing ratios owing to the weak temperature dependence
of the equilibrium constant (Equation 4), while $f{\rm H_{2}}/f{\rm H_{2}O}$
is approximately unity at $\Delta$IW=0.5 at magmatic temperatures.

For C/H$\gtrsim$1, H species can only begin outgassing in earnest once the
melt fraction drops below about 30%, regardless of the redox state. However,
this is contingent on the surface remaining molten to ensure equilibrium
between the interior and atmosphere; this is less likely if the magma ocean
undergoes equilibrium crystallization where crystals and melt freeze together
and a surface lid forms quickly. In comparison, fractional
crystallization—driven by melt–solid separation—can more likely maintain a
molten surface while the deeper mantle crystallizes and exsolves volatiles.
Hence, the transition from a C-rich to H-rich atmosphere depends on the style
of crystallization, and crucially the behavior of the mantle as it transitions
from mostly molten to mostly solid. A large steam-dominated atmosphere can
only form for a relatively oxidized mantle toward the end of mantle
crystallization, although it may persist for a longer duration than the early
carbon-rich atmosphere in the absence of an efficient loss mechanism. An
atmosphere that contains initially more H species by volume mixing ratio
(C/H$\lesssim 0.1$ and particularly for H inventories $>$ 6 oceans) undergoes
less modification during outgassing to reach its final volatile mixing ratios
when the mantle has fully crystallized.

The chemical boundary layer at the surface of the magma ocean facilitates
interior–atmosphere equilibrium, where the equilibrium timescale is rapid
relative to a crystallization timescale of around a million years (Pahlevan et
al., 2019). For low $f{\rm O_{2}}$ or small atmospheres, the crystallization
timescale decreases by no more than two orders of magnitude ($\sim 10^{4}$ yr)
and this decrease would be further mitigated by a greater H2 greenhouse
forcing (Lichtenberg et al., 2021). Other mechanisms to facilitate
equilibrium, such as bubble formation due to volatile supersaturation (Elkins-
Tanton, 2008), would also reduce the equilibrium timescale. Nevertheless, if
crystallization proceeds too rapidly compared to the equilibrium timescale,
then disequilibrium chemistry may arise between the interior and atmosphere.

###  4.2 Hydrogen in the Interior

For otherwise highly soluble volatiles, such as H2O, to exist in the
atmosphere, the magma ocean must have crystallized below 30% melt fraction
before a surface lid could form and persist at the surface. Otherwise,
thermodynamic communication between the interior and atmosphere is broken, and
dissolved volatiles remain trapped in the mantle. Below 30% melt fraction,
continued outgassing of dissolved volatiles requires that significant
solidification occurs through melt percolation and solid compaction before the
rheological front reaches the surface or a quench crust forms (fractional
crystallization). In an end-member scenario, this can enable near-complete
outgassing of all volatiles, which has received the most attention in the
literature. At a given melt fraction, we find that a larger fraction of H can
be retained in the interior for smaller total inventories of H, larger C/H by
weight, and more oxidized interiors. Although C does not influence H
solubility directly, it impacts H retention since it suppresses H outgassing
through its influence on the mean molar mass of the atmosphere (Bower et al.,
2019).

Equilibrium crystallization results in a substantial reservoir of melt and
hence volatiles trapped in the interior when the surface forms a lid. This
later possibility has received less attention in the context of volatile
evolution during the magma ocean stage, but it provides a mechanism to safely
harbor a significant quantity of volatiles in the mantle to protect against
loss from atmospheric escape and impacts. Hence, the high solubility of H2O in
silicate melt may be a crucial property of this life-supporting molecule that
enables it to survive during the violent early years of a terrestrial planet’s
life.

Volatile retention would be further enhanced if a melt–crystal density
crossover enabled the formation of a voluminous basal magma ocean (Labrosse et
al., 2007; Caracas et al., 2019) or if melt is captured as the rheological
front advances through the mantle (Hier-Majumder & Hirschmann, 2017). Trapped
melt, and hence soluble volatiles, could then be sequestered in the deep
mantle by Rayleigh-Taylor instability (Maurice et al., 2017; Miyazaki &
Korenaga, 2021). Moreover, we assume no solubility of H in crystallizing
minerals, whereas experimental data indicate significant quantities of water
may be stored in ringwoodite (Fei & Katsura, 2020). In short, there are
additional processes not included in our models that further conspire against
the complete outgassing of soluble volatiles, namely, water, during the magma
ocean stage. Furthermore, estimates of the volatile concentration in the bulk
silicate Earth are incompatible with complete outgassing (e.g., Hier-Majumder
& Hirschmann, 2017).

Therefore, complete outgassing of soluble volatiles—frequently an outcome of
magma ocean models—only occurs with fractional crystallization if the
aforementioned processes are absent or inefficient. Furthermore, it requires
that a surface lid or quench crust, if present, does not hinder volatile
outgassing as fractional crystallization proceeds. Hence, at the end of the
magma ocean stage it is reasonable to expect that a trapped reservoir of H2O
in the interior will interact with the surface environment and atmosphere.
This could occur owing to post-magma ocean cumulate overturn or owing to
processes operating over geological timescales (millions to billions of
years). The detailed chemical and physical processes governing the formation
and sustenance of cumulate layers and a surface lid in a dynamic magma ocean
require further investigation; external influences such as projectile
bombardment could also stifle lid formation (e.g., Perera et al., 2018).

For an Earth-like planet orbiting a Sun-like star at 1 AU during its magma
ocean stage, we find that atmospheric escape of H2 due to irradiation (energy
and diffusion limited) does not significantly impact the evolution of volatile
reservoirs, in agreement with Hamano et al. (2013) and Katyal et al. (2020).
This is because escape rates are sufficiently low and the magma ocean duration
is sufficiently short (at most a few million years). In future work, a
feedback to probe is that loss of H2 increases $f{\rm O_{2}}$ slightly,
causing more H2 to convert to H2O and thereby partly self-arresting the loss
process. Furthermore, we have not considered photodissociation of H2O, but
this could interplay with the included geochemical reactions that depend on
mantle redox.

###  4.3 Crystallisation Style

Based on energetic considerations of convection versus gravitational settling,
Solomatov & Stevenson (1993a) propose a critical crystal size of 1 cm above
which fractional crystallization is inevitable. Furthermore, they show that
the critical crystal size is weakly dependent on crystal fraction and
therefore depth. Hence, their results are compatible with our fractional
crystallization models with a constant crystal size of 5 cm. In detail,
gravitational settling in our models becomes dominant at the rheological
transition owing to the substantial decrease of turbulence that enabled
crystals to remain mostly in suspension. It is reasonable to assume that
crystals are suspended at high melt fraction (Tonks & Melosh, 1990), although
they could begin to settle and thereby initiate fractional crystallization
before the rheological transition is obtained (Patoc̆ka et al., 2020),
assuming inefficient re-entrainment (Solomatov et al., 1993) or a high
planetary rotation rate (Maas & Hansen, 2019). Settling also depends on
mineral buoyancy, and if the minerals float rather than sink, this could
hinder interior–atmosphere communication. For example, preferential
partitioning of Fe into melt during crystallization continuously changes the
level of neutral buoyancy between crystals and melt (Caracas et al., 2019). In
contrast to fractional crystallization, equilibrium crystallization prevents
differentiation of the mantle and occurs for a smaller crystal size of 1 mm
(Solomatov & Stevenson, 1993b).

Most previous dynamic models of magma ocean cooling do not consider melt–solid
separation (Lebrun et al., 2013; Hamano et al., 2013; Schaefer et al., 2016;
Nikolaou et al., 2019); rather, the interior is assumed to always be
adiabatic. Yet some of these previous models report agreement with the
solidification timescale derived from geochemical models that explicitly
consider fractional crystallization. The previous models implicitly assume
that although the deep mantle reaches the rheological transition first, it
continues cooling as efficiently by convection as the uppermost mantle; this
allows the melt fraction to continue decreasing at a similar rate. However,
this is unlikely given the rheological behavior of a melt–solid aggregate
(e.g., Costa et al., 2009), which at the rheological transition predicts a
large reduction in the convective velocity and rapid melt–solid separation
(Abe, 1993). Hence, somewhat coincidentally, previous models decrease the melt
reservoir and outgas volatiles at a similar rate to geochemical fractional
crystallization models, even though the key ingredient to justify fractional
crystallisation (i.e., melt–solid separation) is not included.

###  4.4 Graphite, Diamond, and Water Precipitation

Our model assumes that all volatile species dissolved at the surface of the
magma ocean continue to remain so, irrespective of pressure (depth of the
magma ocean) and temperature (cooling of the magma ocean). Phase
transformations involving these volatile elements are likely to occur as
pressure and temperature change (e.g., Hirschmann, 2012). However,
experimental data for the speciation and partitioning of H-, C-, and O-bearing
volatile species rarely exceed 7 GPa (Section 2.2), which precludes a
holistic, bottom-up model of magma ocean crystallization. Nevertheless, we now
systematically investigate possible phase transformations of H and C to
identify scenarios under which the pressures of outgassed atmospheres could
diverge from those calculated by our model.

Reducing conditions may induce graphite or diamond precipitation (Hirschmann,
2012; Takahashi et al., 2013; Keppler & Golabek, 2019). The precipitation of
graphite occurs if the CO fugacity of the atmosphere exceeds that defined by
the graphite–carbon monoxide (CCO) buffer:

$\displaystyle\rm{C(s)}+0.5\rm{O_{2}(g)}$ $\displaystyle=\rm{CO(g)},$ (11)
$\displaystyle\log_{10}{\rm K_{eq}}$ $\displaystyle=\frac{6254}{T}+4.334,$

where $\mathrm{K_{eq}}$ as a function of temperature is provided by the JANAF
database (Chase, 1998). Equation 11 indicates that $\log_{10}f$CO in
equilibrium with graphite decreases proportional to $0.5\log_{10}f$O2 and
increases with temperature. Thus, graphite precipitation is favored under more
reducing conditions and lower temperatures (Figure 13). For our cases at
$f\rm{O_{2}}$s of $\Delta$IW=0.5 or more oxidized, graphite precipitation
would not occur because $f$CO never exceeds that defined by the CCO buffer.
This suggests that it is unlikely that a terrestrial magma ocean was saturated
in graphite at high melt fractions ($>30$%). However, it may occur for most
atmospheres as the surface cools below the temperature at which a surface lid
forms around 1650 K (Sossi et al., 2020). At this stage, however, the
atmosphere and mantle are not necessarily in equilibrium and the atmosphere
can evolve as a near-closed system.

For a magma ocean more reduced than IW, the modeled cases at $f\rm{O_{2}}$ of
$\Delta$IW=$-2$ and C/H=5 show that graphite precipitation would occur in a
cooling atmosphere, even accounting for the production of CH4 that buffers
$f$CO (Figure 13). Specifically, higher initial C budgets result in higher
$f$CO, where, for C/H=5, the scenario with 10 oceans crosses the CCO buffer at
$\sim$2300 K, while it does so at $\sim$2000 K for five oceans. As such, the
quoted fugacities of carbon-bearing species for these two scenarios are upper
limits. Cases at C/H=1 and $\Delta$IW=$-2$ do not result in graphite
saturation at any temperature, irrespective of the number of oceans, owing to
the lower $f$CO.

Figure 13: The fugacity of CO in equilibrium with graphite (Equation 11)
assuming ideal gas behavior as a function of temperature at $\Delta$IW=$-2$
(dashed black line) and $\Delta$IW=0.5 (dashed gray line). All cases at
$\Delta$IW=$-2$ and C/H=5 (colored lines) intercept the CCO buffer at
$\Delta$IW=$-2$, indicating graphite precipitation. Only cases that intercept
the CCO buffer are plotted.

We presume that the magma ocean is well mixed, such that it remains
isochemical throughout its depth. However, the saturation of graphite in
equilibrium with silicate melt depends on pressure and temperature (Dasgupta
et al., 2013; Stanley et al., 2014; Chi et al., 2014). As such, chemical
gradients could arise from its precipitation during magma ocean cooling and/or
at depth followed by its isolation from the melt (Hirschmann, 2012). To assess
the potential effect of graphite/diamond precipitation on the calculated
fugacities of carbon-bearing species, we compare the C contents dissolved in
the magma ocean predicted by our models with the graphite/diamond
precipitation curve. Above the IW buffer, C dissolves predominantly as
CO${}_{3}^{2-}$ in mafic and ultramafic melts (Armstrong et al., 2015; Duncan
et al., 2017). Holloway et al. (1992) devised an expression to determine
CO${}_{3}^{2-}$ solubility at graphite saturation. To calculate whether
graphite is expected to precipitate from the silicate melt in our simulations,
we adopt the calibration of Duncan et al. (2017) extrapolated to peridotite
compositions. At constant pressure (1 GPa) and for the temperature range of
cooling considered in our models (1650 to 2500 K), little or no C
precipitation occurs in the magma ocean under any of the scenarios modeled at
IW+0.5 and above. At 1650 K and 1 GPa, the model of Duncan et al. (2017)
predicts that $\sim$350 ppmw CO2 is dissolved in a peridotite melt or 100 ppmw
for a komatiitic melt (expected for an evolved magma ocean) at graphite
saturation, compared to $\sim$300 ppmw in the melt in the most extreme, C/H=5,
10-ocean case.

Below the IW buffer, CO${}_{3}^{2-}$ is no longer the prevailing melt species
of dissolved carbon (instead, carbon likely occurs as some CO-bearing
molecule, Wetzel et al., 2013; Dalou et al., 2019); therefore other
expressions for the prediction of C solubility at graphite saturation need to
be considered. Here we adopt the calibration of Yoshioka et al. (2019) that
links CO dissolved in the melt to its fugacity, $f\rm{CO}$. We fix $f\rm{CO}$
as a function of pressure in equilibrium with graphite at IW$-2$ using
thermodynamic data and the equation of state of Jakobsson & Oskarsson (1994).
At 2 GPa, this yields an $f\rm{CO}$ of $\sim$8.3 GPa, which, using the
expression of Yoshioka et al. (2019), leads to 500 ppmw of dissolved CO. These
estimates far exceed those derived by solely considering CO${}_{3}^{2-}$ as
the dissolved species ($\sim$2 ppmw, Duncan et al., 2017). Given that
dissolved C contents reach 20 ppmw at most in our simulations at IW$-2$, this
implies that graphite and diamond should not precipitate at any stage in the
magma ocean, even under reducing conditions.

Another important phase change can occur in cooling atmospheres, namely, that
between liquid water and steam:

${\rm H_{2}O(l)}={\rm H_{2}O(g)}.$ (12)

Its equilibrium constant as a function of temperature indicates that, at 1400
K, H2O condenses at about 600-bar $f$H2O. This phase transformation is
neglected in our simulations primarily because these pressures and
temperatures exceed the critical point of water, which is 650 K and 220 bars.
Hence at 600-bar $f$H2O the atmosphere is a supercritical fluid rather than an
ideal gas as modeled herein. Moreover, the mixing of CO2 and H2O is nonideal,
such that departures from ideal gas behavior are expected in these solutions
(Duan et al., 1992; Frost & Wood, 1997).

###  4.5 Meteoritic Degassing

Theoretical calculations (Schaefer & Fegley, 2010) and outgassing experiments
(Thompson et al., 2021) predict the composition of gases evolved from the
degassing of meteoritic materials, some of which may be representative of
planetary building blocks. By examining equilibrium between gas and condensed
phases, thermodynamic models find that H2 and CO gases form from reduced
chondrites (such as EH) at high temperature (Schaefer & Fegley, 2010). Only
for the most oxidizing materials (CI chondritic composition) is H2O predicted
to be a major constituent of the gas mixture. Since water is only an important
gaseous species during degassing of one particular class of chondrite, this
calls into question the applicability of steam atmospheres to the evolution of
terrestrial planets at large.

Although instructive, directly correlating gases produced from meteorites to
the composition of planetary atmospheres is complicated for several reasons.
Meteorites are imperfect analogs of the material that went on to form the
terrestrial planets (Sossi, 2021), particularly as regards volatile elements
due to the ubiquitous alteration processes occurring on their parent body
(Brearley, 2006), as well as thermal metamorphism and the likelihood that
planetesimals are already differentiated when they accrete. More importantly,
the likely presence of a deep magma ocean during and immediately after
accretion would have modified the budgets of elements available to degas at
the planetary surface (e.g., Grewal et al., 2020; Dalou et al., 2019). As
shown in our work, equilibrium reactions between silicate melt and volatiles
in the gas phase influence the partitioning of gas species between the
atmosphere and magma ocean according to their solubility. Furthermore, some
volatiles (including H and C) may be permanently sequestered into an iron core
owing to their siderophile nature (e.g., Tagawa et al., 2021).

Therefore, connecting gas mixtures to atmospheric composition necessitates an
assumption that the starting bulk composition—usually derived from
chondrites–adequately captures elemental abundances in the atmosphere during
and after the magma ocean phase. This assumes that volatiles did not
experience significant partitioning due to magma ocean and core formation
processes and that atmosphere formation proceeds by the breakdown of solid
mineral assemblages (i.e., degassing occurs subsolidus). Solid bodies cannot
efficiently replenish surface material from deeper regions, whereas for molten
bodies the entire melt mass can potentially communicate with the atmosphere,
subject to the interior–atmosphere equilibration timescale (e.g., Pahlevan et
al., 2019). Alternatively, partitioning can be permitted if ongoing impact-
induced outgassing of chondritic materials after a surface lid or quench crust
forms (e.g. in a late-veneer scenario or during the late heavy bombardment) is
sufficient to displace or dilute any previously equilibrated atmosphere
(Zahnle et al., 2020).

Extraction of volatiles to solids in the core or mantle locks them away from
outgassing to the atmosphere over a short timescale, and perhaps indefinitely.
By contrast, volatiles partitioned in melt are readily available to outgas
during mantle crystallization as long as the interior and atmosphere remain in
thermodynamic communication (e.g., for fractional crystallization). Accounting
for thermodynamic equilibrium with silicate melt, we demonstrate that the
earliest atmosphere evolves in both mass and composition during magma ocean
cooling. An atmosphere can transition from dominantly CO-rich (”reducing”) to
H2O-rich (”oxidizing”) owing to solubility and redox reactions as the magma
ocean cools (Figure 12). Therefore, the solubility of volatiles in silicate
melt prevents drawing a simple connection between the composition of planetary
building blocks and the earliest atmospheres of terrestrial planets.

##  5 Summary and Conclusions

The investigation of catastrophic outgassing of CO2 and H2O from terrestrial
magma oceans—as well as their influence on the nature of early atmospheres
around rocky planets—was originally motivated by the degassing of hydrated
minerals and oxidized chondritic materials. However, for redox conditions
appropriate for the early Earth, we find that CO is usually the dominant
atmospheric species during the early molten stage, where the melt fraction is
greater than 30% (Figure. 12). This is due to the relatively high solubility
of hydrogen species (H2O) compared to carbon species (CO2), as well as redox
reactions that govern CO–CO2, H2–H2O, and CO2–H2–CH4 equilibria. Only when C/H
by weight is small (C/H$<$0.1) and the H budget large ($>$6 oceans) can H
species constitute more than 50% of the atmosphere, with H2 and H2O in
approximately equal abundance (Figure 4). For more oxidizing conditions, early
atmospheres are dominated by CO2 with H2O a minor species (Figure 8).

For the late molten stage ($<30\%$ melt), more H is retained (relative to the
total budget) for smaller H budgets, larger C/H, and more oxidizing conditions
(Figure 12). At this stage, continued outgassing driven by dissolution
equilibrium requires that the surface remains molten and a persistent lid does
not form; this requirement can be satisfied if a magma ocean undergoes
fractional crystallization. However, for magma oceans that freeze more quickly
by equilibrium crystallization, formation of a surface lid around 30% melt can
prevent subsequent outgassing of highly soluble volatiles such as H2O. This
would leave the atmosphere as carbon/CO dominated by stifling the formation of
a steam atmosphere. Methane only forms for reduced conditions
($\Delta$IW$\lesssim$0.5) when pressures are sufficiently high or temperatures
sufficiently low; otherwise, CH4 is absent. Graphite precipitation from the
atmosphere during magma ocean cooling is expected only for very high C budgets
(higher than those anticipated for Earth) and for $\Delta$IW$\lesssim$0.5.
Nevertheless, for all cases, graphite in the atmosphere may precipitate with
continued cooling.

Complete outgassing of volatiles during magma ocean solidification can arise
as a result of fractional crystallisation. However, by additionally
considering equilibrium crystallization and other processes, we expect that a
substantial reservoir of H could remain in planetary mantles and outgas over
geological timescales, potentially impacting the depth of surface oceans and
the ability to desiccate planets. In this case, hydrogen species would only
play a minor role in determining atmospheric opacity or modulating atmospheric
escape during the early magma ocean; rather, the behavior of carbon species is
dominant. The style of magma ocean crystallization—fractional or
equilibrium—therefore controls the composition and mass of early atmospheres
and the efficiency of volatile delivery to the planetary atmosphere and
surface. Ultimately, the high solubility of H2O in magma oceans may enable its
safe storage during the tumultuous phase of planet formation.

##  6 Data Availability

All of the data generated as part of this study can be obtained by contacting
Dan J. Bower.

DJB acknowledges Swiss National Science Foundation (SNSF) Ambizione Grant
173992. KH is supported by the European Research Council via Consolidator
Grant ERC-2017-CoG-771620-EXOKLEIN. PAS acknowledges Swiss National Science
Foundation (SNSF) Ambizione Grant 180025. PS acknowledges financial support
from the Swiss University Conference and the Swiss Council of Federal
Institutes of Technology through the Platform for Advanced Scientific
Computing (PASC) program. This research was partly inspired by discussions and
interactions within the framework of the National Center for Competence in
Research (NCCR) PlanetS supported by the SNSF. The calculations were performed
on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of
Bern. We thank S. Grimm and D. Kitzmann for computing Rosseland mean
opacities, and A. Wolf, M. Tian, and T. Lichtenberg for discussions and
feedback on the manuscript. Three anonymous referees provided comments that
further enhanced the manuscript.

##  Appendix A Opacities

We computed Rosseland mean opacities at 1 bar between 1700 and 2700 K for H2O,
H2, CO2, CO, and CH4 using Helios-k (Grimm et al., 2021; Grimm & Heng, 2015).
The opacity (Rosseland mean mass absorption coefficient) for H2O using the
full BT2 list (Barber et al., 2006) is $1.5\times 10^{-2}$ m2 kg-1 at 1700 K
and $5.0\times 10^{-4}$ m2 kg-1 at 2700 K; the opacity at 1800 K compares
favorably with $10^{-2}$ m2 kg-1 which has been used extensively in previous
magma ocean studies. For simplicity we also adopt a constant value (i.e.,
independent of temperature) of $10^{-2}$ m2 kg-1. The opacity of H2 is
determined by collision-induced absorption, and using HITRAN (Karman et al.,
2019; Abel et al., 2011) it is essentially constant ($5\times 10^{-5}$ m2
kg-1) from 2700 to 1700 K.

For CO2, HITEMP (Rothman et al., 2010) provided the line list, and we
determined an opacity of $1.5\times 10^{-4}$ m2 kg-1 at 1700 K and $3\times
10^{-5}$ m2 kg-1 at 2700 K. Hence the CO2 opacity at 1700 K is compatible with
$10^{-4}$ m2 kg-1, which is the lowest opacity considered for CO2 in previous
magma ocean models. To be consistent with our selection of the H2O opacity
around 1800 K, we similarly adopt a constant CO2 opacity of $10^{-4}$ m2 kg-1.
For CO we determined an opacity of $1.2\times 10^{-5}$ m2 kg-1 at 1700 K and
$3\times 10^{-6}$ m2 kg-1 at 2700 K using HITEMP (Li et al., 2015a). Again, we
selected an opacity of $10^{-5}$ m2 kg-1 based on the value at 1800 K. We
computed a CH4 opacity of $10^{-2}$ m2 kg-1 at 1800 K based on the line lists
of Yurchenko et al. (2013); Yurchenko & Tennyson (2014). The opacities of
species at 1 bar are reference opacities that are used to compute a pressure-
dependent opacity (Equation A22, Abe & Matsui, 1985).

##  Appendix B Melt–Solid Separation

The fluxes to describe energy transport in a planetary mantle (which is
molten, partially molten, or solid) are described in detail in Bower et al.
(2018); Abe (1995). However, we recap the physics underpinning melt–solid
separation owing to its importance for determining whether a mantle undergoes
fractional or equilibrium crystallization. At high melt fraction, the settling
or flotation of crystals within a magma ocean is determined by Stokes’s law:

$u_{m}-u_{s}=\frac{2a^{2}g(\rho_{m}-\rho_{s})}{9\eta_{m}},$ (B1)

where $u$ is velocity (positive radially outward), $a$ crystal size, $g$
gravity (negative), $\rho$ density, $\eta$ viscosity, and subscripts $m$ and
$s$ denote melt and solid, respectively. Interaction amongst crystals is not
considered, and therefore free settling or flotation of crystals is an upper
limit on the efficiency of melt–solid separation at high melt fraction. At low
melt fraction, separation occurs via Darcy flow:

$u_{m}-u_{s}=\frac{k_{p}g(\rho_{m}-\rho_{s})}{\phi_{p}\eta_{m}},$ (B2)

where $k_{p}$ is permeability and $\phi_{p}$ porosity (synonymous with the
volume fraction of melt). The Rumpf–Gupte permeability law is appropriate at
intermediate porosity (Rumpf & Gupte, 1971),

$k_{rg}=\frac{a^{2}\phi_{p}^{5.5}}{1.4},$ (B3)

and the Blake–Kozeny–Carman permeability law at low porosity,

$k_{bkc}=\frac{a^{2}\phi_{p}^{3}}{1000(1-\phi_{p})^{2}}.$ (B4)

The constants in these permeability laws are constrained by experiments (see
McKenzie, 1984; Abe, 1995, for discussion). The flow laws can be represented
by a single description (Abe, 1995):

$u_{m}-u_{s}=\frac{a^{2}g(\rho_{m}-\rho_{s})}{\eta_{m}}f,$ (B5) $\displaystyle
f=$ $\displaystyle\dfrac{2}{9}$ $\text{for }0.77\leq\phi_{p}$, (B6a)
$\displaystyle f=$ $\displaystyle\dfrac{k_{rg}}{a^{2}\phi_{p}}$ $\text{for
}0.08<\phi_{p}<0.77$, (B6b) $\displaystyle f=$
$\displaystyle\dfrac{k_{bkc}}{a^{2}\phi_{p}}$ $\text{for }\phi_{p}\leq 0.08$,
(B6c)

where the range of applicability of each flow law ensures that $f$ is a
continuous function of porosity. Through consideration of the local
barycentric velocity, the mass flux of melt is

$J_{m}=\rho\phi(1-\phi)(u_{m}-u_{s}),$ (B7)

where $\phi$ is melt fraction (the mass fraction of melt). Equation B7
naturally satisfies the requirement that $J_{m}=0$ for $\phi=0$ and $\phi=1$.
The energy flux of melt–solid separation is therefore

$F_{\rm grav}=J_{m}T_{\rm fus}\Delta S_{\rm fus},$ (B8)

where $T_{\rm fus}$ is the temperature at 50% melt fraction and $\Delta S_{\rm
fus}$ is the entropy of fusion, both of which are functions of pressure
according to the melting curves. Equation B5 reveals the common factors
appearing in all of the flow laws and emphasizes the importance of crystal
size owing to its squared dependence: larger crystals enhance melt–solid
separation. Since we consider a single-component mantle (MgSiO3), everywhere
$\rho_{m}-\rho_{s}<0$ and hence crystals always sink and melt always drains
upward to the surface.

##  Appendix C Magma ocean and atmosphere equilibrium

###  C.1 Volatile mass balance and evolution

We extended Bower et al. (2019) to express the mass balance for a volatile $v$
in a chemically reactive environment with atmospheric escape:

$\displaystyle X_{v}k_{v}M^{s}+X_{v}M^{l}+4\pi
R_{p}^{2}\left(\frac{\mu_{v}}{\bar{\mu}}\right)\frac{p_{v}}{|g|}+m_{v}^{e}+\sum_{w=1}^{r}m_{v}^{w}$
$\displaystyle=X_{v}^{0}M^{m},$ (C1)

where:

  * •

$X_{v}$ is the volatile abundance in the pure melt;

  * •

$k_{v}$ is the distribution coefficient between solid and melt;

  * •

$M^{s}$ is the mantle mass of solid;

  * •

$M^{l}$ is the mantle mass of melt;

  * •

$M^{m}=M^{s}+M^{l}$ is the total mantle mass, which is constant;

  * •

$R_{p}$ is the planetary surface radius;

  * •

$g$ is the surface gravity;

  * •

$\mu_{v}$ is the molar mass of the volatile;

  * •

$\bar{\mu}$ is the mean molar mass of the atmosphere;

  * •

$p_{v}$ is the surface partial pressure of the volatile; and

  * •

$X_{v}^{0}$ is the initial total volatile abundance relative to the total
mantle mass.

A solubility law relates the volatile concentration in the melt ($X_{v}$) to
the partial pressure of the volatile in the atmosphere ($p_{v}$). The first
three terms on the left-hand side of Equation C1 represent the mass of the
volatile stored in the solid mantle, the molten mantle, and the atmosphere.
The final two terms on the left-hand side are fictitious reservoirs used to
account for mass loss due to escape ($m_{v}^{e}$; Appendix C.2), and mass
transfer due to participation of this volatile in reaction $w$ ($m_{v}^{w}$,
Appendix C.3). The time derivative of Equation C1 derives the evolution
equation (cf., Equation A4, Bower et al., 2019):

$\displaystyle\frac{dX_{v}}{dt}\left(k_{v}M^{m}+(1-k_{v})M^{l}\right)+X_{v}(1-k_{v})\frac{dM^{l}}{dt}$
$\displaystyle+\frac{4\pi
R_{p}^{2}\mu_{v}}{|g|}\left[\frac{1}{\bar{\mu}}\frac{dp_{v}}{dt}-\frac{p_{v}}{\bar{\mu}^{2}}\sum_{w=1}^{n}\frac{\mu_{w}}{P_{s}}\left(\frac{dp_{w}}{dt}-\frac{p_{w}}{P_{s}}\frac{dP_{s}}{dt}\right)\right]$
$\displaystyle+\frac{dm_{v}^{e}}{dt}+\sum_{w=1}^{r}\frac{dm_{v}^{w}}{dt}=0,$
(C2)

where $P_{s}$ is the total surface pressure, which is equal to the summation
of the partial pressures of all volatiles according to Dalton’s law. The
derivative of $X_{v}$ with respect to $t$ can be determined from the
solubility law relevant for volatile $v$, which for power-law solubility
(Equation 1) is

$\frac{dX_{v}}{dt}=\frac{dX_{v}}{dp_{v}}\frac{dp_{v}}{dt}=\left(\frac{\alpha_{v}{p_{v}}^{1/\beta_{v}-1}}{\beta_{v}}\right)\frac{dp_{v}}{dt}.$
(C3)

Similarly, $X_{v}$ can be eliminated from Equation C2 using the solubility
law. This results in a system of nonlinear equations that can be solved
simultaneously to determine $dp_{v}/dt$ for each volatile $v$. Solving for
$dp_{v}/dt$ rather than $dX_{v}/dt$ enables us to treat the case of zero
solubility. Oxygen fugacity $f{\rm O_{2}}$ (Equation 7) depends on surface
temperature, so oxygen is not explicitly tracked by mass balance as with the
other volatiles.

###  C.2 Mass loss due to atmospheric escape

The escape flux of H2 (Equation 9) is used to determine a mass-loss rate of H2
(e.g., Katyal et al., 2020):

$\frac{dm_{H_{2}}^{e}}{dt}=4\pi R_{p}^{2}\phi_{\rm H_{2}}\mu_{\rm H_{2}}.$
(C4)

The loss rate depends on the time-dependent volume mixing flux of H2 through
$\phi_{\rm H_{2}}$. In turn, the volume mixing ratio of H2 is computed from
the partial pressure of H2 and the total surface pressure. In the present
work, the mass-loss rate for volatiles other than H2 is set to zero.

###  C.3 Mass exchange due to chemical reactions

During the magma ocean stage, the silicate melt and overlying atmosphere can
accommodate mass exchange via chemical reactions. We consider a general
reaction with two reactants ($A$, $B$) and two products ($C$, $D$):

$aA+bB\rightleftharpoons cC+dD,$ (C5)

where $a$, $b$, $c$, and $d$ are stoichiometric coefficients. The reaction
quotient $Q$ is a function of time $t$:

$Q(t)=\frac{\\{C\\}^{c}\\{D\\}^{d}}{\\{A\\}^{a}\\{B\\}^{b}}\doteq\frac{Q_{p}(t)}{Q_{r}(t)},$
(C6)

where curly brackets $\\{\\}$ denote the activity of the species and the
activity of a gas is its fugacity divided by a reference pressure. At chemical
equilibrium, the reaction quotient is equal to the equilibrium constant $K$,
which is a function of temperature. An appropriate temperature is that of the
magma ocean–atmosphere interface; hence,

$\frac{\\{C\\}^{c}\\{D\\}^{d}}{\\{A\\}^{a}\\{B\\}^{b}}-K(T)=0.$ (C7)

The reaction in Equation C5 can also be expressed in terms of a mass balance
of reactants and products:

$m_{A}+m_{B}\rightleftharpoons m_{C}+m_{D}.$ (C8)

This is because the fractional change of moles corresponds to the fractional
change of mass, so correlating terms in Equation. C5 and C8:

$m_{B}=\left(\frac{b\mu_{B}}{a\mu_{A}}\right)m_{A},\quad
m_{C}=\left(\frac{c\mu_{C}}{a\mu_{A}}\right)m_{A},\quad{\rm etc.}$ (C9)

where $\mu$ is the molar mass. It is convenient in Equation C9 to express the
mass of the reactants and products using the mass of the first reactant
($m_{A}$) multiplied by a constant factor that depends only on molar masses
and reaction stoichiometry; mathematically this choice is arbitrary, and
another mass could equally be used. Hence, each chemical reaction introduces a
single unknown ”reaction mass” ($m_{A}$) and one equation (Equation C7) to the
complete system of equations that we solve. The functional dependencies of all
quantities in Equation C7 are known, since we assume ideal gas behavior to
connect gas activities to volatile partial pressures and the functional
dependence of the equilibrium constant on temperature $T$ is documented in
thermochemical tables (Chase, 1998).

Evolution equations of Equations C7 and C9 are necessary to time-step the
model. The time derivatives of the reactant and product masses are trivial to
compute since stoichiometric coefficients and molar masses are constant. For
example (c.f., Equation C9):

$\frac{dm_{B}}{dt}=\left(\frac{b\mu_{B}}{a\mu_{A}}\right)\frac{dm_{A}}{dt},$
(C10)

and the time derivative of the reactant masses ($m_{C}$, $m_{D}$) can be
similarly computed. The time derivative of the equilibrium condition (Equation
C7) is:

$\frac{1}{Q_{r}}\frac{dQ_{p}}{dt}-\frac{Q_{p}}{Q_{r}^{2}}\frac{dQ_{r}}{dt}-\frac{dK}{dT}\frac{dT}{dt}=0,$
(C11)

where the reaction quotient numerator $Q_{p}=\\{C\\}^{c}\\{D\\}^{d}$ and
denominator $Q_{r}=\\{A\\}^{a}\\{B\\}^{b}$ depend on the volatile partial
pressures and reaction stoichiometry. The equilibrium constant $K$ depends
only on temperature $T$ (Equations 4 and 5) and therefore $dK/dT$ is
straightforward to evaluate. The code also solves for $dT/dt$ in the mantle at
each staggered mesh point, although only the change of surface temperature is
required since this temperature defines the interface of the magma ocean and
atmosphere.

Hence, to include reactions in the time stepper, Equation C11 constitutes an
extra equation per reaction that must be satisfied, where the extra unknowns
we solve for are the change of reaction mass for each reaction ($dm_{A}/dt$).
Equation C10 relates the change of reaction mass for reactant $B$ and products
$C$ and $D$ to the solution quantity ($dm_{A}/dt$). For each volatile, we can
thus evaluate the total rate of change of its mass due to all reactions that
it participates in; this is the final term on the left-hand side of Equation
C2.

###  C.4 Initial volatile abundances

For the initial condition, we solve for the abundance of each volatile in the
melt using a set of coupled equations consisting of all the volatile mass
balances (Equation C1) and $r$ chemical reactions (Equation C7). To do this,
estimates of the initial total volatile abundances $X_{v}^{\rm 0,est}$ are
prescribed, and the initial atmospheric escape masses are set to zero
($m_{v}^{e}=0$). We solve for the partial pressure of each volatile in the
atmosphere and the reaction masses, the latter of which are likely to be
nonzero. Hence, it is convenient—although formally not necessary—to offset the
estimates of the initial total volatile abundances ($X_{v}^{\rm 0,est}$) using
the solved reaction masses to enforce that the initial total volatile
abundances ($X_{v}^{0}$) adhere to chemical equilibrium constraints with zero
reaction mass:

$X_{v}^{0}=X_{v}^{\rm 0,est}-\frac{1}{M^{m}}\sum_{w}^{r}m_{v}^{w}.$ (C12)

Hence, the reaction masses are zero by definition for the initial condition
and can increase or decrease during the time integration.

###  C.5 Alternative using thermodynamic components

The system of volatiles we consider can be represented by two components: C
and H2, where we recall that oxygen is not explicitly tracked. Hence, an
alternative approach to solving the mass balance for each volatile is to
instead solve for the total abundance of C and H2 moles. Conservation laws for
C and H2 abundances can be derived by converting the relevant volatile mass
balances to C or H2 moles and then summing. Without chemical reactions,
conserving species mass or abundance is the same because there is no exchange
of material between species. With reactions, either one can continue to track
species mass by accounting for mass exchange by reactions (as we choose to
do), or alternatively one can instead track the abundance of components.
Considering species mass allows one to retain the familiar mass formulation
derived without reactions, but at the expense of more equations (e.g., five
volatile species plus three reactions). By contrast, considering components
reduces the number of equations (e.g., two components plus three reactions),
but requires a departure from the traditional species mass balance approach.
Ultimately, both approaches are valid and equivalent and will likely find
utility in different applications.

###  C.6 Simplified mass balance and interior depletion

The partitioning of volatiles between a fully molten magma ocean and
atmosphere can be analyzed using a reduced form of the volatile mass balance
(Equation C1). This aids the interpretation of the results—notably
demonstrating their universality—as well as providing analytical expressions
and scaling relations. We consider the mass balance for a single volatile that
is only partitioned between the melt and atmosphere. However, we retain the
mean molar mass of the atmosphere $\bar{\mu}$ to provide insight when the
atmosphere is dominated by the presence of another (or mixture) of volatiles
that are less soluble than the volatile under consideration. Following
Equation C1:

$X_{v}M^{l}+4\pi
R_{p}^{2}\left(\frac{\mu_{v}}{\bar{\mu}}\right)\frac{p_{v}}{|g|}=X_{v}^{0}M^{m}.$
(C13)

For power-law solubility (Equation 1),

$\frac{\phi_{g}}{C}X_{v}+{X_{v}}^{\beta_{v}}=\frac{X_{v}^{0}}{C},$ (C14)

where $\phi_{g}=M^{l}/M^{m}$ is the global melt fraction and

$C=\frac{4\pi R_{p}^{2}\mu_{v}{\alpha_{v}}^{-\beta_{v}}}{M^{m}\bar{\mu}|g|}.$
(C15)

For a single species ($\bar{\mu}=\mu_{v}$), an Earth-sized mantle, and values
of $\alpha_{v}$ and $\beta_{v}$ from the peridotite solubility law (Table 1),
the coefficient $C\approx 17/4$.

####  C.6.1 Volatile depletion for $\beta_{v}=2$

For the peridotite solubility law (as well as other laws), $\beta_{v}=2$ and
Equation C14 is a quadratic equation that can be solved for the volatile
abundance in the melt:

$X_{v}=\frac{-\phi_{g}}{2C}+\frac{1}{2}\sqrt{\frac{\phi_{g}^{2}}{C^{2}}+\frac{4X_{v}^{0}}{C}}.$
(C16)

Hence, comparing to Equation C14, the scaled volatile mass in the liquid is

$m_{v}^{l}=\frac{-\phi_{g}^{2}}{2C^{2}}+\frac{\phi_{g}}{2C}\sqrt{\frac{\phi_{g}^{2}}{C^{2}}+\frac{4X_{v}^{0}}{C}},$
(C17)

the scaled volatile mass in the atmosphere is

$m_{v}^{a}=\frac{X_{v}^{0}}{C}+\frac{\phi_{g}^{2}}{2C^{2}}-\frac{\phi_{g}}{2C}\sqrt{\frac{\phi_{g}^{2}}{C^{2}}+\frac{4X_{v}^{0}}{C}},$
(C18)

and the scaled total volatile mass is

$m_{v}^{t}=\frac{X_{v}^{0}}{C}.$ (C19)

Therefore, the volatile depletion $\mathcal{D}$ expressed as a fraction
relative to the total volatile mass is

$\mathcal{D}=1-\frac{m_{v}^{l}}{m_{v}^{t}}=\frac{m_{v}^{a}}{m_{v}^{t}}=1+\frac{\phi_{g}^{2}}{2CX_{v}^{0}}-\frac{\phi_{g}}{2X_{v}^{0}}\sqrt{\frac{\phi_{g}^{2}}{C^{2}}+\frac{4X_{v}^{0}}{C}}.$
(C20)

Taking the derivative of Equation C20 with respect to the initial volatile
abundance $X_{v}^{0}$, we obtain

$\frac{\partial\mathcal{D}}{\partial
X_{v}^{0}}=\frac{\phi_{g}\left(-\phi_{g}\sqrt{\phi_{g}^{2}+4CX_{v}^{0}}+\phi_{g}^{2}+2CX_{v}^{0}\right)}{2(X_{v}^{0})^{2}C\sqrt{\phi_{g}^{2}+4X_{v}^{0}C}}.$
(C21)

All variables must be positive to be physically meaningful. The denominator is
clearly always positive and the numerator is also positive since the following
condition is always satisfied:

$\phi_{g}^{2}+2CX_{v}^{0}>\phi_{g}\sqrt{\phi_{g}^{2}+4CX_{v}^{0}}\,.$ (C22)

Therefore,

$\frac{\mathcal{\partial D}}{\partial X_{v}^{0}}>0\,.$ (C23)

Hence, at a given melt fraction $\phi_{g}$, the relative depletion of the
interior increases with the initial volatile abundance $X_{v}^{0}$ (Figure 3).

####  C.6.2 Volatile depletion for $\beta_{v}=1$

For a volatile that obeys Henrian solubility behavior ($\beta_{v}=1$), similar
to Appendix C.6.1, we can write:

$X_{v}=\frac{X_{v}^{0}}{C+\phi_{g}},\quad
m_{v}^{l}=\frac{\phi_{g}X_{v}^{0}}{C(C+\phi_{g})},\quad
m_{v}^{a}=\frac{X_{v}^{0}}{C+\phi_{g}}.$ (C24)

Hence, the interior depletion $\mathcal{D}$ is

$\mathcal{D}=\frac{C}{C+\phi_{g}},\qquad\frac{\partial\mathcal{D}}{\partial
X_{v}^{0}}=0.$ (C25)

Notably, interior depletion is independent of the initial volatile abundance
$X_{v}^{0}$ when $\beta_{v}=1$.

##  Appendix D Oxygen fugacity and atmosphere speciation

Figure 4 shows the atmospheric composition at an oxygen fugacity of
$\Delta$IW=0.5. The following figures show the same at $\Delta$IW=$-2$ (Figure
14), $\Delta$IW=2 (Figure 15), and $\Delta$IW=4 (Figure 16).

Figure 14: Atmospheric composition during magma ocean outgassing at $f{\rm
O_{2}}$ of $\Delta$IW=$-2$. See Figure 4 caption. Figure 15: Atmospheric
composition during magma ocean outgassing at $f{\rm O_{2}}$ of $\Delta$IW=2.
See Figure 4 caption. Figure 16: Atmospheric composition during magma ocean
outgassing at $f{\rm O_{2}}$ of $\Delta$IW=4. See Figure 4 caption.

##  Appendix E Comparison with FactSage 8.0

We independently verified the results of our model by comparison to FactSage
calculations, assuming ideal gas behavior (Table E).

Table E: Comparison of Select Outgassed Atmospheres Calculated Using Our Model
Compared to FactSage Calculations

Case | Our model (basr) | FactSage (bars)
---|---|---
H Oceans | C/H |  $f\rm O_{2}$ ($\Delta$IW) | CO | CO2 | H2 | H2O | CH4 | CO | CO2 | H2 | H2O | CH4
3 | 1 | -2.0 | 7.1 | 0.4 | 186 | 23 | 42 | 7.2 | 0.4 | 187.5 | 22.4 | 41
3 | 1 | +0.5 | 51 | 48 | 213 | 460 | 22 | 48.9 | 49.7 | 216.1 | 458.6 | 20.7
1 | 0.1 | +2.0 | 0.7 | 3.7 | 20 | 240 | 0 | 0.6 | 3.7 | 20.2 | 239.9 | 0.00043
1 | 5 | +4.0 | 7 | 368 | 3.7 | 444 | 0 | 6.4 | 367.6 | 3.7 | 445 | 0

Note. — This comparison reveals that the partial pressure of a given volatile
differs by at most a few percent.

## References

  * Abe (1993) Abe, Y. 1993, Thermal Evolution and Chemical Differentiation of the Terrestrial Magma Ocean, ed. E. Takahashi, R. Jeanloz, & D. Rubie, Vol. 74 (AGU, Washington, D.C.), 41–54, doi: https://doi.org/10.1029/GM074p0041
  * Abe (1995) —. 1995, in The Earth’s Central Part: Its Structure and Dynamics, ed. T. Yukutake (Terra Sci. Pub. Com., Tokyo), 215–230.  https://www.terrapub.co.jp/e-library/ecp/pdf/EC0215.PDF
  * Abe & Matsui (1985) Abe, Y., & Matsui, T. 1985, JGRB, 90, C545, doi: https://doi.org/10.1029/JB090iS02p0C545
  * Abel et al. (2011) Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2011, JPCA, 115, 6805, doi: https://doi.org/10.1021/jp109441f
  * Andrault et al. (2011) Andrault, D., Bolfan-Casanova, N., Nigro, G. L., et al. 2011, E&PSL, 304, 251, doi: https://doi.org/10.1016/j.epsl.2011.02.006
  * Andrault et al. (2016) Andrault, D., Monteux, J., Bars, M. L., & Samuel, H. 2016, E&PSL, 443, 195, doi: https://doi.org/10.1016/j.epsl.2016.03.020
  * Ardia et al. (2013) Ardia, P., Hirschmann, M., Withers, A., & Stanley, B. 2013, GeCoA, 114, 52, doi: https://doi.org/10.1016/j.gca.2013.03.028
  * Armstrong et al. (2019) Armstrong, K., Frost, D. J., McCammon, C. A., Rubie, D. C., & Boffa Ballaran, T. 2019, Sci, 365, 903, doi: https://doi.org/10.1126/science.aax8376
  * Armstrong et al. (2015) Armstrong, L. S., Hirschmann, M. M., Stanley, B. D., Falksen, E. G., & Jacobsen, S. D. 2015, GeCoA, 171, 283, doi: https://doi.org/10.1016/j.gca.2015.07.007
  * Bale et al. (2009) Bale, C., Bélisle, E., Chartrand, P., et al. 2009, Calphad: Computer Coupling of Phase Diagrams and Thermochemistry, 33, 295, doi: https://doi.org/10.1016/j.calphad.2008.09.009
  * Barber et al. (2006) Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087, doi: https://doi.org/10.1111/j.1365-2966.2006.10184.x
  * Barth et al. (2021) Barth, P., Carone, L., Barnes, R., et al. 2021, AsBio, doi: https://doi.org/10.1089/ast.2020.2277
  * Benedikt et al. (2020) Benedikt, M., Scherf, M., Lammer, H., et al. 2020, Icarus, 347, 113772, doi: https://doi.org/10.1016/j.icarus.2020.113772
  * Berndt et al. (2002) Berndt, J., Liebske, C., Holtz, F., et al. 2002, AmMin, 87, 1717, doi: https://doi.org/10.2138/am-2002-11-1222
  * Bonati et al. (2019) Bonati, I., Lichtenberg, T., Bower, D. J., Timpe, M., & Quanz, S. 2019, A&A, 621, doi: https://doi.org/10.1051/0004-6361/201833158
  * Bower et al. (2019) Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, A&A, 631, doi: https://doi.org/10.1051/0004-6361/201935710
  * Bower et al. (2018) Bower, D. J., Sanan, P., & Wolf, A. S. 2018, PEPI, 274, 49, doi: https://doi.org/10.1016/j.pepi.2017.11.004
  * Bower et al. (2021) —. 2021, SPIDER: Simulating Planetary Interior Dynamics with Extreme Rheology, Version 0.2.1, Zenodo, doi: https://doi.org/10.5281/zenodo.5682523
  * Brearley (2006) Brearley, A. J. 2006, The Action of Water, ed. D. S. Lauretta & H. Y. McSween (University of Arizona Press), 587–624
  * Caracas et al. (2019) Caracas, R., Hirose, K., Nomura, R., & Ballmer, M. D. 2019, E&PSL, 516, 202, doi: https://doi.org/10.1016/j.epsl.2019.03.031
  * Charnoz et al. (2021) Charnoz, S., Sossi, P. A., Lee, Y.-N., et al. 2021, Icarus, 364, 114451, doi: https://doi.org/10.1016/j.icarus.2021.114451
  * Chase (1998) Chase, M. 1998, NIST-JANAF Thermochemical Tables, 4th Edition (American Institute of Physics), doi: https://doi.org/10.18434/T42S31
  * Chi et al. (2014) Chi, H., Dasgupta, R., Duncan, M. S., & Shimizu, N. 2014, GeCoA, 139, 447, doi: https://doi.org/10.1016/j.gca.2014.04.046
  * Corgne et al. (2008) Corgne, A., Wood, B. J., & Fei, Y. 2008, GeCoA, 72, 2409, doi: https://doi.org/10.1016/j.gca.2008.03.001
  * Costa et al. (2009) Costa, A., Caricchi, L., & Bagdassarov, N. 2009, GGG, 10, doi: https://doi.org/10.1029/2008GC002138
  * Dalou et al. (2019) Dalou, C., Hirschmann, M. M., Jacobsen, S. D., & Le Losq, C. 2019, GeCoA, 265, 32, doi: https://doi.org/10.1016/j.gca.2019.08.029
  * Dasgupta (2013) Dasgupta, R. 2013, RvMG, 75, 183, doi: https://doi.org/10.2138/rmg.2013.75.7
  * Dasgupta et al. (2013) Dasgupta, R., Chi, H., Shimizu, N., Buono, A. S., & Walker, D. 2013, GeCoA, 102, 191, doi: https://doi.org/10.1016/j.gca.2012.10.011
  * Dixon et al. (1995) Dixon, J. E., Stolper, E. M., & Holloway, J. R. 1995, JPet, 36, 1607, doi: https://doi.org/10.1093/oxfordjournals.petrology.a037267
  * Duan et al. (1992) Duan, Z., Møller, N., & Weare, J. H. 1992, GeCoA, 56, 2605, doi: https://doi.org/10.1016/0016-7037(92)90347-L
  * Duncan et al. (2017) Duncan, M. S., Dasgupta, R., & Tsuno, K. 2017, E&PSL, 466, 115, doi: https://doi.org/10.1016/j.epsl.2017.03.008
  * Elkins-Tanton (2008) Elkins-Tanton, L. 2008, E&PSL, 271, 181, doi: https://doi.org/10.1016/j.epsl.2008.03.062
  * Fei & Katsura (2020) Fei, H., & Katsura, T. 2020, E&PSL, 531, 115987, doi: https://doi.org/10.1016/j.epsl.2019.115987
  * Frost & Wood (1997) Frost, D. J., & Wood, B. J. 1997, GeCoA, 61, 3301, doi: https://doi.org/10.1016/S0016-7037(97)00168-3
  * Gaillard et al. (2021) Gaillard, F., Bouhifd, M. A., Füri, E., et al. 2021, SSRv, 217, 22, doi: https://doi.org/10.1007/s11214-021-00802-1
  * Gooding & Muenow (1977) Gooding, J. L., & Muenow, D. W. 1977, Metic, 12, 401, doi: https://doi.org/10.1111/j.1945-5100.1977.tb00458.x
  * Grewal et al. (2020) Grewal, D. S., Dasgupta, R., & Farnell, A. 2020, GeCoA, 280, 281, doi: https://doi.org/10.1016/j.gca.2020.04.023
  * Grimm & Heng (2015) Grimm, S. L., & Heng, K. 2015, ApJ, 808, 182, doi: https://doi.org/10.1088/0004-637X/808/2/182
  * Grimm et al. (2021) Grimm, S. L., Malik, M., Kitzmann, D., et al. 2021, ApJS, 253, 30, doi: https://doi.org/10.3847/1538-4365/abd773
  * Gronoff et al. (2020) Gronoff, G., Arras, P., Baraka, S. M., et al. 2020, JGRA, 125, e2019JA027639, doi: https://doi.org/10.1029/2019JA027639
  * Hamano et al. (2013) Hamano, K., Abe, Y., & Genda, H. 2013, Natur, 497, 607, doi: https://doi.org/10.1038/nature12163
  * Hamano et al. (2015) Hamano, K., Kawahara, H., Abe, Y., Onishi, M., & Hashimoto, G. L. 2015, ApJ, 806, 216, doi: https://doi.org/10.1088/0004-637X/806/2/216
  * Hamilton & Oxtoby (1986) Hamilton, D., & Oxtoby, S. 1986, JG, 94, 626
  * Hier-Majumder & Hirschmann (2017) Hier-Majumder, S., & Hirschmann, M. M. 2017, GGG, 18, doi: https://doi.org/10.1002/2017GC006937
  * Hirschmann et al. (2012) Hirschmann, M., Withers, A., Ardia, P., & Foley, N. 2012, E&PSL, 345–348, 38, doi: https://doi.org/10.1016/j.epsl.2012.06.031
  * Hirschmann (2012) Hirschmann, M. M. 2012, E&PSL, 341–344, 48, doi: https://doi.org/10.1016/j.epsl.2012.06.015
  * Hirschmann (2016) —. 2016, AmMin, 101, 540, doi: https://doi.org/10.2138/am-2016-5452
  * Hirschmann (2018) —. 2018, E&PSL, doi: https://doi.org/10.1016/j.epsl.2018.08.023
  * Holloway & Blank (1994) Holloway, J. R., & Blank, J. G. 1994, in Volatiles in magmas, ed. M. R. Carroll & J. R. Holloway, Vol. 30 (Mineralogical Society of America), 187–230, doi: https://doi.org/10.1515/9781501509674-012
  * Holloway et al. (1992) Holloway, J. R., Pan, V., & Gudmundsson, G. 1992, EJMin, 4, 105
  * Iacono-Marziano et al. (2012) Iacono-Marziano, G., Morizet, Y., Trong, E. L., & Gaillard, F. 2012, GeCoA, 97, 1, doi: https://doi.org/10.1016/j.gca.2012.08.035
  * Jakobsson & Oskarsson (1994) Jakobsson, S., & Oskarsson, N. 1994, GeCoA, 58, 9
  * Karman et al. (2019) Karman, T., Gordon, I. E., van der Avoird, A., et al. 2019, Icarus, 328, 160, doi: https://doi.org/10.1016/j.icarus.2019.02.034
  * Katyal et al. (2020) Katyal, N., Ortenzi, G., Grenfell, J. L., et al. 2020, A&A, 643, doi: https://doi.org/10.1051/0004-6361/202038779
  * Keppler & Golabek (2019) Keppler, H., & Golabek, G. 2019, Geochemical Perspectives Letters, 11, doi: https://doi.org/10.7185/geochemlet.1918
  * Labrosse et al. (2007) Labrosse, S., Hernlund, J. W., & Coltice, N. 2007, Natur, 450, 866, doi: https://doi.org/10.1038/nature06355
  * Lange & Ahrens (1982) Lange, M. A., & Ahrens, T. J. 1982, Icarus, 51, 96, doi: https://doi.org/10.1016/0019-1035(82)90031-8
  * Lebrun et al. (2013) Lebrun, T., Massol, H., Chassefière, E., et al. 2013, JGRE, 118, 1155, doi: https://doi.org/10.1002/jgre.20068
  * Lécuyer et al. (1998) Lécuyer, C., Gillet, P., & Robert, F. 1998, ChGeo, 145, 249, doi: https://doi.org/10.1016/S0009-2541(97)00146-0
  * Li et al. (2015a) Li, G., Gordon, I. E., Rothman, L. S., et al. 2015a, ApJS, 216, 15, doi: https://doi.org/10.1088/0067-0049/216/1/15
  * Li et al. (2015b) Li, Y., Dasgupta, R., & Tsuno, K. 2015b, E&PSL, 415, 54, doi: https://doi.org/10.1016/j.epsl.2015.01.017
  * Lichtenberg et al. (2021) Lichtenberg, T., Bower, D. J., Hammond, M., et al. 2021, JGRE, e2020JE006711, doi: https://doi.org/10.1029/2020JE006711
  * Maas & Hansen (2019) Maas, C., & Hansen, U. 2019, E&PSL, 513, 81, doi: https://doi.org/10.1016/j.epsl.2019.02.016
  * Marty (2012) Marty, B. 2012, E&PSL, 313-314, 56, doi: https://doi.org/10.1016/j.epsl.2011.10.040
  * Maurice et al. (2017) Maurice, M., Tosi, N., Samuel, H., et al. 2017, JGRE, 122, 577, doi: https://10.1002/2016JE005250
  * McKenzie (1984) McKenzie, D. 1984, JPet, 25, 713, doi: https://doi.org/10.1093/petrology/25.3.713
  * McMillan (1994) McMillan, P. F. 1994, in Volatiles in Magmas, ed. M. R. Carroll & J. R. Holloway, Vol. 30 (Mineralogical Society of America), 131–156, doi: https://doi.org/10.1515/9781501509674-010
  * Miyazaki & Korenaga (2021) Miyazaki, Y., & Korenaga, J. 2021, Does detecting water vapors on rocky planets indicate the presence of oceans?: An insight from self-consistent mantle degassing models.  https://arxiv.org/abs/2108.03759
  * Mosenfelder et al. (2009) Mosenfelder, J. L., Asimow, P. D., Frost, D. J., Rubie, D. C., & Ahrens, T. J. 2009, JGRB, 114, doi: https://doi.org/10.1029/2008JB005900
  * Mysen et al. (2011) Mysen, B. O., Kumamoto, K., Cody, G. D., & Fogel, M. L. 2011, GeCoA, 75, 6183, doi: 10.1016/j.gca.2011.07.035
  * Newcombe et al. (2017) Newcombe, M., Brett, A., Beckett, J., et al. 2017, GeCoA, 200, 330, doi: https://doi.org/10.1016/j.gca.2016.12.026
  * Nikolaou et al. (2019) Nikolaou, A., Katyal, N., Tosi, N., et al. 2019, ApJ, 875, 11, doi: https://doi.org/10.3847/1538-4357/ab08ed
  * Olson & Sharp (2019) Olson, P. L., & Sharp, Z. D. 2019, PEPI, 294, 106294, doi: https://doi.org/10.1016/j.pepi.2019.106294
  * O’Neill & Eggins (2002) O’Neill, H. S., & Eggins, S. M. 2002, ChGeo, 186, 151, doi: https://doi.org/10.1016/S0009-2541(01)00414-4
  * Pahlevan et al. (2019) Pahlevan, K., Schaefer, L., & Hirschmann, M. M. 2019, E&PSL, 526, 115770, doi: https://doi.org/10.1016/j.epsl.2019.115770
  * Pan et al. (1991) Pan, V., Holloway, J. R., & Hervig, R. L. 1991, GeCoA, 55, 1587, doi: https://doi.org/10.1016/0016-7037(91)90130-W
  * Patoc̆ka et al. (2020) Patoc̆ka, V., Calzavarini, E., & Tosi, N. 2020, PhRvF, 5, doi: https://doi.org/10.1103/PhysRevFluids.5.114304
  * Perera et al. (2018) Perera, V., Jackson, A. P., Elkins-Tanton, L. T., & Asphaug, E. 2018, JGRP, 123, 1168, doi: https://doi.org/10.1029/2017JE005512
  * Putirka & Rarick (2019) Putirka, K. D., & Rarick, J. C. 2019, AmMin, 104, 817, doi: https://doi.org/10.2138/am-2019-6787
  * Raymond & Izidoro (2017) Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134, doi: https://doi.org/10.1016/j.icarus.2017.06.030
  * Raymond et al. (2007) Raymond, S. N., Quinn, T., & Lunine, J. I. 2007, AsBio, 7, 66, doi: https://doi.org/10.1089/ast.2006.06-0126
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, JQSRT, 111, 2139, doi: https://doi.org/10.1016/j.jqsrt.2010.05.001
  * Rubie et al. (2015) Rubie, D., Nimmo, F., & Melosh, H. 2015, in Treatise on Geophysics, 2nd edn., ed. G. Schubert, Vol. 9 (Elsevier), 43–79, doi: https://doi.org/10.1016/B978-0-444-53802-4.00154-8
  * Rumpf & Gupte (1971) Rumpf, H. C. H., & Gupte, A. R. 1971, Chemie Ingenieur Technik, 43, 367, doi: https://doi.org/10.1002/cite.330430610
  * Salvador et al. (2017) Salvador, A., Massol, H., Davaille, A., et al. 2017, JGRE, 122, 1458, doi: https://doi.org/10.1002/2017JE005286
  * Scaillet & Gaillard (2011) Scaillet, B., & Gaillard, F. 2011, Natur, 480, 48, doi: https://doi.org/10.1038/480048a
  * Schaefer & Fegley (2010) Schaefer, L., & Fegley, B. 2010, Icarus, 208, 438, doi: https://doi.org/10.1016/j.icarus.2010.01.026
  * Schaefer & Fegley (2017) —. 2017, ApJ, 843, 120, doi: https://doi.org/10.3847/1538-4357/aa784f
  * Schaefer et al. (2016) Schaefer, L., Wordsworth, R. D., Berta-Thompson, Z., & Sasselov, D. 2016, ApJ, 829, 63, doi: https://doi.org/10.3847/0004-637X/829/2/63
  * Solomatov (2000) Solomatov, V. S. 2000, in Origin of the Earth and Moon, ed. R. M. Canup & K. Righter, Space Science Series (The University of Arizona Press), 323–338
  * Solomatov et al. (1993) Solomatov, V. S., Olson, P., & Stevenson, D. J. 1993, E&PSL, 120, 387, doi: https://doi.org/10.1016/0012-821X(93)90252-5
  * Solomatov & Stevenson (1993a) Solomatov, V. S., & Stevenson, D. J. 1993a, JGRE, 98, 5375, doi: https://doi.org/10.1029/92JE02948
  * Solomatov & Stevenson (1993b) —. 1993b, JGRE, 98, 5391, doi: https://doi.org/10.1029/92JE02579
  * Sossi (2021) Sossi, P. A. 2021, NatAs, doi: https://doi.org/10.1038/s41550-021-01353-9
  * Sossi et al. (2020) Sossi, P. A., Burnham, A. D., Badro, J., et al. 2020, SciA, 6, doi: https://doi.org/10.1126/sciadv.abd1387
  * Sossi et al. (2022) Sossi, P. A., Tollan, P., Badro, J., & Bower, D. J. 2022, E&PSL
  * Stanley et al. (2014) Stanley, B. D., Hirschmann, M. M., & Withers, A. C. 2014, GeCoA, 129, 54, doi: https://doi.org/10.1016/j.gca.2013.12.013
  * Stolper (1982) Stolper, E. 1982, GeCoA, 46, 2609, doi: https://doi.org/10.1016/0016-7037(82)90381-7
  * Stolper & Holloway (1988) Stolper, E., & Holloway, J. R. 1988, E&PSL, 87, 397, doi: https://doi.org/10.1016/0012-821X(88)90004-0
  * Tagawa et al. (2021) Tagawa, S., Sakamoto, N., Hirose, K., et al. 2021, NatCo, 12, 2588, doi: https://doi.org/10.1038/s41467-021-22035-0
  * Takahashi et al. (2013) Takahashi, S., Ohtani, E., Terasaki, H., et al. 2013, PCM, 40, 647, doi: https://doi.org/10.1007/s00269-013-0600-x
  * Thompson et al. (2021) Thompson, M. A., Telus, M., Schaefer, L., et al. 2021, NatAs, doi: https://doi.org/10.1038/s41550-021-01338-8
  * Tonks & Melosh (1990) Tonks, W., & Melosh, H. 1990, in Origin of the Earth, ed. H. E. Newsom & J. H. Jones (Oxford University Press, New York), 151–174
  * Tu et al. (2015) Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3, doi: https://doi.org/10.1051/0004-6361/201526146
  * Wetzel et al. (2013) Wetzel, D. T., Rutherford, M. J., Jacobsen, S. D., Hauri, E. H., & Saal, A. E. 2013, PNAS, 110, 8010, doi: https://doi.org/10.1073/pnas.1219266110
  * Wilson & Head III (1981) Wilson, L., & Head III, J. W. 1981, JGRB, 86, 2971, doi: https://doi.org/10.1029/JB086iB04p02971
  * Wolf & Bower (2018) Wolf, A. S., & Bower, D. J. 2018, PEPI, 278, 59, doi: https://doi.org/10.1016/j.pepi.2018.02.004
  * Yoshioka et al. (2019) Yoshioka, T., Nakashima, D., & Keppler, H. 2019, GeCoA, 259, 129, doi: https://doi.org/10.1016/j.gca.2019.06.007
  * Yurchenko & Tennyson (2014) Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649, doi: https://doi.org/10.1093/mnras/stu326
  * Yurchenko et al. (2013) Yurchenko, S. N., Tennyson, J., Barber, R. J., & Thiel, W. 2013, JMoSp, 291, 69, doi: https://doi.org/10.1016/j.jms.2013.05.014
  * Zahnle et al. (2019) Zahnle, K. J., Gacesa, M., & Catling, D. C. 2019, GeCoA, 244, 56, doi: https://doi.org/10.1016/j.gca.2018.09.017
  * Zahnle et al. (2020) Zahnle, K. J., Lupu, R., Catling, D. C., & Wogan, N. 2020, PSJ, 1, 11, doi: https://doi.org/10.3847/psj/ab7e2c