# EVOLUTION OF MERCURY’S EARLIEST ATMOSPHERE

Noah Jäggi Physics Institute, University of Bern, Sidlerstrasse 5, 3012 Bern,
Switzerland  Diana Gamborino Physics Institute, University of Bern,
Sidlerstrasse 5, 3012 Bern, Switzerland  Dan J. Bower 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  Aaron S. Wolf Earth and Environmental Sciences, University of
Michigan, 1100 North University Avenue, Ann Arbor, MI 48109-1005, USA  Apurva
V. Oza Jet Propulsion Laboratory, California Institute of Technology,
Pasadena, USA  Physics Institute, University of Bern, Sidlerstrasse 5, 3012
Bern, Switzerland  Audrey Vorburger Physics Institute, University of Bern,
Sidlerstrasse 5, 3012 Bern, Switzerland  André Galli Physics Institute,
University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland  Peter Wurz
Physics Institute, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

(Received March 26, 2021 ; Revised September 2, 2021; Accepted October 8,
2021)

###### Abstract

MESSENGER observations suggest a magma ocean formed on proto-Mercury, during
which evaporation of metals and outgassing of C- and H-bearing volatiles
produced an early atmosphere. Atmospheric escape subsequently occurred by
plasma heating, photoevaporation, Jeans escape, and photoionization. To
quantify atmospheric loss, we combine constraints on the lifetime of surficial
melt, melt composition, and atmospheric composition. Consideration of two
initial Mercury sizes and four magma ocean compositions determine the
atmospheric speciation at a given surface temperature. A coupled
interior–atmosphere model determines the cooling rate and therefore the
lifetime of surficial melt. Combining the melt lifetime and escape flux
calculations provide estimates for the total mass loss from early Mercury.
Loss rates by Jeans escape are negligible. Plasma heating and photoionization
are limited by homopause diffusion rates of $\sim 10^{6}$ kg/s. Loss by
photoevaporation depends on the timing of Mercury formation and assumed
heating efficiency and ranges from $\sim 10^{6.6}$ to $\sim 10^{9.6}$ kg/s.
The material for photoevaporation is sourced from below the homopause and is
therefore energy-limited rather than diffusion-limited. The timescale for
efficient interior–atmosphere chemical exchange is less than ten thousand
years. Therefore, escape processes only account for an equivalent loss of less
than 2.3 km of crust ($0.3\%$ of Mercury’s mass). Accordingly, $\leq 0.02\%$
of the total mass of H2O and Na is lost. Therefore, cumulative loss cannot
significantly modify Mercury’s bulk mantle composition during the magma ocean
stage. Mercury’s high core:mantle ratio and volatile-rich surface may instead
reflect chemical variations in its building blocks resulting from its solar-
proximal accretion environment.

Solar system, Solar system astronomy, Mercury (planet), Planetary surface,
Planetary atmospheres, Atmospheric science

††journal: The Planetary Science Journal††software: SPIDER (Bower et al.,
2019, 2018), VapoRock (Wolf et al., 2021), VULCAN (Tsai et al., 2017), DISHOOM
(Oza et al., 2019; Gebek & Oza, 2020), E-MC (Wurz & Lammer, 2003)

##  1 Introduction

MESSENGER data from X-ray, gamma ray, and neutron spectrometers constrain the
composition of Mercury’s surface, and motivate theories and models to
understand Mercury’s bulk composition, formation, and evolution. The surface
composition and geology of Mercury is compatible with partial melting of
cumulates that were originally formed by magma ocean crystallization (McCoy et
al., 2018). Subsequent impact excavation exposed the cumulates at the surface
(McCoy et al., 2018; Charlier et al., 2013). The low oxygen fugacity (_f_ O_2)
of the uppermost layer of Mercury’s regolith, together with Mercury’s large
core size, suggest a reduced mantle where nominally lithophile elements such
as Ca, Mn, Cr, and Ti are present in sulfides rather than silicates (Vander
Kaaden & McCubbin, 2016). Relative to basaltic rocks exposed at the surface of
other terrestrial planets, a large amount of the moderately volatile element
Na (3–5 wt%) is detected on Mercury’s surface (Peplowski et al., 2014).
Observations of Na variation in Mercury’s exosphere may relate to night-side
deposit formation and dawn re-emission (e.g., Cassidy et al., 2016). Hence, it
remains an open question how moderately volatile elements such as Na may have
accumulated on the surface—whether from an extant or now extinct process—and
how their abundance compares to Mercury’s bulk composition.

Magma oceans are pivotal in determining the initial conditions and subsequent
evolution and chemical differentiation of terrestrial planets in the solar
system (e.g., Elkins-Tanton, 2012; Chao et al., 2021). Radiometric dating
reveals that magmatic iron meteorites, which represent planetesimal cores,
formed within 2 Myr of solar system formation (Kruijer et al., 2014). The
rocky planet whose mass is most similar to that of Mercury, and for which
samples are available, Mars, likely accreted, formed an iron core, and
underwent complete solidification of its magma ocean ocean within about 20 Myr
of solar system formation (Bouvier et al., 2018). Crucially, rapid core
formation in terrestrial planets requires a magma ocean to enable efficient
metal segregation (Stevenson, 1990). By analogy, and given its solar-proximal
location, extensive melting is therefore expected to have occurred on proto-
Mercury (Vander Kaaden & McCubbin, 2016; Brown & Elkins-Tanton, 2009).
Following its crystallisation, partial melting of magma ocean cumulates have
been invoked to explain Mercury’s contemporary surface composition (McCoy et
al., 2018).

Energy from accretion and radiogenic heat (e.g. from 26Al) may have driven the
differentiation of Mercury if it formed sufficiently early in solar system
history (Bhatia & Sahijpal, 2017; Siegfried & Solomon, 1974). Following a
phase of rapid growth, the subsequent reduction of impactor flux would have
enabled Mercury’s magma ocean to cool and crystallize without additional
large-scale remelting. During this time, the mantle is expected to have
stratified into a basal layer of olivine and a plagioclase and clinopyroxene
dominated crust which is now observed on Mercury’s surface (Brown & Elkins-
Tanton, 2009). During the cooling of the magma ocean when the surface remains
mostly molten, chemical species readily exchange between the interior,
atmosphere, and exosphere—as occurred for other terrestrial planets in the
solar system (e.g., Elkins-Tanton, 2008).

Fegley & Cameron (1987) addressed the hypothesis that the anomalously high
bulk density of Mercury (owing to a high core/mantle ratio) is the result of
evaporation of silicate melt components from the surface of a Hermean magma
ocean. They presumed atmospheric loss was sufficiently slow that the
atmosphere remained in equilibrium with the magma ocean. In their model, vapor
was removed in a step-wise fashion and the composition of the magma ocean
evolved accordingly. In reality, however, evaporated species are transported,
mixed, and lost from the atmosphere and exosphere, with the flux at which loss
occurs integrated over the magma ocean lifetime ultimately dictating the total
mass loss. Therefore, consideration of interior, atmospheric, and exospheric
processes are necessary to assess whether significant quantities of rock-
derived atmospheres can be lost during the Hermean magma ocean stage.

Based on observations of solar-mass stars, the early solar extreme ultraviolet
(EUV) and X-ray fluxes were likely 400 times larger than they are today. This
would have made photoionization a highly efficient non-thermal, and
photoevaporation a highly efficient thermal, atmospheric escape mechanism
(Johnstone et al., 2015; Tu et al., 2015). Other loss mechanisms of potential
importance include atmospheric sputtering and kinetic escape (e.g. Jeans
escape) that occur over the lifetime of the magma ocean. Non-thermal loss
rates can be constrained by known plasma pressures at proto-Mercury due to the
incoming solar wind, as well as EUV luminosities of the early Sun estimated
from population studies of nearby Sun-like stars (Ribas et al., 2014; Tu et
al., 2015).

In this paper we establish the extent of element evaporation and loss from
Mercury during its early magma ocean phase. Models are constructed of the (1)
coupled evolution of the magma ocean and atmosphere, (2) evaporation of metals
and metal-oxide species from the Hermean magma ocean, (3) the mixing ratios
and abundances of molecular species throughout the atmosphere and at the
exobase, and finally (4) loss rates of these species from the upper
atmosphere. We discuss these results in the context of the chemical evolution
of Mercury’s surface environment, bulk composition, and present-day
observations.

##  2 METHOD

###  2.1 Overview

Our combined modeling strategy provides insight into the initial composition
and evolution of Mercury’s exosphere by considering (1–2) energy and mass
exchange between the interior and atmosphere, (3) speciation in the
atmosphere, (4-5) loss from the exosphere:

  1. 1.

SPIDER (Bower et al., 2021, 2019, 2018) is a coupled interior–atmosphere model
used to determine the surface temperature and lifetime of melt at the surface,
as well as the pressure–temperature structure of the atmosphere. Magma ocean
cooling is regulated by the atmospheric opacity, which depends on the pressure
of atmospheric species and hence the solubility of species in melt.

  2. 2.

VapoRock calculates the equilibrium partial pressures of metal-bearing gas
species of the elements Si, Al, Mg, Ca, Na, Fe and K above the magma ocean
surface (Wolf et al., 2021). This determines the metal-bearing composition of
the atmosphere as a function of temperature and the bulk composition of the
magma ocean. It utilises ENKI’s ThermoEngine (http://enki-portal.org) and
combines estimates for element activities in silicate melts with thermodynamic
data for metal and metal oxide vapor species (Lamoreaux et al., 1987;
Lamoreaux & Hildenbrand, 1984).

  3. 3.

VULCAN (Tsai et al., 2017, 2021) solves for the equilibrium chemistry of the
atmosphere as a function of altitude by using element abundances for metals
(output by VapoRock), volatile abundances (output by SPIDER), and the
atmospheric pressure–temperature structure (also output from SPIDER). This
provides the mixing ratios of atmospheric species, which are required to
calculate escape at the exobase.

  4. 4.

DISHOOM (Oza et al., 2019; Gebek & Oza, 2020) is an atmospheric evolution
model which computes the total mass loss of gaseous species to space due to
ultraviolet (EUV) heating, surface heating, as well as plasma heating.

  5. 5.

Exospheric Monte Carlo (E-MC) model (Gamborino et al., 2019; Vorburger et al.,
2015; Wurz & Lammer, 2003) determines the rate of exospheric escape of
particles due to Jeans escape and photoionization. It tracks particle
trajectories using a thermal energy distribution that depends on the
temperature at the exobase.

###  2.2 Cooling of the magma ocean

Previous thermal modeling of Mercury’s interior has focused either on the
accretion phase (Bhatia & Sahijpal, 2017) or its long-term evolution over
billions of years (e.g., Tosi et al., 2013; Grott et al., 2011; Spohn, 1991;
Stevenson et al., 1983). Here, we model the thermal evolution of Mercury’s
magma ocean at the end of its accretion phase. At this time, the final magma
ocean cools and crystallizes on a timescale short enough such that there is
negligible disruption so our results remain independent of its accretion
history. We model the thermal evolution of Mercury’s magma ocean using SPIDER
(Bower et al., 2018, 2019) to constrain the duration of melt at the surface as
it cools from 2400 K to 1500 K. This is necessary to compute the evaporation
of metals and metal oxides at the planetary surface, prior to the formation of
a surface lid around 1500 K. Heating by the decay of radiogenic isotopes 26Al,
40K, 232Th, 235U, 238U is included in our model and the model starts from
solar system time zero to obtain an upper estimate of the surface cooling
time. The main parameters are provided in Table 6 (Appendix A)and are guided
by the parameters and results from previous models of Mercury (Bhatia &
Sahijpal, 2017; Tosi et al., 2013).

Cases prefixed by ‘S’ (‘Small Mercury’, Table 1) have a planetary radius of
2440 km, which is the present-day radius of Mercury. Cases prefixed by ‘L’
have a radius of 3290 km, which assumes Mercury was larger than at present
day, perhaps due to mantle stripping driven by an impactor (Chau et al., 2018;
Asphaug & Reufer, 2014; Benz et al., 2008). Cases with ‘V’ (volatile) consider
the partitioning of carbon and hydrogen species, here termed volatiles,
between the melt and atmosphere as well as redox reactions (Bower et al.,
2021). By contrast, cases with ‘N’ (non-volatile) do not consider volatiles
but rather assume SiO is the only IR absorbing species; a suffix of ‘5’
denotes a low SiO opacity ($10^{-5}$ m2 kg-1) and a suffix of ‘3’ denotes a
large SiO opacity ($10^{-3}$ m2 kg-1), both at a reference pressure of
$3\times 10^{-6}$ bar (Fig. 2, Semenov et al., 2003). For non-volatile cases,
the surface pressure of SiO is imposed in SPIDER as a function of surface
temperature as determined by VapoRock calculations. Then the atmospheric
opacity and hence magma ocean cooling rate can be determined.

Table 1:  Parameters for the magma ocean cases.

Case | $R_{P}$ | $g$ |  H_2O |  H_2 |  CO_2 |  CO |  SiO
---|---|---|---|---|---|---|---
| km | m/s2 | Pressure (bar) at 2000 K
SN5 | 2440 | 3.7 | — | — | — | — | $1.4E-4$
SN3 | 2440 | 3.7 | — | — | — | — | $1.4E-4$
SV | 2440 | 3.7 | 0.7 | 3.2 | 0.05 | 1.1 | —
LN5 | 3290 | 4.0 | — | — | — | — | $1.4E-4$
LN3 | 3290 | 4.0 | — | — | — | — | $1.4E-4$
LV | 3290 | 4.0 | 1.2 | 5.8 | 0.2 | 4.9 | —

Note. — Small (S) and large (L) Mercury models are inspired by models M1 and
M4 from Bhatia & Sahijpal (2017), respectively. Second letter of the case name
denotes with volatiles (V) and no volatiles (N). Non-volatile cases have an
additional number of either 5 or 3, to denote small ($10^{-5}$ m2 kg-1) and
large ($10^{-3}$ m2 kg-1) SiO opacity, respectively (Semenov et al., 2003).

For cases SV and LV, carbon and hydrogen can either exist as reduced (CO, H_2)
or oxidised (CO_2, H_2O) species, where the _f_ O_2 is constrained to one log
unit below the iron-wustite buffer ($\Delta$IW = $-1$, IW-1 hereafter). This
is marginally higher than the most recent estimates for the _f_ O_2 of
Mercury’s mantle (Cartier & Wood, 2019). For cases SV and LV, we determine the
total H and C abundances by calculating the ppmw necessary for an Earth-size
planet to produce a 100 bar CO_2 (i.e. Venus-like atmosphere) and 270 bar H_2O
(i.e. one Earth ocean mass) atmosphere. The abundances of H and C are
equivalent to 330 ppmw of H_2O and 120 ppmw of CO_2, respectively. The mass of
large Mercury’s mantle is about a factor of 5 larger than small Mercury’s
mantle, resulting in a 5 times increase in the total volatile budget.

###  2.3 Evaporation from the magma ocean

At the high surface temperatures that characterize a magma ocean ($>$1500 K),
the partial pressures of the vapor species of the major rock-forming oxides
(e.g. SiO_2, NaO0.5, KO0.5) can become significant (Sossi & Fegley, 2018;
Sossi et al., 2019; Visscher & Fegley, 2013). Gas-liquid equilibria for these
elements are described by congruent evaporation, generalised as:

$\mathrm{M}^{x+n}\mathrm{O}_{(x+n)/2}(l)=\mathrm{M}^{x}\mathrm{O}_{x/2}(g)+\frac{n}{4}\mathrm{\ce{O_{2}}}(g),$
(1)

where M is the metal, $x$ is the oxidation state of the metal in its gaseous
state, and $n$ the number of electrons exchanged in the reaction. Both $x$ and
$n$ are integer values that may be $\geq{0}$ or $\leq{0}$. At equilibrium, the
partial pressure of any given metal or metal-oxide species in an ideal gas is
given by:

$p(\mathrm{M}^{x}\mathrm{O}_{x/2})=\frac{K_{(1)}X(\mathrm{M}^{x+n}\mathrm{O}_{(x+n)/2})\gamma(\mathrm{M}^{x+n}\mathrm{O}_{(x+n)/2})}{(f\mathrm{\ce{O_{2}}})^{(n/4)}},$
(2)

where $K_{(1)}$ is the equilibrium constant of the reaction (Eq. 1), $X$ is
the mole fraction, and $\gamma$ is the activity coefficient of the metal oxide
melt species, $\mathrm{M}^{x+n}\mathrm{O}_{(x+n)/2}$. Equilibrium constants
involving 31 gas species (Table 7, Appendix B) are calculated according to
their thermodynamic properties given in Lamoreaux et al. (1987); Lamoreaux &
Hildenbrand (1984). Evident from Eq. 2 is that estimates for the composition
of the silicate melt in addition to the activity coefficients of its
constituent components are required to correctly predict partial pressures. To
this end, likely compositions representative of Mercury’s crust, mantle and
possible precursors are shown in Table 2. The MELTS algorithm is used to
estimate activity coefficients of melt oxide species (Ghiorso & Sack, 1995).
The _f_ O_2 is constrained to lie one log unit below the IW buffer (IW-1),
which is calculated according to O’Neill & Eggins (2002):

$\log
f\mathrm{\ce{O_{2}}}(\mathrm{IW})=2\frac{-244118+115.559\,T-8.474\,T\,\mathrm{ln}(T)}{\mathrm{ln}(10)R\,T},$
(3)

where $T$ is temperature and $R$ the gas constant. These ingredients together
comprise the VapoRock code and permit calculation of equilibrium partial
pressures over a range of temperatures, _f_ O_2, and silicate melt
compositions (Wolf et al., 2021).

###  2.4 Magma ocean composition

The composition of the early Hermean mantle is uncertain. To address how
variability in the surface composition affects the evolved partial pressures
of metal and metal oxide gas species, four compositions are investigated
(Table 2): (1) Enstatite Chondrites (EH4), (2) Bencubbin Chondrites (CB), (3)
Northern Smooth Plain (NSP) lava, and (4) NSP source. The first composition
assumes that core–mantle differentiation was sluggish, with the composition of
the magma ocean being approximated by Enstatite Chondrites (EH4), often cited
as appropriate starting compositions for Mercury due to their high bulk iron
content, strongly reduced nature, and the resemblance of partial melts thereof
to Hermean surface compositions (Nittler et al., 2011). We assume that all FeO
is extracted in the form of metallic iron to form Mercury’s core, which
results in high SiO2 and MgO contents in the complementary silicate fraction.

The second composition is based on chondrules found in the Bencubbin-class of
Carbonaceous Chondrites (CB), which best reproduce Mercury’s surface
composition based on MESSENGER data (McCoy et al., 2018; Brown & Elkins-
Tanton, 2009). Spectrometric measurements of Mercury’s surface show a crust
rich in Na and S and poor in Fe relative to other basaltic rocks (Nittler &
Weider, 2019). Presuming these abundances are representative of bulk Mercury,
two other compositions are investigated; the Northern Smooth Plain lava (NSP
melt) that represents a volatile-rich composition observed on the surface, and
its inferred mantle source (NSP source) (Nittler & Weider, 2019). These two
compositions represent Mercury’s crust and mantle, respectively. Although it
is not anticipated that magma ocean crystallisation produced the NSP melt
composition directly, it is included to define an end-member opposing the CB
composition that is Na and K poor, and comparably Fe-rich (Table 2).

Table 2: Magma ocean surface composition.

Oxide | EH4 | CB | NSP | NSP
---|---|---|---|---
(wt%) |  |  | source | lava
SiO2 | 62.73 | 50.70 | 53.67 | 58.70 |  |
Al2O3 | 2.58 | 4.60 | 4.75 | 13.80 |  |
MgO | 30.24 | 36.90 | 36.89 | 13.90 |  |
CaO | 1.99 | 3.30 | 2.26 | 5.81 |  |
FeO | 0.00 | 3.50 | 0.02 | 0.04 |  |
Na2O | 1.71 | 0.19 | 1.97 | 7.00 |  |
K2O | 0.20 | 0.05 | 0.05 | 0.20 |  |
Total | 99.45 | 99.24 | 99.61 | 99.45 |  |

Note. — Compositions based on Enstatite chondrites (EH4) (Wiik, 1956); CB
chondrite chondrule data with bulk CB Na and K mass balanced for chondrules to
fit bulk meteorite iron-silicate ratio (Lauretta et al., 2007; Weisberg et
al., 2000, 1990); and northern smooth plains (NSP) composition for the lava
and source (Nittler & Weider, 2019; Namur et al., 2016)

###  2.5 Atmospheric structure

The atmospheric structure is constrained by the temperature at the magma
ocean–atmosphere interface, the planetary equilibrium temperature, and the
atmospheric composition and pressure (Appendix E). For non-volatile cases, the
calculated vapor pressures of Si, Na, K, Fe, Mg, Al, and Ca oxide species in
equilibrium with the magma ocean (Fig. 2) are used directly in the exospheric
loss model (Section 2.6). SiO vapour pressures reported in Table 1 at
$T_{\mathrm{surf}}=2000$ K are not strongly affected by the magma ocean
composition and range from $10^{-4\pm 0.1}$ bar. For cases SV and LV in which
Mercury’s atmosphere contains outgassed H- and C-bearing gases, the partial
pressures of H_2, H_2O, CO, and CO_2 are calculated by SPIDER according to
volatile solubility and _f_ O_2 buffered by the magma ocean at IW-1. A
modified version of the VULCAN code is then used to compute the equilibrium
chemical speciation in the atmosphere that contains both the metal-bearing
gases and H and C volatiles (Tsai et al., 2017, 2021). VULCAN computes the
atmospheric mixing ratios using the pressure–temperature ($P$–$T$) structure
of the atmosphere.

VULCAN by default includes about 300 reactions for C, H, O, and N to which we
added reactions involving Si, Mg, Ca, Fe, Na and K to obtain their equilibrium
speciation (Table 8, Appendix C). Table 3 shows the initial element-to-
hydrogen ratios used in the VULCAN calculations. Surface vapor pressures of Ca
and Al-bearing species did not exceed $10^{-6}$ bar in the magma ocean
temperature range investigated, and are thus excluded. For the remaining
species we used a case dependent $P$–$T$ profile from SPIDER at a surface
temperature of 2000 K to determine the mixing ratio of species in the
atmosphere as a function of altitude. The $P$–$T$ profile may imply
condensation of certain elements initially present in the vapor, which may
rain out of the atmosphere prior to escape. To assess this possibility, Gibbs
Free Energy minimisation of the atmospheric composition were preformed
throughout the atmospheric column using FactSage 7.3. (Bale et al., 2016).

Table 3: Surface element ratios at $T_{\mathrm{surf}}=2000$ K.

Composition | Element Ratios
---|---
| C/H | O/H | Mg/H | Si/H | Na/H | K/H | Fe/H
SV Case
EH4 | 1.532E-01 | 2.472E-01 | 1.251E-06 | 1.469E-05 | 1.254E-04 | 1.121E-05 | 1.053E-08
CB | 1.680E-06 | 1.039E-05 | 7.155E-05 | 5.785E-06 | 9.388E-06
NSP source | 1.740E-06 | 1.020E-05 | 2.282E-04 | 5.254E-06 | 7.155E-07
NSP melt | 9.642E-07 | 1.161E-05 | 3.480E-04 | 7.139E-06 | 1.351E-06
LV Case
EH4 | 3.608E-01 | 4.641E-01 | 6.806E-07 | 7.997E-06 | 6.825E-05 | 6.104E-06 | 5.732E-09
CB | 9.144E-07 | 5.653E-06 | 3.894E-05 | 3.148E-06 | 5.109E-06
NSP source | 9.469E-07 | 5.551E-06 | 1.242E-04 | 2.859E-06 | 3.894E-07
NSP melt | 5.248E-07 | 6.321E-06 | 1.894E-04 | 3.886E-06 | 7.354E-07

Note. — Ratios are based on SPIDER and VapoRock results for magma ocean
compositions given in Table 2. We neglect Al and Ca as their vapor pressure do
not exceed $10^{-9}$ bar at 2000 K for any composition.

SPIDER, and for volatile cases also VULCAN, provide descriptions of the
atmospheric structure and composition needed to determine the altitude of the
homopause and exobase. The homopause is the altitude at which molecular
diffusion exceeds diffusion by eddies and thus separates the well-mixed lower
atmosphere from the mass-separated upper atmosphere. The exobase is the
altitude at which gas is loosely bound to the planet and is collisionless
(Knudsen number = 1) resulting in efficient escape.

####  2.5.1 Homopause level and diffusion

To determine the homopause level $z_{\mathrm{hom}}$ we require the particle
density at the homopause $n_{\mathrm{hom}}$. For a steady-state homopause
height, the molecular coefficient $D_{ik}$ is equal to the eddy diffusion
coefficient $K_{zz}$, allowing us to solve for the particle density
$n_{\mathrm{hom}}$. The diffusion coefficient $D_{ik}$ (m2/s) within the
homosphere is calculated for each major species using the Chapman-Enskog
relation (Chapman & Cowling, 1970). It determines the binary diffusion rate of
a gaseous species $i$ with mass $m_{i}$ within a gas of average mass $m_{k}$:

$D_{ik}=\frac{3\sqrt{k_{B}T_{\mathrm{skin}}}}{8n_{\mathrm{hom}}\sigma^{2}_{ik}\Omega_{ik}}\sqrt{\frac{1}{2\pi}\frac{m_{i}+m_{k}}{m_{i}m_{k}}},$
(4)

where, $k_{B}$ is the Boltzmann constant and $T_{\mathrm{skin}}$ the absolute
temperature at the homopause (skin temperature). Homopause pressures
($P_{\mathrm{hom}}=n_{\mathrm{hom}}k_{B}T_{\mathrm{skin}}$) for each species
are thereby about $10^{-6}$ bar for all cases. We approximate the
intermolecular distance $\sigma_{ik}=\frac{\sigma_{i}+\sigma_{k}}{2}$ with the
radius of the species relative to the mean species diameter weighted by the
mixing ratio. The dimensionless collision integral $\Omega_{ik}$ is assumed to
be unity.

The eddy velocity is often approximated by the atmospheric species thermal
speed $v_{\mathrm{th}}$, and the characteristic eddy length scale $L_{eddy}$
is approximated by the atmospheric scale height $H$ (e.g. Atreya et al.,
1986). Values for $K(z)=v_{\mathrm{eddy}}\,L_{\mathrm{eddy}}$ that are
calculated based on this assumption exceed the suggested eddy diffusion
coefficient upper limit of 320 m2/s by several orders of magnitude (Vlasov &
Kelley, 2015). Hence we use this upper limit in the volatile cases to
determine $n_{\mathrm{homo}}$, which is based on the energy dissipation rate
within the Earth’s atmosphere.

To compute the $z_{\mathrm{hom}}$ and $T_{\mathrm{hom}}$ for the volatile
cases we determine the altitude at which the $P$–$T$ profile reaches a number
density $n_{\mathrm{homo}}(K(z))$. As the number density at the homopause only
depends on $K_{zz}$ with the same order of magnitude for both volatile and
non-volatile cases, $n_{\mathrm{homo}}(K(z))$ is approximately
$10^{18}-10^{19}\,\mathrm{at}\,\mathrm{m}^{-3}$. For the non-volatile cases we
do not obtain a P(z)-T profile from VULCAN for the non-volatile cases due to
the absence of H based species. We therefore use the barometric formula with
gravity as a function of height to compute $z_{\mathrm{hom}}$.

####  2.5.2 Exobase level

Due to the large difference in number density between the homopause and the
exobase, the barometric formula is not applicable assuming an isothermal upper
atmosphere with height-dependent gravity. We approximate the exobase height of
early Mercury, which is subject to extensive loss, by finding an exobase
height that results in a loss rate that is in equilibrium with the homopause
diffusion rate. The loss from the exobase is proportional to the exobase
height (increasing surface area) whereas diffusion from the homopause to the
exobase is inversely proportional to the exobase height (decreasing density
gradient). The exobase altitude $z_{\mathrm{exo}}$ of each species is
determined for all cases by setting the homopause diffusion rate
$\dot{M}_{\mathrm{diff},i}$ equal to the largest, diffusion limited mass loss
rate of photoionization $\dot{M}_{\mathrm{ion}}$ (Eq. 9, Section 2.6).

The homopause diffusion rate $\dot{M}_{\mathrm{diff}}$ in kg/s of a species
$i$ is obtained by multiplying the diffusion coefficient $D_{ik}$ with the
species number density gradient, the species mixing ratio at the homopause
$n_{i}/n_{hom}$, and the homopause surface area $A_{\mathrm{hom}}$:

$\dot{M}_{\mathrm{diff},i}=-D_{ik}\frac{\Delta n_{\mathrm{exo-hom}}}{\Delta
z_{\mathrm{exo-hom}}}\,m_{i}\frac{n_{i}}{n_{\mathrm{hom}}}\,A_{\mathrm{hom}}.$
(5)

The number density of particles at the exobase $n_{\mathrm{exo}}$ is necessary
to determine the number density gradient between the homopause and exobase and
ultimately $z_{\mathrm{exo}}$. For a single species $i$, the exobase is
defined at a altitude at which the particle free path
($\lambda_{col}=1/n_{\mathrm{exo}}\sigma_{\mathrm{col}}$) is equal to the
exospheric scale height ($H=k_{B}T/mg$), therefore (i.e., Gronoff et al.,
2020):

$n_{i,\mathrm{exo}}=\frac{1}{H_{i}\sigma_{i}}=\frac{m_{i}\,g(h)}{k_{B}\,T_{\mathrm{skin}}\,\sigma_{i}},$
(6)

with the collision cross section $\sigma_{i}$, the skin temperature
$T_{\mathrm{skin}}$, and the acceleration of gravity $g(h)$ at the exobase
altitude $z_{\mathrm{exo}}$. In a multi-species atmosphere, each species has a
specific mass and collision cross section, leading to a species-specific scale
height and exobase density and altitude. The collision cross sections (CCSs)
of each species are approximated as their respective atom or molecule size
(Table 9, Appendix D). Typical values for $n_{\mathrm{exo}}$ are around
$10^{12}-10^{13}\,\mathrm{atom}\,\mathrm{m}^{-3}$ which coincides with $\sim
10^{-12}$ bar. The skin temperature used for determining $z_{\mathrm{exo}}$
and $z_{\mathrm{hom}}$ is derived from the atmosphere model (Appendix E):

$T_{\mathrm{skin}}=\left[\frac{\epsilon(T_{\mathrm{surf}}^{4}-T_{\infty}^{4})}{2}+T_{\infty}^{4}\right]^{1/4},$
(7)

with the magma ocean surface temperature $T_{\mathrm{surf}}$, emissivity
$\epsilon$ (depends on optical depth and hence atmospheric composition and
pressure), and equilibrium temperature $T_{\infty}$.

###  2.6 Exospheric loss

The E-MC escape model focuses on Jeans escape, and photoionization and
photodissociation to investigate the loss of proto-Mercury’s exosphere. These
mechanisms compete for importance; Jeans escape acts at high exospheric
temperatures, whereas photoionization and photodissociation act at large solar
EUV and X-ray fluxes present during early times, respectively. The E-MC model
simulates escape by tracking $\sim 10^{5}$ exospheric particles, with
trajectories initiated at the exobase with a random angle and energy selected
from a Maxwellian velocity distribution function (Vorburger et al., 2015). The
initial energy of the exospheric particles depends on the exobase temperature,
which decreases with time due to magma ocean and atmospheric cooling.

The loss processes in the E-MC model are calculated on a particle by particle
basis. As soon as a particle reaches Mercury’s Hill radius (Table 4), it is
assumed to have escaped Mercury’s gravitational attraction and is subsequently
removed from the simulation. Another loss process is through interaction with
photons. At each altitude step starting from the exobase and moving away from
the planet, the E-MC model calculates the probability of a particle being
photo-dissociated or photoionized. If the particle is photodissociated, the
code calculates the corresponding trajectories of the fragments and assesses
the chance of escaping the gravitational well and the potential for subsequent
photoionization. Ionized particles are considered lost from the exosphere,
assuming they are picked up by the electromagnetic forces of the solar wind
plasma or Mercury’s magnetospheric plasma. Photoionization and
photodissociation rates are scaled for each dominant species using the EUV
flux of the early Sun. The EUV flux and mass loss are dependent on the
rotational evolution of the Sun, with a fast rotator being much more active
than a slow rotator (Johnstone et al., 2015; Tu et al., 2015). We consider a
moderately fast rotating Sun, where the EUV luminosity $L_{\mathrm{EUV}}$
(J/s) is:

$L_{\mathrm{EUV}}=(4.7\times 10^{25})\ t^{-1.18},$ (8)

where $t$ (Ma) is the time since the formation of the solar system. Typical
values for the incident EUV fluxes at Mercury at 1 Ma and 5 Ma are thereby
$10^{3.0}$ and $10^{2.2}$ J/s m-2 respectively.

The loss is equal to the sum of the exospheric particles that have either been
photoionized or lost through gravitational escape. The loss of a given species
from the exosphere at a given time is then calculated using its mixing ratio
at the exobase. For a given exobase temperature, the loss rate
($\dot{M}=dm/dt$) from the exosphere by photoionization
$\dot{M}_{\mathrm{ion}}$ is:

$\dot{M}_{\mathrm{ion}}=\sum_{i=0}\,\psi^{\mathrm{source}}_{i}\,\xi_{\mathrm{i}}\,m_{i}A_{\mathrm{exo}},$
(9)

where $A_{\mathrm{exo}}$ is the surface area of the exobase and
$0\leq\xi_{i}\leq 1$ is the fraction of lost particles of a species $i$ with
mass $m_{i}$, thermal speed of $v_{i}^{\mathrm{therm}}$, and particle flux
leaving the exobase of
$\psi^{\mathrm{source}}_{\mathrm{i}}=n^{exo}_{i}\,v_{i}^{therm}$. The area of
the exobase is equal to sum of the $R_{P}$ and $z_{\mathrm{exo}}$ (Section
2.5). The total loss is determined by integrating the loss flux over the
lifetime of surficial melt.

###  2.7 Atmospheric loss and surface evaporation

Atmospheric loss by thermal processes (photoevaporation) and non-thermal
processes (plasma-heating) are determined using DISHOOM (Oza et al., 2019),
using Eq. 10 and Eq. 11, respectively. Preliminary calculation of Jeans-escape
using DISHOOM demonstrated negligible loss due to surface heating compared to
all other mechanisms. The escape parameters appropriate to proto-Mercury are
summarized in Table 4.

Table 4:  Escape parameters $\lambda_{0}$ for small Mercury and large Mercury
cases at $T_{\mathrm{surf}}=2000$ K.

Species | $\lambda_{0}$ | $\mu_{\mathrm{atm}}$ | $T_{\mathrm{exo}}$ | $z_{\mathrm{exo}}$ | Case
---|---|---|---|---|---
|  | [amu] | [K] | [km] |
Small Mercury, $R_{H}$ = 72 $R_{P}$
Na, K, Fe | 7.8 | 24.1 | 1613 | 2670 | SN5, SN3
H, C, O | 7.6 | 14.3 | 1021 | 2370 | SV
Large Mercury, $R_{H}$ = 90 $R_{P}$
Na, K, Fe | 15.7 | 24.1 | 1615 | 1890 | LN5, LN3
H, C, O | 15.3 | 14.3 | 893 | 2160 | LV

Note. — The escape parameter $\lambda_{0}$ with the respective mean molecular
weight of the upper atmosphere $\mu_{\mathrm{atm}}$, as well as the exobase
temperature $T_{\mathrm{exo}}$ and altitude $z_{\mathrm{exo}}$. The Hill
radius $R_{H}$ in Mercury radii $R_{P}$ describes the gravitational field of
influence of Mercury in each case.

Irradiation from the impinging solar wind plasma and high energy photons may
heat the atmosphere and drive escape at a level that is significantly larger
than surface heating and photoionization. Upper-atmospheric heating
(photoevaporation) is caused by incoming X-ray and EUV photons that deposit
heat into a neutral medium via molecular absorption (Watson et al., 1981) or
photoelectric heating (Murray-Clay et al., 2009). This expands the atmospheric
envelope beyond the gravitational influence of the body ($R_{H}$ in Table 4).
The heating can be estimated by energy-limited escape driven by EUV photons
(Watson et al., 1981), which is a reasonable approximation to thermally-driven
hydrodynamic escape (Krenn et al., 2021; Volkov & Johnson, 2013):

$\begin{split}\dot{M}_{U}&=\frac{\eta_{\mathrm{EUV}}\,L_{\mathrm{EUV}}\,2\pi\,z_{\mathrm{abs}}^{2}}{U_{\mathrm{env}}},\end{split}$
(10)

where, $z_{\mathrm{abs}}$ is the absorption altitude, generally taken to be
1.25 planetary radii for an outgassing atmosphere (e.g., Johnson et al., 2015)
where the X-ray and EUV photons can absorb and thereby deposit heat into the
atmospheric molecules. We note that we use 1.25 Mercury radii $R_{P}$, as a
conservative lower-limit in absorption altitude as the homopause situated at
$\sim$ 1.4 $R_{P}$ represents an upper-limit. The efficiency at which the
atmosphere is heated, $\eta_{EUV}$, is uncertain so we use $10^{-3}$ as a
conservative lower estimate (Ito & Ikoma, 2021) and $10^{-1}$ as an upper
limit (Mordasini, 2020). Both efficiencies that we used were previously
applied to atmospheres that use vastly different planet parameters but similar
enough as they consider a metal oxide (non-volatile) or a H/He (volatile)
atmosphere respectively. Hot Jupiter H/He envelopes as well as volcanic
atmospheres suggest $\eta_{EUV}$ may be as large as $0.35$ (Murray-Clay et
al., 2009; Lellouch et al., 1992).

Mass loss due to plasma-heating is observed at Jupiter’s moon Io and is
fundamentally driven by plasma ram pressure and magnetic pressure interacting
with the atmosphere (e.g., Johnson, 1990). Therefore, we estimate the
atmospheric loss from an impinging plasma on Mercury by scaling to the plasma
pressures measured at the Galilean satellites (Johnson, 2004). Following Oza
et al. (2019); Gebek & Oza (2020), the mass loss of a species $i$ by plasma
heating at proto-Mercury is:

$\begin{split}\dot{M}_{P,i}&=\frac{x_{i}}{\dot{M}_{i,\mathrm{Io}}}\,\hat{P}\,\hat{U}\,\hat{v}_{\mathrm{ion}}\,\left(\hat{h}_{\mathrm{exo}}\right)^{2}\end{split}$
(11)

where $x_{i}$ is the element fraction of the species $i$ in the atmosphere,
and $\dot{M}_{i,\mathrm{Io}}$ its atmospheric sputtering loss rate at Io. The
total plasma pressure $\hat{P}$, gravitational binding energy $\hat{U}$, and
ion velocity $\hat{v}_{\mathrm{ion}}$ of Mercury are expressed as non-
dimensional values that are scaled to Io’s corresponding values. The total
plasma pressure is additive where
$P_{\mathrm{tot}}=P_{\mathrm{mag}}+P_{\mathrm{ram}}$. For the calculations we
use a magnetic pressure $P_{mag}=1.7$ nPa based on an estimation of Mercury’s
magnetic moment of 2.76 $\times 10^{12}$ T m3 at the magnetopause stand-off
distance of 1.4 $R_{P}$. The ram pressure due to the solar wind varies from
$\sim$ 10–30 nPa (Korth et al., 2012), yielding a total pressure of $\sim$
12–32 nPa.

Evaporation from the magma ocean and atmospheric loss have to be equal to
retain the atmospheric pressure and thus a steady state. The evaporation rate
of a species is approximated by the Hertz-Knudsen-Langmuir equation. The
evaporation rate of a species $i$ with molar mass $M_{i}$ over the surface of
Mercury in mol/s is given by:

$\dot{M}_{\mathrm{evap},i}=M_{i}\frac{dn_{i}}{dt}=-4\pi
R_{P}^{2}M_{i}\frac{\gamma_{ev}P_{i,eq}-\gamma_{en}P_{i,s}}{\sqrt{2\pi
M_{i}RT}},$ (12)

with the evaporation and condensation coefficients $\gamma$ (set to unity for
a liquid), surface pressure $P_{s}$, equilibrium pressure $P_{eq}$, and the
Mercury radius $R_{P}$. By setting the homopause diffusion rate
$\dot{M}_{\mathrm{diff},i}$ equal to the evaporation rate
$\dot{M}_{\mathrm{evap},i}$, the equation is solved for the ratio of surface
to equilibrium pressure, $p_{i,s}/p_{i,eq}$ for each species at each
temperature step. For homopause diffusion rates, the ratio lies $>$ 0.99,
therefore the atmosphere up to the homopause is considered to be in
equilibrium.

##  3 RESULTS

###  3.1 Surficial melt lifetime and atmospheric structure

The surface temperature of the Hermean magma ocean cools from 2400 to 1500 K
in around 400 – 9000 years, depending on the planetary size (i.e. mantle mass)
and efficiency of radiative energy loss to space (Fig. 1). The cooling rate is
inversely proportional to $R_{P}$ since it depends on the ratio of the
planetary surface area to mantle mass. Hence, a large Mercury takes longer to
cool than a small Mercury for otherwise identical parameters.

####  3.1.1 Non-volatile cases

The cooling trajectory of non-volatile cases is characterised by two episodes.
First, when the surface temperature is high (early time), the pressure and
hence opacity of SiO is large and therefore cooling is slow. The second
episode of cooling is rapid, since even a small decrease in surface
temperature produces a drastic fall in both SiO pressure and opacity, driving
the planet towards cooling like an ideal black body. Therefore, the cooling
timescale for N5 cases is only marginally greater than for an ideal black
body, which bounds the minimum cooling time to 470 years. Increasing SiO
opacity by 2 orders of magnitude increases the minimum cooling time to around
1000 years (N3 cases), marginally affected by the chosen magma ocean
composition.

Non-volatile cases at 2400 K have thin atmospheres of $<0.1$ bar and comprise
$>99.8\%$ gaseous Na, SiO, Fe, K, and Mg. The major constituents are Na and
SiO at high temperatures whereas the mixing ratio of SiO rapidly decreases
below 2400 K (Fig. 2). The partial pressure of Mg behaves similar to SiO but
does not exceed $1\%$ of the atmospheric mixing ratio for any composition. The
mixing ratios of Fe and K, however, are more variable because they depend on
the assumed magma ocean composition and reach their highest mixing ratios of
$9\%$ and $2\%$, respectively, for CB and EH4 compositions with high FeO and
K2O (Table 2, Fig. 2). Refractory components—AlO and Ca—have negligible
partial pressures ($<10^{-6}$ bar at 2400 K) and are thus ignored in further
calculations. The highest total surface pressure of metal-bearing species at
low temperatures is obtained with the NSP melt composition. At 2400 K, the
total surface pressure is $6.16\times 10^{-2}$ bar, which decreases to
$8.52\times 10^{-6}$ bar at 1500 K.

Figure 3 shows the homopause levels of Na for the non-volatile cases as a
function of time for CB and NSP melt compositions. Exobase levels lie within a
few 100s of kilometers of the homopause and are omitted in the log-log plot as
a result. During the magma ocean phase, the levels evolve within 470–660 years
for N5 cases, and 1100–1480 years for N3 cases. The homopause and exobase
locations are only weakly sensitive to the planet size and gravity. The magma
ocean composition, however, exerts a strong influence on the atmospheric
structure. For the CB case, the homopause lies at 685 km whereas for the NSP
melt composition it lies around 1258 km at a magma ocean surface temperature
of 2400 K. The early inflation of an atmosphere above a cooling magma ocean is
due to increasing $T_{\mathrm{skin}}$ caused by decreasing IR-opacity as the
partial pressure of SiO decreases (Eq. 7). Following this stage, the homopause
altitude falls to 83 km (CB) and 439 km (NSP melt) at 1500 K. The exobase
density and location is further dependent on the mean cross section of
atmospheric species, which is tied to the composition-dependent vapor
pressures (Eq. 6). The high vapor pressure of Na in the NSP melt relative to
the CB composition lowers the mean molecular weight and the mean collision
cross section of the atmosphere, which both increase its extent.

Figure 1: Evolution of the surface temperature of the Hermean magma ocean. For
non-volatile cases, the upper bound of cooling is provided by EH4 composition
and the lower bound by Cb composition. See Table 1 for case parameters.

####  3.1.2 Volatile cases

Volatile-bearing cases result in cooling times of 3400 years (Case SV) and
8900 years (Case LV). Both small and large proto-Mercury have the same initial
volatile abundances of C and H by ppmw, but this manifests in a larger total
reservoir size of volatiles for a large proto-Mercury compared to a small one.
The mass of volatiles in the atmosphere defines the surface atmospheric
pressure, which in turn determines the optical thickness of the atmosphere and
hence the efficiency of radiative cooling.

Figure 2:  Oxide partial vapor pressures calculated using the CB, EH4 and NSP
source and NSP melt composition. The mixing ratios are thereby reflected by
the activities of the elements in the given melt composition. The high-Fe
composition CB shows a large Fe Partial pressure, comparable with K in the
high-K compositions EH4 and NSP melt.  Figure 3: Homopause altitudes for
sodium derived for the non-volatile cooling scenarios for CB and NSP melt.
Lines terminate when the surface temperature reaches 1500 K. The equilibrium
exobase levels are not plotted as they lie close to the homopause levels.
Homopause levels for other species show the same trends and reach comparable
elevations.

Atmospheres of volatile cases around a small and large proto-Mercury reach
surface pressures of about 5 and 12 bar at a magma ocean surface temperature
$T_{\mathrm{surf}}$ of 2000 K, respectively (Fig. 4). This result is
independent of the partial pressures of the metal-bearing species, as their
contribution is $\leq 0.1$ bar at $T_{\mathrm{surf}}$ = 2000 K. Thus, it is
the outgassed hydrogen and carbon species (which depends on their
solubilities) that dictates the surface pressure. VULCAN is then used to
compute the equilibrium chemistry of the atmosphere accounting for the
outgassed volatiles as well as the metals and oxides. For both small and large
Mercury, the atmosphere below the homopause is dominated by H2 and CO with
about 60 and 27 vol.%, respectively, at $T_{\mathrm{surf}}$ = 2000 K. Between
the homopause and the exobase the dominant H, C and O based species dissociate
to monoatomic gases.

For all compositions, Na and K are the dominant metallic elements at
$T_{\mathrm{surf}}=2000$ K. At the surface their hydroxide forms NaOH and KOH
are fairly abundant, but dissociate towards the homopause (Table 4) Sodium
hydride (NaH) is also present at the surface at pressures one order of
magnitude lower than NaOH and remains about constant throughout the
atmosphere, reaching similar mixing ratios to K. Potassium hydride (KH) is
ignored as no rate constant exists in the NIST kinetics database
(kinetics.nist.gov).

Exobase levels are situated up to 2910 and 2590 km for small and large Mercury
with homopause levels down to 2360 and 2160 km, respectively. The $P$–$T$
structure of the atmosphere gives $T_{\mathrm{skin}}$ values for the homopause
and exobase of about 1021 and 893 K (Table 4), respectively, for small and
large cases, which lies well below the $T_{\mathrm{skin}}$ calculated for the
non-volatile cases ($T_{\mathrm{skin}}=1615$ K). The exobase levels of the
volatile cases are thus comparable to the non-volatile cases. Unlike the non-
volatile cases, the planet size has a large impact on the atmospheric
structure, which is solely due to the difference in the total volatile
reservoir (Section 3.1).

Figure 4: Major and minor element composition of the atmosphere for NSP melt
composition with a surface temperature of 2000 K for (a) SV and (b) LV Mercury
models. The surface pressure that normalises the y-axis is 5.0 bar for SV and
12.1 bar for LV. The homopause level of H is plotted as a black, horizontal
line.

###  3.2 Atmospheric loss

We find that the major atmospheric escape mechanism is photoevaporation (Eq.
10) with a lower limit of photoevaporation constrained to 106.6 kg/s and an
upper-limit of 109.6 kg/s for both non-volatile (Na) and volatile cases (H, C,
and O). The upper limit is thereby roughly three orders of magnitude larger
than the photoionization of the major atmospheric species. For the high
heating efficiency ($\eta_{\mathrm{EUV}}=10^{-1}$) case, loss rates become
evaporation limited when reaching 1600 K as the surface to equilibrium
pressure ratio approaches zero. For low heating efficiencies
($\eta_{\mathrm{EUV}}=10^{-3}$), the ratio of surface to equilibrium pressure
$p_{i,s}/p_{i,eq}$ remains at $>$0.93 (Eq. 12). Photoevaporation as an
approximation of thermally-driven hydrodynamic escape (Eq. 10) therefore
expresses the highest uncertainty on the stability of the atmosphere.

Figure 5 shows the integrated mass loss over the most extensive surficial melt
lifetime of 8900 years. The photoevaporative erosion $dR$ of the surface can
be estimated by assuming mass conservation where:

$\frac{dM}{dR}=4\pi R_{P}^{2}\rho_{mantle}.$ (13)

Using a mantle density of $\rho_{mantle}$ = 3.5 g/cm3 and assuming a high EUV
heating efficiency of $10^{-1}$ at a large Mercury size allows for $\sim$ 2.3
km loss of crust over the 8900 years of the volatile case surficial melt
lifetime. We have shown using Eq. 12 that this case becomes evaporation
limited due to the high photoevaporation rates. The total integrated loss when
assuming an EUV heating efficiency of $10^{-1}$ is therefore lower than shown
in Figure 5 but not significantly. This is due to most loss occurring during
the early magma ocean stage, when the surface temperatures are high and
evaporation is not limiting the potentially high photoevaporation rates.

Figure 5: Integrated mass loss rates of the bulk atmosphere through
photoevaporation from 2400 K to 1500 K for large Mercury cases. The shaded
areas represent the uncertainty of loss rates based on initial conditions. For
the N3 case only the upper and lower limit are shown as dashed lines as it
lies within both V and N5 areas. The highest loss rates assume upper limits
for EUV luminosity only 1 Ma after Sun formation, LEUV(1 Ma), and heating
efficiency of $10^{-1}$, and the lowest assume LEUV(5 Ma), and heating
efficiency of $10^{-3}$. Using a mantle density of $\rho_{\mathrm{mantle}}$ =
3.5 g/cm3 results in a maximum of crust material being lost to space ranging
between 2.3 km and 16 cm, depending on the degree of XUV intensity and heating
efficiency.

Plasma-driven escape is diffusion-limited as a supply is required at the
exobase (Eq. 5). The $\dot{M}_{P}$ calculated lie orders of magnitudes below
$\dot{M}_{\mathrm{ion}}$ as given in (Table 5) and therefore do not affect the
atmospheric structure. Using an intermediate total pressure of 25 nPa, we find
about $\dot{M}_{P}$ = 103.4 kg/s for small and large, non-volatile Mercury
cases and about $\dot{M}_{P}$ = 103.1 kg/s for volatile Mercury cases,
respectively. The loss rates are thereby comparable to the plasma-driven
escape observed on Jupiter’s moon Io in SO_2 $\sim$ 103 kg/s (Thomas et al.,
2004).

Time-averaged mass loss rates by photoionization are given in (Table 5). Like
non-thermal plasma-driven escape, non-thermal escape due to photoionization at
the exobase (Eq. 9) is diffusion-limited. As photoionization rates of non-
volatile cases follow the same trends independent of planet size, we report
results focusing on a small proto-Mercury only, omitting the large Mercury
photoionization rates, which are mostly within a factor two for the dominant
species (Table 5). The results of non-volatile N5 cases are also not reported,
as an almost isothermal atmosphere results in less than a factor two larger
loss rates at high temperatures.

Regarding the presence of metal oxide derived gaseous species, Gibbs Free
Energy minimisation of the vapor phase (using FactSage) along a case-dependent
atmospheric $P$–$T$ profile (Figure 9, Appendix E) indicates that Mg and SiO
condense into clinopyroxene (1900 K) and then into olivine ($\approx$ 1700 K)
during cooling, by which temperature their fraction remaining in the gas is
negligible. Iron persists in the vapor to lower temperatures, condensing
partially into olivine before iron metal condenses at 1350 K. Therefore, while
Mg, Si (and Ca and Al) all condense prior to reaching the exobase ($T\approx$
1680 K), Fe is likely to partially reside in the vapor phase. Sodium never
fully condenses (Nepheline, its major host mineral, condenses in very minor
proportions below 1500 K), while K remains entirely in the vapor phase down to
at least 950 K.

####  3.2.1 Loss from non-volatile atmospheres

$\dot{M}_{\mathrm{ion}}$ in non-volatile cases is sensitive to the chosen
initial composition (Fig. 6). The loss fluxes of the non-volatile species of
interest —SiO, Na and K—are proportional to their mixing ratios in the
atmosphere (Fig. 2). In cases with high initial SiO partial pressures the
mixing ratios and hence the diffusion-limited loss rates of Na and K increase
during initial cooling as SiO becomes less abundant. This is most evident in
the CB loss flux with an initially increasing loss rate despite decreasing
temperatures (Fig. 6).

Figure 6:  Mass loss fluxes of exospheric species from 2400 K to 1550 K
plotted for the small Mercury N3 cases. The loss fluxes for large Mercury
follow the same trends and are about a factor of two smaller. The N5 cases
with isothermal atmospheres show a factor of two larger loss rates at high
temperatures but the same species related trends.

In the low-Na and low-K composition CB, where SiO is the dominant metal oxide
at high temperatures, loss fluxes of SiO reach up to $9.6\times 10^{5}$ kg/s
at $T_{\mathrm{surface}}$ = 2400 K. The vapor pressure of SiO declines with
respect to other dominant gas species (Na and K), thereby reducing its mixing
ratio rapidly with decreasing temperature. The lower mixing ratio of SiO
causes loss rates to drop to $9.4\times 10^{2}$ kg/s at 1500 K. For the same
composition, loss rates for Na and K are around $2.0\times 10^{5}$ kg/s and
$3.1\times 10^{3}$ kg/s at $T_{\mathrm{surf}}$ = 2400 K, increasing to
$3.2\times 10^{5}$ kg/s and $4.2\times 10^{3}$ kg/s at $T_{\mathrm{surf}}$ =
1500 K, respectively. For the high-Na end-member composition of NSP melt, the
diffusion-limited loss rates for Na and K are $5.9\times 10^{5}$ kg/s and
$1.2\times 10^{4}$ kg/s respectively when at $T_{\mathrm{surf}}$ = 2400 K, and
decrease to $3.9\times 10^{5}$ kg/s and $4.3\times 10^{3}$ kg/s respectively
when $T_{\mathrm{surf}}$ = 1500 K. The ratio of Na to K diffusion rates is
$\sim$ 100 and hence about one order of magnitude higher than their mixing
ratios in the atmosphere. Integrated diffusion-limited losses over the
surficial melt lifetimes are shown in Figure 7.

Figure 7: Integrated mass loss rates of exospheric oxides and elements through
photoionization as the magma ocean cools from a surface temperature of 2400 K
to 1500 K. SiO loss is shown although it is unlikely to persist in the upper
atmosphere as it is consistently below its highest condensation temperature of
$T=1900$.

If we assume based on FactSage results, that SiO, and Mg are absent and only
Na, K and Fe remain in the atmosphere, the diffusion-limited loss rates for
all cases peak in the small Mercury NSPm case with a fairly temperature
independent Na loss rate of about $10^{6}$ kg/s. Loss rates for K and Fe are
thereby about three orders of magnitude lower than that of Na and do not
significantly contribute to the total loss. The difference in change of
$\dot{M}_{\mathrm{diff}}$ and thus $\dot{M}_{\mathrm{ion}}$ between species
with continued magma ocean cooling results from the constantly dropping total
surface pressure simultaneously to shifting partial pressures. At low magma
ocean surface temperatures, Na exerts most of the pressure, whereas at high
temperatures condensing species, such as SiO, are the predominant contributors
to the total pressure and mean molecular mass of the atmosphere (Fig. 2). A
lower surface pressure results in a lower homopause level and leads to a
smaller diffusion area and rate, however the rapid decline of the SiO partial
pressure in the atmosphere ($n_{\mathrm{SiO}}/n_{\mathrm{hom}}$) results in
the sharp drop of SiO homopause diffusion rates (i.e. Eq 5) but increases the
$T_{\mathrm{skin}}$ and therefore the extent of the atmosphere (Fig. 3).
Similarly, the loss rates through ionization at the exobase for a relatively
low mean molecular weight (‘light’) Na, K, and Fe based atmosphere are up to a
factor five higher than the $2\times 10^{5}$ kg/s loss of Na from a ’heavy’
atmosphere, which includes Mg and notably SiO.

These limits are significant for atmospheric escape estimates by plasma
heating, photoionization or Jeans escape. All three of those processes are
calculated from the exobase, which is for plasma heating and photoionization
where ions and photons can access a rarefied neutral atmosphere. This is
therefore diffusion-limited as a high flux is required to source the neutral
species experiencing a momentum transfer from the plasma. Photoevaporation is
not necessarily diffusion-limited so long as a sufficiently large column
exists at the altitude where EUV photons are able to absorb on to infrared
emitting molecules ($z_{\mathrm{abs}}$, Fig. 8). In N2/CH4 atmospheres (e.g.
Kuiper Belt Objects) the critical column density is estimated to be
$\gtrsim\,$1018/cm2 (Johnson et al., 2015) which is easily achieved at a
fiducial absorption altitude of $z_{\mathrm{abs}}$ situated at $\sim$ 1.25
$R_{P}$, where the column density is equivalent to $\sim$1021/cm2 for an
isothermal scale height of H $\approx$ 150 km at a $T_{\mathrm{surf}}$ = 2000
K. For a magma–silicate atmosphere, as studied here, SiO or a similar species
would be able to re-emit in the infrared resulting in upper atmospheric
expansion, and Roche-lobe overflow to space.

####  3.2.2 Loss from volatile atmosphere

In the volatile cases, assuming a speciation as encountered at the H
homopause, the diffusion limited loss fluxes of primary species lie between
$10^{4}$ kg/s and $10^{5}$ kg/s for all major species (H2, CO, H2O and CO2) in
the small and large proto-Mercury cases.

The loss flux of minor species Na, and K are several orders of magnitude lower
than those of the major species, at $10^{1}$, and $10^{0}$ kg/s respectively.
In the non-volatile cases, loss fluxes of Na and K are directly proportional
to their thermodynamic activities in the melt. Sodium activity increases by
about a factor 4.5 from the CB to NSP melt composition and K by a factor two
between NSP source and EH4, respectively. Relative to non-volatile cases, loss
fluxes for Na and K are several orders of magnitude lower in the high
pressure, volatile-rich atmosphere at $T_{\mathrm{surf}}$ = 2000 K.

If we assume the loss fluxes of the dominant H, C, and O based species at
$T_{\mathrm{surf}}$ = 2000 K to be constant over the lifetime of the molten
surface (Section 3.1) and integrate them for small and large Mercury volatile
cases we obtain a total mass loss by photoionization of $4.1\times 10^{16}$ kg
and $1.8\times 10^{17}$ kg, respectively. This exceeds the total
photoionization mass loss from the low absorbing N5 non-volatile case by only
about one order of magnitude (Fig. 7). The mass loss of Na in the volatile
cases, however, only contributes about $10^{12}$ kg of the total, which is
about four orders of magnitude below the total mass loss of the non-volatile
cases. Again, this assumes that the Na loss flux is constant in the volatile
case. This is deemed appropriate, because Na is only a minor component of such
atmospheres, is lost at slow rates that represent an insignificant fraction of
its total budget, and does not condense before reaching
$T_{\mathrm{surf}}=1500$ K.

##  4 Discussion

###  4.1 Mass loss of proto-Mercury

Table 5 tabulates the total atmospheric loss rates due to the following escape
mechanisms: ionization $\dot{M}_{\mathrm{ion}}$, photoevaporation
$\dot{M}_{U}$, and plasma-heating $\dot{M}_{P}$. Figure 8 illustrates the
atmospheric level from where the degassed magma ocean atmosphere is escaping.

Table 5: Time-averaged mass loss rates in kg/s.

Mass loss [log10]
---
Process | Size | Emissivity
|  | N5, N3 | V
$\dot{M}_{P}$ | S | 3.4 | 3.1
L | 3.4 | 3.2
$\dot{M}_{\mathrm{ion}}$ ($\dot{M}_{\mathrm{diff},i}$) | S | 5.6 | 5.6
L | 5.2 | 5.8
$\dot{M}_{U}$ |  | LEUV
|  $t=1$ Myr |  $t=5$ Myr
$\eta_{\mathrm{EUV}}(10^{-3})$ | S | 7.5 | 6.6
L | 7.6 | 6.7
$\eta_{\mathrm{EUV}}(10^{-1})$ | S | 9.5 | 8.6
L | 9.6 | 8.7

Note. — Loss rates of plasma heating $\dot{M}_{P}$, photoionization
$\dot{M}_{\mathrm{ion}}$, and photoevaporation $\dot{M}_{U}$. Photoevaporation
is insensitive to the atmosphere’s emissivity and composition but depends on
the EUV flux and the EUV heating efficiency (end members of $10^{-3}$ and
$10^{-1}$). The EUV flux is a function of the age of the solar system (Eq. 8,
after Johnstone et al., 2015; Tu et al., 2015).

Figure 8: Mass loss processes and their rates demonstrate the coupling between
various atmospheric layers. $R_{P}$ is the planet radius, $z_{\mathrm{exo}}$
is the exobase altitude, and $z_{\mathrm{hom}}$ is the homopause altitude
which governs exospheric loss processes of Jeans escape, plasma heating
$\dot{M}_{P}$, and photo-ionization $\dot{M}_{\mathrm{ion}}$.
$z_{\mathrm{abs}}$ is the absorption altitude where upper-atmospheric heating
(photoevaporation) $\dot{M}_{U}$ commences. The absorption altitude is assumed
to lie below the homopause and therefore photoevaporation is not limited by
homopause diffusion.

We find that the loss fluxes from the exobase caused by photoionization and
atmospheric sputtering are supply-limited (in all volatile and non-volatile
cases) by homopause diffusion $\dot{M}_{\mathrm{diff}}$, which dictates the
exobase elevation in order to remain in steady state. The maximum
$\dot{M}_{\mathrm{ion}}$ of volatile and non-volatile cases are comparable,
even though their atmospheres are comprised of different major species (non-
volatile case: Na, and SiO, volatile case: H2, H2O, CO and CO2). This
similarity is attributed to the higher $T_{\mathrm{skin}}$ of the non-volatile
cases at $T_{\mathrm{surf}}=2000$ K, caused by IR opacity which is tied to the
mixing ratio of SiO. SiO pressure rapidly decreases with decreasing
temperatures, which is contrary to CO_2 and H_2O in the volatile cases. A
higher skin temperature in the non-volatile cases hence compensates for the
lower atmospheric pressures.

The photoevaporation rate is limited by the degree of upper atmospheric
heating efficiency, $\eta_{EUV}$, rather than by $\dot{M}_{\mathrm{diff}}$ ,
as well as by the supply of gases from surface evaporation (eq. 12). The
assumption of a constant photoevaporation rate is thus only valid for
evaporation at high temperatures above 1600 K at which evaporation rates are
fast enough for supply to be sustained, or for moderate mass loss rates of
about 107 kg/s. The majority of mass loss occurs at high temperatures (Fig. 5)
when surface evaporation rates are high compared to photoevaporative loss
rates, the later of which are independent of temperature and depend instead on
the EUV flux.

The diffusion-limited loss rates by photoionization of the four major volatile
species: H2, H2O, CO and CO2 from a volatile-rich atmosphere total $\sim
10^{5}$ kg/s. The loss of Na from a thick, volatile rich atmosphere is
inhibited by its low mixing ratios at the homopause and exobase. Therefore,
diffusion-limited loss of Na is most efficient when the atmosphere is thin,
reaching a few $10^{5}$ kg/s, which coincides with the total mass loss rates
from volatile cases. The total integrated mass loss by photoionization from
Mercury’s exosphere is low for small volatile and non-volatile cases with
$\leq\,4.1\times 10^{16}$ kg and $\leq\,3.0\times 10^{16}$ kg, respectively.

The mass loss of single species is negligible compared to the total inventory
of the magma ocean reservoir. For example, 0.033 wt% H2O and a low estimate of
0.1 wt% Na in a total mass of M${}_{\mathrm{MO}}\,\approx\,\times 10^{23}$ kg
corresponds to a reduction of the total H2O reservoir mass (volatile cases)
and Na (non-volatile cases) by $\leq 0.02\%$. Assuming a well-mixed mantle
reservoir, the bulk composition of Mercury would not significantly change even
for species with low abundance in the reservoir and large loss rates such as
H2O, CO2, and Na. Energy-limited escape via photoevaporation, however, can
erode up to 2.3 km of Mercury’s crust which coincides with 0.3% of small
Mercury (Figure 5). Assuming small heating efficiencies as well as a lower EUV
flux, leads to integrated photoevaporation losses and eroded crust thicknesses
that are reduced by four orders of magnitude.

Physical segregation between crystal and liquid during magma ocean cooling
will induce chemical fractionation of element abundances with respect to those
of the bulk mantle. Namely, the incompatible lithophile elements (Na, K, Al,
and Ca) become enriched in late-stage liquids of a Hermean magma ocean. This
effect is simulated by considering the composition of the NSP melt as a
surface magma ocean analogue relative to that of its inferred source. These
differences notwithstanding, the partial pressures of metal-bearing gas
species vary only marginally among EH4, CB and NSP compositions. This is due
to two factors: 1) vapor pressures of different elements vary by orders of
magnitude among one another (e.g. compare Na with AlO) whereas abundances of
these major elements vary only by a factor of 2–3 in most cases; and 2) higher
mole fractions of Na and K in the NSP melt are partially compensated by their
lower activity coefficients relative to the NSP source or EH4 composition. As
we show in Section 3, all elements other than Na and K, and potentially Fe
condense before reaching the exobase. We can therefore conclude that the
atmospheric pressure and speciation around proto-Mercury only depends on the
abundance of volatile- and moderately volatile elements.

####  4.1.1 Early origin of surface Na

In order to determine the potential impact of a magma ocean-generated
atmosphere on the surface composition of a small proto-Mercury, we calculate
the total mass of Na in the atmosphere. We consider a hypothetical scenario in
which the atmosphere collapses as soon as the first crust forms, coinciding
with the termination of the magma ocean stage when a surface temperature of
1500 K is reached. To obtain a result which is consistent with the notion of
Na-poor building blocks (e.g., Humayun & Cassen, 2000), we use the low-Na CB
composition (Table 2) and its H-, and C-absent, pure Na atmosphere
composition. The resulting total amount of Na integrated over the whole
atmosphere yields about $10^{11}$ kg.

Whether it precipitates as Na metal or as another compound depends on the
composition of the atmosphere that exists. Although not considered in our
model, such an atmosphere would contain significant quantities of other
moderately volatile elements that could combine with Na to form complex
molecules, namely, F, Cl and S. The species NaCl is inferred to be stable
among volcanic gases (e.g., Aiuppa et al., 2003) and has been directly
observed in Io’s atmosphere (Lellouch et al., 2003; Moullet et al., 2010), and
may therefore be a potential candidate to form surficial deposits. Sodium
chloride is also observed as a stable precipitate from experimentally-
generated volcanic gas analogues (Renggli & Klemme, 2020), and is therefore
likely to occur as a Na-bearing phase on the Mercurian surface. This is
supported by the coinciding distribution pattern of Na and Cl from more recent
volcanic deposits found in gamma ray spectrometer data (Evans et al., 2015).

Here we consider a simplified case, for which the mass of sodium in the
atmosphere is uniformly distributed over the surface of Mercury as pure, low
density, Na metal, results in a layer less than 1 mm for the CB case. Using a
more Na-rich composition like EH4 combined with an increased atmospheric
reservoir size of large proto-Mercury would lead to a factor four thicker Na
layer, but still less than 1 mm. In the volatile cases (i.e., with CO2 and
H2O), the amount of Na in the atmosphere is identical to the volatile-free
cases, as, in our model, the partial pressure of Na is independent of the
presence of volatiles. The small dissolved quantities of CO2 and H2O in the
silicate melt ($\leq$ 1000 ppm) should thereby not influence the activity
coefficients of the major rock-forming species. This hypothetical Na metal
layer would not outlast meteorite impacts, which are assumed to have removed
50 m to 10 km of early crust (Hyodo et al., 2021). For the enrichment to be
preserved, the atmospheric sodium would have to be incorporated into a layer
with a thickness exceeding the removed crust. However, for a minimum layer of
50 m we obtain a total Na wt% increase of merely $\sim$1 ppm and $\sim$10 ppb
for small and large proto-Mercury cases, respectively. We thus conclude, that
the collapse of an early Na-rich atmosphere would not contribute to a notable
increase of Na in the surface.

###  4.2 Controls on mass loss

The mean column density at the exobase depends on the weighted average of the
dominant species’ collision cross sections (CCSs). Loss rates are directly
related to the exobase density. However, using CCSs from Kim & Desclaux (2002)
that are about one order of magnitude smaller would reduce the homopause and
exobase levels by a few tens of kilometers and decreases the homopause
diffusion limited loss by $\leq 2\%$. The sensitivity of mass loss to the
chosen CCSs is therefore weak. In non-volatile cases, if we consider that all
species except Na, K, and Fe condense (FactSage in Section 3), then mean
molecular mass and the collision cross section of the atmosphere decreases
which enhances molecular diffusion (Eq. 4). This pushes the homopause and
therefore the exobase further from the planet surface, increasing the
atmospheric surface area and therefore loss. Furthermore, the absence of SiO
leads to a hotter skin temperature as the atmosphere becomes IR transparent,
further enhancing loss. The difference of the ionization mass loss rate at the
exobase between a Na, K, and Fe atmosphere and an atmosphere where SiO is a
major component at high temperature is thereby about a factor three larger for
all cases.

For rocky exoplanets on short orbits, the atmospheric temperature around our
homopause levels ($10^{-7}$ bar) can be as high as 3800 K for a surface
temperature of 2400 K (Ito et al., 2015). Mercury posesses different planet
parameters (1 $M\mathrm{Earth}$ and 0.02 AU vs. 0.055 $M\mathrm{Earth}$ and
0.3 AU for Mercury), however, the more intense early UV flux experienced by
Mercury could similarly boost the temperature at the homopause. Calculations
with an increased skin temperature of 3800 K at the homopause resulted in
about a factor two higher photoionization loss rates for all cases.

We used photoevaporation as a proxy for thermally driven hydrodynamic escape.
Krenn et al. (2021) has shown for a large range of parameters that
photoevaporation can underestimate hydrodynamic escape especially at low EUV
fluxes. Given our large incident EUV fluxes of about $10^{2}-10^{3}$ Js-1m-2
and our escape parameters (Table 4) we expect to be within one order of
magnitude of thermally-driven hydrodynamic escape rates (compare EUV fluxes
and escape parameters to Fig. 4 in Krenn et al., 2021, although our escape
parameters for small Mercury lie just below the shown range).

###  4.3 Atmospheric evolution and structure

Figure 8 illustrates the transport of mass away from different levels in the
atmosphere. In our model, atmospheric escape can be either energy-limited
(e.g. $\dot{M}_{U}$) or diffusion-limited (Jeans escape, $\dot{M}_{P}$,
$\dot{M}_{\mathrm{ion}}$). Below we describe the role of enhanced atmospheric
heating or cooling on diffusion and energy-limited escape.

Diffusion rates are tied to the homopause density, which determines the
homopause altitude. The eddy diffusion coefficient ($K_{zz}$) needed to
determine $n_{\mathrm{hom}}$ bears large uncertainties, however. For all
cases, a $K_{zz}$ larger than the Earth-derived upper limit of $3.2\times
10^{6}$ cm2/s would most likely be adequate to accommodate proto-Mercury’s
increased atmospheric temperature, increasing $z_{\mathrm{hom}}$ and lowering
$n_{\mathrm{hom}}$, leading to a slightly larger diffusion and therefore loss
rate. Even if we assume a larger $K_{zz}$, however, homopause diffusion will
remain the limiting factor for mass loss. We find for volatile cases, that
even if $K_{zz}$ is three orders of magnitude larger, the total loss for
volatile cases increases by a factor of less than two. The sensitivity of
$\dot{M}_{\mathrm{diff}}$ and therefore $\dot{M}_{\mathrm{ion}}$ to the eddy
diffusion coefficient is therefore weak.

Ionization could further increase the exobase temperature, and hence the
reported diffusion-limited loss fluxes could be a lower limit. Whether mass
loss occurs from the exobase surface, or whether it is the result of an
advective outflow, is canonically assessed by the escape parameter (e.g. Genda
& Abe, 2003). If the escape parameter $\lambda_{0}\leq 3$ the atmosphere
experiences mass outflow due to its non-zero net velocity, and escape occurs
inward of the exobase at the sonic point (where the thermal velocity exceeds
the sound speed). If $\lambda_{0}\gg 3$ the atmosphere escapes because the
mean free path is longer than the scale height, and Jeans escape prevails.
Table 4 shows how the exospheric escape parameters are all
$4\leq\lambda_{0}\leq 15$, a near-transitional escape regime between Jeans-
and hydrodynamic end members, which was recently determined to be relevant for
the putative magma ocean on the Moon (Tucker et al., 2021). These authors
demonstrated via DSMC simulations (Bird, 1994) that cooling due to escape is
important for $\lambda_{0}\lesssim 15$. Therefore in Table 4 based on our
escape parameters, it would appear that although ionization may further
enhance escape, cooling may temper this loss. In addition, the significant
ionization rates of $\dot{M}_{\mathrm{ion}}\simeq$ 106 kg/s at the semi-major
axis of proto-Mercury promotes the generation of an ionosphere that is
modulated by the planetary magnetic field. Simulations on an early Mars
analogue have demonstrated that ion escape is efficient at removing material
(Egan et al., 2019). Therefore, it is possible that we are underestimating
escape by not considering magnetic interactions.

In the concurrent ’energy-limited’ regime it would appear that if EUV photons
are able to absorb on to a sufficiently high flux of molecules (Section
3.2.1), heating would overwhelm cooling. However, based on the escape
parameters in Table 4 it appears that cooling associated with escape may be
important, arresting loss. At the same time, the study of low-mass, close-in
exoplanets orbiting Sun-like stars has posited the idea that low-mass planets
are nevertheless born with hydrogen/helium (H/He) envelopes although these are
rapidly lost due to photoevaporation (Mordasini, 2020). For a H/He envelope
equivalent to $1\%$ the mass of proto-Mercury, we find that our upper-limit on
photoevaporation results in the dissipation of a H/He envelope in 104.4 years,
which is larger than the lifetime of the molten surface. The possibility of a
H/He envelope to persist during the molten surface lifetime is therefore non-
trivial, and could result in significant heating which could not only enhance
escape but also elongate the melt lifetime past the 104 years we study here.

###  4.4 Origin and evolution of Mercury

The elevated core:mantle ratio, coupled with an Na- and S-rich surface,
distinguish Mercury from the other terrestrial planets. Two key hypotheses
exist to account for these characteristics; (1) the preferential loss of
silicate material, either by evaporation (Fegley & Cameron, 1987) or by
collisional stripping (Benz et al., 1988) and (2) equilibrium condensation and
sorting of metal from silicate in the solar nebula (Lewis, 1972;
Weidenschilling, 1978).
In evaluating hypothesis (1), Fegley & Cameron (1987) concluded that $\sim$
75–79% of silicate material would need to be lost during a fractional
vaporisation hypothesis to reproduce the core:mantle ratio of present-day
Mercury. In this work, we show that such high fractions of loss of silicate
material are untenable, be it from a small or a large proto-Mercury (total
mass losses are below 0.3%). The principal reason is that atmospheric cooling
timescales are too rapid with respect to evaporation and escape timescales,
meaning that integrated loss rates over $\sim$ 104 yr are small with respect
to the mass of proto-Mercury. Moreover, substantial amounts of atmospheric- or
collisional escape of Mercury’s crust is not represented in the high K/U ratio
of its surface (McCubbin et al., 2012), as preferential loss of silicate
material will predominantly deplete its incompatible lithophile element budget
(O’Neill & Palme, 2008).

There are several caveats to our conclusions, namely, that our results are
valid for dry or C-, and H-bearing atmospheres, but do not consider the effect
of other minor volatiles (Cl, S, F) on the volatility behaviour of metals.
Metal chlorides and metal sulfides may be important gaseous species under
moderate temperatures ($\sim$ 1000 K, Renggli et al., 2017), increasing their
volatility. Secondly, conditions on the surface of Mercury may have been
considerably more reduced than modelled herein (IW-5; Cartier & Wood, 2019).
Because the partial pressures of most metal-bearing species increase with
decreasing $f$O2 (Eq. 1) vaporisation rates for alkali metals may be an order
of magnitude higher (considering that their exponent, $n=1/4$; Sossi et al.,
2019).

However, these faster evaporation rates may be offset by the presence of a
surficial graphite layer on the magma ocean (Keppler & Golabek, 2019). Such a
layer is promoted under reducing conditions as the solubility of C in silicate
melt decreases from $\sim$ 360 ppm at the IW buffer to 1 ppm at IW-4 (Duncan
et al., 2017; Keppler & Golabek, 2019). The extent of a graphite layer
therefore depends on the C content of Mercury and its $f$O2, both of which are
poorly known. A surficial lid would additionally delay cooling of the mantle,
unless the lid is regularly broken as possibly occurred for the flotation
crust on the moon (Perera et al., 2018). Nevertheless, the net effect of a
graphite lid on Mercury’s magma ocean would be to reduce the extent of
degassing calculated herein. Therefore, we conclude that the physico-chemical
characteristics of Mercury cannot have been produced during a magma ocean
stage on a near fully-grown planet.

These obstacles are ameliorated when considering vapour loss from planetary
building blocks. Should Mercury have accreted from smaller, km-size
planetesimals, then melting and vaporisation on the precursor bodies would
have led to more efficient mass loss (e.g., Hin et al., 2017). Thus,
vaporisation may still be a physically viable mechanism to explain Mercury’s
composition, provided it occurred on its precursor bodies. However, another
problem arises because moderately volatile elements, such as Na, S and K, are
always more volatile (i.e., their partial pressures are higher for a given
activity) than the major mantle components, such as Mg and Si (Sossi et al.,
2019). Moreover, as demonstrated herein, Na is more easily lost with respect
to Mg and Si due to its lower molar mass and higher tendency to remain in the
gas phase in an adiabatically-expanding atmosphere (3.2.1). As such, appealing
to evaporative loss of Mg and Si to increase the core:mantle ratio while
retaining Na and K is inconsistent with evaporation from a silicate melt on
small planetary bodies. Therefore, other hypotheses should be considered.

##  5 Conclusions

We combined chemical and thermodynamic equilibrium models of the thermal
evolution of Mercury’s magma ocean and gaseous species derived thereof, to
model the thermochemical evolution of an early atmosphere on Mercury. For an
initially large Hermean mantle with initial C and H budgets comparable to
those of other rocky planets, namely Earth (‘volatile cases’), the lifetime of
surficial melt may have reached almost 104 years. Compared to a present-day
sized proto-Mercury without a greenhouse atmosphere, this lifetime is an order
of magnitude larger and thereby may enable early atmospheric mass loss to
occur over an extended duration. Cases with C and H show that Mercury could
have started with a 5–12 bar atmosphere. By contrast, excluding the presence
of C and H species results in a thin, short-lived metal- and metal oxide-
bearing atmosphere. The upper atmospheres of volatile cases are dominated by
H2, and CO whereas non-volatile cases are mostly Na and SiO.

Photoionization is a minor exospheric loss mechanism, limited by homopause
diffusion ($\dot{M}_{\mathrm{diff}}$) up to a maximum of a few $10^{5}$ kg/s.
If C and H volatiles are absent from the atmosphere, the
$\dot{M}_{\mathrm{diff}}$ limit applies to SiO and Na. Mass loss rates via
photoevaporation, $\dot{M}_{U}\leq$109.5 kg/s, exceed those from all other
known mechanisms due to the high EUV luminosity of the early Sun. This could
in the best case scenario erode an equivalent thickness of up to 2.3 km of
proto-Mercury’s crust when assuming high EUV heating efficiencies of
$10^{-1}$. Atmospheric sputtering $\dot{M}_{U}\sim$ 103.4 kg/s (also limited
by $\dot{M}_{\mathrm{diff}}$) occurs at the exobase, knocking-off neutral gas
molecules due to the ram pressure of the solar wind.

By integrating atmospheric loss rates over surficial melt lifetimes, we
bracket the expected total mass loss from Mercury’s early atmosphere. Based on
photoionization, Jeans escape, and plasma heating, the evaporation and loss of
the magma ocean of proto-Mercury did not significantly modify its bulk
composition. This is because magma ocean cooling times are too short to drive
substantial total loss for the determined atmospheric loss fluxes.
Photoevaporation can remove an equivalent crustal thickness of up to 2.3 km in
about 10,000 years, which is approximately $\sim$ 1020 kg of material.
Integrated losses of even the most volatile elements considered here, Na and
K, are insignificant with respect to their total budgets when escape is
diffusion-limited ($\leq 0.02\%$ decrease of the initial Na composition which
would be a difference of $3\times 10^{-4}$ wt%). Hence the present Na-rich
surface composition may indicate that catastrophic volatile loss during the
magma ocean stage did not occur, and that Mercury’s peculiar composition is
inherited from that of the solar-proximal region of the nebula from which it
accreted.

Financial support has been provided by the Swiss National Science Foundation
(SNSF) Fund (200021L182771/1). D. J. Bower acknowledges SNSF Ambizione Grant
173992. P. A. Sossi was supported by SNSF Ambizione grant 180025 and A. Wolf
by the National Science Foundation EAR 1725025 as well as the Turner
Postdoctoral Fellowship. Thanks to S. Suriano, P. Saxena, A. Heays, S.-M.
Tsai, and N Ligterink for discussions relating to this work. Part of this work
was conducted at the Jet Propulsion Laboratory, California Institute of
Technology, under contract with NASA.

##  Appendix A Magma ocean model

The evolving surface temperature of the Hermean magma ocean is calculated
using the SPIDER code, which is described in detail in Bower et al. (2018,
2019, 2021). Table 6 shows the parameters used to model proto-Mercury. The
mass absorption coefficients of H and C volatile species are determined at
1.01 bar, and the coefficients of SiO at $3\times 10^{-6}$ bar.

Table 6: Standard parameters for magma ocean cases.

Parameter | Value | Units
---|---|---
Core heat capacity | 850 | J kg-1 K-1
Core density | 7200 | kg m-3
Core radius | 2000 | km
Equilibrium temperature, $T_{\infty}$ | 440 | K
Gravity, $g$ | Table 1 | m s-2
Planetary radius, $R_{P}$ | Table 1 | km
Boundary layer scaling, b | $10^{-7}$ | K-2
Al abundance†† | 19500 | ppmw
26Al/Al (zero time) | $5.25\times 10^{-5}$ | —
K abundance† | 403 | ppmw
40K/K (present time) | $1.17\times 10^{-4}$ | —
Th abundance† | 49 | ppbw
232Th/Th (present time) | 1 | —
U abundance† | 28 | ppbw
235U/U (present time) | 0.007 | —
238U/U (present time) | 0.993 | —
H_2 mass absorption (CIA) | $5\times 10^{-5}$ | m2 kg-1
H_2 solubility law | ∗
H_2O mass absorption | $10^{-2}$ | m2 kg-1
H_2O solubility law | ∗∗
CO mass absorption | $10^{-5}$ | m2 kg-1
CO solubility law | ∗
CO_2 mass absorption | $10^{-4}$ | m2 kg-1
CO_2 solubility law | ∗∗
SiO mass absorption (large) | $10^{-3}$ | m2 kg-1
SiO mass absorption (small) | $10^{-5}$ | m2 kg-1
Initial surface temp | 2400‡ | K

${\dagger}{\dagger}$ Average Al abundance based on the composition of EH4 and
NSP source (Table 2).
${\dagger}$ Average current estimates for bulk heat source from Tosi et al.
(2013) and natural abundances from Ruedas (2017).
${\ddagger}$ Similar to maximum temperature estimate of Mercury’s surface
during accretion and differentiation (Bhatia & Sahijpal, 2017).
$*$ Lichtenberg et al. (2021).
${**}$ Bower et al. (2019).

##  Appendix B VapoRock species

The species included in VapoRock are given in Table 7.

Table 7: Species included in VapoRock (Wolf et al., 2021). Species
---
Al |  AlO |  AlO2 |  Al2 |  Al2O |  Al2O2
Si |  SiO |  SiO2 |  Si2 |  Si2O2 |  Si3
K |  KO |  KO2 |  K2 |  K2O |
Na |  NaO |  Na2 |  Na2O |  |
Mg |  MgO |  Mg2 |  |  |
Ca |  CaO |  Ca2 |  |  |
Fe |  FeO |  |  |  |
O |  O2 |  |  |  |

##  Appendix C Modified VULCAN

We incorporated Na, Si, Mg, K, Fe and their derivatives into VULCAN by adding
12 reactions (from kinetics.nist.gov) to the preexisting chemistry network
based on C, H and O (Table 8). We initially added more reactions, but removed
those that had a negligible impact on the resulting atmospheric speciation
when omitted.

Table 8: Key reactions added to VULCAN.

Reaction
---
OH + SiO |  -¿ |  SiO2 + H
OH + Si |  -¿ |  SiO + H
Si + O2 |  -¿ |  SiO + O
NaO + O |  -¿ |  Na + O2
Na + H2O |  -¿ |  NaOH + H
H2O + NaO |  -¿ |  NaOH + OH
H2 + NaO |  -¿ |  NaOH + H
HCO + Na |  -¿ |  CO + NaH
Mg + O2 |  -¿ |  MgO + O
H2O + KO |  -¿ |  KOH + OH
CO2 + Fe |  -¿ |  CO + FeO
3-body reactions
OH + K + M |  -¿ |  KOH + M
Na + O2 + M |  -¿ |  NaO2 + M
NaOH + M |  -¿ |  OH + Na + M
FeO + H2O + M |  -¿ |  Fe(OH)2 + M

Note. — The reactions given affect the speciation of Si, Mg, Fe, Na and K and
Si between $T$ = 2000 - 873 K and $P$ = 11.7 - 10-7 bar.

##  Appendix D Collision cross sections

The CCSs are shown in Table 9, which were approximated by the circular area of
radius equal to the atom or bond length. Furthermore, all bonds were
approximated to be covalent.

Table 9: Species cross sections (CS) used

Species | CS [$\text{\AA}^{2}$] |  | Species | CS [$\text{\AA}^{2}$]
---|---|---|---|---
H | 0.88 |  | CO2 | 10.3 |
H2 | 1.29 |  | KOH | 19.6 |
H2O | 2.84 |  | Na | 11.3 |
O | 1.13 |  | K | 18.6 |
O2 | 4.08 |  | SiO | 8.45 |
C | 1.41 |  | Mg | 6.61 |
CO | 4.01 |  | Fe | 7.65 |

Note. — CS values are based on sizes of atomic, single, double, and triple
bond data (Clementi et al., 1967; Pyykkö & Atsumi, 2009a, b; Pyykkö et al.,
2005).

##  Appendix E Atmospheric P–T profile

SPIDER determines an atmospheric pressure–temperature profile through an
analytical solution to the radiative transfer equations (see Appendix in Abe &
Matsui (1985) and Sect. 3.7.2 in Andrews (2010)). The solution gives rise to
the skin temperature equation (Eq. 7). Fig. 9 shows the volatile (V) and non-
volatile (N3, N5) atmospheric pressure–temperature profiles that are used for
FactSage and VULCAN calculations. An unphysical outcome of assuming only
radiative equilibrium (no convection) is a temperature discontinuity between
the base of the atmosphere and the surface of the magma ocean, which is
visually more evident for the non-volatile cases that have a small optical
depth. Nevertheless, for all cases the surface temperature is 2000 K.

Figure 9: Atmospheric $P$–$T$ profiles of the volatile (V) and non-volatile
(N3, N5) atmospheres surrounding a small (S) and large Mercury. The LN3 and
LN5 profiles are omitted as they are visually indistinguishable from the
respective S case profiles.

## References

  * Abe & Matsui (1985) Abe, Y., & Matsui, T. 1985, J. Geophys. Res. Solid Earth, 90, C545, doi: 10.1029/JB090iS02p0C545
  * Aiuppa et al. (2003) Aiuppa, A., Dongarrà, G., Valenza, M., Federico, C., & Pecoraino, G. 2003, in Geophysical Monograph Series, Vol. 139 (Blackwell Publishing Ltd), 41–54, doi: 10.1029/139GM03
  * Andrews (2010) Andrews, D. G. 2010, An Introduction to Atmospheric Physics, 2nd edn. (Cambridge University Press)
  * Asphaug & Reufer (2014) Asphaug, E., & Reufer, A. 2014, Nature Geoscience, 7, 564, doi: 10.1038/ngeo2189
  * Atreya et al. (1986) Atreya, S. K., Sandel, B. R., Romani, P. N., et al. 1986, Bull. Am. Astron. Soc., 18, 758
  * Bale et al. (2016) Bale, C. W., Bélisle, E., Chartrand, P., et al. 2016, Calphad: Computer Coupling of Phase Diagrams and Thermochemistry, 55, 1, doi: 10.1016/j.calphad.2016.07.004
  * Benz et al. (2008) Benz, W., Anic, A., Horner, J., & Whitby, J. A. 2008, in Mercury (Springer), 7–20
  * Benz et al. (1988) Benz, W., Slattery, W. L., & Cameron, A. 1988, Icarus, 74, 516
  * Bhatia & Sahijpal (2017) Bhatia, G. K., & Sahijpal, S. 2017, Meteorit. Planet. Sci., 52, 295, doi: 10.1111/maps.12789
  * Bird (1994) Bird, G. A. 1994, Molecular Gas Dynamics And The Direct Simulation Of Gas Flows (Springer)
  * Bouvier et al. (2018) Bouvier, L. C., Costa, M. M., Connelly, J. N., et al. 2018, Nature, 558, 586, doi: 10.1038/s41586-018-0222-z
  * Bower et al. (2021) Bower, D. J., Hakim, K., Sossi, P. A., & Sanan, P. 2021, Planet. Sci. J.
  * Bower et al. (2019) Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, Astron. Astrophys., 631, A103, doi: 10.1051/0004-6361/201935710
  * Bower et al. (2019) Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, Astron. Astrophys., 631, A103, doi: 10.1051/0004-6361/201935710
  * Bower et al. (2018) Bower, D. J., Sanan, P., & Wolf, A. S. 2018, Phys. Earth Planet. Inter., 274, 49, doi: 10.1016/j.pepi.2017.11.004
  * Brown & Elkins-Tanton (2009) Brown, S. M., & Elkins-Tanton, L. T. 2009, Earth Planet. Sci. Lett., 286, 446, doi: 10.1016/j.epsl.2009.07.010
  * Cartier & Wood (2019) Cartier, C., & Wood, B. J. 2019, Elements, 15, 39, doi: 10.2138/gselements.15.1.39
  * Cassidy et al. (2016) Cassidy, T. A., McClintock, W. E., Killen, R. M., et al. 2016, Geophys. Res. Lett., 43, 11,121, doi: 10.1002/2016GL071071
  * Chao et al. (2021) Chao, K.-H., deGraffenried, R., Lach, M., et al. 2021, Geochemistry, 125735
  * Chapman & Cowling (1970) Chapman, S., & Cowling, T. G. 1970, Press, Cambridge, England
  * Charlier et al. (2013) Charlier, B., Grove, T. L., & Zuber, M. T. 2013, Earth Planet. Sci. Lett., 363, 50, doi: 10.1016/j.epsl.2012.12.021
  * Chau et al. (2018) Chau, A., Reinhardt, C., Helled, R., & Stadel, J. 2018, Atrophys. J., 865, 35
  * Clementi et al. (1967) Clementi, E., Raimondi, D. L., & Reinhardt, W. P. 1967, The Journal of Chemical Physics, 47, 1300, doi: 10.1063/1.1712084
  * Duncan et al. (2017) Duncan, M. S., Dasgupta, R., & Tsuno, K. 2017, Earth and Planetary Science Letters, 466, 115, doi: 10.1016/j.epsl.2017.03.008
  * Egan et al. (2019) Egan, H., Jarvinen, R., Ma, Y., & Brain, D. 2019, MNRAS, 488, 2108, doi: 10.1093/mnras/stz1819
  * Elkins-Tanton (2008) Elkins-Tanton, L. 2008, Earth Planet. Sci. Lett., 271, 181, doi: 10.1016/j.epsl.2008.03.062
  * Elkins-Tanton (2012) Elkins-Tanton, L. T. 2012, Ann. Rev. Earth Planet. Sci., 40, 113, doi: 10.1146/annurev-earth-042711-105503
  * Evans et al. (2015) Evans, L. G., Peplowski, P. N., McCubbin, F. M., et al. 2015, Icarus, 257, 417, doi: 10.1016/J.ICARUS.2015.04.039
  * Fegley & Cameron (1987) Fegley, B., & Cameron, A. G. W. 1987, Earth Planet. Sci. Lett., 82, 207, doi: 10.1016/0012-821X(87)90196-8
  * Gamborino et al. (2019) Gamborino, D., Vorburger, A., & Wurz, P. 2019, Ann. Geophys., 37, 455, doi: 10.5194/angeo-37-455-2019
  * Gebek & Oza (2020) Gebek, A., & Oza, A. V. 2020, MNRAS, 497, 5271, doi: 10.1093/mnras/staa2193
  * Genda & Abe (2003) Genda, H., & Abe, Y. 2003, Icarus, 164, 149, doi: 10.1016/S0019-1035(03)00101-5
  * Ghiorso & Sack (1995) Ghiorso, M. S., & Sack, R. O. 1995, Contributions to Mineralogy and Petrology, 119, 197, doi: 10.1007/BF00307281
  * Gronoff et al. (2020) Gronoff, G., Arras, P., Baraka, S., et al. 2020, Journal of Geophysical Research: Space Physics, 125, doi: 10.1029/2019JA027639
  * Grott et al. (2011) Grott, M., Breuer, D., & Laneuville, M. 2011, Earth Planet. Sci. Lett., 307, 135, doi: 10.1016/j.epsl.2011.04.040
  * Hin et al. (2017) Hin, R. C., Coath, C. D., Carter, P. J., et al. 2017, Nature, 549, 511, doi: 10.1038/nature23899
  * Humayun & Cassen (2000) Humayun, M., & Cassen, P. 2000, Origin of the Earth and Moon, 3
  * Hyodo et al. (2021) Hyodo, R., Genda, H., & Brasser, R. 2021, Icarus, 354, 114064, doi: 10.1016/j.icarus.2020.114064
  * Ito & Ikoma (2021) Ito, Y., & Ikoma, M. 2021, Monthly Notices of the Royal Astronomical Society, 502, 750, doi: 10.1093/mnras/staa3962
  * Ito et al. (2015) Ito, Y., Ikoma, M., Kawahara, H., et al. 2015, The Astrophysical Journal, 801, 144, doi: 10.1088/0004-637X/801/2/144
  * Johnson (1990) Johnson, R. E. 1990, Energetic Charged-Particle Interactions with Atmospheres and Surfaces (Springer)
  * Johnson (2004) —. 2004, ApJ, 609, L99, doi: 10.1086/422912
  * Johnson et al. (2015) Johnson, R. E., Oza, A., Young, L. A., Volkov, A. N., & Schmidt, C. 2015, ApJ, 809, 43, doi: 10.1088/0004-637X/809/1/43
  * Johnstone et al. (2015) Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015, Astronomy and Astrophysics, 577, A28, doi: 10.1051/0004-6361/201425301
  * Keppler & Golabek (2019) Keppler, H., & Golabek, G. 2019, Geochemical Perspectives Letters, 12, doi: 10.7185/geochemlet.1918
  * Kim & Desclaux (2002) Kim, Y. K., & Desclaux, J. P. 2002, Physical Review A - Atomic, Molecular, and Optical Physics, 66, 127081, doi: 10.1103/PhysRevA.66.012708
  * Korth et al. (2012) Korth, H., Anderson, B. J., Johnson, C. L., et al. 2012, Journal of Geophysical Research: Space Physics, 117, n/a, doi: 10.1029/2012JA018052
  * Krenn et al. (2021) Krenn, A. F., Fossati, L., Kubyshkina, D., & Lammer, H. 2021, Astronomy & Astrophysics, 650, A94, doi: 10.1051/0004-6361/202140437
  * Kruijer et al. (2014) Kruijer, T. S., Touboul, M., Fischer-Gödde, M., et al. 2014, Science, 344, 1150, doi: 10.1126/science.1251766
  * Lamoreaux & Hildenbrand (1984) Lamoreaux, R. H., & Hildenbrand, D. L. 1984, Journal of Physical and Chemical Reference Data, 13, 151, doi: 10.1063/1.555706
  * Lamoreaux et al. (1987) Lamoreaux, R. H., Hildenbrand, D. L., & Brewer, L. 1987, Journal of Physical and Chemical Reference Data, 16, 419, doi: 10.1063/1.555799
  * Lauretta et al. (2007) Lauretta, D. S., Goreva, J. S., Hill, D. H., & Killgore, M. 2007, in Lunar and Planetary Science Conference, Vol. 38, 2236
  * Lellouch et al. (1992) Lellouch, E., Belton, M., De Pater, I., et al. 1992, Icarus, 98, 271, doi: 10.1016/0019-1035(92)90095-O
  * Lellouch et al. (2003) Lellouch, E., Paubert, G., Moses, J. I., Schneider, N. M., & Strobel, D. F. 2003, Nature, 421, 45, doi: 10.1038/nature01292
  * Lewis (1972) Lewis, J. S. 1972, Earth and Planetary Science Letters, 15, 286, doi: 10.1016/0012-821X(72)90174-4
  * Lichtenberg et al. (2021) Lichtenberg, T., Bower, D. J., Hammond, M., et al. 2021, J. Geophys. R. Planets, e2020JE006711, doi: https://doi.org/10.1029/2020JE006711
  * McCoy et al. (2018) McCoy, T. J., Peplowski, P. N., McCubbin, F. M., & Weider, S. Z. 2018, in Mercury: The View after MESSENGER, ed. S. C. Solomon, L. R. Nittler, & B. J. Anderson (Cambridge University Press), 176–190, doi: 10.1017/9781316650684
  * McCubbin et al. (2012) McCubbin, F. M., Riner, M. A., Vander Kaaden, K. E., & Burkemper, L. K. 2012, Geophysical Research Letters, 39
  * Mordasini (2020) Mordasini, C. 2020, Astronomy and Astrophysics, 638, A52, doi: 10.1051/0004-6361/201935541
  * Moullet et al. (2010) Moullet, A., Gurwell, M. A., Lellouch, E., & Moreno, R. 2010, Icarus, 208, 353, doi: 10.1016/j.icarus.2010.02.009
  * Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, Astrophysical Journal, 693, 23, doi: 10.1088/0004-637X/693/1/23
  * Namur et al. (2016) Namur, O., Collinet, M., Charlier, B., et al. 2016, Earth and Planetary Science Letters, 439, 117, doi: 10.1016/j.epsl.2016.01.030
  * Nittler & Weider (2019) Nittler, L. R., & Weider, S. Z. 2019, Elements, 15, 33, doi: 10.2138/gselements.15.1.33
  * Nittler et al. (2011) Nittler, L. R., Starr, R. D., Weider, S. Z., et al. 2011, Science, 333, 1847, doi: 10.1126/science.1211567
  * O’Neill & Eggins (2002) O’Neill, H. S. C., & Eggins, S. M. 2002, Chemical Geology, 186, 151, doi: 10.1016/S0009-2541(01)00414-4
  * O’Neill & Palme (2008) O’Neill, H. S. C., & Palme, H. 2008, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 366, 4205
  * Oza et al. (2019) Oza, A. V., Johnson, R. E., Lellouch, E., et al. 2019, Astrophys. J., 885, 168, doi: 10.3847/1538-4357/ab40cc
  * Peplowski et al. (2014) Peplowski, P. N., Evans, L. G., Stockstill-Cahill, K. R., et al. 2014, Icarus, 228, 86, doi: 10.1016/J.ICARUS.2013.09.007
  * Perera et al. (2018) Perera, V., Jackson, A. P., Elkins-Tanton, L. T., & Asphaug, E. 2018, J. Geophys. Res.-Planet., 123, 1168, doi: https://doi.org/10.1029/2017JE005512
  * Pyykkö & Atsumi (2009a) Pyykkö, P., & Atsumi, M. 2009a, Chemistry - A European Journal, 15, 12770, doi: 10.1002/chem.200901472
  * Pyykkö & Atsumi (2009b) —. 2009b, Chemistry - A European Journal, 15, 186, doi: 10.1002/chem.200800987
  * Pyykkö et al. (2005) Pyykkö, P., Riedel, S., & Patzschke, M. 2005, Chemistry - A European Journal, 11, 3511, doi: 10.1002/chem.200401299
  * Renggli et al. (2017) Renggli, C. J., King, P. L., Henley, R. W., & Norman, M. D. 2017, Geochimica et Cosmochimica Acta, 206, 296, doi: 10.1016/j.gca.2017.03.012
  * Renggli & Klemme (2020) Renggli, C. J., & Klemme, S. 2020, Journal of Volcanology and Geothermal Research, 400, 106929, doi: 10.1016/j.jvolgeores.2020.106929
  * Ribas et al. (2014) Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54, doi: 10.1051/0004-6361/201322597
  * Ruedas (2017) Ruedas, T. 2017, Geochem. Geophy. Geosys., 18, 3530, doi: 10.1002/2017GC006997
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, Astronomy & Astrophysics, 410, 611, doi: 10.1051/0004-6361:20031279
  * Siegfried & Solomon (1974) Siegfried, R. W., & Solomon, S. C. 1974, Icarus, 23, 192 , doi: https://doi.org/10.1016/0019-1035(74)90005-0
  * Sossi & Fegley (2018) Sossi, P. A., & Fegley, B. 2018, Reviews in Mineralogy and Geochemistry, 84, 393, doi: 10.2138/rmg.2018.84.11
  * Sossi et al. (2019) Sossi, P. A., Klemme, S., O’Neill, H. S., Berndt, J., & Moynier, F. 2019, Geochimica et Cosmochimica Acta, 260, 204, doi: 10.1016/J.GCA.2019.06.021
  * Spohn (1991) Spohn, T. 1991, Icarus, 90, 222 , doi: 10.1016/0019-1035(91)90103-Z
  * Stevenson (1990) Stevenson, D. J. 1990, in Origin of the Earth, ed. H. E. Newsom & J. H. Jones (Oxford University Press, New York), 231–249
  * Stevenson et al. (1983) Stevenson, D. J., Spohn, T., & Schubert, G. 1983, Icarus, 54, 466, doi: 10.1016/0019-1035(83)90241-5
  * Thomas et al. (2004) Thomas, N., Bagenal, F., Hill, T. W., & Wilson, J. K. 2004, in Jupiter. the planet, satellites and magnetosphere, ed. F. Bagenal, T. E. Dowling, & W. B. McKinnon, Vol. 1 (New York: Cambridge University Press), 561–591.  https://ui.adsabs.harvard.edu/abs/2004jpsm.book..561T
  * Tosi et al. (2013) Tosi, N., Grott, M., Plesa, A.-C., & Breuer, D. 2013, J. Geophys. Res.-Planet., 118, 2474, doi: 10.1002/jgre.20168
  * Tsai et al. (2017) Tsai, S.-M., Lyons, J. R., Grosheintz, L., et al. 2017, Astrophys. J. Sup., 228, 20
  * Tsai et al. (2021) Tsai, S.-M., Malik, M., Kitzmann, D., et al. 2021.  https://arxiv.org/abs/2108.01790
  * Tu et al. (2015) Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, Astron. Astrophys., 577, L3, doi: 10.1051/0004-6361/201526146
  * Tucker et al. (2021) Tucker, O. J., Killen, R. M., Johnson, R. E., & Saxena, P. 2021, Icarus, 359, 114304, doi: 10.1016/j.icarus.2021.114304
  * Vander Kaaden & McCubbin (2016) Vander Kaaden, K. E., & McCubbin, F. M. 2016, Geochim. Cosmochim. Ac., 173, 246, doi: 10.1016/J.GCA.2015.10.016
  * Visscher & Fegley (2013) Visscher, C., & Fegley, B. 2013, The Astrophysical Journal Letters, 767, 12, doi: 10.1088/2041-8205/767/1/L12
  * Vlasov & Kelley (2015) Vlasov, M. N., & Kelley, M. C. 2015, Ann. Geophys., 33, 857, doi: 10.5194/angeo-33-857-2015
  * Volkov & Johnson (2013) Volkov, A. N., & Johnson, R. E. 2013, Astrophysical Journal, 765, 90, doi: 10.1088/0004-637X/765/2/90
  * Vorburger et al. (2015) Vorburger, A., Wurz, P., Lammer, H., Barabash, S., & Mousis, O. 2015, Icarus, 262, 14, doi: 10.1016/j.icarus.2015.07.035
  * Watson et al. (1981) Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150, doi: 10.1016/0019-1035(81)90101-9
  * Weidenschilling (1978) Weidenschilling, S. J. 1978, Icarus, 35, 99, doi: 10.1016/0019-1035(78)90064-7
  * Weisberg et al. (2000) Weisberg, M. K., Prinz, M., Clayton, R. N., et al. 2000, Meteoritics and Planetary Science, 36, 401, doi: 10.1111/j.1945-5100.2001.tb01882.x
  * Weisberg et al. (1990) Weisberg, M. K., Prinz, M., & Nehru, C. E. 1990, Meteoritics, 25, 269, doi: 10.1111/j.1945-5100.1990.tb00711.x
  * Wiik (1956) Wiik, H. B. 1956, Geochim. Cosmochim. Ac., 9, 279, doi: 10.1016/0016-7037(56)90028-X
  * Wolf et al. (2021) Wolf, A. S., Jäggi, N., Sossi, P. A., Bower, D. J., & Ghiorso, M. S. 2021, VapoRock: Thermodynamics of vaporized silicate rocks & melts for modeling magma ocean atmospheres and stellar nebula, v0.1, Zenodo, doi: 10.5281/zenodo.4594226
  * Wurz & Lammer (2003) Wurz, P., & Lammer, H. 2003, Icarus, 164, 1, doi: 10.1016/S0019-1035(03)00123-4