# The Extent, Nature, and Origin of K and Rb Depletions and Isotopic Fractionations in Earth, the Moon, and Other Planetary Bodies

Nicolas Dauphas, Nicole X. Nie, Marc Blanchard, Zhe J. Zhang, Hao Zeng, Justin Y. Hu, Merlin Meheut, Channon Visscher, Robin Canup, Timo Hopp

► To cite this version:

Nicolas Dauphas, Nicole X. Nie, Marc Blanchard, Zhe J. Zhang, Hao Zeng, et al.. The Extent, Nature, and Origin of K and Rb Depletions and Isotopic Fractionations in Earth, the Moon, and Other Planetary Bodies. The Planetary Science Journal, 2022, 3, 10.3847/PSJ/ac2e09. insu-03661489

HAL Id: insu-03661489

https://insu.hal.science/insu-03661489

Submitted on 7 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Extent, Nature, and Origin of K and Rb Depletions and Isotopic Fractionations in Earth, the Moon, and Other Planetary Bodies

Nicolas Dauphas<fnref n="1" />,<fnref n="2" /> ID, Nicole X. Nie<fnref n="1" />,<fnref n="2" /> ID, Marc Blanchard<fnref n="3" /> ID, Zhe J. Zhang<fnref n="1" /> ID, Hao Zeng<fnref n="1" />,<fnref n="4" />, Justin Y. Hu<fnref n="1" /> ID, Merlin Meheut<fnref n="3" /> ID, Channon Visscher<fnref n="5" /> ID, Robin Canup<fnref n="6" /> ID, and Timo Hopp<fnref n="1" /> ID

Received 2021 April 6; revised 2021 October 2; accepted 2021 October 5; published 2022 February 3

## Abstract

Moderately volatile elements (MVEs) are depleted and isotopically fractionated in the Moon relative to Earth. To understand how the composition of the Moon was established, we calculate the equilibrium and kinetic isotopic fractionation factors associated with evaporation and condensation processes. We also reassess the levels of depletions of K and Rb in planetary bodies. Highly incompatible element ratios are often assumed to be minimally affected by magmatic processes, but we show that this view is not fully warranted, and we develop approaches to mitigate this issue. The K/U weight ratios of Earth and the Moon are estimated to be 9704 and 2448, respectively. The ⁸⁷Rb/⁸⁶Sr atomic ratios of Earth and the Moon are estimated to be 0.072 5 and 0.015 4, respectively. We show that the depletions and heavy isotopic compositions of most MVEs in the Moon are best explained by evaporation in 99%-saturated vapor. At 99% saturation in the protolunar disk, Na and K would have been depleted to levels like those encountered in the Moon on timescales of ~40–400 days at 3500–4500 K, which agrees with model expectations. In contrast, at the same saturation but a temperature of 1600–1800 K relevant to hydrodynamic escape from the lunar magma ocean, Na and K depletions would have taken 0.1–103 Myr, which far exceeds the 1000 yr time span until plagioclase flotation hinders evaporation from the magma ocean. We conclude that the protolunar disk is a much more likely setting for the depletion of MVEs than the lunar magma ocean.

*Unified Astronomy Thesaurus concepts: Earth-moon system (436); Isotopic abundances (867); Lunar composition (948); Chemical abundances (224)*

*Supporting material: tar.gz file*

## 1. Introduction

Models of the formation of the Moon are constrained by the dynamical, chemical, and isotopic characteristics of the Moon (Canup & Asphaug 2001; Canup 2004a, 2012; Ćuk & Stewart 2012; Reufer et al. 2012; Rufu et al. 2017; Citron et al. 2018; Lock et al. 2018; Hosono et al. 2019; Wissing & Hobbs 2020; Asphaug et al. 2021). Tremendous progress in establishing those characteristics has been achieved by applying new methodologies and instrumentation to study lunar meteorites and lunar samples recovered by the Apollo, Luna, and Chang'E-5 missions (Saal et al. 2008; Qin et al. 2010; Sharp et al. 2010; Armytage et al. 2012; Fitoussi & Bourdon 2012; Paniello et al. 2012; Zhang et al. 2012; Hui et al. 2013; Herwartz et al. 2014; Kato et al. 2015; Kruijer et al. 2015; Touboul et al. 2015; Wang & Jacobsen 2016; Young et al. 2016; Pringle & Moynier 2017; Greenwood et al. 2018; Mougel et al. 2018; Nie & Dauphas 2019; Wang et al. 2019; Cano et al. 2020; Tian et al. 2020; Che et al. 2021; Li et al. 2021). Missions such as Clementine (Giguere et al. 2000), Lunar Prospector (Feldman et al. 1998; Lawrence et al. 1998), Kaguya (Kobayashi et al. 2010; Yamashita et al. 2010;

Kobayashi et al. 2012), GRAIL (Wieczorek et al. 2013), Chang'E-1 (Zhu et al. 2010), and Chang'E-4 (Li et al. 2019) have also provided new insights into the composition of the lunar surface and its internal structure. This wealth of data has shaken the paradigm of formation of the Moon by a giant impact involving the collision of the proto-Earth with a Mars-size impactor at a velocity just above the mutual escape velocity and an impact angle of 45° (often referred to as the “canonical model” of lunar formation). While this model can readily explain the angular momentum of the Earth-Moon system and the small size of the lunar core (Canup & Asphaug 2001; Canup 2004a), the isotopic similarities in O (Wiechert et al. 2001; Spicuzza et al. 2007; Herwartz et al. 2014; Young et al. 2016; Greenwood et al. 2018; Cano et al. 2020), Ti (Zhang et al. 2012), Cr (Qin et al. 2010; Mougel et al. 2018), and W (Kruijer et al. 2015; Touboul et al. 2015) between lunar and terrestrial rocks and the large depletions (LSPET 1969; Tera et al. 1970; Larimer 1971; Ganapathy & Anders 1974; Morgan et al. 1978; Wänke 1981; Taylor 1982; O'Neill 1991; McDonough et al. 1992; Warren 2005; Taylor & Wieczorek 2014; Albarède et al. 2015) and significant isotopic fractionations (Paniello et al. 2012; Kato et al. 2015; Wang & Jacobsen 2016; Kato & Moynier 2017; Pringle & Moynier 2017; Nie & Dauphas 2019; Wang et al. 2019; Tian et al. 2020) in moderately volatile elements (MVEs) of lunar rocks are more difficult to explain. For those reasons, more

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Origins Laboratory, Department of the Geophysical Sciences and Enrico Fermi Institute, The University of Chicago, 5734 South Ellis Avenue, Chicago, IL 60637, USA; dauphas@uchicago.edu</footnote>

<footnote n="2">² Earth and Planets Laboratory, Carnegie Institution for Science, Washington, DC 20015, USA</footnote>

<footnote n="3">³ Géosciences Environnement Toulouse (GET), Université de Toulouse, CNRS, IRD, CNES, UPS, Toulouse, France</footnote>

<footnote n="4">⁴ Institute of Molecular Engineering and Department of Chemistry, The University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="5">⁵ Department of Chemistry, Dordt College, Sioux Center, IA 51250, USA</footnote>

<footnote n="6">⁶ Planetary Science Directorate, Southwest Research Institute, Boulder, CO 80302, USA</footnote>
complex models of the Moon-forming giant impact or its aftermath have been considered (Pahlevan & Stevenson 2007; Canup 2012; Ćuk & Stewart 2012; Reufer et al. 2012; Rufu et al. 2017; Citron et al. 2018; Lock et al. 2018; Asphaug et al. 2021).

An important constraint on lunar formation scenarios is the depletion in MVEs of the Moon relative to Earth (Table 1, Figure 1; see Section 1 for details). In the present study, we examine the origin of the depletions of the Moon in K, Cu, Rb, Zn, Ga, and Sn, with a focus on K and Rb, as these have been extensively studied in lunar rocks and their strongly lithophile and incompatible nature makes it easier to assess the bulk lunar composition. Earth is depleted in K and Rb relative to CI (Ivuna-type) carbonaceous chondrites by factors of ~7 and ~12, respectively. Lunar rocks show further depletions in K and Rb relative to CI of ~29 and ~57, respectively. Those chemical depletions are accompanied by isotopic fractionation (Table 2; Wang & Jacobsen 2016; Pringle & Moynier 2017; Nie & Dauphas 2019; Tian et al. 2020).

Early studies investigating the isotopic composition of K (expressed here as $\delta^{41}K(\%e) = [(^{41}K/^{39}K)_{SMP}/(^{41}K/^{39}K)_{STD} - 1] \times 1000$, where SMP is sample and STD is the SRM3141a reference material) focused on lunar soils and breccias, where they found large isotopic variations (Garner et al. 1975; Church et al. 1976; Humayun & Clayton 1995a). Sodium and potassium are important constituents of the lunar atmosphere, which is constantly replenished by vaporization from exposed rocks (Potter & Morgan 1988; Tyler et al. 1988). Vaporization, transport at the lunar surface, recondensation, and loss to space are responsible for the K isotopic variations measured in lunar soils and breccias. Humayun and Clayton (1995b) measured lunar samples that were not affected by these processes and found terrestrial-like isotopic compositions at a precision of ~±0.5‰. Wang & Jacobsen (2016) and Tian et al. (2020) improved the precision of K isotopic analyses and found that mare basalts, taken as proxies for the bulk Moon, had heavy K isotopic composition relative to terrestrial basalts (+0.41 ± 0.09‰). Tian et al. (2020) found heavy K isotopic composition in KREEP-bearing lunar samples relative to normal mare basalts, which they interpreted to reflect volatilization of K from the magma ocean late in the differentiation history of the Moon, as has been suggested for other elements (Dhaliwal et al. 2018; Day et al. 2020). These KREEP-bearing samples are, however, not representative of the silicate Moon composition.

Early studies investigating the isotopic composition of Rb (expressed as $\delta^{87}Rb (\%e) = [(^{87}Rb/^{85}Rb)_{SMP}/(^{87}Rb/^{85}Rb)_{STD} - 1] \times 1000$, where SMP is sample and STD is the SRM984 reference material) focused on lunar soils and breccias, where no fractionation was found but the error bars were ±1.5‰ (Garner et al. 1975). The focus has only recently shifted to establishing the Rb isotopic composition of the Moon. Nebel et al. (2011) measured some lunar samples, but their precision was still limited, and they failed to detect any difference with terrestrial rocks. Pringle & Moynier (2017) and Nie & Dauphas (2019) measured lunar samples with improved precision and found that lunar samples were slightly enriched in the heavy isotope of Rb by +0.16 ± 0.04‰.

Several other MVEs are isotopically fractionated in lunar rocks relative to terrestrial rocks (Table 2), notably Zn (+1.1 ± 0.5 ‰ for the $^{66}\text{Zn}/^{64}\text{Zn}$ ratio; Paniello et al. 2012; Kato et al. 2015), Cu (+0.5 ± 0.1 ‰ for the $^{65}\text{Cu}/^{63}\text{Cu}$ ratio; Herzog et al. 2009), Ga (+0.3 ± 0.2 ‰ for the $^{71}\text{Ga}/^{69}\text{Ga}$

ratio; Kato & Moynier 2017; Wimpenny et al. 2022), and Sn (−0.48 ± 0.15‰ for the $^{124}\text{Sn}/^{16}\text{Sn}$ ratio; Wang et al. 2019). For all these elements except Sn (Wang et al. 2019), lunar rocks were found to have heavy isotopic compositions. Some of these elements show significant isotopic variations among different lunar lithologies, which makes it difficult to estimate the isotopic composition of the bulk Moon. Other elements such as Cr have complex behaviors during core formation, mantle melting, and magmatic differentiation, which can be associated with significant isotopic fractionation, further complicating interpretations. Sossi et al. (2018) observed that many lunar rocks had light Cr isotopic compositions relative to terrestrial igneous rocks. They explained this difference by Cr volatilization from the proto-Moon under oxidizing equilibrium conditions (near the FMQ buffer) at low temperature (1700 K). This interpretation was questioned, however, by Shen et al. (2020), who argued that after correction for Cr isotopic fractionation induced by magmatic processes (Bonnand et al. 2016; Xia et al. 2017; Shen et al. 2018), Earth and the Moon have similar Cr isotopic compositions. Because of these uncertainties and the fact that Cr is not obviously depleted in the Moon relative to Earth (Cr ~ 2190–2463 ppm in the lunar mantle; Seifert & Ringwood 1988; compared to ~2625 ppm in the terrestrial mantle; McDonough & Sun 1995), we do not use this element to constrain the conditions of MVE depletion in the Moon.

Developing a quantitative understanding of the origin of the depletion in MVEs of the Moon has proven to be challenging (Canup et al. 2015; Charnoz & Michaut 2015; Lock et al. 2018), but the measured isotopic fractionations for MVEs in lunar rocks may provide new insights into the conditions that prevailed when these elements were depleted. Interpretation of available isotopic data is, however, not straightforward. For example, Wang et al. (2019) and Sossi et al. (2020) argued that the isotopic compositions of K, Cu, Zn, Rb, and Sn could be explained by equilibrium isotopic fractionation between vapor and liquid, while Nie & Dauphas (2019) made the case that the same data (except Sn) could be explained by kinetic isotopic fractionation induced by evaporation in near-saturated conditions. The critical questions that remain to be answered about the depletion and isotopic fractionation of MVEs in the Moon are the following: (i) Do the measured isotopic fractionations in lunar rocks reflect equilibrium or kinetic isotopic fractionation? (ii) If the isotopic fractionations are equilibrium in nature, what do they tell us about the conditions of lunar formation, notably temperature? (iii) If the isotopic fractionations are kinetic in nature, what was the undersaturation level during MVE loss, and can it help constrain lunar formation models?

We revisit the question of the origin of the depletion in MVEs of the Moon by carefully examining these three questions. While working on this manuscript, we realized that improvements could be brought to estimates of the depletions of K and Rb in the Moon and other planetary bodies. In Section 2, we therefore start by reassessing the depletions in K and Rb in planetary bodies by using literature data, focusing on the effect of magmatic differentiation. In Section 3, we present ab initio calculations to estimate equilibrium isotopic fractionation between vapor and liquid to assess whether this could be an explanation for the isotopic compositions of MVEs in lunar rocks. In Section 4, we evaluate whether the measured isotopic fractionations in lunar rocks can be explained by kinetic effects during vaporization in an undersaturated medium, and we
<table><caption>**Table 1**
K and Rb Depletions and Enrichments in Selected Planetary Materials (see Text for Details)</caption><thead><tr><th></th><th>K<br>g g<sup>-1</sup></th><th>Ba<br>g g<sup>-1</sup></th><th>Th<br>g g<sup>-1</sup></th><th>U<br>g g<sup>-1</sup></th><th>La<br>g g<sup>-1</sup></th><th>Rb<br>g g<sup>-1</sup></th><th>Sr<br>g g<sup>-1</sup></th><th>K/U<br>g g<sup>-1</sup></th><th>Rb/Ba<br>g g<sup>-1</sup></th><th>±</th><th>K/Rb<br>g g<sup>-1</sup></th><th>Rb/Sr<br>g g<sup>-1</sup></th><th>±</th><th><sup>87</sup>Rb/<sup>86</sup>Sr<br>at.</th><th>±</th><th>K/U × CI</th><th>±</th><th>Rb/Sr × CI</th><th>±</th><th>K/Rb × CI</th><th>±</th></tr></thead><tbody><tr><td>CI</td><td>5.50E-04</td><td>2.46E-06</td><td>2.83E-09</td><td>7.70E-08</td><td>2.35E-07</td><td>2.33E-06</td><td>7.73E-06</td><td>71429</td><td>0.9472</td><td></td><td>236</td><td>0.3014</td><td></td><td>0.8722</td><td></td><td>1.00</td><td></td><td>1.00000</td><td></td><td>1.00</td><td></td></tr><tr><td>CM</td><td>4.00E-04</td><td>3.30E-06</td><td>4.00E-08</td><td>1.10E-08</td><td>3.17E-07</td><td>1.70E-06</td><td>1.01E-05</td><td>36364</td><td>0.5152</td><td></td><td>235</td><td>0.1683</td><td></td><td>0.4870</td><td></td><td>0.51</td><td></td><td>0.55841</td><td></td><td>1.00</td><td></td></tr><tr><td>CO</td><td>3.45E-04</td><td>4.29E-06</td><td>4.50E-08</td><td>1.30E-08</td><td>3.87E-07</td><td>1.45E-06</td><td>1.27E-05</td><td>26538</td><td>0.3380</td><td></td><td>238</td><td>0.1142</td><td></td><td>0.3304</td><td></td><td>0.37</td><td></td><td>0.37878</td><td></td><td>1.01</td><td></td></tr><tr><td>CV</td><td>3.10E-04</td><td>4.90E-06</td><td>6.00E-08</td><td>1.70E-08</td><td>4.86E-07</td><td>1.25E-06</td><td>1.53E-05</td><td>18235</td><td>0.2551</td><td></td><td>248</td><td>0.0817</td><td></td><td>0.2364</td><td></td><td>0.26</td><td></td><td>0.27105</td><td></td><td>1.05</td><td></td></tr><tr><td>H</td><td>7.80E-04</td><td>4.20E-06</td><td>4.20E-08</td><td>1.20E-08</td><td>2.95E-07</td><td>2.90E-06</td><td>1.00E-05</td><td>65000</td><td>0.6905</td><td></td><td>269</td><td>0.2900</td><td></td><td>0.8391</td><td></td><td>0.91</td><td></td><td>0.96210</td><td></td><td>1.14</td><td></td></tr><tr><td>L</td><td>8.25E-04</td><td>3.70E-06</td><td>4.30E-08</td><td>1.30E-08</td><td>3.10E-07</td><td>3.10E-06</td><td>1.11E-05</td><td>63462</td><td>0.8378</td><td></td><td>266</td><td>0.2793</td><td></td><td>0.8081</td><td></td><td>0.89</td><td></td><td>0.92654</td><td></td><td>1.13</td><td></td></tr><tr><td>LL</td><td>7.90E-04</td><td>4.80E-06</td><td>4.30E-08</td><td>1.30E-08</td><td>3.15E-07</td><td>3.10E-06</td><td>1.11E-05</td><td>60769</td><td>0.6458</td><td></td><td>255</td><td>0.2793</td><td></td><td>0.8081</td><td></td><td>0.85</td><td></td><td>0.92654</td><td></td><td>1.08</td><td></td></tr><tr><td>EH</td><td>8.00E-04</td><td>2.60E-06</td><td>3.00E-08</td><td>9.00E-09</td><td>2.35E-07</td><td>2.60E-06</td><td>7.20E-06</td><td>88889</td><td>1.0000</td><td></td><td>308</td><td>0.3611</td><td></td><td>1.0449</td><td></td><td>1.24</td><td></td><td>1.19802</td><td></td><td>1.30</td><td></td></tr><tr><td>EL</td><td>7.35E-04</td><td>3.03E-08</td><td>3.50E-08</td><td>1.00E-08</td><td>1.90E-07</td><td>2.50E-06</td><td>8.20E-06</td><td>73500</td><td>0.8242</td><td></td><td>294</td><td>0.3049</td><td></td><td>0.8822</td><td></td><td>1.03</td><td></td><td>1.01146</td><td></td><td>1.25</td><td></td></tr><tr><td>Bulk silicate Earth</td><td>1.97E-04</td><td>6.49E-06</td><td>7.46E-08</td><td>2.03E-08</td><td>6.20E-07</td><td>5.10E-05</td><td>2.04E-05</td><td>9704</td><td>2054</td><td>0.0787</td><td>0.0157</td><td>386</td><td>112</td><td>0.0250</td><td>0.0050</td><td>0.0725</td><td>0.0145</td><td>0.13586</td><td>0.02876</td><td>0.08309</td><td>0.01658</td><td>1.64</td><td>0.48</td></tr><tr><td>Bulk silicate Moon</td><td>4.40E-05</td><td>5.74E-06</td><td>6.60E-08</td><td>1.80E-08</td><td>5.48E-07</td><td>9.58E-08</td><td>1.80E-05</td><td>2448</td><td>140</td><td>0.0167</td><td>0.0018</td><td>459</td><td>56</td><td>0.0053</td><td>0.0006</td><td>0.0154</td><td>0.0017</td><td>0.03427</td><td>0.00196</td><td>0.01763</td><td>0.00190</td><td>1.94</td><td>0.24</td></tr><tr><td>Bulk silicate Mars</td><td>2.77E-04</td><td>5.11E-06</td><td>5.88E-08</td><td>1.60E-08</td><td>4.88E-07</td><td>9.30E-07</td><td>1.61E-05</td><td>17331</td><td>1771</td><td>0.182</td><td>0.036</td><td>298</td><td>66</td><td>0.0579</td><td>0.0115</td><td>0.168</td><td>0.033</td><td>0.24263</td><td>0.02479</td><td>0.19215</td><td>0.03801</td><td>1.26</td><td>0.28</td></tr><tr><td>Bulk silicate Vesta</td><td>5.75E-05</td><td>4.44E-06</td><td>5.11E-08</td><td>1.39E-08</td><td>4.24E-07</td><td>3.97E-08</td><td>1.40E-05</td><td>4140</td><td>805</td><td>0.00895</td><td>0.00172</td><td>1448</td><td>396</td><td>0.00285</td><td>0.00055</td><td>0.0082</td><td>0.0016</td><td>0.05796</td><td>0.01127</td><td>0.00945</td><td>0.00182</td><td>6.13</td><td>1.68</td></tr><tr><td>Angrite (ADOR)</td><td>1.76E-05</td><td></td><td></td><td>1.36E-07</td><td></td><td>3.94E-08</td><td>1.37E-04</td><td>129</td><td>74</td><td></td><td></td><td>447</td><td>131</td><td>0.000287</td><td>0.000142</td><td>0.00083</td><td>0.00041</td><td>0.00181</td><td>0.00104</td><td>0.00095</td><td>0.00047</td><td>1.89</td><td>0.55</td></tr><tr><td>Mercury (surface)</td><td>1.29E-03</td><td></td><td></td><td>9.00E-08</td><td></td><td></td><td></td><td>14311</td><td>4108</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.20035</td><td>0.05751</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note. The chondrite values are from Wasson & Kallmeyn (1988), except for the CI values, which are from Barrat et al. (2012).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Planetary depletions in K and Rb in planetary materials (see text and Table 1 for details). Because Sr and U are two refractory lithophile elements, K/U and Rb/Sr ratios normalized to CI (Ivuna-type) chondrite are direct measures of the depletions in MVEs K and Rb (Halliday & Porcelli 2001; Davis 2006). (a) K/U vs. Rb/Sr in the angrite parent body (APB), the Moon, Vesta (HED meteorites), Earth, Mercury, Mars (SNC meteorites), and CV, CO, CM, CI, LL, L, H, EL, and EH chondrites (listed from the most volatile-element-depleted to the most volatile-element-enriched objects). The diagonal line corresponds to equal depletions in K and Rb. (b) Relative volatilities of K and Rb. While to first order the two elements show similar depletions (panel (a)), in detail Rb appears to be more volatile than K such that the most volatile-element-depleted objects display K/Rb ratios that are higher than CI.</figcaption></figure>

discuss possible settings where such fractionation could have been produced.

## 2. A Reassessment of K and Rb Depletions in Earth, the Moon, Mars, Mercury, Vesta, and the Angrite Parent Body

The depletions of planetary bodies in K and Rb are fundamental constraints for studies aimed at understanding the setting under which MVE depletion occurred. Elemental ratios involving lithophile (rock-loving) highly incompatible (they prefer silicate melts as opposed to minerals) elements (e.g., K/U, K/La, and Rb/Ba) are often used to estimate the extents of those depletions. However, magmatic processes can fractionate these ratios, a complication that is often not well accounted for. The focus of this section is therefore to unravel magmatic processes to better evaluate the depletions of Earth and the Moon in K and Rb. We also discuss the depletions of K and Rb in Vesta, Mars, and the angrite parent body, as these provide insights into the degree of MVE depletion that one might consider for the Moon-forming impactor Theia. The results from this section are summarized in Table 1.

Potassium and rubidium are two lithophile, incompatible, and moderately volatile elements. To quantify their depletions, they are typically normalized to refractory, lithophile, and incompatible elements (e.g., U and Sr; Halliday & Porcelli 2001; Davis 2006). Uranium is a suitable element for K-normalization because these two elements have similar behaviors during magmatic differentiation and K, U, and Th can be measured remotely by spacecraft using gamma-ray spectrometers, allowing one to (i) assess the compositions of planetary bodies that have not yet been sampled and (ii) map the distribution of these elements on planetary surfaces (the Moon, Feldman et al. 1998; Lawrence et al. 1998; Prettyman et al. 2006; Kobayashi et al. 2010; Yamashita et al. 2010; Zhu et al. 2010; Kobayashi et al. 2012; Mars, Taylor et al. 2006; Boynton et al. 2007; Taylor et al. 2010; Mercury, Peplowski et al. 2011; Vesta, Prettyman et al. 2015). Strontium is used to normalize Rb (Davis 2006) because in some circumstances the Rb/Sr ratio can be constrained precisely from ⁸⁷Sr isotopic analyses (⁸⁷Rb decays into ⁸⁵Sr with

a half-life of 49.2 Gyr; DePaolo & Wasserburg 1976; O’Nions et al. 1977; Allègre et al. 1979). Unless otherwise stated, all elemental ratios are given as weight ratios (ppm.ppm⁻¹ or g.g⁻¹). We use the CI chondrite concentrations (taken as a proxy for solar composition) of selected elements for normalization purposes, and these are K = 550 ppm, Rb = 2.33 ppm, Sr = 7.73 ppm, Ba = 2.46 ppm, U = 0.0077 ppm, Th = 0.0283 ppm, and La = 0.235 ppm (Barrat et al. 2012). The reason for adopting the values from Barrat et al. (2012) is that we are primarily interested in elemental ratios and adopting concentrations measured from the same digestion aliquots minimizes potential bias. For comparison, Palme et al. (2014) reported CI chondrite concentration values of K = 546 ppm, Rb = 2.32 ppm, Sr = 7.79 ppm, Ba = 2.42 ppm, U = 0.0081 ppm, Th = 0.0300 ppm, and La = 0.2414 ppm. Lodders (2021) reported CI chondrite concentration values of K = 539 ppm, Rb = 2.22 ppm, Sr = 7.79 ppm, Ba = 2.39 ppm, U = 0.00816 ppm, Th = 0.0298 ppm, and La = 0.244 ppm. The CI K/Rb, K/U, and Rb/Sr weight ratios (ppm/ppm) are 236, 71,429, and 0.3014, respectively. The present best estimates of the solar photosphere K/Rb and Rb/Sr weight ratios are 170 ± 34 and 0.43 ± 0.09 g.g⁻¹, respectively (Scott et al. 2015; Grevesse et al. 2015).

The approach that we use to estimate the depletions of K and Rb is to tie their behaviors with those of elements that share the same incompatible and lithophile character but are highly refractory. The rationale for doing so is that refractory lithophile elements can reasonably be assumed to be in chondritic proportions in planetary mantles, so they provide a means of assessing how igneous processes modified elemental ratios involving MVEs such as K/U or Rb/Ba ratios. The elements La, K, and U have similar behaviors during melting and planetary differentiation (Heier & Rogers 1963; Wasserburg et al. 1964; Jochum et al. 1983; Wänke & Dreibus 1988; McDonough et al. 1992; Arevalo et al. 2009; Palme & O'Neill 2014). The elements Ba and Rb also behave very similarly during magmatic processes (Hofmann & White 1983).
<table><caption>**Table 2**
Lunar Depletion, Isotopic Fraction, and Kinetic and Equilibrium Fractionation Factors of MVEs</caption><thead><tr><th rowspan="2">Element</th><th colspan="2">50% Te<sup>a</sup></th><th rowspan="2">Depletion</th><th rowspan="2">References</th><th rowspan="2">Isotopic ratio</th><th rowspan="2">ΔMoon -BSE %</th><th rowspan="2">±</th><th rowspan="2">References</th><th colspan="2">Isotopic</th><th rowspan="2">ΔKin(%<sub>c</sub>)</th><th rowspan="2">⟨F⟩<sub>t</sub></th><th rowspan="2">±</th><th rowspan="2">⟨F⟩<sub>t</sub></th><th rowspan="2">±</th><th rowspan="2">⟨F⟩ > <sub>v</sub><sup>c</sup></th><th rowspan="2">B<sup>b</sup></th><th colspan="2">Δ Eq liq-vap (%<sub>c</sub>) versus T (K)</th></tr><tr><th>K</th><th>Moon/BSE</th><th>masses</th><th>silicate N m<sup>-1</sup></th><th>metal N m<sup>-1</sup></th><th>N m<sup>-1</sup></th><th>silicate-vapor 1700</th><th>metal-vapor 1700</th></tr></thead><tbody><tr><td>K</td><td>1006</td><td>0.252</td><td>This study</td><td></td><td><sup>41</sup>K/<sup>39</sup>K</td><td>0.41</td><td>0.07</td><td>Wang &amp; Jacobsen 2016; Tian et al. 2020</td><td>40.96</td><td>38.96</td><td>-24.69</td><td>52</td><td>26</td><td></td><td></td><td>0</td><td>5498</td><td>0.099</td><td>0.023</td></tr><tr><td>Cu</td><td>1037</td><td>0.047</td><td>O'Neill 1991; McDonough &amp; Sun 1995</td><td></td><td><sup>65</sup>Cu/<sup>63</sup>Cu</td><td>0.50</td><td>0.10</td><td>Herzog et al. 2009</td><td>64.93</td><td>62.93</td><td>-15.51</td><td>54</td><td>27</td><td>94</td><td>47</td><td>0</td><td>2148</td><td>0.040</td><td>0.009</td><td>0.070</td><td>0.016</td></tr><tr><td>Zn</td><td>726</td><td>0.005</td><td>Albarède et al. 2015</td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td>1.12</td><td>0.50</td><td>Paniello et al. 2012; Kato et al. 2015</td><td>65.93</td><td>63.93</td><td>-15.26</td><td>122</td><td>61</td><td>62</td><td>31</td><td>0</td><td>2081</td><td>0.088</td><td>0.021</td><td>0.045</td><td>0.011</td></tr><tr><td>Ga</td><td>968</td><td>0.106</td><td>O'Neill 1991; McDonough &amp; Sun 1995</td><td></td><td><sup>71</sup>Ga/<sup>69</sup>Ga</td><td>0.32</td><td>0.20</td><td>Kato &amp; Moynier 2017; Wimpenny et al. 2022</td><td>70.92</td><td>68.93</td><td>-14.19</td><td>368</td><td>184</td><td>90</td><td>45</td><td>0</td><td>1796</td><td>0.229</td><td>0.054</td><td>0.056</td><td>0.013</td></tr><tr><td>Rb</td><td>800</td><td>0.212</td><td>This study</td><td></td><td><sup>87</sup>Rb/<sup>85</sup>Rb</td><td>0.16</td><td>0.04</td><td>Pringle &amp; Moynier 2017; Nie &amp; Dauphas 2019</td><td>86.91</td><td>84.91</td><td>-11.56</td><td>52</td><td>26</td><td></td><td></td><td>0</td><td>1189</td><td>0.021</td><td>0.005</td><td></td><td></td></tr><tr><td>Sn</td><td>704</td><td>0.031</td><td>O'Neill 1991; McDonough &amp; Sun 1995</td><td></td><td><sup>124</sup>Sn/<sup>116</sup>Sn</td><td>-0.48</td><td>0.15</td><td>Wang et al. 2019</td><td>123.91</td><td>115.90</td><td>-29.02<sup>d</sup></td><td>254</td><td>127</td><td>155</td><td>77</td><td>179</td><td>2448</td><td>0.064</td><td>0.015</td><td>-0.020</td><td>-0.005</td></tr></tbody></table>

**Notes.**

<sup>a</sup> 50% condensation temperature under solar nebula conditions (Lodders 2003).

<sup>b</sup> 1000lnβ = B<F>/T<sup>2</sup>.

<sup>c</sup> The thermodynamically favored tin vapor species is SnO. Other MVEs exist as monoatomic gas species.

<sup>d</sup> The mass of SnO is used, as it is the stable vapor species.
## 2.1. Earth

Although terrestrial rocks are readily accessible and the major, minor, and trace element compositions of tens of thousands of bulk igneous rocks have been measured, Earth is a geologically active planet, which makes estimating its bulk composition challenging. Incompatible elements are minimally fractionated in melt during partial melting, but the residual mantle can have highly fractionated incompatible element ratios. Extraction of the continental crust left behind a highly depleted mantle that is continuously molten by decompression melting at spreading centers to make mid-ocean ridge basalts (MORBs), which inherit incompatible element fractionation from their depleted source (Hofmann 1988; Sun & McDonough 1989; Salters & Stracke 2004; Gale et al. 2013). Recycling of basaltic oceanic crust and sediments at subduction zones introduces chemical and mineralogical heterogeneities in the mantle that are sampled in MORBs and oceanic island basalts (OIBs; Hofmann & White 1982; Hirschmann & Stolper 1996; Stracke et al. 2003; Pertermann et al. 2004; Sobolev et al. 2007; Stracke & Bourdon 2009; Nielsen 2010; Yang et al. 2020). During subduction, saline hydrous fluids enriched in fluid-mobile elements can be transferred to the overlying mantle wedge, imparting a peculiar signature to arc magmas and subducted slabs (Rustioni et al. 2019). These fluid-mobile incompatible elements can also be fractionated from fluid-immobile incompatible elements during weathering and riverine transport (Gaillardet et al. 2003). All these processes and others have shaped the relative abundances of incompatible elements, and no single terrestrial sample can be taken as representative of the bulk silicate Earth (BSE) composition.

The incompatible element ratio K/U has received considerable attention because (i) both K and U are important heat-producing elements in Earth (Wasserburg et al. 1964; Korenaga 2008; Arevalo et al. 2009; Huang et al. 2013; Farcy et al. 2020; O'Neill et al. 2020) and (ii) the 40 K-40Ar, 235 U-207Pb, and 238 U-206Pb decay systems are used for estimating the timing and pace of Earth's degassing (Ozima 1975; Allègre et al. 1983; Bender et al. 2008), as well as dating the formation of Earth and its core (Allègre et al. 1995; Wood & Halliday 2005; Rudge et al. 2010). Although K and U have similar behaviors during magmatic differentiation, the K/U ratio varies significantly among mantle-derived magmas and continental crust, which makes it difficult to estimate the BSE value. Recent studies have estimated the BSE K/U ratio based on mass-balance calculations between the mantle source of oceanic basalts (MORBs and OIBs) and the continental crust. A difficulty with this approach is that all those reservoirs have fractionated K/U ratios and their K and U contents are uncertain.

We use here a different approach that focuses on samples that show the least fractionated incompatible element ratios (Figure 2). To assess whether incompatible elements are fractionated or not, we rely on strongly lithophile highly refractory elements, as these can reasonably be assumed to be in chondritic proportions in the silicate Earth. We focus on La, U, Th, and Ba, as these bracket K and Rb in terms of incompatibility during peridotite melting (Sun & McDonough 1989) and include elements that have been used previously as proxies for K (La and U) and Rb (Ba) (Heier & Rogers 1963; Wasserburg et al. 1964; Hofmann & White 1983; Jochum et al. 1983; Wänke & Dreibus 1988; McDonough et al. 1992; Arevalo et al. 2009; Palme & O'Neill 2014). Highly

incompatible element behavior can be affected by many processes other than peridotite melting (Hofmann & White 1982; Hirschmann & Stolper 1996; Gaillardet et al. 2003; Stracke et al. 2003; Pertermann et al. 2004; Sobolev et al. 2007; Stracke & Bourdon 2009; Rustioni et al. 2019; Yang et al. 2020), and the elements selected provide a good cross section of incompatible element behaviors in global geochemical cycles. We did not include high field strength elements Nb and Ta in our analysis because all accessible terrestrial reservoirs have a low Nb/Ta ratio relative to chondrites (Münker et al. 2003), possibly due to partitioning of Nb in Earth's core (Wade & Wood 2001) or fractionation in rutile during dehydration/melting of subducted slabs, and storage of the residual high Nb/Ta eclogite in a deep hidden mantle reservoir (Kamber & Collerson 2000; Rudnick et al. 2000). In the following, we assume that refractory lithophile elements are in CI-like chondrite relative proportions in the BSE. Terrestrial rocks, however, have a -4.5% abundance anomaly in refractory lithophile Rare Earth Element (REE) Tm relative to CI chondrites and neighboring REEs (Dauphas & Pourmand 2015; Barrat et al. 2016). This is most likely due to the presence of a component in meteorites that is fractionated owing to nebular processes. This component could have a composition like group II meteoritic refractory inclusions, which have large positive Tm anomalies and fractionated abundances of other refractory lithophile elements. Subtracting this component from CI chondrites to account for the negative Tm anomalies in terrestrial rocks yields corrected CI-normalized concentrations for Ba, Th, U, and La of 0.98, 1, 0.99, and 0.95, respectively (CI* from Table 6 of Dauphas & Pourmand 2015). The assumption that Ba, Th, U, and La are in CI proportions within Earth is thus likely to be valid within ~5%.

In a sample, incompatible elements can be variably depleted or enriched, but what matters to us is whether their relative abundances are chondritic or not. The metric that we use to quantify the flatness of highly incompatible element abundances relative to CI is the relative standard deviation (RSD),

$$RSD = \frac{\sqrt{\frac{n\sum_i([E_i]^* - \sum_i[E_i]^*/n)^2}}{\sum_i[E_i]^*}, \quad (1)$$

where $[E_i]^*$ is the CI chondrite normalized concentration of element $E_i$ in the sample, $n$ is the number of elements considered (here 4), and the sum runs over the highly incompatible refractory lithophile elements Ba, Th, U, and La. To select the samples with the flattest abundance pattern for these four elements, we use three databases: the MORB whole rock compilation of Gale et al. (2013), the continental crust compilation (including island arcs) of Ptáček et al. (2020), and an oceanic island whole rock composition database compiled from PetDB. From these databases, we select the samples for which concentration data are available for K, Rb, Ba, Th, U, Ta, Nb, and La (1639 samples for MORB, 23,576 for continental crust and island arcs, and 1020 for ocean islands). For all the samples, we calculate the RSD metric and select the 655 samples whose RSD is below the 2.5th percentile of the population (i.e., the 2.5% samples with the flattest incompatible refractory lithophile element abundance pattern; this cutoff was selected because it is relatively stringent but retains enough samples to calculate statistically meaningful quantities)
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Estimates of the BSE K/U and Rb/Ba ratios. Three databases were used to estimate those ratios: mid-oceanic ridge basalts (MORB; Gale et al. 2013), Oceanic Islands (OI; PetDB query), and Continental Crust and a small proportion of modern Island Arc samples (CC+IA; Ptáček et al. 2020). We select from those databases the samples that contain concentration data for K, Rb, Ba, Th, U, Ta, Nb, and La. We then normalize the concentrations of refractory lithophile incompatible elements Ba, Th, U, and La to CI chondrite composition (Barrat et al. 2012) and calculate the relative standard deviation (RSD) of those four normalized concentrations. Samples containing Ba, Th, U, and La in near-CI proportions will have RSDs (Equation (1)) near zero (the incompatible element pattern is labeled as “flat” if the samples are among the 2.5% with the lowest RSD), while those with fractionated or imprecise concentrations will have high RSD. Refractory lithophile elements are expected to be in chondritic proportions in the BSE. Panels (a) and (b) show how K/U and Rb/Ba weight ratios converge to values of 9845 ± 2054 and 0.0787 ± 0.0157, respectively, as defined by the 2.5% samples with the lowest RSD (corresponding to a K/Rb ratio of 386 ± 109). Panels (c) and (d) show the joint distributions of K/U-La/U and Rb/Ba-Th/Ba (like U, La is often taken as a proxy to K during mantle melting; Th is highly incompatible like Ba and Rb). The red filled circles are the 2.5% samples with the lowest RSD, and the red filled diamond is our best estimate for the BSE. As shown, the CI-normalized La/U and Th/Ba ratios of the low-RSD samples are evenly distributed around unity, and no correlation is seen between K/U-La/U or Rb/Ba-Th/Ba, so there is no indication that our estimate could be biased. Panel (e) shows which samples belong to the group with flat incompatible Ba-Th-U-La pattern. Although we do not select against differentiated igneous rocks, almost all the samples in the “flat” subgroup are basalts. Panel (f) shows that the K/U and K/Rb ratios of samples with flat Ba-Th-U-La pattern correlate. The slightly lower K/Rb and K/U of OIBs might be due to recycling of altered oceanic crust.</figcaption></figure>
(Figures 2(a) and (b)). Although we applied no criterion to exclude differentiated igneous rocks, most of the samples with low RSD are basaltic (Figure 2(e)). Given that Ba, Th, U, and La present a range of incompatible element behaviors that bracket K and Rb, it is reasonable to assume that in that low-RSD sample set the K/U and Rb/Ba ratios are representative of the BSE.

In the samples with flat patterns, we find no correlation between K/U and La/U ratios (Figure 2(c)). The median K/U weight ratio of the 2.5% (655) flattest samples from all three subsets is 9845, while the average is $10,238 \pm 402$ ($\pm 2\sigma_{\text{mean}}$). We favor the median value over the average because it is less influenced by outliers. The 655 flattest samples comprise 95 MORB, 494 CC+IA, and 66 OI. The median K/U weight ratios of the flat subset of MORB, CC+IA, and OI are $10,600 \pm 482$, $9940 \pm 506$, and $8191 \pm 957$ ($\pm 2\sigma_{\text{mean}}$), respectively. Overall, there is good agreement between K/U ratios in the samples with flat patterns, except for the oceanic island subset, which gives a lower K/U ratio but with a large uncertainty. The flat sample subsets show a correlation between K/U and K/Rb ratios, with OIs also characterized by lower K/Rb ratios than both MORBs and CC+IA (Figure 2(f)). This correlation may be due to the presence of recycled altered oceanic crust with low K/Rb (Hart 1969; Staudigel et al. 1981; Hart & Staudigel 1982; Staudigel & Hart 1983) and possibly K/U ratios (Hofmann 1997; Lassiter 2004; Nielsen 2010; Hanyu et al. 2011) in the source of OIBs.

Lanthanum is often used as another proxy for K. The median K/La weight ratio of the 655 samples with flat patterns is 313. Assuming a CI La/U for the BSE of 30.52 (Barrat et al. 2012), we estimate the BSE K/U weight ratio to be 9563. Taking the K/La in the flat subsets of MORB, CC+IA, and OI, we get K/U ratios of 10,031, 9947, and 7584, respectively. These values are very similar to those inferred using K/U ratios directly. We take the average of the values given by K/U and K/La ratios to estimate a BSE K/U ratio of 9704 and ascribe it an uncertainty of $\pm 2054$ by calculating $2\sigma$ for all the estimates given above (9845, 10,600, 9940, 8191, 9563, 10,031, 9947, 7584).

Among the samples with flat patterns, we find no correlation between Rb/Ba and Th/U ratios (Figure 2(d)). Overall, Rb/Ba ratios agree well between the three sample subsets (MORB, CC+IA, and OI). The median Rb/Ba weight ratio of the 2.5% (655) flattest samples from all three subsets is $0.0787 \pm 0.0048$ (Ba/Rb = 12.7 $\pm$ 0.8), which we take as our preferred BSE value. For comparison, the median Rb/Ba weight ratios of the MORB, CC+IA, and OI 2.5% flattest samples are $0.0860 \pm 0.0020$, $0.0768 \pm 0.0063$, and $0.0715 \pm 0.0063$, respectively. While we can only use Ba as a proxy element for Rb, we can use both La and U as proxy elements for K, which allows us to better assess the limitations of our approach ($\pm 20\%$ uncertainty). Rather than using the $2\sigma$ of 0.0120 for the four Rb/Ba BSE estimates given above as the uncertainty, we therefore ascribe the same relative uncertainty to the BSE Rb/Ba ratio as what we have obtained for the K/U ratio ($\pm 0.0157$). We can combine the BSE Rb/Ba ratio of $0.0787 \pm 0.0157$ with a CI Ba/U ratio of 319.5 (Barrat et al. 2012) and the BSE K/U ratio of 9704 $\pm$ 2054 to estimate the BSE K/Rb weight ratio: $(\text{K/U})_{\text{BSE}}/(\text{Rb/Ba})_{\text{BSE}}/(\text{Ba/U})_{\text{CI}} = 386 \pm 112$. Using our Rb/Ba estimate, we calculate an Rb/Sr weight ratio of $0.0250 \pm 0.0050$ (i.e., an $^{87}$Rb/$^{86}$Sr atomic ratio of $0.0725 \pm 0.0145$) for the BSE assuming a chondritic CI-like Ba/Sr ratio of 0.318 (Barrat et al. 2012).

To assess the accuracy of our approach based on examination of igneous rock samples with flat CI-normalized refractory element patterns, we apply it to ratios of refractory elements whose relative abundances should be chondritic by excluding the element of interest from the flatness criterion. For example, we apply the approach outlined above to the U/Th ratio by taking the median value of the 2.5% samples with the flattest refractory lithophile element pattern while excluding U from the RSD flatness criterion (i.e., $i = Ba, Th, La$ and $n = 3$). Doing so, we obtain the following CI-normalized ratios (the element in square brackets is the one that was excluded from the criterion to select the 2.5% flattest patterns; asterisk denotes normalized to CI): $([\text{Th}]/\text{Ba})^* = 0.99$, $([\text{U}]/\text{Th})^* = 0.99$, $([\text{Th}]/\text{U})^* = 1.01$, $([\text{U}]/\text{La})^* = 0.95$. The 95% confidence interval of these ratios, considering the student t-factor, is $\sim 0.08$. This shows that this approach can reliably estimate BSE incompatible element ratios involving U and Th with an accuracy of $\sim 8\%$, which is well within the $\pm 20\%$ uncertainty quoted above.

Normalized to CI (Barrat et al. 2012), we calculate $(\text{K/U})^* = 0.136 \pm 0.029$ (i.e., a depletion in K of a factor of 7.4 $\pm$ 1.6 relative to CI) and $(\text{K/Rb})^* = 1.64 \pm 0.46$. Using a concentration for U of $20.3 \pm 4$ ng g$^{-1}$ in the BSE (McDonough & Sun 1995) and our estimate of the K/U ratio of 9704 $\pm$ 2054, we calculate a K concentration of 197 $\pm$ 55 ppm. Assuming a CI Ba/U ratio for the BSE and using our estimate of the Rb/Ba weight ratio of $0.0787 \pm 0.0157$, we calculate an Rb concentration of $0.51 \pm 0.14$ ppm. Lyubetskaya & Korenaga (2007) proposed a BSE composition that contains less refractory lithophile elements, with U and Ba concentrations given as 17.3 and 5.08 ppb, respectively. Combining these concentrations with our K/U and Rb/Ba ratios would yield K and Rb concentrations of 170 and 0.40 ppm, respectively.

Our estimated K/U ratio of 9704 $\pm$ 2054 for the BSE is within the range of values proposed in previous studies, with most values within 8000–14,000 (Figure 3) (Gast 1960; Heier & Rogers 1963; Wasserburg et al. 1964; Hurley 1968; Larimer 1971; Ganapathy & Anders 1974; Shaw 1972; Smith 1977; O’Nions et al. 1978; Jacobsen & Wasserburg 1979; Jagoutz et al. 1979; Armstrong 1981; Sun 1982; Anderson 1983; Jochum et al. 1983; McDonough & Sun 1995; Allègre et al. 2001; Lassiter 2004; Lyubetskaya & Korenaga 2007; Arevalo et al. 2009; Palme & O’Neill 2014; Farcy et al. 2020; O’Neill et al. 2020). The most recent estimates for Earth’s K/U ratio are 13,800 (Arevalo et al. 2009), 14,050 (Palme & O’Neill 2014), 12,100 (Farcy et al. 2020), and 11,820 (O’Neill et al. 2020). These values overlap with our estimate but seem to favor higher K/U ratios. Arevalo et al. (2009), Farcy et al. (2020), and O’Neill et al. (2020) rely on a global mass balance between continental crust, MORB, and OIB sources, which is difficult to establish because K/U is fractionated in these reservoirs and the K content in each reservoir is uncertain. As discussed above, our approach recovers a CI U/La ratio when U is not used in the flatness criterion. The validity of other approaches could be assessed by similarly verifying whether they can recover a near-Cl U/La ratio for the BSE.

Our estimated Rb/Sr ratio of $0.025 \pm 0.005$ for the BSE is within the range of values proposed in previous studies of $\sim 0.02$ to 0.035 (Figure 3; Gast 1960; Hurley 1968; Larimer 1971; Shaw 1972; Ganapathy & Anders 1974;
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison between this study (Table 1; red filled circles with error bars) and literature estimates of K/U and Rb/Sr ratios in the BSE and bulk silicate Moon. Our estimates rely on geochemical similarities between Rb-Ba-Th and K-U-La by focusing on samples that have near-chondritic proportions of refractory lithophile elements Ba-Th-U-La (low-RSD samples; Figures 2, 4) or for the Moon by interpolating Rb/Ba-Th/Ba and K/U-La/U correlations to chondritic Th/Ba and La/U ratios (Figure 4). (a) K/U ratio of BSE; mostly estimated assuming constant K/U or K/La ratios (Gast 1960; Heier & Rogers 1963; Wasserburg et al. 1964; Hurley 1968; Larimer 1971; Shaw 1972; Ganapathy & Anders 1974; Smith 1977; O’Nions et al. 1978; Jacobsen & Wasserburg 1979; Jagoutz et al. 1979; Armstrong 1981; Sun 1982; Anderson 1983; Jochum et al. 1983; McDonough & Sun 1995; Allègre et al. 2001; Lassiter 2004; Lyubetskaya & Korenaga 2007; Arevalo et al. 2009; Palme & O’Neill 2014; Farcy et al. 2020; O’Neill et al. 2020). (b) K/U ratio of the bulk Moon; mostly estimated assuming constant K/U or K/La ratios (LSPET 1969; Tera et al. 1970; Larimer 1971; Ganapathy & Anders 1974; Morgan et al. 1978; Wänke 1981; Taylor 1982; O’Neill 1991; McDonough et al. 1992; Warren 2005; Taylor & Wieczorek 2014; Albarède et al. 2015). (c) Rb/Sr ratio of BSE; after 1976 often estimated from $^{87}$Sr-$^{143}$Nd systematics (Gast 1960; Hurley 1968; Larimer 1971; Shaw 1972; Ganapathy & Anders 1974; DePaolo & Wasserburg 1976; Jahn & Nyquist 1976; O’Nions et al. 1977; Smith 1977; O’Nions et al. 1978; Allègre et al. 1979; Jacobsen & Wasserburg 1979; Jagoutz et al. 1979; Armstrong 1981; Sun 1982; Anderson 1983; Hofmann & White 1983; McDonough et al. 1992; McDonough & Sun 1995; Allègre et al. 2001; Lyubetskaya & Korenaga 2007; McDonough & Arevalo 2008; Palme & O’Neill 2014). (d) Rb/Sr ratio of bulk silicate Moon; mostly estimated assuming constant Rb/Ba ratio (Larimer 1971; Ganapathy & Anders 1974; Morgan et al. 1978; Taylor 1982; O’Neill 1991; McDonough et al. 1992; Taylor & Wieczorek 2014; Albarède et al. 2015).</figcaption></figure>

Jahn & Nyquist 1976; DePaolo & Wasserburg 1976; O’Nions et al. 1977; Smith 1977; O’Nions et al. 1978; Allègre et al. 1979; Jacobsen & Wasserburg 1979; Jagoutz et al. 1979; Armstrong 1981; Sun 1982; Anderson 1983; Hofmann & White 1983; McDonough et al. 1992; McDonough & Sun 1995; Allègre et al. 2001; Lyubetskaya & Korenaga 2007; Palme & O’Neill 2014; McDonough & Arevalo 2008). Assuming an initial solar system $^{87}$Sr/$^{86}$Sr ratio of ~0.69879 (0.69876 ± 0.00002 for ALL, Gray et al. 1973; 0.69883 ± 0.00002 for ADOR, Wasserburg et al. 1977) and assuming a time of Rb/Sr fractionation of 4.56 Gyr, our BSE Rb/Sr ratio would correspond to a present-day BSE $^{87}$Sr/$^{86}$Sr ratio of 0.7035 ± 0.0010. The most precise estimates of the BSE Rb/Sr ratio either use the Rb/Ba ratio and assume a chondritic Ba/Sr ratio (Hofmann & White 1983) or interpolate the correlation observed in mantle-derived samples between $^{87}$Sr/$^{86}$Sr and $^{143}$Nd/$^{144}$Nd ratios to a chondritic $^{143}$Nd/$^{144}$Nd ratio (DePaolo & Wasserburg 1976; O’Nions et al. 1977; Allègre et al. 1979). Both approaches yield an Rb/Sr ratio of ~0.03, which is identical within error to the value that we obtain. As discussed by Hofmann & White (1983), the isotopic approach based on the Sr-Nd is not intrinsically superior to using trace element abundances, as significant scatter is present in the Sr-Nd mantle array and the assumption that the BSE should fall on the best-fit line to the mantle array is not necessarily justified (Armstrong 1981; White & Hofmann 1982; Hofmann & White 1983).

## 2.2. The Moon

Given the high incompatibility of K, Rb, and U, the assumption is often made that one can take the average K/U and K/Rb ratios of lunar rocks as direct measures of the composition of the bulk Moon. A potential difficulty with this approach is that lunar magma ocean (LMO) crystallization involved the formation of a plagioclase-rich flotation crust and dense mafic cumulates in which nonnegligible amounts of those elements could have partitioned. For example, terrestrial anorthosites and anorthositic gabbros are known to have high K/Rb ratios, due to the preferential partitioning into plagioclase of K relative to Rb (Philpotts & Schnetzler 1970).

To assess the effect of LMO differentiation on K/U, Rb/Ba, and K/Rb ratios, we have compiled the concentrations of K, Rb, Ba, Th, U, and La in 112 mare basalts and eight KREEP-rich Apollo 14 samples (14148–14149–14156, 14163, 14006, 14141, 14305, 14301, 14303, and 14321) from the Lunar Sample Compendium (LSC; Meyer 2012). Obvious outlier concentrations were eliminated before calculation of sample averages. The data compiled in the LSC were measured using diverse techniques, often on different fragments of the same rock. To evaluate whether averaging such disparate data can be a cause of inaccuracy or imprecision, we have compiled Ba, Th, U, La, and Rb concentration data of mare basalts measured by Inductively Coupled Plasma Mass Spectrometry (ICPMS) on homogenized powder by Neal (2001). ICPMS is the modern gold standard for precise and accurate concentration analysis of a variety of trace elements on single sample aliquots. As
<figure><img src="image_7.png" /></figure>

**Figure 4.** Correlations between K/U-La/U, Rb/Ba-Th/Ba, and K/Rb-U/Ba in lunar basalts (blue circles) and KREEP-rich soils and breccias from Apollo 14 (yellow circles). All these elements (K, Rb, Ba, Th, U, and La) are highly incompatible during partial melting and magmatic differentiation. The purple diamonds are calculated compositions for a model of LMO crystallization corresponding to (1) the KREEP liquid after 95% LMO crystallization and (2) a mixture of 98% mafic cumulate+1% plagioclase+1% KREEP liquid that would be representative of the source of mare basalts (Snyder et al. 1992; see text for details). As shown, the calculated KREEP and mare basalt compositions fall within the range of compositions measured in lunar samples, supporting the view that the trends depicted in this diagram reflect LMO differentiation and more specifically plagioclase flotation. Ba, Th, La, and U are highly refractory elements, and the La/U, Th/Ba, and U/Ba ratios are expected to be approximately chondritic (CI-like) in bulk planetary objects. Extrapolating the K/U vs. La/U relationship (a) to a CI-like La/U ratio yields a bulk lunar K/U weight ratio of 2448 ± 140 g g⁻¹. Similarly, extrapolating the Rb/Ba vs. Th/Ba relationship (b) to a CI-like Th/Ba ratio yields a bulk lunar Rb/Ba weight ratio of 0.0167 ± 0.0018 g g⁻¹. Combining the bulk Moon K/U and Rb/Ba ratios yields a bulk Moon K/Rb ratio of 459 ± 56 g g⁻¹, consistent with the values measured in lunar samples with U/Ba ratios close to the CI value. The estimated bulk lunar weight ratios are displayed as red squares. Another approach to estimating bulk Moon K/U, Rb/Ba, and K/Rb ratios is to focus on samples that have Ba, Th, U, and La in near-chondritic proportions. The degree to which the samples depart from CI proportions is quantified using the RSD of the CI-normalized abundances of these elements (100 RSD; Equation (1); the samples highlighted with thick contours are the 25% with the lowest RSD). As shown, the samples with low RSD have K/U (D), Rb/Ba (E), and K/Rb (F) ratios that are scattered around the bulk Moon values inferred from the regressions depicted in panels (a) and (b). The data were compiled from the LSC (Meyer 2012). The mare basalt samples are all those for which K, Rb, Ba, La, and U data are reported in the compendium (n = 113). The Apollo 14 KREEP-rich samples are those used by Warren & Wasson (1979) to define the KREEP end-member composition and for which bulk K, Rb, Ba, and U data are reported in the compendium (n = 8; samples 14148–14149–14156, 14163, 14006, 14141, 14305, 14301, 14303, 14321).

discussed below, we use K/U, Rb/Ba, K/Rb, La/U, Th/Ba, and U/Ba ratios to estimate the degrees of depletion in K and Rb of the Moon (Figure 4). In Figure 5, we compare the ratios from the ICPMS data set of Neal (2001) (K is from the LSC) with those for the same samples compiled from the LSC (Meyer 2012) but excluding the values from Neal (2001). This comparison shows that the ICPMS data (Figures 5(a)–(c)) show similar dispersion to the non-ICPMS data, and the values generally agree (Figures 5(d)–(i)). Because our LSC compilation contains approximately twice as many samples as those reported by Neal (2001) and ICPMS data show similar scatter to the non-ICPMS data, we use the full LSC compilation (including the data from Neal 2001) in our estimate of the composition of the bulk Moon.

The La/U, U/Ba, and Th/Ba ratios in the bulk Moon should be approximately chondritic given that these elements are refractory and lithophile. These elements could, however, have been fractionated by LMO differentiation. In Figures 4(a)–(c), we plot K/U weight ratio versus CI-normalized La/U ratio,

K/Rb weight ratio versus CI-normalized U/Ba ratio, and Rb/Ba weight ratio versus CI-normalized Th/Ba ratio for the compiled lunar samples, and we find relatively well-defined correlations: $(K/U)_{wt} = (1983 \pm 282)(La/U)_{CI} + (466 \pm 366)$, $(Rb/Ba)_{wt} = (0.00793 \pm 0.00449)(Th/Ba)_{CI} + (0.00881 \pm 0.00414)$, and $(K/Rb)_{wt} = (-387 \pm 119)(U/Ba)_{CI} + (898 \pm 113)$. By regressing the lunar data to CI La/U, Th/Ba, and U/Ba ratios, we estimate the bulk Moon to have the following composition: $K/U = 2448 \pm 140$, $K/Rb = 512 \pm 43$, and $Rb/Ba = 0.0167 \pm 0.0018$. Assuming that the bulk Moon Ba/U ratio is like CI, we can also combine the bulk Moon K/U and Rb/Ba ratios estimated from the regressions to calculate a $K/Rb$ ratio of $459 \pm 56$, which agrees within error with the value obtained by regressing $K/Rb$ versus U/Ba but should be more accurate.

Another approach for calculating the bulk Moon composition is to focus on the samples that have refractory lithophile incompatible element abundances closest to CI. As was done for terrestrial samples, we calculated the RSD flatness metric
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between La-U-Th-Ba-Rb data from the LSC (Meyer 2012) and high-precision ICPMS data from Neal (2001). There are 120 entries in our mare basalt+KREEP compilation from the LSC that were measured by a variety of techniques (e.g., OES, NAA, IDMS, XRF, radiation counting, ICPMS, wet chemical analysis, AAS), including some that have relatively low precision and high detection limit. Neal (2001) reported analyses of the compositions of 65 mare basalts by the technique of ICPMS, which has high precision and low detection limit. K₂O was not measured in Neal (2001), and the values used here are all averages from the LSC. Panels (a)–(c) are like panels (a)–(c) from Figure 4, but using only the ICPMS data from Neal (2001). The red lines are the regressions through the whole LSC data set, and the red squares are our best estimates of bulk silicate Moon values (Figures 4(a)–(c)). As shown, the data from Neal (2001) are consistent with the LSC compilation, and data scatter is not improved by using only ICPMS data. Panels (d)–(i) compare elemental ratios obtained by ICPMS from single aliquots (Neal 2001) (x-axis) with ratios for the same rocks calculated by taking averages of concentrations reported in various studies compiled in the LSC (obvious outliers were excluded from sample averages) but excluding ICPMS data from Neal (2001) (y-axis). The dashed lines are 1:1. As shown, non-ICPMS data compiled from the LSC agree with ICPMS data from Neal (2001).</figcaption></figure>

for lunar basalts and KREEP samples based on CI-normalized Ba, Th, U, and La concentrations. For terrestrial samples, we only considered the 2.5% samples with the flattest abundance profile. The number of basalts and KREEP samples in our database is only 120, so we relax this criterion and consider the 25% samples with the most chondritic compositions (33 samples in total; Figures 4(d)–(f)). We calculate the following median weight ratios for that subset of samples: K/U = 2443 ± 189, K/Rb = 443 ± 64, and Rb/Ba = 0.0165 ± 0.0041. Assuming a CI Ba/U ratio, the K/U and Rb/Ba ratios of the bulk silicate Moon yield a K/Rb ratio of 464 ± 122. These values are very similar to the ones calculated using the

regression approach. Therefore, we adopt the regressed values as our best estimates for the bulk Moon (K/U = 2448 ± 140, K/Rb = 459 ± 56, and Rb/Ba = 0.0167 ± 0.0018).

We compare in Figure 3 our estimates of bulk Moon K/U and Rb/Sr ratios with previously reported values from the literature. Our Rb/Sr ratio of 0.0053 ± 0.0006 is within the range of previously reported values (Larimer 1971; Ganapathy & Anders 1974; Morgan et al. 1978; Taylor 1982; O'Neill 1991; McDonough et al. 1992; Taylor & Wieczorek 2014; Albarède et al. 2015). Our K/U ratio of 2448 ± 140 is also within the range of previously reported values (LSPET 1969; Tera et al. 1970; Larimer 1971; Ganapathy & Anders 1974; Morgan et al. 1978;
<table><caption>**Table 3**
Partition Coefficients Used in Modeling Trace Element Behavior during LMO Crystallization</caption><thead><tr><th>D<sub>mineral-melt</sub></th><th>Plagioclase</th><th>Clinopyroxene</th><th>Orthopyroxene (+Pigeonite)</th><th>Olivine</th><th>Ilmenite</th><th>References</th></tr></thead><tbody><tr><td>K</td><td>0.39</td><td>0.03</td><td>0.01</td><td>0.007</td><td>0</td><td>1 (pl), 2 (cxp, opx/pig,ol)</td></tr><tr><td>Ba</td><td>0.087</td><td>0.04</td><td>0.01</td><td>0.01</td><td>0</td><td>1 (pl), 2 (cxp, opx/pig,ol)</td></tr><tr><td>Th</td><td>0.005</td><td>0.005</td><td>0.003</td><td>0.0000058</td><td>0.001</td><td>3</td></tr><tr><td>U</td><td>0.00042</td><td>0.0064</td><td>0.0031</td><td>0.00001</td><td>0.0039</td><td>3</td></tr><tr><td>La</td><td>0.04<sup>a</sup></td><td>0.054</td><td>0.002</td><td>0.0004</td><td>0.0072</td><td>4 (cxp, opx/pig,ol), 5 (ilm)</td></tr><tr><td>Rb</td><td>0.084</td><td>0.02</td><td>0.02</td><td>0.01</td><td>0</td><td>1 (pl), 2 (cxp, opx/pig,ol)</td></tr></tbody></table>

**References.** (1) Sun et al. 2017; (2) Philpotts & Schnetzler 1970; (3) Fonseca et al. 2014; (4) McKenzie & O’Nions 1991; (5) Nielsen et al. 1992.

Wänke 1981; Taylor 1982; O’Neill 1991; McDonough et al. 1992; Warren 2005; Taylor & Wieczorek 2014; Albarède et al. 2015) but is significantly higher than recent estimates that seem to cluster around ~2000 (O’Neill 1991; McDonough et al. 1992; Warren 2005; Taylor & Wieczorek 2014; Albarède et al. 2015). An approach often used to estimate bulk Moon K/U or K/La ratios involves plotting the concentrations of K against La or U, and taking the slope of the correlation forced through the origin (Wänke 1981; O’Neill 1991; Taylor & Wieczorek 2014). In such diagrams, the KREEP samples rich in incompatible elements have more leverage to define the slope. The La/U ratio of KREEP is, however, lower than chondritic, and the correlation between K/U and La/U ratios in lunar rocks (Figure 4) suggests that the K/U ratio of KREEP is also lower than the bulk Moon value owing to removal of a high K/U plagioclase flotation crust. This could explain why our bulk silicate Moon K/U ratio is higher than some recent estimates. The K/Th ratio of the lunar surface was measured remotely by gamma-ray spectroscopy, yielding values of 359–368 (Prettyman et al. 2006) and 346 ± 60 (Kobayashi et al. 2010). Assuming a chondritic U/Th ratio, these values correspond to K/U ratios of 1319–1352 and 1272 ± 221, respectively, which is significantly lower than the value calculated here. The reason for this discrepancy is unknown.

To calculate the concentrations of K, U, and Th in the Moon, we use the bulk silicate Moon Th concentration of 66 ppb reported in Warren (2005), a Th/U ratio of 3.67 (Pourmand & Dauphas 2010; Barrat et al. 2012), and the K/U ratio of 2448 reported here. We thus have Th = 66 ppb, U = 18.0 ppb, and K = 44.0 ppm. Using a similar approach, we obtain Rb = 0.0958 ppm, Sr = 18 ppm, and Ba = 5.74 ppm (Table 1).

The correlations highlighted above between ratios involving K and Rb, on the one hand, and refractory lithophile incompatible elements Ba, Th, U, and La, on the other hand, must reflect magmatic processes. A likely culprit is LMO crystallization. To test this idea, we have quantitatively modeled trace element behavior during LMO crystallization by implementing the crystallization model of Snyder et al. (1992) in a Mathematica program (provided as Supplementary Online Material). We do not expect more recent and realistic models of LMO crystallization (Elkins-Tanton et al. 2011; Lin et al. 2017; Charlier et al. 2018) to yield drastically different magmatic evolutions for the highly incompatible elements that we are primarily interested in. In the Snyder et al. (1992) model, the crystallization of the LMO is divided into five stages. Stage 1 (0%–40% crystallization) involves batch equilibrium crystallization of olivine. Stage 2 (40%–78%) also involves batch equilibrium crystallization. It keeps olivine at 40%, and the fraction of orthopyroxene steadily increases until 78% crystallization (at the end of stage 2, olivine represents

40/78 = 51% of the total solid, while orthopyroxene represents 49% of the solid). Stage 3 (78%–86%) involves fractional crystallization of 53% plagioclase, 25% olivine, and 22% pigeonite. Stage 4 (86%–95%) involves fractional crystallization of 38% clinopyroxene, 36% plagioclase, and 26% pigeonite. Stage 5 (95–99.5%) involves fractional crystallization of 24% clinopyroxene, 31% plagioclase, 34% pigeonite, and 11% ilmenite. To calculate how elemental ratios would evolve during differentiation of the LMO, one must know the mineral/melt partition coefficients of all the minerals involved. We take as initial LMO composition the bulk Moon concentrations used in the calculation are summarized in Table 3, and more details are provided hereafter.

We assume that the partition coefficients for pigeonite are identical to those of orthopyroxene. The partition coefficients for U in plagioclase, clinopyroxene, orthopyroxene (and pigeonite), olivine, and ilmenite are $D_{plg-melt}^U = 0.00042 D_{cxp-melt}^U = 0.0064$, $D_{opx-pig-melt}^U = 0.0031$, $D_{ol-melt}^U = 0.00001$, and $D_{ilm-melt}^U = 0.0039$ (Fonseca et al. 2014). The partition coefficients for Th in plagioclase, clinopyroxene, orthopyroxene (and pigeonite), olivine, and ilmenite are $D_{plg-melt}^{Th} = 0.005$, $D_{cxp-melt}^{Th} = 0.005$, $D_{opx, pig-melt}^{Th} = 0.003$, $D_{ol-melt}^{Th} = 0.0000058$, and $D_{ilm-melt}^{Th} = 0.001$ (Fonseca et al. 2014). For K, Rb, and Ba in plagioclase, we use the following values: $D_{plg-melt}^K = 0.39$, $D_{plg-melt}^Rb = 0.084$, and $D_{plg-melt}^{Ba} = 0.086$ (Sun et al. 2017). For La in plagioclase, we use a value of 0.04, which is higher than the value proposed by Sun et al. (2017) for plagioclase in the LMO but is still acceptable given present uncertainties. The reason for taking a higher value is that otherwise we cannot explain (i) the coupling between K and La in KREEP and mare basalts and (ii) the subchondritic La/U ratio of KREEP. The other partition coefficients for K, Rb, and Ba are from Philpotts & Schnetzler (1970): $D_{ol-melt}^K = 0.007$, $D_{ol-melt}^Rb = 0.01$, $D_{ol-melt}^{Ba} = 0.01$, $D_{opx, pig-melt}^K = 0.01$, $D_{opx, pig-melt}^Rb = 0.02$, $D_{opx, pig-melt}^{Ba} = 0.01$, $D_{cxp-melt}^K = 0.03$, $D_{cxp-melt}^Rb = 0.02$, $D_{cxp-melt}^{Ba} = 0.04$, $D_{cxp-melt}^K = 0.03$, $D_{cxp-melt}^Rb = 0.02$, $D_{cxp-melt}^{Ba} = 0.04$. The partition coefficients for La in orthopyroxene, clinopyroxene, and olivine are from McKenzie and O’Nions (1991): $D_{cxp-melt}^{La} = 0.054$, $D_{opx, pig-melt}^{La} = 0.002$, and $D_{ol-melt}^{La} = 0.0004$. The La partition coefficient in Ilmenite is from Nielsen et al. (1992): $D_{ilm-melt}^{La} = 0.0072$. To our knowledge, there are no partition data for K, Rb, and Ba in ilmenite because these elements are highly incompatible in this mineral, so we set their D values to 0.

Using these partition coefficients, our estimated bulk Moon composition, and the proportions of the minerals crystallizing

<footnote><sup>a</sup> The partition coefficient of La in plagioclase is higher than that given by Sun et al. (2017) for reasons explained in the text.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Maps of K/Th ratio normalized to the BSE value (top panel) and K concentration normalized to the maximum value (bottom panel). Far-side K/Th values show more noise owing to lower K and Th concentrations, but the K/Th ratio is the same on the near and far sides. K depletion is a global feature of the Moon. Data from Lunar Prospector (Feldman et al. 1998; Lawrence et al. 1998).</figcaption></figure>

from Snyder et al. (1992), we calculated the concentrations of K, Rb, Ba, Th, U, and La in the minerals and melt that are expected to have formed during the crystallization of the LMO. KREEP is the liquid remaining after ~95%-99% crystallization. In the K/U versus La/U, Rb/Ba versus Th/Ba, and K/Rb versus U/Ba diagrams, the calculated KREEP composition plots toward the field defined by KREEP-rich samples but with a less fractionated composition (Figure 4). The source of mare basalts is more complex, as it involves mixing of (i) mafic cumulate minerals formed after plagioclase crystallization and flotation, (ii) pure plagioclase, and (iii) trapped liquid. We follow Snyder et al. (1992) and calculate the composition corresponding to mixing between 98% mafic cumulate (itself a mixture of 80% mafic cumulate at 86% LMO crystallization and 20% mafic cumulate at 95% crystallization), 1% plagioclase (at 95% LMO crystallization), and 1% KREEP liquid (at 95% crystallization). In the K/U versus La/U, Rb/Ba versus Th/Ba, and K/Rb versus U/Ba diagrams, the calculated mare basalt source composition plots within the field of mare basalt samples (Figure 4). Mare basalts have high K/U and La/U ratios, and low Th/Ba, Rb/Ba, U/Ba, and K/Rb ratios because they contain significant amounts of the plagioclase and mafic cumulate components. KREEP-rich

samples have low K/U and La/U ratios, and high Th/Ba, Rb/Ba, U/Ba, and K/Rb ratios because KREEP is the residual liquid produced by crystallization and sequestration of mafic cumulate and plagioclase. We conclude that LMO crystallization is most likely responsible for the chemical fractionation of incompatible elements K, Rb, Ba, Th, U, and La documented in our compilation of lunar rock compositions.

To assess whether K depletion was a global feature or was biased by the sampling of the Procellarum region by the Apollo mission as was recently suggested by Tartèse et al. (2021), we mapped the K/Th ratio of the Moon using data from Lunar Prospector (Lawrence et al. 1998; Feldman et al. 1998), and we compiled K, La, and U data of lunar meteorites from the lunar meteorite compendium. As shown in Figure 6, the K/Th ratio is the same on the near and far sides. There is more noise in the K/Th ratio on the far side owing to lower concentrations of K and Th, except for a small region around the South Pole-Aitken (Moriarty et al. 2021; Figure 6, bottom panel; Figures 7(a), (b)), but the K/Th ratio is the same. Lunar meteorites, which should be less biased than Apollo missions in their sampling, also show uniformly low K/U ratios relative to BSE (Figure 7(c)). The relatively small variations in this ratio among lunar meteorites could be due to magmatic
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Potassium depletion at the lunar surface. (a) K/Th ratio as a function of great-circle distance from the center of the Procellarum basin (data from Lunar Prospector; Lawrence et al. 1998; Feldman et al. 1998). As shown, the depletion in K (normalized to refractory lithophile incompatible element Th) is uniform across the lunar surface. (b) K concentration as a function of great-circle distance from the center of the Procellarum basin. The KREEP region and its antipode in the South Pole–Aitken basin show less dispersion owing to higher K and Th concentrations (Figure 6 bottom panel). (c) K/U vs. La/U ratio diagram for lunar meteorites, which sample a more representative area of the Moon, including the far side (data compiled from the lunar meteorite compendium). Lunar meteorites show the same level of K depletion as what is seen in Apollo mare basalts, demonstrating that the latter are not biased by their belonging to the Procellarum basin.</figcaption></figure>

differentiation and terrestrial weathering during residence at Earth's surface. To summarize, available evidence points to the fact that the depletion in K is laterally uniform at the surface of the Moon. Mare basalts are thought to have formed by melting of the lunar mantle at depths of ~200–400 km (Kesson & Lindsley 1976; Grove & Krawczynski 2009 and references therein). All mare basalts are depleted in K, and as discussed above, variations in K/U or K/La can be explained by magmatic differentiation in the LMO. Available evidence thus suggests that the depletion in K extends to a depth of at least ~400 km. The degree of depletion for K and Rb calculated here (Figure 4, Table 1) is thus relevant to at least the upper 400 km of the Moon. Tartèse et al. (2021) argued for a preferential depletion of MVEs in the near side by a large impact in the Procellarum region, which is not supported by available data. Furthermore, the GRAIL results show that Procellarum is most likely not an impact structure (Andrews-Hanna et al. 2014).

## 2.3. Mars and Vesta

The compositions of Vesta and Mars are estimated based on analyses of HED (Howardite-Eucrite-Diogenite) and SNC (Martian; Shergottite-Nakhlite-Chassignite) meteorites, respectively. There is no clear difference between the compositions of meteorite finds and falls for these two meteorite groups, so finds and falls are considered together. We use primarily the compilations of Kitts & Lodders (1998) for eucrites and of Lodders (1998) for Martian meteorites.

Martian (SNC) meteorites reveal a correlation between the K/U and La/U ratios, which interpolates to a K/U weight ratio of 17.331 ± 1771 at a CI La/U ratio (Figure 8(a)). We also find a weak correlation between Rb/Ba and Th/Ba ratios that interpolates to an Rb/Ba weight ratio of 0.182 ± 0.036 at a CI Th/Ba ratio. (Figure 8(b)) Combining these estimates of the Martian K/U and Rb/Ba ratios with a CI U/Ba ratio, we estimate the bulk Mars K/Rb weight ratio to be 298 ± 66, which is like the values measured in SNC meteorites (Figure 8(c)). Assuming a CI Ba/Sr ratio, the Rb/Ba ratio of 0.182 ± 0.036 gives a Rb/Sr weight ratio for bulk Mars of 0.0579 ± 0.0115.

Using gamma-ray spectroscopy, Taylor et al. (2006) estimated that the K/Th global ratio at the surface of Mars

was 5330 ± 220, corresponding to a K/U ratio of 19.589 ± 809, which is within error identical to the value inferred from Martian meteorites. Yoshizaki & McDonough (2020) used Mars Odyssey gamma-ray spectroscopy data and $^{87}$Sr-$^{143}$Nd isotopic systematics in SNC meteorites to estimate bulk Mars ratios of K/U = 20,000, Rb/Ba = 0.22, and K/Rb = 300, which are broadly consistent with the values that we derive. To calculate the concentrations of heat-producing elements in the Martian mantle, we use the U concentration of Taylor (2013) of 16 ppb and the elemental ratios derived here. We thus have U = 16 ppb, Th = 59 ppb, and K = 27.7 ppm for bulk silicate Mars (Table 1). Adopting the U concentration of 18 ppb proposed by Yoshizaki & McDonough (2020) would not change the results much.

Vesta's HED meteorites comprise three groups: howardite, eucrite, and diogenite. Among these, eucrites are basaltic in composition and have a simpler history than both howardite and diogenite meteorites. The K/U ratio in eucrites correlates with the La/U ratio (Figure 9(a)). Similarly, the Rb/Ba ratio correlates with the Th/Ba ratio (Figure 9(b)). Interpolating to CI La/U and Th/Ba ratios, we calculate a K/U weight ratio of 4140 ± 805 and an Rb/Ba weight ratio of 0.00895 ± 0.00172 for bulk Vesta. Combining these values with a CI Ba/U ratio, we calculate a K/Rb weight ratio of 1448 ± 396. This value is within the range of measured K/Rb ratios in HED meteorites (Figure 9(c)). The high-precision measurements of Tera et al. (1970) give a similar high K/Rb ratio, suggesting that this estimate is reliable. The K/Th ratio of Vesta's regolith measured remotely by gamma-ray spectroscopy is 900 ± 400 (Prettyman et al. 2015), corresponding to a K/U ratio of 3308 ± 1470. This is in good agreement with the value estimated here based on eucrite measurements. Our estimated Rb/Ba ratio (0.00895 ± 0.00172) is also close to the value (0.0132) obtained by Mittlefehldt et al. (2013). Assuming a CI Ba/Sr ratio, the Rb/Ba ratio of 0.00895 ± 0.00172 gives a Rb/Sr weight ratio for bulk Vesta of 0.00285 ± 0.00055. To calculate the concentrations of elements in the mantle of Vesta, we use the inferred Al$_2$O$_3$ concentration from Toplis et al. (2013) (2.70 wt%) together with CI chondrite U/Al and Th/Al ratios (Barrat et al. 2012) and our estimated K/U ratio. We thus calculate U = 13.9 ppb, Th = 51.1 ppb, and K = 57.5 ppm for bulk silicate Vesta (Table 1).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Evaluation of the degree of K and Rb depletions in Mars. The data points represent the compositions of SNC meteorites (Lodders 1998). The K/U and Rb/Ba ratios of Mars are calculated (red squares) by interpolating the data to CI chondrite La/U (panel (a)) and Th/Ba (panel (b)) ratios. The K/Rb ratio calculated from those interpolations and assuming a CI U/Ba ratio is consistent with the K/Rb ratio measured in SNC meteorites (panel (c)).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Evaluation of the degree of K and Rb depletions in Vesta. The data points represent the compositions of eucrites (Kitts & Lodders 1998). The K/U and Rb/Ba ratios of Vesta are calculated (red squares) by interpolating the data to CI chondrite La/U and Th/Ba ratios. The K/Rb ratio calculated from those interpolations and assuming a CI U/Ba ratio is ~6 times higher than CI but is consistent with the K/Rb ratio measured in HED meteorites (panel (c)).</figcaption></figure>

## 2.4. The Angrite Parent Body

Angrites are the most volatile-element-depleted/refractory-element-enriched meteorites available. Their extreme depletions in K and Rb make determination of the concentrations of these elements particularly challenging. Another difficulty is that all angrites except Angra dos Reis (AdoR hereafter) are meteorite finds, meaning that they have been exposed to terrestrial weathering for extended periods of time. Earth's crust is very rich in K and Rb, so the concentrations measured

in most angrites are unreliable owing to terrestrial contamination. To estimate the degree of K and Rb depletions in the angrite parent body, we therefore rely solely on AdoR. Great attention was paid in previous studies to measuring the concentration of Rb in AdoR (and mitigating any terrestrial contamination), as the Sr isotopic composition of this meteorite is used to define the solar system initial $^{87}$Sr/$^{86}$Sr ratio. The K/Rb weight ratio of AdoR is 447 ± 131, which is the average of the following values: 508, 419 (Wassergburg et al. 1977), and 415 (Tera et al. 1970). The Rb/Sr weight
ratio of AdoR is 0.000 287 ± 0.000 142, which is the average of the following values: 0.000 502, 0.000 284, 0.000 471, 0.000 298 (Wasserburg et al. 1977), 0.000 094 4, 0.000 217, and 0.000 146 (Lugmair & Galer 1992). Application of the Chauvenet criterion flagged the value of 0.000 854 from Nyquist et al. (1994) as a possible outlier, so it was not included in the average. The Sr concentration of AdoR is 137 ± 8 ppm, which is the average of the following values: 166, 133, 135, 136 (Wasserburg et al. 1977), 133 (Tera et al. 1970), 142.0, 139.3, 126.0 (Lugmair & Galer 1992), 134 (Nyquist et al. 1994), and 123 ppm (Riches et al. 2012). If we assume that the angrite parent body has chondritic U/Sr ratio, the U concentration should be 0.136 ppm in AdoR. The measured U concentrations in AdoR are 0.238 (Tissot et al. 2017), 0.20 (Riches et al. 2012), and 0.164 ppm (Wasserburg et al. 1977). There are variations in reported U concentrations, but the value of Wasserburg et al. (1977) is close to that calculated assuming a CI U/Sr ratio. We calculate the K/U weight ratio in AdoR as $(K/U)_{AdoR} = (K/Rb)_{AdoR} \times (Rb/Sr)_{AdoR} \times (Sr/U)_{CI} = 129 \pm 74$. The CI-normalized K/Rb and K/U ratios are 1.89 ± 0.55 and 0.001 81 ± 0.00104. The angrite parent body is thus depleted by a factor of ~554 in K relative to CI. The concentrations given above for AdoR are for this specific meteorite and not the parent body, but given the incompatible element nature of K, Rb, U, and Sr, it is reasonable to assume that the K/U, K/Rb, and Rb/Sr that we report in Table 1 are representative of the angrite parent body as a whole.

## 2.5. Mercury

The MESSENGER mission measured the K, U, and Th concentrations of rocks at the surface of Mercury and obtained the following values: $K = 1288 \pm 234$ ppm, $U = 0.090 \pm 0.020$ ppm, and $Th = 0.155 \pm 0.054$ ppm (Prettyman et al. 2015; Nittler et al. 2018). The extent to which these values are representative of bulk Mercury is unknown but taken at face value; the Mercury K/U weight ratio is $14311 \pm 4108$. The CI-normalized K/U ratio is thus $0.200 \pm 0.058$, corresponding to a depletion in K of a factor of $5.0 \pm 1.4$ relative to CI.

# 3. Calculation of Equilibrium and Kinetic Isotopic Fractionation during Evaporation and Condensation

## 3.1. Equilibrium Fractionation

The extent of equilibrium isotopic fractionation between the condensed phase (liquid for the relatively high pressure conditions relevant to the Moon-forming impact and solid for low-pressure nebular conditions) and vapor will depend on temperature and speciation in the gas and liquid/solid. Both metal and silicate liquid could have coexisted in the aftermath of the Moon-forming impact and in the growing Moon. While K and Rb are strongly lithophile, Cu, Zn, Ga, and Sn can partition in metal (Capobianco et al. 1999; Righter et al. 2010; Mahan et al. 2017; Righter 2019). However, judging from the size of the lunar core (~1% of the lunar mass; Williams et al. 2014), little metal would have been present in the protolunar disk or the growing Moon. Furthermore, any metal present would have segregated from silicate in the Moon or moonlet cores, and in the disk midplane, such that only silicate liquid would have been in contact with vapor. It is thus likely that silicate played a key role in controlling equilibrium isotopic fractionation between liquid and vapor for Cu, Zn, Ga, and Sn.

That said, we allow for the possibility that these elements were in either silicate or metal, and we calculate equilibrium isotopic fractionation for both. We do not consider the role of sulfide, as it would probably have a subordinate role, except possibly for Cu during differentiation of the Moon (Xia et al. 2019).

Wang et al. (2019) estimated equilibrium isotopic fractionation factors between vapor and liquid in the aftermath of the giant impact for K, Zn, and Sn. For K they used an ionic bond model (Young et al. 2015), an approach that was also used by Tang & Young (2020). For Zn, they calculated the isotopic fractionation between vapor and solid zincite (ZnO₂) based on previous ab initio calculations (Ducher et al. 2016). For Sn, they used the fractionation factor for metallic Sn based on a previous Mössbauer study (Polyakov et al. 2005). In the present contribution, we examine the effect of adopting different model species for K, Zn, and Sn. More importantly, we extend the study to Rb, Cu, and Ga. We rely on the literature and new first-principles density functional theory (DFT) calculations to compute equilibrium fractionation factors between vapor and silicate. Liquids are notoriously difficult to model using ab initio techniques, so instead we calculate isotopic fractionation factors (β-factors or reduced partition function ratios) of model minerals, bearing in mind that little is known about the structure of melts under conditions relevant to the Moon-forming impact. While spectroscopic studies provide insights into coordination environments in melts, the focus is almost always on magmas relevant to igneous petrology such as basalts. The coordination environments of minor and trace elements in peridotitic melts relevant to Moon formation are largely unknown, which is a caveat to the calculations presented below. Given these uncertainties, the DFT calculations should be seen as first-order estimates of equilibrium isotopic fractionation under Moon-forming conditions. They may not be accurate to better than ~50%, but this is sufficient for our purpose of evaluating scenarios of MVE depletion in the Moon.

The equilibrium isotopic fractionation factor of an element X between phases a and b ($\alpha_{a-b,X} = R_{a,X}/R_{b,X}$, where $R_{a,X}$ and $R_{b,X}$ are the isotopic ratios in phases a and b at equilibrium) can be calculated from the reduced partition function ratios ($\beta$-factor) of phases a and b (Bigeleisen & Mayer 1947),

$$ \alpha_{a-b,X} = \frac{\beta_{a,X}}{\beta_{b,X}}, \qquad (2) $$

or in $\Delta$ notation

$$ \begin{aligned} \Delta_{a-b,eq}^{X} (\%c) &= 1000(\alpha_{a-b,X} - 1) \\ &= 1000\left(\frac{\beta_{a,X}}{\beta_{b,X}} - 1\right) \approx 1000(\ln \beta_{a,X} - \ln \beta_{b,X}), \end{aligned} \qquad (3) $$

where $\beta_{a,X}$ is the reduced partition function ratio for element X between phase a and perfect monoatomic gas. The β-factor of solid can be calculated from vibrational frequencies using the harmonic approximation (Blanchard et al. 2017),

$$ \beta_{a,X} = \left( \prod_{i=1}^{3N} \prod_{q} \frac{v_i}{v_i'} \times \frac{e^{-\frac{h\nu_{q,i}}{kT}}}{1 - e^{-\frac{h\nu_{q,i}}{kT}}} \times \frac{1 - e^{-\frac{h\nu_{q,i}}{kT}}}{e^{-\frac{h\nu_{q,i}}{2kT}}} \right)^{\frac{1}{n \cdot N_q}}, \quad (4) $$

where N is the total number of atoms in the unit cell, $v_{q,i}$ and $v'_{q,i}$ are the frequencies of vibrational mode i for two isotopes at the
wavevector $q$, $N_q$ is the total number of $q$-vectors, $n$ is the number of isotopic sites in the unit cell, $h$ and $k$ are Planck's and Boltzmann's constants, respectively, and $T$ is temperature in K. The vibrational frequencies can also be used to calculate the phonon density of states (PDOS) $g(E)$, whose second moment yields the mean force constant $\langle F \rangle$ of the bonds of the element of interest with coordinating atoms (Dauphas et al. 2012),

$$ \langle F \rangle = \frac{M}{\hbar^2} \int_0^{+\infty} E^2 g(E) dE, \qquad (5) $$

where $E$ stands for vibrational energy, $M$ is the atomic mass, and $\hbar$ is the reduced Planck constant. A Bernoulli expansion of the reduced partition function ratio in powers of the frequencies of the substituted isotopes (Bigeleisen 1958; Dauphas et al. 2012) or an expansion of the kinetic energy in powers of the inverse of the temperature (Polyakov & Kharlashina 1995) gives the following expression for the $\beta$-factors:

$$ 1000 \ln \beta \approx A_1 x + A_2 x^2 + A_3 x^3, \qquad (6) $$

with $x = 1/T^2$. The higher-order terms $A_2x^2$ and $A_3x^3$ are only significant at low temperature. At the elevated temperatures that are relevant to (i) condensation in the solar nebula, and (ii) planetary accretion and differentiation, 1000 $\ln \beta$ scales as the inverse of the square of the temperature (Herzfeld & Teller 1938; Bigeleisen & Mayer 1947; Dauphas et al. 2012),

$$ 1000 \ln \beta(i/j) \approx 1000 \left( \frac{1}{m_j} - \frac{1}{m_i} \right) \frac{\hbar^2}{8k^2 T^2} \langle F \rangle, \qquad (7) $$

where $m_i$ and $m_j$ are the atomic masses of the two isotopes $i$ and $j$. The formula is sometimes written as $1000 \ln \beta(i/j) \approx 1000(1/m_j - 1/m_i)\hbar^2 A/(24k^2 T^2)$, the factor of 3 difference arising from the fact that the mean force constant $\langle F \rangle$ in Equation (7) is averaged along each direction, while the latter formula uses the sum $A$ of the three restoring force constants corresponding to displacement in three perpendicular directions. For the elements that we are interested in, this equation yields

$$ 1000 \ln \beta(^{41}K / ^{39}K) \approx 5498 \langle F \rangle / T^2, \qquad (8) $$

$$ 1000 \ln \beta(^{87}\text{Rb} / ^{85}\text{Rb}) \approx 1189 \langle F \rangle / T^2, \qquad (9) $$

$$ 1000 \ln \beta(^{65}\text{Cu} / ^{63}\text{Cu}) \approx 2148 \langle F \rangle / T^2, \qquad (10) $$

$$ 1000 \ln \beta(^{66}\text{Zn} / ^{64}\text{Zn}) \approx 2081 \langle F \rangle / T^2, \qquad (11) $$

$$ 1000 \ln \beta(^{71}\text{Ga} / ^{69}\text{Ga}) \approx 1796 \langle F \rangle / T^2, \qquad (12) $$

$$ 1000 \ln \beta(^{124}\text{Sn} / ^{116}\text{Sn}) \approx 2448 \langle F \rangle / T^2. \qquad (13) $$

We used DFT to compute phonon frequencies in selected minerals (Hohenberg & Kohn 1964; Kohn & Sham 1965), from which we calculated the mean force constant of the chemical bonds of the element of interest with coordinating atoms (Table 4). Our DFT calculations used the exchange-correlation functional of Perdew, Burke, and Ernzerhof; a plane-wave basis set; and atomic pseudopotentials as implemented in the Quantum Espresso package (Perdew et al. 1996; Giannozzi et al. 2009, 2017). The cutoff energies for wave functions ($\epsilon_{cut}$), charge density, and a uniform (Monkhorst-Pack) k-point mesh were determined such that the total

energies converged within 15 meV atom⁻¹ (Monkhorst & Pack 1976). Although this is outside of the scope of the present contribution, some of the minerals investigated below would be relevant to solar nebula condensation. During cooling of solar nebula gas, Lodders (2003) predicts condensation of K, Rb in feldspar, Cu in Fe alloy, Zn in forsterite and enstatite, Ga in Fe alloy and feldspar, and Sn in Fe alloy. Wood et al. (2019) predict condensation of Zn in troilite and Ga in Fe alloy.

### 3.1.1. Potassium and Rubidium

The structural position of alkali elements in silicate melts depends on their proportions relative to Al (Lacy 1963; Mysen et al. 1981; Matson et al. 1983; Mysen 1983). In the BSE (possibly representative of the protolunar disk composition before MVE loss), the atomic ratio (K+Na)/Al is equal to 0.12. With such a low ratio (<1), we would expect K and Rb to charge-compensate tetrahedral Al³⁺. We have therefore adopted feldspar as a model mineral for K and Rb calculations because Al is fourfold coordinated and DFT calculations have already been published for that mineral structure (Li et al. 2019; Zeng et al. 2019), allowing us to compare our results. Zeng et al. (2019) calculated the force constant of K in orthoclase, microcline, and anorthite, and they examined the effect of Na-K substitution in the albite-microcline solid solution. They also calculated the force constant of Rb in orthoclase. Li et al. (2019) focused on the effect of Na-K substitution in alkali feldspars. For the same compositions, the K force constants calculated by Li et al. (2019) are almost a factor of two higher than those calculated by Zeng et al. (2019). The discrepancy between these studies stems from the use of a one-valence pseudopotential by Li et al. (2019) compared to a nine-valence pseudopotential by Zeng et al. (2019). We follow the latter study and use a nine-valence pseudopotential. We focus on the effect of Al-Si disorder on equilibrium isotopic fractionation factors of K in feldspar, as this had not been carefully investigated before.

Depending on temperature, K-feldspars (and Rb-feldspars) can show various degrees of disorder in the site occupancy of Al and Si. Below 400°C, microcline is the triclinic ordered stable structure of K-feldspars; between ~500 and 900°C, more disordered monoclinic orthoclase is the stable structure; above ~900°C, highly disordered monoclinic sanidine is the stable structure. This disorder corresponds to possible site occupancy for aluminum, which occupies positions #2–3 in microcline, #1–4 (with 50% probability for each site) in orthoclase, and #1–8 (with 25% probability for each site) in sanidine (the site positions refer to Figure 10(a)).

In practice, a total of seven configurations were investigated by placing Al at different possible sites (#1#2, #2#3, #2#4, #6#7, #5#6, #3#6, #1#6). The results of these simulations were combined to calculate the $\beta$-factors of the different feldspars according to the probabilities of site occupancy. For microcline, we use #2–3 alone. For orthoclase the $\beta$-values are given by (#1#2+#2#3+#2#4)/3. For sanidine, not all possible site configurations were investigated, but we still weigh the $\beta$-values according to the probabilities of site occupancy, (#1#2+#2#3+#2#4+#5#6+#6#7)/9 + (#1#6+#3#6)/2/9. Rb-bearing feldspars were modeled by replacing K with Rb and placing Al atoms at positions #1#2 and #6#7, respectively. In addition, Rb-bearing
<table><caption>**Table 4**

Reduced Partition Function Ratios and Lattice Parameters of K, Rb, and Zn in Selected Minerals and Gas Species</caption><thead><tr><th>Species</th><th>Configuration, Weighing</th><th>Isotopic Ratio</th><th>a<br>Calc<sup>b</sup><br>(Å)</th><th>a<br>Meas<sup>b</sup><br>(Å)</th><th>b<br>Calc<sup>b</sup><br>(Å)</th><th>b<br>Meas<sup>b</sup><br>(Å)</th><th>c<br>Calc<sup>a</sup><br>(Å)</th><th>c<br>Meas<sup>b</sup><br>(Å)</th><th>⟨F⟩ (N m<sup>-1</sup>)</th><th>A<sub>1</sub></th><th>A<sub>2</sub></th><th>A<sub>3</sub></th></tr></thead><tbody><tr><td>K-feldspar</td><td>#2#3</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>8.7485</td><td>8.5714</td><td>13.1306</td><td>12.9646</td><td>7.3351</td><td>7.2217</td><td>26</td><td>1.432E-01</td><td>-3.060E-04</td><td>4.690E-06</td></tr><tr><td>K-feldspar</td><td>#1#2</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.885</td><td>7.792</td><td>7.885</td><td>7.792</td><td>7.289</td><td>7.192</td><td>25</td><td>1.391E-01</td><td>-3.250E-04</td><td>5.240E-06</td></tr><tr><td>K-feldspar</td><td>#2#3</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.732</td><td>7.792</td><td>8.036</td><td>7.792</td><td>7.328</td><td>7.192</td><td>25</td><td>1.368E-01</td><td>-3.150E-04</td><td>5.300E-06</td></tr><tr><td>K-feldspar</td><td>#2#4</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.869</td><td>7.792</td><td>7.876</td><td>7.792</td><td>7.149</td><td>7.192</td><td>25</td><td>1.399E-01</td><td>-3.230E-04</td><td>5.310E-06</td></tr><tr><td>K-feldspar</td><td>#6#7</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.946</td><td>7.782</td><td>7.975</td><td>7.782</td><td>7.234</td><td>7.194</td><td>21</td><td>1.128E-01</td><td>-3.080E-04</td><td>5.300E-06</td></tr><tr><td>K-feldspar</td><td>#5#6</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.971</td><td>7.782</td><td>7.971</td><td>7.782</td><td>7.238</td><td>7.194</td><td>22</td><td>1.219E-01</td><td>-3.130E-04</td><td>5.250E-06</td></tr><tr><td>K-feldspar</td><td>#3#6</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.843</td><td>7.782</td><td>8.016</td><td>7.782</td><td>7.265</td><td>7.194</td><td>25</td><td>1.378E-01</td><td>-3.230E-04</td><td>5.240E-06</td></tr><tr><td>K-feldspar</td><td>#1#6</td><td><sup>41</sup>K/<sup>39</sup>K</td><td>7.949</td><td>7.782</td><td>7.949</td><td>7.782</td><td>7.217</td><td>7.194</td><td>23</td><td>1.262E-01</td><td>-3.150E-04</td><td>5.320E-06</td></tr><tr><td>Microcline<sup>a</sup></td><td>#2#3</td><td><sup>41</sup>K/<sup>39</sup>K</td><td></td><td></td><td></td><td></td><td></td><td></td><td>25</td><td>1.400E-01</td><td>-3.105E-04</td><td>4.995E-06</td></tr><tr><td>Orthoclase<sup>a</sup></td><td>0.33#12+0.33#23<br>+0.33#24</td><td><sup>41</sup>K/<sup>39</sup>K</td><td></td><td></td><td></td><td></td><td></td><td></td><td>25</td><td>1.397E-01</td><td>-3.195E-04</td><td>5.182E-06</td></tr><tr><td>Sanidine<sup>a</sup></td><td>0.11(#12+#23+#24+#56<br>+#67)+0.22(#16+#36)</td><td><sup>41</sup>K/<sup>39</sup>K</td><td></td><td></td><td></td><td></td><td></td><td></td><td>24</td><td>1.295E-01</td><td>-3.164E-04</td><td>5.252E-06</td></tr><tr><td>Average K-feldspars</td><td></td><td><sup>41</sup>K/<sup>39</sup>K</td><td></td><td></td><td></td><td></td><td></td><td></td><td>25</td><td>1.364E-01</td><td>-3.155E-04</td><td>5.143E-06</td></tr><tr><td>KRb-microcline</td><td>#2#3</td><td><sup>87</sup>Rb/<sup>85</sup>Rb</td><td>8.8871</td><td></td><td>13.1331</td><td></td><td>7.3488</td><td></td><td>33</td><td>3.939E-02</td><td>-9.190E-05</td><td>2.200E-06</td></tr><tr><td>Rb-microcline</td><td>#2#3</td><td><sup>87</sup>Rb/<sup>85</sup>Rb</td><td>9.0037</td><td></td><td>13.1244</td><td></td><td>7.3638</td><td></td><td>28</td><td>3.282E-02</td><td>-4.720E-05</td><td>5.740E-07</td></tr><tr><td>Rb-feldspar</td><td>#1#2</td><td><sup>87</sup>Rb/<sup>85</sup>Rb</td><td>7.9489</td><td></td><td>7.9489</td><td></td><td>7.14932</td><td></td><td>22</td><td>2.621E-02</td><td>-5.680E-05</td><td>1.010E-06</td></tr><tr><td>Rb-feldspar</td><td>#6#7</td><td><sup>87</sup>Rb/<sup>85</sup>Rb</td><td>8.0177</td><td></td><td>8.04762</td><td></td><td>7.25355</td><td></td><td>22</td><td>2.602E-02</td><td>-5.590E-05</td><td>1.030E-06</td></tr><tr><td>Average Rb-<br>feldspars</td><td></td><td><sup>87</sup>Rb/<sup>85</sup>Rb</td><td></td><td></td><td></td><td></td><td></td><td></td><td>26</td><td>3.111E-02</td><td>-6.295E-05</td><td>1.204E-06</td></tr><tr><td>Zn<sub>2</sub>SiO<sub>4</sub>, M1 site</td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td>4.8437</td><td></td><td>10.4594</td><td></td><td>6.0874</td><td></td><td>125</td><td>2.596E-01</td><td>-3.371E-05</td><td>6.061E-09</td></tr><tr><td>Zn<sub>2</sub>SiO<sub>4</sub>, M2 site</td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td>4.8437</td><td></td><td>10.4594</td><td></td><td>6.0874</td><td></td><td>107</td><td>2.224E-01</td><td>-2.499E-05</td><td>3.693E-09</td></tr><tr><td>(Zn<sub>0.06</sub>Mg<sub>1.94</sub>)SiO<sub>4</sub>, Zn in M<sub>1</sub> site</td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td>4.8024</td><td>4.753</td><td>10.3205</td><td>10.19</td><td>6.0452</td><td>5.978</td><td>133</td><td>2.763E-01</td><td>-3.818E-05</td><td>6.944E-09</td></tr><tr><td>(Zn<sub>0.06</sub>Mg<sub>1.94</sub>)SiO<sub>4</sub>, Zn in M<sub>2</sub> site</td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td>4.8024</td><td>4.753</td><td>10.3205</td><td>10.19</td><td>6.0452</td><td>5.978</td><td>110</td><td>2.295E-01</td><td>-2.633E-05</td><td>4.185E-09</td></tr><tr><td>Zn dissolved in<br>forsterite<sup>a</sup><br>0.5M<sub>1</sub>+0.5M<sub>2</sub></td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td></td><td></td><td></td><td></td><td></td><td></td><td>122</td><td>2.529E-01</td><td>-3.226E-05</td><td>5.564E-09</td></tr><tr><td>hcp-Zn<sup>d</sup></td><td></td><td><sup>66</sup>Zn/<sup>64</sup>Zn</td><td></td><td></td><td></td><td></td><td></td><td></td><td>62</td><td>1.290E-01</td><td></td><td></td></tr><tr><td>Cu-Albite, (Cu<sub>0.25</sub>Na<sub>0.75</sub>)AlSi<sub>3</sub>O<sub>8</sub></td><td></td><td><sup>65</sup>Cu/<sup>63</sup>Cu</td><td>8.292</td><td>8.142</td><td>12.886</td><td>12.785</td><td>7.188</td><td>7.159</td><td>54</td><td>1.160E-01</td><td>-3.46E-04</td><td>3.41E-06</td></tr><tr><td>Cu-Albite, (Cu<sub>0.0625</sub>Na<sub>0.9375</sub>)AlSi<sub>3</sub>O<sub>8</sub></td><td></td><td><sup>65</sup>Cu/<sup>63</sup>Cu</td><td>8.281</td><td>8.142</td><td>12.93</td><td>12.785</td><td>7.233</td><td>7.159</td><td>41</td><td>8.807E-02<sup>c</sup></td><td></td><td></td></tr><tr><td>Ga-Albite, Na(Ga<sub>0.25</sub>Al<sub>0.75</sub>)Si<sub>3</sub>O<sub>8</sub></td><td></td><td><sup>71</sup>Ga/<sup>69</sup>Ga</td><td>8.287</td><td>8.142</td><td>12.969</td><td>12.785</td><td>7.267</td><td>7.159</td><td>368</td><td>6.606E-01</td><td>-5.993E-03</td><td>7.778E-05</td></tr><tr><td>α-Ga<sup>d</sup></td><td></td><td><sup>71</sup>Ga/<sup>69</sup>Ga</td><td></td><td></td><td></td><td></td><td></td><td></td><td>90</td><td>1.616E-01</td><td></td><td></td></tr><tr><td>Sn-Albite, (Sn<sub>0.25</sub>Na<sub>0.75</sub>)Al(Al<sub>0.083</sub>Si<sub>0.917</sub>)3O<sub>8</sub></td><td></td><td><sup>124</sup>Sn/<sup>116</sup>Sn</td><td>8.369</td><td>8.142</td><td>13.009</td><td>12.785</td><td>7.255</td><td>7.159</td><td>80</td><td>1.957E-01</td><td>-9.665E-04</td><td>1.471E-05</td></tr><tr><td>KO gas<sup>e</sup></td><td></td><td><sup>41</sup>K/<sup>39</sup>K</td><td></td><td></td><td></td><td></td><td></td><td></td><td>27</td><td>1.483E-01</td><td>-2.065E-04</td><td>4.109E-07</td></tr></tbody></table>

Notes. 1000 ln $\beta = A_1x + A_2x^2 + A^3x^3$, with $x = 10^6/T^2$, where T is in K. The force constants are directionally averaged, meaning that they should be multiplied by a factor of 3 to calculate the sum of the three restoring force constants corresponding to displacement in three perpendicular directions.

a The values for K-feldspars were weighted according to the probabilities of site occupancy of Al and Si in tetrahedral sites (Figure (10)). Note that for sanidine not all possible configurations were investigated. The value for Zn dissolved in forsterite took into account the fact that at high temperature Zn is partitioned with equal probabilities in M1 and M2 sites.

b “Calc” and “Meas” stand for calculated and measured, respectively.

c For this large supercell containing 208 atoms, only the Cu force constant was calculated in order to assess the dilution effect. The full frequency calculation has not been performed.

d Calculated from previously published PDOS.

e Calculated from the measured wavenumber and Equation (14).
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Model structures of the minerals investigated. (a) Unit cell of K- and Rb-feldspars (K,Rb)AlSi₃O₈. Aluminum (here in pink) can occupy two tetrahedral sites out of eight (the remaining six are occupied by Si, here in blue). In microcline Al occupies positions #2–3. In orthoclase, Al may occupy #1–4 with 50% probability for each site. In sanidine, Al can occupy #1–8 with 25% probability for each site. K is in green, and O is in red. (b) Unit cell model Zn₂SiO₄ olivine structure (Zn is in gray, Si in blue, and O in red). (c) Supercell of model Mg₂SiO₄ (olivine “forsterite” structure, 2 × 1 × 2) with 1 Zn atom in M1 site (Zn in gray, Mg in orange, Si in blue, and O in red). The proportions of Zn and Mg in the supercell are 1/32 and 31/32. (d) Same as panel (c), but with Zn in the M2 rather than M1 site. It is estimated that at 300, 1000, and 2000 K, ~90%, ~66%, and ~58% of Zn is in M1, respectively (the rest is in M2; Umemoto et al. 2011). (e) Unit cells of Cu-feldspar (Cu,Na)AlSi₃O₈ (Cu in green, Na in yellow, Si and Al in dark and light blue, and O in red). Cu⁺ substitutes for Na⁺ and forms an O-Cu-O linear configuration. (f) Unit cell of Ga-feldspar Na(Ga,Al)Si₃O₈, where Ga³⁺ substitutes for Al³⁺ (Ga in green, Na in yellow, Si and Al in dark and light blue, and O in red). (g) Unit cell of Sn-feldspar (Sn,Na)Al(Al,Si)₃O₈, where we modeled a coupled substitution of (Sn²⁺, Al³⁺) for (Na⁺, Si⁴⁺) (Sn in gray, Na in yellow, Si and Al in dark and light blue, and O in red).</figcaption></figure>

microcline was modeled by replacing one (KRb-microcline) or two K (Rb-microcline) with Rb.

The pseudopotentials used come from the ONCV library (Schlipf & Gygi 2015) with a $ϵ$ cut set at 85 Ry. Other pseudopotentials were tested. For Al, Si, and O, those from Méheut & Schauble (2014) were used in combination with those from the PSlibrary (Dal Corso 2014) for K and Rb, with a $ϵ$ cut set at 80 Ry. Both sets of pseudopotential lead to a difference of only 1 N m⁻¹ (i.e., 25 vs. 26 N m⁻¹) for the K force constant in microcline. The electronic structure integration is performed by sampling the first Brillouin zone with a uniform 2 × 2 × 2 k-point grid. Increasing the size of the q-point grid from 2 × 2 × 2 to 3 × 3 × 3 does not change the ln β value by more than 0.1‰.

We find little variation in β-factor from one Al-Si configuration to another. The total range in $\langle F \rangle$ values for the various configurations is only 5 N m⁻¹. As shown in this study, the effect of Al-Si ordering is negligible, which is a new result. Our best estimate for the value of $\langle F \rangle$ for K-feldspar is 25 N m⁻¹, which agrees with the results reported by Zeng et al. (2019). The force constants of K calculated previously for muscovite, orthoclase, sylvite, phlogopite, illite, microcline, albite, and anorthite (with similar pseudopotentials and parameters) range between 22 and 69 N m⁻¹ (Zeng et al. 2019). The stiffest bonds are found for K at low concentration in albite and anorthite. Given the high atomic Ca/(Na+K) ratio

of the BSE of ~5, we take the force constant 52 N m⁻¹ of K in microcline-albite at 12.5% K dilution (Zeng et al. 2019) as representative of K in silicate liquid. To account for uncertainties in melt structure and discrepancies between existing studies (Li et al. 2019; Zeng et al. 2019), we conservatively ascribe an uncertainty of 50% to this value (52 ± 26 N m⁻¹). For comparison, the K force constants for K-feldspar calculated by Wang et al. (2019) and Tang & Young (2020) using an ionic bond model are 72 and 28 N m⁻¹, respectively. The average force constant of Rb in feldspar calculated here is 26 N m⁻¹, consistent with previous observations that K and Rb display very similar force constants (Zeng et al. 2019). We therefore adopt the same force constant for Rb as that used for K in silicate liquid (52 ± 26 N m⁻¹).

The dominant vapor species for potassium and rubidium in the protolunar disk (Visscher & Fegley 2013) and solar nebula (Fegley 1993) are atomic K and Rb. The value of 1000 ln β of perfect monoatomic K and Rb vapor is 0. We have also calculated the value of 1000 ln β of KO and K₂O, two minor gaseous species of K that thermodynamic calculations predict should be present at a small level in vapor in equilibrium with hot silicate moon composition (Visscher & Fegley 2013). These gaseous species were modeled by placing KO and K₂O molecules in a cubic box with length of 20 bohr, and the boundary conditions were applied using the Martyna–Tuckerman method (Martyna & Tuckerman 1999). The computed
values of $\langle F \rangle$ for KO and K₂O gas are 69 and 72 N m⁻¹, respectively. We can also calculate the force constant of KO using the measured wavenumber $\tilde{\nu}$ of KO gas of 300–400 cm⁻¹ (Hirota 1995), as well as the relationship that relates force constant (A in N m⁻¹) to wavenumber for a diatomic molecule, $A = \mu(2\pi c\tilde{\nu})^2$, where $\mu = m_1m_2/(m_1 + m_2)$ is the reduced mass and c is the speed of light. For $\tilde{\nu} = 350 \pm 50$ cm⁻¹ (Hirota 1995), we thus calculate $A = 82 \pm 23$ N m⁻¹, or directionally averaged $\langle F \rangle = 27 \pm 8$ N m⁻¹. Our DFT calculation gives a significantly higher wavenumber $\tilde{\nu}$ of ~556 cm⁻¹, which explains why our $\langle F \rangle$ value is also much higher (72 N m⁻¹). For comparison, spectroscopy gives vibrational frequencies for ⁸⁵RbO and ⁸⁷RbO of 387.22 and 386.52 cm⁻¹, respectively (Yamada & Hirota 1999). This corresponds to a force constant for RbO of $A = 119$ or $\langle F \rangle = 40$ N m⁻¹, which is close to the value estimated from spectroscopy for KO but is again lower than our DFT estimate. Most compounds seem to have similar force constants for K and Rb (Zeng et al. 2019), suggesting that our DFT estimate for gas molecules may be inaccurate. We therefore use the experimentally determined wavenumber for KO of 350 cm⁻¹ to calculate its β-factor using the following formula (see the Appendix for its derivation; Table 4):

$$1000 \ln \beta_{41/39} \approx 1000 \left( \frac{1}{m_{39\text{K}}} - \frac{1}{m_{41\text{K}}} \right) \mu \times \left[ \frac{1}{24} \left( \frac{hc\tilde{\nu}}{kT} \right)^2 - \frac{1}{1440} \left( \frac{hc\tilde{\nu}}{kT} \right)^4 + \frac{1}{60480} \left( \frac{hc\tilde{\nu}}{kT} \right)^6 \right], \quad (14)$$

where $m_{39\text{K}}$ and $m_{41\text{K}}$ are the masses of K isotopes; $h$, $k$, and $c$ are the Planck constant, Boltzmann constant, and speed of light, respectively; $\mu = 1/(1/m^{39\text{K}} + 1/m^{16\text{O}})$ is the reduced mass of KO; $\tilde{\nu}$ is the wavenumber; and $T$ is the temperature. The main speciation of K in the vapor is as monoatomic gas, so the uncertainty in the β-factor of KO is inconsequential.

Using a force constant of $52 \pm 26$ N m⁻¹ for K and Rb in silicate melt and a value of 0 for monoatomic gas, we calculate the following temperature dependence of K and Rb equilibrium isotopic fractionation between liquid and vapor: $\Delta^{41}\text{K}/^{39}\text{K}_{\text{f-v,eq}} \approx (286,000 \pm 143,000)/T^2$ and $\Delta^{87}\text{Rb}/^{85}\text{Rb}_{\text{f-v,eq}} \approx (62,000 \pm 31,000)/T^2$.

### 3.1.2. Zinc

Partitioning of Zn²⁺ among olivine, orthopyroxene, and basaltic melt is like Fe²⁺ (Le Roux et al. 2010). In basaltic glass, available data suggest that Zn is in predominantly fourfold coordination, but some sixfold coordination may be present (Kuzmin et al. 2017). Zinc coordination seems to be influenced by the alkali element content, with alkali-rich compositions favoring fourfold coordination and alkali-poor compositions favoring sixfold coordination (Galoisy 2006). The coordination of Fe²⁺ in geological materials has been studied more extensively, and recent XANES studies indicate that Fe²⁺ in basaltic melt is in fivefold average coordination with oxygen (Wilke et al. 2004; Jackson et al. 2005). X-ray diffraction data also indicate that Fe²⁺ in fayalite melt is in roughly fivefold coordination (Sanloup et al. 2013). Ab initio molecular dynamics simulations find that Fe²⁺ in ultramafic melt is in a low coordination of ~3–4 (Solomatova & Caracas 2019; Ghosh & Karki 2020), but this value may not be reliable, as those simulations predict a higher coordination

number for Mg than for Fe, which is opposite to expectations (Solomatova & Caracas 2019). To summarize, considerable uncertainties remain on the exact nature of the coordination environment of Zn in mafic and utramafic magmas, with coordinations of 4 or 5 more likely, but higher values also possible in alkali-poor compositions. Ducher et al. (2016) investigated several minerals with Zn in fourfold (franklinite, gahnite, hemimorphite, zincite, sphalerite, wurtzite), mixed five- and sixfold (adamite), mixed four- and sixfold (hydro-zincite), and sixfold (smithsonite, gunningite) coordinations. Of all of these, only hemimorphite would be relevant to silicate, but even then, it is a hydrated mineral that may not be a good proxy for Zn in mafic and ultramafic magmas. To complement that work, we have therefore examined Zn-substituted olivine Zn₂SiO₄, with Zn in sixfold coordination (Figures 10(b)–(d)). Zinc in olivine is approached in two ways. We modeled the virtual Zn end-member of olivine (Zn₂SiO₄), relaxing the cell parameters and atomic positions. We also considered the substitution of one Zn for one Mg in a 2 × 1 × 2 forsterite supercell (each of the two nonequivalent cation sites, M1 and M2, were considered). In that case, atomic positions were relaxed, while the cell parameters were kept fixed to the optimized values for pure forsterite.

The computational parameters are the same as in Ducher et al. (2016). Ultrasoft pseudopotentials are from the GBRV library (Garrity et al. 2014). The wave function and charge-density cutoffs were set at 80 and 720 Ry, respectively. The electronic structure integration is performed by sampling the first Brillouin zone with a 2 × 2 × 2 shifted k-points grid for the olivine unit cell and one shifted k-point for the supercell. Zinc force constants (Table 4) were calculated using a shifted q-point grid of 4 × 2 × 4 for the unit cell and one shifted q-point for the supercell.

There is little difference between Zn substituting for all Mg atoms and Zn in dilution in a supercell. A larger difference is found between the force constant of Zn in the M1 and M2 sites. At high temperature, Zn is expected to partition approximately equally between those two sites, and the corresponding average force constant of Zn in olivine is 122 N m⁻¹. Wang et al. (2019) used the value for solid zincite (ZnO) of 148 N m⁻¹ as a proxy for Zn in silicate melt. This value is close to the one calculated here for olivine. Ducher et al. (2016) calculated the force constant of Zn in a variety of minerals. They examined a single silicate, hemimorphite, and obtained a force constant of 166 N m⁻¹. Other minerals had force constants that ranged between 99 N m⁻¹ (wurtzite; ZnS) and 177 N m⁻¹ (gahnite; ZnAl₂O₄). For Zn in silicate melt, we adopt a force constant of 122 ± 61 N m⁻¹, which corresponds to the force constant that we calculated for Zn-substituted olivine. We ascribe to this value an uncertainty of ±50% based on the range of values calculated previously in minerals with different coordination geometries and atoms (Ducher et al. 2016). The force constant of hemimorphite of 166 N m⁻¹, a silicate mineral that contains Zn in fourfold coordination, is well within the error of our proposed BSE melt value.

Black et al. (2011) used the previously published PDOS (Chesser & Axe 1974) of hcp-Zn to calculate a force constant of ~66 N m⁻¹ in this metal. We digitized the PDOS g(E) of hcp-Zn calculated by Chesser & Axe (1974, their Figure 8) and Li & Tse (2000, their Figure 2), from which we computed the force constant using the formula $\langle F \rangle = M \int_0^{+\infty} E^2 g(E) dE / \hbar^2$ (Dauphas et al. 2012; Hu et al. 2013), where M is the atomic
mass of Zn, *E* stands for energy, and *ħ* is the reduced Planck constant. For hcp-Zn, we calculate force constant values of 66 N m⁻¹ from the PDOS of Chesser & Axe (1974) (in agreement with the value of Black et al. 2011) and 62 N m⁻¹ from the PDOS of Li & Tse (2000). We use the latter as a proxy for Zn⁰ dissolved in metal in the protolunar disk or moonlet/Moon cores (we ascribe an uncertainty of ±31 N m⁻¹ to this value). Zinc in the vapor is expected to be primarily in monoatomic form (Visscher & Fegley 2013), meaning that its force constant should be zero. Zinc equilibrium isotopic fractionation between liquid and vapor is therefore Δ⁶⁶Zn/⁶⁴Znᵢ silicate-v,eq ≃ (254,000 ± 127,000)/T² for Zn in silicate and Δ⁶⁶Zn/⁶⁴Znᵢ metal-v,eq ≃ (129,000 ± 64,000)/T² for Zn in metal. Bridgestock et al. (2014), Mahan et al. (2017), and Xia et al. (2019) experimentally studied Zn equilibrium isotopic fractionation between metal and silicate at temperatures above 1450 K. They did not detect any systematic difference in δ⁶⁶Zn within ±0.1‰. Our force constant estimates give Δ⁶⁶Zn/⁶⁴Znᵢ met-sil,eq ≃ -125,000/T², corresponding to a difference of -0.06‰ at 1450 K, which would not have been detected in these experimental studies.

### 3.1.3. Copper

At high oxygen fugacity (above ~ΔFMQ+1), Cu can exist as both Cu⁺ and Cu²⁺ in silicate melt, but under Moon-forming conditions, we expect it to be present as Cu⁺ (Ripley & Brophy 1995; Holzheid & Lodders 2001; Liu et al. 2014). Laboratory experiments in mafic systems show that the partitioning of Cu between clinopyroxene and melt increases with increasing Na content of clinopyroxene, suggesting that Cu⁺ could substitute for Na⁺ in this mineral (Liu et al. 2014). Zajacz et al. (2013) examined the solubility behavior of Cu alloy in silicate melt and found that it correlated with peraluminosity, suggesting that Cu⁺ was charge balancing Al³⁺ tetrahedra. The speciation of Cu in basaltic glasses was studied directly by XAS (Maurizio et al. 2000; Lanzirotti et al. 2019). The XAS spectra that provide the best matches are those of aqueous Cu⁺ coordinated in linear complexes with Cl, S, or O ligands. At the low Cl and S contents of the protolunar disk, we would expect Cu to be complexed primarily with O. As a first pass on this question, we have decided to calculate the force constant of Cu⁺ substituting for K and Na in a feldspar silicate structure (Table 4, Figure 10(e)). As discussed below, Cu⁺ in this substitution forms a linear structure with O, consistent with XAS studies.

We chose the low albite structure (NaAlSi₃O₈) and investigated the dilution effect by substituting one Cu atom for one Na atom in the albite unit cell (Cu/(Cu+Na) = 1/4) and in the 2 × 1 × 2 albite supercell (Cu/(Cu+Na) = 1/16). The pseudopotentials come from the ONCV library (Schlipf & Gygi 2015). The wave functions and charge-density cutoffs were set at 85 and 510 Ry, respectively. The electronic structure integration is performed by sampling the first Brillouin zone with a 2 × 2 × 2 shifted *k*-points grid for the unit cell and one *k*-point for the supercell.

The force constant of Cu⁺ in the feldspar structure varies slightly with the copper concentration. Its value decreases from 54 to 41 N m⁻¹ when the Cu/(Cu+Na) ratio decreases from 1/4 to 1/16. In both optimized structures, Cu displays an almost linear O-Cu-O configuration (Figure 10(e)), which agrees well with XAS findings (Lanzirotti et al. 2019; Maurizio et al. 2000). The calculated Cu-O bond lengths are 2.05 and 2.10 Å

in the unit cell and supercell, respectively, while Maurizio et al. (2000) found a value of 1.87 ± 0.04 Å. We adopt a force constant of 54 ± 27 N m⁻¹ for Cu in silicate melt, with the caveat that a slightly higher force constant might be expected in natural samples given our overestimate of the bond length, but we expect the true value to be within the error of the value given here. Liu et al. (2021) calculated the β-factors of a variety of Cu-bearing minerals. For Cu⁺-bearing minerals, they reported force constants between 54 and 135 N m⁻¹, but all those minerals are either sulfides or oxides. They modeled silicate minerals liebauite and enstatite, but these contain copper as Cu²⁺, which is irrelevant to lunar formation. For native copper (fcc), they obtain a force constant of 94 N m⁻¹, which is the value that we adopt for Cu⁰ dissolved in liquid Fe-Ni, with an uncertainty of ±47 N m⁻¹.

Copper in the vapor is expected to be primarily in monoatomic form, meaning that its force constant should be zero. Copper equilibrium isotopic fractionation between liquid and vapor is therefore Δ⁶⁵Cu/⁶³Cuᵢ silicate-v,eq ≃ (116,000 ± 58,000)/T² if Cu is in silicate and Δ⁶⁵Cu/⁶³Cuᵢ metal-v,eq ≃ (202,000 ± 101,000)/T² if Cu is in metal. Experimental studies found a heavy Cu isotope enrichment in metal relative to silicate (Savage et al. 2015; Xia et al. 2019), but significant scatter was present and no correlation with temperature was observed, suggesting that the measured heavy Cu isotopic composition of metal relative to silicate may not reflect equilibrium conditions. Our force constant estimates give Δ⁶⁵Cu/⁶³Cuᵢ met-sil,eq ≃ +86,000/T², corresponding to a difference of +0.04‰ at 1550 K.

### 3.1.4. Gallium

Gallium in silicate melts exists as Ga³⁺. Several lines of evidence point to its substitution for Al³⁺ in silicates (Henderson et al. 1985; Baker 1995). As a melt proxy, we therefore calculated reduced partition function ratios for albite with one substitution of Ga for Al in the low albite unit cell, corresponding to a composition Na(Ga₀.₂₅Al₀.₇₅)Si₃O₈ (Figure 10(f)). We did not investigate the effect of Si-Ga disorder. We used the same pseudopotentials and computational parameters as those used to calculate the reduced partition function ratio of K in microcline and Cu in albite. The lattice parameters and force constant are compiled in Table 4. We calculate a force constant for Ga in albite of 368 N m⁻¹. We are not aware of any study that examined Ga equilibrium isotopic fractionation factors. Conservatively, we ascribe a 50% uncertainty to this value (368 ± 184 N m⁻¹). To our knowledge, the β-factor of metallic Ga has not been calculated before. We digitized the PDOS *g*(*E*) of α-Ga calculated by Remsing et al. (2018, their Figure 2), from which we computed a force constant (Dauphas et al. 2012; Hu et al. 2013) of 90 N m⁻¹ for α-Ga, which we use as a proxy for Ga in metal in the protolunar disk or moonlet/Moon cores (we ascribe an uncertainty of ±45 N m⁻¹ to this value). α-Ga is a metallic molecular crystal with strong Ga dimer covalent bonds and weaker intermolecular forces (Gong et al. 1991; Remsing et al. 2018), so it is uncertain to what extent it is representative of Ga dissolved in Fe-Ni alloy. We did not calculate the vapor speciation of Ga in the protolunar disk, but based on a comparison of the free energies of oxidation of Ga with those of other elements for which we did a speciation calculation, Ga is likely to be present in the vapor in a monoatomic form, and its force constant should be zero. Equilibrium isotopic
fractionation between liquid and vapor is therefore $\Delta^{71}\text{Ga}/^{69}\text{Ga}_{\text{silicate-v,eq}} \simeq (661,000 \pm 330,000)/T^2$ if Ga is in silicate and $\Delta^{71}\text{Ga}/^{69}\text{Ga}_{\text{metal-v,eq}} \simeq (162,000 \pm 81,000)/T^2$ if Ga is in metal.

### 3.1.5. Tin

In relatively reducing conditions, Sn is predominantly present as $Sn^{2+}$ in silicate liquids, while $Sn^{4+}$ is stabilized under oxidizing conditions (above ~FMQ; see Figure 11 of Dauphas et al. 2018). In the protolunar disk, we expect tin to exist as $Sn^{2+}$ in silicate or to be dissolved in Fe-Ni alloy. XANES spectroscopy shows that in silicate melts, $Sn^{2+}$ behaves as a network modifier like $Na^{+}$, compensating the charge of $Al^{3+}$ tetrahedra (Farges et al. 2006). For silicate, we have therefore decided to use an albite model structure by substituting $Sn^{2+}$ for $Na^{+}$ and $Al^{3+}$ for $Si^{4+}$ in a neighboring site (Figure 10(g)). The model composition is $(Sn_{0.25}Na_{0.75})Al(Al_{0.083}Si_{0.917})_3O_8$. The pseudopentials and parameters are the same as those used to calculate the force constant of K in microcline and Ga and Cu in low albite. The lattice parameters and force constant are compiled in Table 4. We calculate a force constant for $Sn^{2+}$ substituting for $Na^{+}$ in albite of $80 \text{ N m}^{-1}$. Tin possesses a Mössbauer isotope ($^{119}\text{Sn}$), and its force constant can be measured using nuclear resonant inelastic X-ray scattering spectroscopy (NRIXS; Dauphas et al. 2018). Roskosz et al. (2020) measured the force constant of tin in various synthetic glasses at different oxygen fugacities. The force constants for $Sn^{2+}$ in rhyolite, basalt, and anorthite composition glasses are 250, 254, and $204 \text{ N m}^{-1}$, respectively. These values are significantly higher than the value calculated here. As for Cu substituted for Na in albite, the calculated bond lengths may be overestimated in the present model, and the absolute value of the force constant could be underestimated. However, this uncertainty cannot explain the observed discrepancy. In the compositions studied by Roskosz et al. (2020), Sn is more diluted than in our model. However, as shown for Zn-olivine and Cu-feldspar, the effect of dilution on the force constant is small. Therefore, it is likely that our simple substitution model in albite does not reflect the main structural environment of Sn in liquid silicate. We therefore adopt the force constant measured in basalt by Roskosz et al. (2020) and ascribe a 50% uncertainty to this value ($254 \pm 127 \text{ N m}^{-1}$). Based on a thermodynamical calculation, Wang et al. (2019) argued that at the oxygen fugacity set by a liquid of bulk silicate Moon composition, Sn would be present predominantly as metal at the high temperature relevant to the protolunar disk. Their calculation, however, only considers pure compounds ($Sn^0$, $SnO$, $SnO_2$) and does not consider the activities of tin dissolved in metal and silicate, so it may not be accurate. For metal, we use the force constant of $155 \text{ N m}^{-1}$ reported by Roskosz et al. (2020) based on NRIXS measurements (we ascribe an uncertainty of $\pm 77 \text{ N m}^{-1}$ to this value). It is much higher than the force constant of $81 \text{ N m}^{-1}$ adopted by Wang et al. (2019), but that value was derived from conventional Mössbauer and is presumably less reliable. As discussed by Wang et al. (2019), the dominant gas species in protolunar disk conditions is likely to be $SnO$. They estimated the directionally averaged force constant of $SnO$ bonds based on its measured vibration frequency and obtained a value of $179 \text{ N m}^{-1}$ (the uncertainty on the wavenumber measurement is ~1%, translating to a $\pm 4 \text{ N m}^{-1}$ uncertainty on the force constant).

We calculate the following temperature dependence of Sn equilibrium isotopic fractionation between liquid and vapor: $\Delta^{124}\text{Sn}/^{116}\text{Sn}_{\text{silicate-v,eq}} \simeq 2448(254 - 179)/T^2 = (184,000 \pm 381,000)/T^2$ for Sn in silicate and $\Delta^{124}\text{Sn}/^{116}\text{Sn}_{\text{metal-v,eq}} \simeq 2448(155-179)/T^2 = -(59,000 \pm 290,000)/T^2$ for Sn in metal. The predicted metal-silicate equilibrium is $\Delta^{124}\text{Sn}/^{116}\text{Sn}_{\text{met-sil,eq}} \simeq -243,000/T^2$. Wang et al. (2019) found light Sn isotopic compositions in lunar samples, which cannot be explained by evaporation in a kinetic regime. Using a force constant of $81 \text{ N m}^{-1}$ for Sn in the liquid, they made the case that at equilibrium the vapor could have heavy Sn isotopic composition relative to the liquid. Our calculation and the results of Roskosz et al. (2020) show that this may not be true, as Sn in silicate or metal has a force constant similar to, or even higher than, that of gas $SnO$.

### 3.1.6. Temperature Dependence of Isotopic Fractionation Factors

Calculations of the temperature at which MVEs were fractionated in the aftermath of the Moon-forming impact give values of ~3500 K (Canup et al. 2015; Lock et al. 2018). The force constants measured at room temperature by NRIXS, or calculated at 0 K using DFT, can be extrapolated to high temperature only if (i) the harmonic approximation holds and (ii) the melt structure does not change. Anharmonicity would probably be associated with a weakening of the bonds at high temperature. Close to the critical point, the melt is expected to adopt a structure that resembles vapor. In both cases, our calculated equilibrium isotopic fractionation factors would represent conservative upper limits.

## 3.2. Kinetic Isotopic Fractionation Factors

The Hertz-Knudsen equation gives the degree of kinetic isotopic fractionation associated with pure evaporation and condensation (Richter et al. 2007),

$$ \alpha_{\text{Kin}}^{i/j} = \frac{\gamma_i}{\gamma_j} \sqrt{\frac{m_j}{m_i}}, \qquad (15) $$

where i and j are two isotopes of the element of interest, $\gamma_i$ and $\gamma_j$ are the evaporation or condensation coefficients of the two isotopes considered, and $m_i$ and $m_j$ are the masses of the isotopes (or isotopologues for gas molecules). An alternate formula sometimes used to account for departures from ideality in evaporation or condensation experiments is (Richter et al. 2011)

$$ \alpha_{\text{Kin}}^{i/j} = \left( \frac{m_j}{m_i} \right)^b, \qquad (16) $$

or in $\Delta$ notation,

$$ \Delta_{\text{Kin}}^{i/j} = 1000(\alpha_{\text{Kin}}^{i/j} - 1) = 1000 \left[ \left( \frac{m_j}{m_i} \right)^b - 1 \right], \qquad (17) $$

where b is a constant. No experiment has been performed to document the value of b during condensation from highly supersaturated vapor, but unless the condensation coefficients depend strongly on isotopic mass, one would expect $b \simeq 0.5$. Experiments of evaporation from silicate melts give $b \simeq 0.5$ for both K (Yu et al. 2003; Richter et al. 2011; Zhang et al. 2021) and Rb (Zhang et al. 2021). We are unaware of vacuum
evaporation experiments for Cu, Zn, Ga, and Sn from silicate melt. Previous experiments performed for other elements in different compositions yielded b exponents that never departed much from 0.5 (Richter et al. 2009), so this is the value that we adopt for all elements. The elements K, Rb, Cu, Zn, and Ga are expected to be in monoatomic form in the vapor, so the masses used in the equation above are those of the isotopes. Tin likely exists as SnO in the vapor, so we use the masses of the ¹⁶Sn¹⁶O and ¹²⁴Sn¹⁶O isotopologues. The kinetic isotopic fractionation factors Δᵢⱼₖₘ for the ⁴¹K/³⁹K, ⁸⁷Rb/⁸⁵Rb, ⁶⁵Cu/⁶³Cu, ⁶⁶Zn/⁶⁴Zn, ⁷¹Ga/⁶⁹Ga, and ¹²⁴Sn¹⁶O/¹¹⁶Sn¹⁶O are thus -24.7, -11.6, -15.5, -15.3, -14.2, and -29.0‰, respectively (Table 2).

### 3.3. Combining Equilibrium and Kinetic Effects

At equilibrium, the flux of vapor atoms impinging the surface of the condensate is exactly equal to the flux of atoms evaporated, and the isotopic fractionation between the liquid and the vapor is the equilibrium value. If the vapor pressure of an element ($P_{Vap}$) is higher than the equilibrium vapor pressure ($P_{Vap,Eq}$), then the flux of atoms impinging the surface of the liquid will exceed the evaporative flux, which corresponds to net condensation. Conversely, if the vapor pressure is lower than the equilibrium vapor pressure, the condensate will experience net evaporation. For a net evaporative flux, we define the liquid-vapor isotopic fractionation for the isotopic ratio i/j as

$$ \alpha_{Evap}^{i/j} = \frac{N_{i,l}/N_{j,l}}{J_{i,l \to v}/J_{j,l \to v}}, \quad (18) $$

where $N_{i,l}/N_{j,l}$ is the isotopic ratio i/j in the liquid and $J_{i,l \to v}/J_{j,l \to v}$ is the isotopic ratio in the escaping vapor. For a net condensing flux, we define the liquid-vapor isotopic fractionation as

$$ \alpha_{Cond}^{i/j} = \frac{J_{i,v \to l}/J_{j,v \to l}}{N_{i,v}/N_{j,v}}, \quad (19) $$

where $N_{i,v}/N_{j,v}$ is the isotopic ratio i/j in the vapor and $J_{i,v \to l}/J_{j,v \to l}$ is the isotopic ratio of the vapor atoms that are condensing at the surface. During both evaporation and condensation, kinetic isotopic fractionation can be present if the vapor is supersaturated (condensation) or undersaturated (evaporation) (Richter et al. 2009).

During condensation, the kinetic theory of gases predicts that the flux of atoms impinging the surface of the liquid should be enriched in the light isotopes relative to the vapor. During evaporation, the same theory predicts that the flux of atoms evaporated from the surface of the liquid should be enriched in the light isotopes relative to the liquid. In general, the degree of isotopic fractionation during either evaporation or condensation is an interplay between equilibrium and kinetic processes. The factor that controls the overall fractionation is the degree of supersaturation (for condensation) or undersaturation (for evaporation),

$$ S = P_{Vap}/P_{Vap,Eq}. \quad (20) $$

Following the pioneering work of Jouzel & Merlivat (1984), who investigated H and O isotopic fractionation during snow formation, several studies have established the theoretical

foundations for describing stable isotopic fractionation during evaporation and condensation in planetary materials. Using the Hertz-Knudsen equation, Richter et al. (2002) derived the equations that give the isotopic fractionation for evaporation for any value of S when there is no equilibrium fractionation. Richter et al. (2007) gave the equation for evaporation when S = 0 and there is equilibrium fractionation. Richter (2004) gave the equation for condensation for any value of S when there is no equilibrium fractionation. Dauphas & Rouxel (2006) gave the equations for evaporation and condensation for any value of S when there is no equilibrium fractionation. Simon & DePaolo (2010) adapted the formula of Jouzel & Merlivat (1984) to describe isotopic fractionation during condensation when both equilibrium and kinetic effects are present. Dauphas et al. (2015) gave the formulae for both evaporation and condensation when equilibrium isotopic fractionation is present and for any value of S. During evaporation, the isotopic fractionation between liquid and escaping vapor ($\Delta_{l \to l \to v}^{i/j} = \delta_{Liquid}^{i/j} - \delta_{Vap Flux}^{i/j}$; l stands for liquid, and $l \to v$ stands for evaporation flux from liquid to vapor) for the ratio i/j is given by Dauphas et al. (2015) (the isotopic composition of the vapor is assumed to be set by the net evaporative flux from the liquid to the vapor),

$$ \Delta_{l \to l \to v}^{i/j} = \Delta_{Eq}^{i/j} - (1 - S)\Delta_{Kin}^{i/j}. \quad (21) $$

where $\Delta_{Eq}^{i/j} = (\alpha_{Eq, l \to v}^{i/j} - 1) \times 1000$ is the equilibrium fractionation between liquid and vapor, and $\Delta_{Kin}^{i/j} = (\alpha_{Kin}^{i/j} - 1) \times 1000$ is defined in Equation (17). During condensation, the isotopic fractionation between vapor and condensing liquid ($\Delta_{v \to l \to v}^{i/j} = \delta_{Cond Flux}^{i/j} - \delta_{Vap}^{i/j}$; v stands for l stands for condensation flux from vapor to liquid, and v stands for vapor) takes the form (Dauphas et al. 2015; the isotopic composition of the liquid is assumed to be set by the net condensation flux from the vapor to the liquid)

$$ \Delta_{v \to l \to v}^{i/j} = \Delta_{Eq}^{i/j}/S + (1 - 1/S)\Delta_{Kin}^{i/j}. \quad (22) $$

Note that this formula (Dauphas et al. 2015) is mathematically equivalent to that given by Simon et al. (2017).

## 4. Discussion

In Section 2, we reassessed the degrees of depletion of K and Rb. We find that K and Rb are depleted by factors of 3.96 ± 0.87 and 4.71 ± 1.07 in the Moon relative to Earth. If we take the estimates from O'Neill (1991) for the Moon and Earth, the depletion factors for Cu, Zn, Ga, and Sn are 10, 26, 6, and 30, respectively. If we take the estimates of O'Neill (1991) for the Moon and McDonough & Sun (1995) for Earth, the depletion factors for Cu, Zn, Ga, and Sn are 21, 42, 9, and 32, respectively. Albarède et al. (2015) concluded, based on Zn/Fe analyses, that the Moon was depleted in Zn by a factor of ~200. Using the most recent estimates and assuming that the BSE is representative of the material that made the Moon, we estimate that 74.8% ± 5.5% of K, 78.8% ± 4.8% of Rb, 95.3% of Cu, 99.5% of Zn, 89.4% of Ga, and 96.9% of Sn must have been lost in the aftermath of the Moon-forming impact.

Earth and the Moon are not the most volatile depleted objects in the solar system. Several groups of differentiated meteorites show larger depletions, notably IVB iron meteorites,
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Main models of depletion of MVEs in the Moon. (a) Incomplete condensation model of Canup et al. (2015). In that model, a proto-Moon representing ~40% of the final lunar mass is rapidly formed beyond the Roche limit. The remaining 60% is accreted over 100–1000 yr by viscous spreading of MVE-depleted material from inside the Roche limit. (b) Incomplete condensation synestia model of Lock et al. (2018). In this scenario the Moon is born in a synestia structure produced by a high-energy impact. MVE depletion is due to incomplete condensation in the synestia envelope. (c) MRI-driven vapor drainage model of Charnoz & Michaut (2015) and Nie & Dauphas (2019). In this model, MRI in the protolunar disk (also see Carballido et al. 2016; Gammie et al. 2016; Mullen & Gammie 2020) drives a large fraction of the MVE-rich vapor layer to flow inward to be accreted by Earth. MVEs lost onto Earth are constantly replenished by evaporation from the liquid layer. (d) Tidally driven hydrodynamic escape of MVEs from the LMO (Sossi et al. 2019; Charnoz et al. 2021). Escape of MVEs from the LMO is a two-step process involving evaporation from silicate melt, followed by hydrodynamic escape to space facilitated by tidal pull from Earth. In their simplest forms, A and B are partial condensation models, while C and D are fractional evaporation models.</figcaption></figure>

HEDs, and angrites (e.g., Figure 1). The canonical impact model of lunar formation predicts that the fraction of impactor material in the disk in the aftermath of the giant impact should be around 0.73–0.89 (Canup 2004b). If the impactor material was depleted in MVEs to a similar level to angrites, the 11%–27% BSE contribution would have yielded a uniform three- to ninefold depletion for MVEs in the Moon relative to Earth. A volatile-element-depleted impactor can thus explain the roughly fourfold depletions in K and Rb. It cannot explain, however, the larger depletions documented for other elements. At face value the large depletions measured for Cu, Zn, and Sn rule out a scenario involving solely an impactor depleted in volatile elements.

Several scenarios have been proposed to explain the MVE depletion of the Moon (Figure 11).

(1) Canup et al. (2015) consider the canonical giant impact scenario (Figure 11(a)). In a representative model outcome, approximately 40% of the disk material is originally ejected beyond the Roche limit (at ~3 $R_⊕$), where it forms a proto-Moon in a few weeks to a year. The remaining 60% is accreted from material that was originally ejected within the Roche limit and viscously spreads outward on a timescale of 100–1000 yr. In the inner disk, the effective viscosity arises from the tidal disruption by Earth of short-lived clumps of liquid. Approximately 10% of the material ejected beyond the Roche limit is in the form of vapor, but this vapor would most likely be accreted by the Moon. In the Canup et al. (2015) scenario, the Moon inherits its depletion in MVEs from the MVE-depleted component that is transferred from within the Roche limit to the Moon-forming region. Canup et al. (2015) argue that the lunar mantle may be stratified so that only the late accreted material would be sampled at the lunar surface. To first order, this is a partial condensation scenario whereby the Moon is

depleted in volatile elements because they stay behind in the gas when the liquid condensate is transferred beyond the Roche limit to be accreted by the growing Moon. MVEs in the vapor are eventually accreted by Earth.

(2) Lock et al. (2018) consider a scenario whereby the Moon is born in a synestia structure produced by a high-energy impact that leaves the Earth–Moon system in a state characterized by a high angular momentum (Figure 11(b)). In this model, vapor plumes raising adiabatically in the synestia structure see their temperature and pressure decrease until condensates form. When an ascending parcel of vapor becomes optically thin, it radiates away its energy, leading to rapid and extensive condensation. The volatile-element-depleted condensates thus formed fall back in the structure, where they are captured by a seed moonlet that formed rapidly beyond the Roche limit. The condensates can also be partially or completely revaporized. The growing moonlet can continue exchanging/equilibrating with the surrounding gas running past it. Lock et al. (2018) modeled the depletion in volatile elements of the Moon as a partial condensation model. They could reproduce the abundance pattern of MVEs in the Moon, notably Ge, Zn, Na, Cu, and K, for pressures above 20 bar and temperatures of ~3500–4000 K.

(3) Sossi et al. (2019) and Charnoz et al. (2021) considered a model involving loss of MVE from the LMO by hydrodynamic escape facilitated by the tidal pull of Earth (Figure 11(c)). In that scenario, MVEs are lost incrementally at a temperature of ~1600–1800 K. Tang & Young (2020) made the case that the vapor generated would have the opportunity to equilibrate with the surface of the magma ocean, thereby limiting escape efficiency and isotopic fractionation. Charnoz et al. (2021) showed, however, that tidal pull from Earth would have been able to sustain a vigorous hydrodynamic escape flow that could
have driven significant undersaturation in the vapor in contact with the LMO, allowing the production of significant kinetic isotope effects needed to explain available isotopic data (Nie & Dauphas 2019). According to Charnoz et al. (2021), this hydrodynamic wind would have been sufficient to drive the loss of Na, K, and Zn in less than ~1000 yr before the LMO was capped by a solid flotation crust. This model cannot explain simultaneously the depletions of Na, K, and Zn. The reason for the discrepancy between Tang & Young (2020) and Charnoz et al. (2021) is that Tang & Young (2020) prescribed a diffusive boundary layer at the vapor/magma interface that limits escape efficiency, while Charnoz et al. (2021) describe the atmosphere as a single layer with advective transport throughout.

(4) Charnoz & Michaut (2015) and Nie & Dauphas (2019) consider a disk with a vapor layer that is abnormally viscous owing to magnetorotational instability (MRI; also see Carbalido et al. 2016; Gammie et al. 2016; Mullen & Gammie 2020; Figure 11(d)). A large fraction of the vapor layer flows inward through viscous spreading so that volatile elements are lost by accretion onto Earth. The volatile elements lost onto Earth are constantly replenished by vaporization from the liquid layer. In this scenario of fractional evaporation, vapor is constantly removed from the system, and once accreted onto Earth, it has no chance to back react with the liquid layer.

(5) Machida & Abe (2004) and Desch & Taylor (2012) made the case that volatiles could have been lost from the disk by hydrodynamic escape. Nakajima & Stevenson (2018) showed, however, that hydrogen loss would have been limited by diffusion through heavy atoms and molecules (e.g., SiO and O). Because hydrogen escape was weak, escape of MVEs would have been negligible, so this model is not discussed further below.

In essence, the models of Canup et al. (2015) and Lock et al. (2018) are partial condensation models. Depending on the degree of supersaturation of the vapor, they could involve kinetic or equilibrium processes. The models of Charnoz & Michaut (2015) and Nie & Dauphas (2019) are fractional evaporation models. Nie & Dauphas (2019) invoke a small degree of vapor undersaturation (~99%), which produces kinetic isotope effects. The model of Sossi et al. (2019) and Charnoz et al. (2021) is a fractional evaporation model, involving equilibrium (Sossi et al. 2019; Tang & Young 2020) and possibly kinetic isotopic fractionation, depending on the strength of the hydrodynamic escape (Charnoz et al. 2021). Elemental abundances provide no direct clues on whether MVE depletion in lunar rocks involved partial condensation or evaporation, and whether equilibrium or kinetic processes were at play. In contrast, isotopic compositions are extremely sensitive to the conditions of vapor-liquid exchanges and can help us test scenarios of MVE depletion in the Moon.

## 4.1. High-temperature Fractional Kinetic Evaporation from the Protolunar Disk

The Charnoz & Michaut (2015) and Nie & Dauphas (2019) models are primarily evaporation models involving continuous removal of the vapor phase by viscous accretion onto Earth. Equilibrium thermodynamics predicts that a significant fraction of vapor Na and K could have been ionized in the protolunar disk (Visscher & Fegley 2013). The disk could have been embedded in Earth's dipole field, and it is conceivable that MRI in the vapor layer would develop

(Charnoz & Michaut 2015; Carballido et al. 2016; Gammie et al. 2016; Mullen & Gammie 2020). Magnetic fields in partly ionized disks characterized by slower rotational velocity in more distant layers, like in a Keplerian disk, can produce significant turbulence and viscosity. The underlying mechanism is the coupling through magnetic field lines of parcels of ionized gas on different orbits. As a result of viscosity, most of the mass is transported inward, and the associated loss of angular momentum is balanced by the outward transport of some material. Vapor in the disk could have been accreted onto Earth by this process but would have been constantly replenished by vaporization from the liquid layer, driving MVE depletion in the residual liquid, which eventually formed the Moon.

If the vapor is constantly removed, the evaporative process can be modeled using a Rayleigh distillation, where the vapor need not be in equilibrium with the liquid (Equation (21)),

$$ \delta_l \approx \delta_0 - [\Delta_{\text{Eq}} - (1-S)\Delta_{\text{Kin}}] \ln f, \quad (23) $$

where f is the fraction of the denominator isotope (e.g., $^{39}$K) left in the liquid after evaporation. If $S \sim 1$, at each increment of vapor removal the vapor has time to equilibrate with the liquid and the fractionation is entirely equilibrium. If $S \ll 1$, the vapor is removed at a faster rate than the timescale needed for equilibration and the isotopic fractionation is almost entirely kinetic.

Wang et al. (2019) found light Sn isotopic composition in lunar rocks relative to BSE. This cannot be explained by evaporation in a kinetic regime, as this would leave the residue (the Moon) enriched in the heavy isotopes, which is opposite to observations. Wang et al. (2019) therefore invoked equilibrium isotopic fractionation in a distillation model to explain the Sn data. The temperature that they considered is 2500 K, which is low, as the temperature in the protolunar disk in the immediate aftermath of the giant impact would have been higher, and MVE element depletion patterns point to temperatures of ~3500 K (Canup et al. 2015; Lock et al. 2018). We used the known depletions of MVEs in the Moon, together with equilibrium and kinetic isotopic fractionation factors (Table 2), to calculate the isotopic fractionation of elements remaining in the liquid after fractional evaporation at 2500 K, where each increment of vapor lost is in equilibrium with the liquid (Figure 12). As shown, equilibrium isotopic fractionation ($S=1$) at 2500 K can only explain the isotopic composition of Ga in the Moon. Unlike Wang et al. (2019), we cannot explain the light isotopic composition of Sn because we use a higher force constant for Sn in melt than the value that they used (for both Sn dissolved in metal or silicate liquid). The force constants that we used were measured by NRIXS and should be more accurate than those used by Wang et al. (2019). Based on available evidence, we can therefore exclude this scenario.

Nie & Dauphas (2019) examined quantitatively the implications of a model of kinetic isotopic fractionation associated with evaporation of MVEs by MRI-powered viscous drainage onto Earth. This model can potentially explain why elements of different volatilities experienced vaporization under similar undersaturations. The degree of saturation in the vapor layer is given by Nie & Dauphas (2019) (for $S$ close to 1),

$$ S_i \approx 1 - \frac{3\pi\nu}{2\gamma_i} \sqrt{\frac{m_i/\bar{m}}{RGM_\oplus}}, \quad (24) $$
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison between measured Moon–Earth isotopic differences (y-axis) and predicted fractionations in the residual liquid for fractional equilibrium evaporation (Table 2 and Equation (23) and (35)). The temperature of 1000 K was considered by Sossi et al. (2019), but this is below the basalt solidus. The temperatures of 1250 and 1550 K are the approximate low-pressure basalt solidus (incipient melting) and liquidus (complete melting). The liquidus temperature could be more relevant to hydrodynamic escape from the LMO (Tang & Young 2020). The temperature of 2500 K was proposed by Wang et al. (2019) to explain the light Sn isotopic composition of the Moon by equilibrium evaporative loss from the protolunar disk. Metal in the LMO could have segregated into a core, so we only plot the results for isotopic fractionation between vapor and silicate at 1000, 1250, and 1550 K. In the disk at 2500 K, metal could have been present, so we show the predictions for both metal (blue filled circles) and silicate (red filled circles).</figcaption></figure>

where $v$ is the kinematic viscosity (in $m^2 s^{-1}$), $\gamma_i$ is the evaporation coefficient, $m_i$ and $\bar{m}$ are the mass of the vapor species and mean molecular mass of the bulk vapor, $R$ is the distance from Earth's center, $G$ is the gravitation constant, and $M_⊕$ is the mass of Earth. The effective kinematic viscosity in the vapor layer $v$ needed to sustain an undersaturation level of 0.99 is $10^{-2}-10^{-3}$ (Nie & Dauphas 2019), which agrees with values estimated by Carballido et al. (2016) for an MRI-active protolunar disk where alkali elements are partly ionized.

Two parameters can vary between elements in Equation (24) ($\gamma_i$ and $m_i$), which could potentially produce differences in saturation $S_i$ (Nie & Dauphas 2019). The square roots of the masses of MVEs are, however, all very similar (within a factor of ~2; for example, we have $\sqrt{39} \approx 6$ and $\sqrt{85.5} \approx 9$ for K and Rb, respectively). The only parameter in this equation that could vary significantly from one element to another is the evaporation coefficient $\gamma_i$, which is the ratio of the actual

evaporation flux to that given by the kinetic theory of gases assuming that at equilibrium all gas atoms impinging the surface are incorporated in the liquid and evaporation balances that flux. Determination of evaporation coefficients is fraught with difficulties, as it is most often based on experiments of evaporation under vacuum and measurement or calculation of equilibrium vapor pressures. Evaporation coefficients for silicate melts at ~1200–1400°C are highly uncertain, with values of 0.07–0.3 for Na, K, and Rb reported in the literature (Zhang et al. 2021, and references therein). The kinetics of evaporation of Zn, Ga, Sn, and Cu are poorly known because their thermodynamics are less well constrained and, to our knowledge, no vacuum evaporation experiment has been performed for these MVEs. The temperature that can reproduce the pattern of MVE depletion in lunar rocks is ~3500 K (Canup et al. 2015; Lock et al. 2018), with a large uncertainty attached to it owing to difficulties in extrapolating thermodynamic data well beyond their calibration range.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Comparison between measured Moon–Earth isotopic differences (y-axis) and predicted fractionations in the residual liquid for fractional evaporation in different undersaturation levels and temperatures of 3500 K (panels (a)–(c)) and 1600 K (panels (d)–(f)) (Table 2 and Equation (23)). The isotopic composition of the Moon for K can only be reproduced by evaporation with $S \simeq 0.99$ at 3500 K and $S \simeq 0.995$ at 1600 K (see text and Nie & Dauphas 2019 for details). The δ-notations refer to $^{41}$K/$^{39}$K, $^{87}$Rb/$^{85}$Rb, $^{65}$Cu/$^{63}$Cu, $^{66}$Zn/$^{64}$Zn, $^{71}$Ga/$^{69}$Ga, and $^{124}$Sn/$^{160}$Sn ratios. The 3500 K calculation would be relevant to the protolunar disk, when metal and silicate could have coexisted, so we show the predicted isotopic fractionations for the two liquids (metal: blue filled circles; silicate: red filled circles).</figcaption></figure>

Kobsch & Caracas (2020) ran first-principles molecular dynamics simulations to calculate the critical temperature of Na and K-feldspar composition and to study the speciation of these elements in liquid and vapor. The critical temperatures that they calculate are in the range 5000–5500 K. Expectedly, at the high temperatures relevant to MVE depletion, the melts become de-polymerized, and the species present are more similar to those encountered in the vapor. It is thus conceivable that, given structural changes in the melt, the evaporation coefficients could converge to values near unity (near the apex of the liquid-vapor dome, liquid and vapor are nearly indistinguishable and the evaporation/sticking coefficient should be 1). It is thus reasonable to assume that $\gamma_i \simeq 1$ and vapor saturation given by Equation (14) could be very similar for all species.

In Figures 13(a)–(c), we compare measured isotopic differences between Moon and Earth with predictions from the Nie & Dauphas (2019) model. Given structural changes in the melt at >3500 K, the force constants that we calculate are likely to be inaccurate, but this is inconsequential, as equilibrium fractionation factors decrease with the inverse of the square of the temperature and are entirely negligible at 3500 K. We plot predictions for different degrees of saturation: $S=1$ (equilibrium), 0.99, and 0.98. With $S=1$, the predicted fractionations are much smaller than the values measured, while $S=0.98$ produces larger shifts than measured. As shown, kinetic effects at an undersaturation of $S=0.99$ can reproduce very well the data for K, Rb, Cu, Zn, and Ga, a

conclusion that had been reached previously by Nie & Dauphas (2019). The isotopic composition of tin cannot be reproduced, and we have no explanation for this departure, other than assuming that it was fractionated isotopically by large-scale magmatic differentiation or core formation on Earth or the Moon.

Van Kooten et al. (2020) calculated a lower saturation of 0.9 using the same approach as Nie & Dauphas (2019). The difference stems from the fact that they used reduced masses of coordinating atoms to calculate the kinetic isotopic fractionation factor (e.g., $40 \times 16/(40+16)$ for K coordinated with O). However, using the reduced mass for the kinetic isotopic fractionation factor during evaporation violates Dalton’s law and the kinetic theory of gases. At equilibrium, the kinetic theory of gases states that the flux of atoms impinging a surface must scale as the square root of the mass of the atom/molecule in the gas phase. Because at equilibrium the flux *in* must be equal to the flux *out*, the flux of atoms leaving the surface must also scale as the square root of the mass, not the reduced mass. There could be some mass dependence on the evaporation coefficients, but this is considered in Equation (16). There is also no evidence from vacuum evaporation experiments for the involvement of reduced masses of coordinating atoms in the liquid. For example, experiments of evaporation in vacuum of basaltic melt yield a kinetic isotopic fractionation factor for K of $(38.964/40.962)^{0.48} \approx 0.97628$ (Yu et al. 2003; Richter et al. 2011; Zhang et al. 2021), which is very close to the ideal
value of (38.964/40.962)$^{0.5}$ = 0.97531. Using the reduced mass that van Kooten et al. (2020) propose would yield a kinetic fractionation factor of [38.964 × 16/(38.964 + 16)]$^{0.5}$/ [40.962 × 16/(40.962 + 16)]$^{0.5}$ = 0.99288 = (38.964/40.962)$^{0.14}$, which disagrees with experiments. Reduced mass is sometimes involved in modeling evaporation when transport is limited by diffusive transport in a vapor of different mean atomic weight than the element evaporated (Richter et al. 2011; Bourdon & Fitoussi 2020; Sossi et al. 2020). This is not the model envisioned by van Kooten et al. (2020), and it is also not relevant to the measurements of tektites and nuclear fallout samples that they invoke, as transport in the vapor is unlikely to be limited by diffusion in those settings. To summarize, there is no theoretical or experimental support for using the reduced mass of coordinating atoms to calculate kinetic isotopic fractionation factors, and isotopic fractionation of MVEs in the bulk Moon is consistent with evaporation in a medium that was ~99% saturated (Nie & Dauphas 2019).

Below, we examine more closely the timescale required for evaporation of alkali elements, as this provides a straightforward test to the model of MVE loss through viscous drainage in the protolunar disk. A caveat to this calculation is that thermodynamic data needed to calculate equilibrium vapor pressures are largely unconstrained at the temperatures considered for the Moon-forming impact, which is the reason why we focus on Na and K that have been more extensively studied than trace elements Rb, Cu, Zn, Ga, and Sn. The net flux of alkali element *i* from the liquid layer to the vapor layer per unit time integrated over the whole surface *s* of the disk is

$$J_i \times s = \frac{(1 - S_i) \gamma_i P_{eq,i}}{\sqrt{2\pi m_i RT}} s \quad (25)$$

The equilibrium vapor pressure can be calculated from the standard Gibbs free energy $\Delta G$ of the reaction $AO_{0.5,l} \leftrightarrow A_v + 0.25O_{2,v}$ (A = Na or K). We thus have for the flux (the factor of $10^5$ accounts for the fact that $P_{eq,i}$ in Equation (25) is in Pa, while the law of mass action uses a standard pressure of 1 bar in the definition of fugacity)

$$\frac{dN_i}{dt} = -10^5 \frac{(1 - S_i) \gamma_i \Gamma_i e^{-\Delta G_i / RT}}{\sqrt{2\pi m_i RT} f_{O_2}^{1/4}} X_i s \quad (26)$$

where $f_{O_2}$ is the oxygen fugacity (the $O_2$ partial pressure in bar), $\Gamma_i$ is the activity coefficient of $AO_{0.5,l}$, $X_i$ is the mole fraction of $AO_{0.5,l}$ (the main species considered are $SiO_2$, $TiO_2$, $Al_2O_3$, $Cr_2O_3$, $MnO$, $FeO$, $NiO$, $MgO$, $CaO$, $NaO_{0.5}$, and $KO_{0.5}$) in the melt layer of volume $v$, and $N_i$ is the total amount of $AO_{0.5,l}$ (A = Na or K) in the layer ($N_i = X_i N_{total}$; $N_{total} = \rho_l v / M_l$), with $M_l$ and $\rho_l$ the molar mass and density of the liquid layer, respectively). This equation can be rewritten as

$$\frac{d\ln N_i}{dt} = -10^5 \frac{(1 - S_i) \gamma_i \Gamma_i e^{-\Delta G_i / RT} M_l}{\sqrt{2\pi m_i RT} f_{O_2}^{1/4} \rho_l v / s} \quad (27)$$

which can be integrated to calculate the time to achieve a given degree of depletion,

$$t_i = -10^{-5} \frac{\sqrt{2\pi m_i RT} f_{O_2}^{1/4} \rho_l v / s}{(1 - S_i) \gamma_i \Gamma_i e^{-\Delta G_i / RT} M_l} \ln f_i \quad (28)$$

Injecting Equation (23) into this equation, we can eliminate the saturation by relating it to the measured isotopic fractionation,

$$t_i = -10^{-5} \frac{\sqrt{2\pi m_i RT} f_{O_2}^{1/4} \rho_l \Delta_{\text{Kin},i} v / s}{(\delta_{l,i} - \delta_0 + \Delta_{\text{Eq},i} \ln f_i) \gamma_i \Gamma_i e^{-\Delta G_i / RT} M_l} (\ln f_i)^2 \quad (29)$$

At high temperature when equilibrium isotopic fractionation is negligible, the saturation required for explaining the heavy isotope enrichment of MVEs is $S = 0.99$ (see text above and Nie & Dauphas 2019). For Na and K, the degrees of depletions $f_i$ in the Moon are 0.15 (Dauphas et al. 2014) and 0.252 (Table 2). The temperature relevant to depletion of MVEs in the protolunar disk is ~3500 K (Canup et al. 2015; Lock et al. 2018). The molar masses $m_i$ of atomic K and Na are 0.039 and 0.023 kg mol$^{-1}$ respectively. The volume-to-surface ratio $v/s$ of the liquid layer is its half-height $h/2$. A lunar-mass disk of density $\rho_l = 3000$ kg m$^{-3}$ spread between Earth and the Roche radius (~3 Earth radius) would have a half-height of 12,000 m (12 km). For comparison, Pahlevan & Stevenson (2007) examined a 2 Moon mass disk extending to 5 Earth radii, corresponding to a half-height (or volume-to-surface ratio) of 8 km. The molar mass of the liquid for a BSE composition is $M_l = 0.052$ kg mol$^{-1}$. The oxygen fugacity ($f_{O_2}$ in bar) is the one calculated for vapor in equilibrium with BSE liquid (Visscher & Fegley 2013; T is the temperature in K; the formula is valid between 1500 and 4000 K),

$$f_{O_2} = 10^{1.35 \times 10^7 / T^2 - 28179 / T + 6.5137}, \quad (30)$$

the impact on the evaporation timescale of not knowing the oxygen fugacity precisely is small. The unitless product of evaporation and activity coefficients $\gamma_i \Gamma_i$ can be expressed as (Zhang et al. 2021)

$$\gamma_{\text{Na}} \Gamma_{\text{Na}} = e^{-\frac{22239}{T} + 2.4979}, \quad (31)$$

$$\gamma_{\text{K}} \Gamma_{\text{K}} = e^{-\frac{29265}{T} + 3.4593}. \quad (32)$$

The product $\gamma_{\text{Na}} \Gamma_{\text{Na}}$ increases from $1.1 \times 10^{-5}$ at 1600 K to 0.02 at 3500 K. The product $\gamma_{\text{K}} \Gamma_{\text{K}}$ increases from $3.6 \times 10^{-7}$ at 1600 K to 0.007 at 3500 K. These increases primarily reflect increases in the activity coefficients $\Gamma$. The standard Gibbs free energies $\Delta G_i$ of the evaporation reaction $NaO_{0.5,l} \leftrightarrow Na_v + 0.25O_{2,v}$ and $KO_{0.5,l} \leftrightarrow K_v + 0.25O_{2,v}$ were calculated using the thermodynamic data from Knacke et al. (1991) for $NaO_{0.5,l}$, $Na_v$, $K_v$, and $O_{2,v}$ and from Lamoreaux & Hildenbrand (1984) for $KO_{0.5,l}$. We fitted the standard Gibbs free energies of the reactions using second-order polynomials,

$$\Delta G_{\text{Na}} = 0.0054 \times T^2 - 138.34 \times T + 279829, \quad (33)$$

$$\Delta G_{\text{K}} = 0.007 \times T^2 - 143.54 \times T + 246083. \quad (34)$$

In Figure 14, we plot the disk evaporation timescales for Na and K as a function of temperature. As shown, at 3500 K, evaporation of Na and K in a 99% saturated vapor would have reached Moon-like levels of depletion for Na and K in ~1 yr. The lifetime of the disk is inferred to have been on the order of ~100 yr (Canup 2004a), meaning that there would have been ample time for MVEs to be lost by volatilization from the disk and be accreted onto Earth.

Magnetic fields could have also dramatically influenced the evolution of the disk. Gammie et al. (2016) and
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Evaporation timescales for reaching Moon-like depletions for Na and K in the protolunar disk (panel (a)) and the LMO (panel (b)) at 99% saturation (Equation (28); Figure 13). The lifetimes of the protolunar disk stage (~100 yr; panel (a)) and exposed LMO (~1000 yr until a plagioclase flotation crust forms; panel (b)) are indicated as dashed gray lines. For the protolunar disk, we assume an oxygen fugacity buffered by BSE liquid (Equation (30); Visscher & Fegley 2013; Canup et al. 2015). For the LMO, we consider two oxygen fugacities: one at BSE (Equation (30)), and the other at IW-1.5. Loss of Na and K in the protolunar disk at ~3500 K would take ~1 yr, meaning that there would be ample time to reach the levels of depletions seen in the Moon. Loss of Na and K from the LMO at ~1700 K would take 0.5–18 Myr depending on the oxygen fugacity, meaning that a flotation crust would have formed before Na and K had a chance to escape the Moon. The two reasons for the difference in timescales between the two settings are (i) the lower equilibrium vapor pressures and escape flux in the LMO (1700 K) compared to the protolunar disk (3500 K) and (ii) the higher surface-to-volume ratio of a thin disk compared to a sphere.</figcaption></figure>

Mullen & Gammie (2020) made the case that MRI could have heated the disk rapidly, leading to its complete vaporization and spreading to 10 Earth radii within ~4 days. In that scenario, MVEs could have been lost very early, as the disk was being vaporized by MRI and the vapor was being accreted by Earth or possibly lost by a magnetized wind (Gammie et al. 2016). Once all the disk material was vaporized, mass loss by accretion onto Earth would still have taken place, but this would not have induced any chemical or isotopic fractionation because at that point the vapor would have had the composition of the bulk Moon and transport through MRI-driven viscous spreading is not expected to fractionate elements or isotopes. As the

MVE-depleted disk spread outward and cooled, it would have quantitatively condensed, imparting an MVE-depleted and isotopically fractionated signature to the Moon. It is thus conceivable that the heating and vaporization episode took place in a few days following the giant impact. We calculated above an evaporation timescale of ~1 yr for MVE loss from the protolunar disk at 3500 K. The temperature could have been higher and the evaporation timescale shorter. At 4500 K, we calculate, for example, an evaporation timescale of ~40 days, which is an order of magnitude larger than the duration of the disk heating and spreading phase calculated by Gammie et al. (2016) and Mullen & Gammie (2020). Given all sources of uncertainties, the calculated evaporation timescales for Na and K depletions are in reasonable agreement with model predictions for vapor loss from the protolunar disk.

In the Appendix, we show how the model of viscous drainage of volatiles in the protolunar disk is part of a larger class of models that involve advective removal of volatile species from the vapor in proportion to their partial vapor pressures, resulting in evaporation under similar undersaturations of elements with vastly different volatilities. Further work is needed to evaluate whether this condition can be achieved in a different setting than that envisioned by Nie & Dauphas (2019) and the one proposed here, for example, in the synestia model of Lock et al. (2018).

## 4.2. Low-temperature Fractional Evaporation from the Lunar Magma Ocean

Several lines of evidence support the view that the Moon experienced a stage of magma ocean (Smith et al. 1970; Wood et al. 1970; Warren 1985; Charlier et al. 2018). At ~80% crystallization, a flotation crust of plagioclase would have formed a conductive lid, thereby reducing the rate of heat loss to space (Elkins-Tanton et al. 2011). Before that, it is likely that any proto-crust that formed foundered back in the magma ocean and that a large fraction of the lunar surface was covered by magma. Volatile species in such a free-radiating magma ocean could have been volatilized and subsequently lost to space. The Jeans escape parameter for an Na atmosphere above the LMO is $\lambda_e \approx 5$ (Tang & Young 2020, Appendix), corresponding to a slow hydrodynamic or Jeans escape regime (Volkov et al. 2011). Saxena et al. (2017) modeled the evolution of the primordial lunar atmosphere heated by radiation from the post-impact Earth. They examined Jeans escape, which happens through loss at the exobase of atoms, in the tail of the Maxwellian energy distribution, that have vertical velocities in excess of the escape velocity of the Moon. They concluded that Na loss by this process would have been small and could not explain the degree of MVE depletion documented in lunar rocks. We also quantified Na loss from the Moon by Jeans escape, and we concur with this assessment (see the Appendix for details). Tang & Young (2020) made the case that the proximity of the newly formed Moon to Earth would have promoted hydrodynamic over Jeans escape, whereby volatile elements could have been lost through an organized flow (Volkov et al. 2011). They calculated that during hydrodynamic escape each increment of the escaping vapor would have been in near isotopic equilibrium with the surface of the LMO. Sossi et al. (2019) also considered hydrodynamic escape with the escaping vapor in equilibrium with the LMO. They, concluded that equilibrium isotopic fractionation was a viable explanation for the heavy isotopic
composition of the Moon for K and Zn. However, the temperature required of only 1000 K is unrealistic, as it is below the solidus of basalt of ~1300 K. Tang & Young (2020) considered a temperature of 1550 K, and they concluded that loss from the LMO could not account for the heavy K isotope enrichment of the Moon.

We have calculated equilibrium isotopic fractionation factors of K, Rb, Cu, Zn, Ga, and Sn between vapor and liquid (Section 3.1). We can thus reassess the idea put forward by Sossi et al. (2019) that tidal pull of Earth could have stripped the Moon of its volatiles through hydrodynamic escape with each increment of volatile removed in isotopic equilibrium with the LMO. The equation governing such removal would be that of a Rayleigh distillation,

$$ \delta_l \approx \delta_0 - \Delta_{l \to l_v} \ln f, \quad (35) $$

where $\Delta_{l \to l_v} = \Delta_{Eq}$ ($S=1$ in Equation (21)) and $f$ is the fraction of the element remaining in the LMO. In Figure 12, we plot the measured Moon–BSE isotopic differences for K, Rb, Zn, Cu, Ga, and Sn against the predicted difference ($\delta_l - \delta_0$) using our reassessments of the degrees of depletion and equilibrium isotopic fractionation factors (Table 2). We only consider melt species dissolved in silicate liquid (i.e., not metal) because in an LMO setting metal would have segregated into the core. For hydrodynamic escape from the LMO involving equilibrium between liquid and vapor removed, we consider three temperatures: (i) 1000 K, which Sossi et al. (2019) showed could theoretically explain K isotopic data, although at this temperature the lunar surface would not be covered with mafic magma; (ii) 1250 K, which is the solidus temperature of basalt at low pressure; and (iii) 1550 K, which was used by Tang & Young (2020) in their study and is the liquidus temperature of basalt. The degree of depletion and isotopic fractionation of K and Rb are better known than for other MVEs, so we focus on those two elements.

At 1000 K, we calculate Moon–Earth differences of $+0.39\% \pm 0.20\%$ for $\delta^{41}$K and $+0.10\% \pm 0.05\%$ for $\delta^{87}$Rb, which compares well with the measured differences of $+0.41\% \pm 0.07\%$ for $\delta^{41}$K (Wang & Jacobsen 2016; Tian et al. 2020) and $+0.16\% \pm 0.04\%$ for $\delta^{87}$Rb (Pringle & Moynier 2017; Nie & Dauphas 2019). However, at 1000 K, the lunar surface would be solid, and it is an impossible scenario. At 1250 K, the temperature of incipient melting of basalt, we calculate isotopic fractionations of $+0.25\% \pm 0.13\%$ for $\delta^{41}$K and $+0.06\% \pm 0.03\%$ for $\delta^{87}$Rb, which are both lower than the values measured, especially for Rb. A temperature of 1250 K is also not realistic, as, at that temperature, the lunar surface would still be primarily solid. A more realistic minimum temperature is the low-pressure liquidus temperature of basalt of ~1550 K (Tang & Young 2020). At 1550 K, we calculate isotopic fractionations of $+0.16 \pm 0.08\%$ for $\delta^{41}$K and $+0.04 \pm 0.02\%$ for $\delta^{87}$Rb, which are clearly too low and inconsistent with the data. At ~1900–2000 K, the liquidus temperature of BSE (the solidus is at ~1400 K; Hirschmann 2000; Andrault et al. 2011), the isotopic fractionations for K and Rb would be only $+0.07\% \pm 0.04\%$ for $\delta^{41}$K and $+0.02\% \pm 0.01\%$ for $\delta^{87}$Rb. Fractional equilibrium evaporation from the LMO also misses the mark for other elements even at low temperature, most notably Ga. We ascribed very conservative error bars to our calculated force constants ($\pm 50\%$), and the temperatures of ~1550 K should be within

the range of applicability of the harmonic approximation. The discrepancy between predictions and measurements is significant, and we concur with Tang & Young (2020) that if near-equilibrium melt-vapor prevailed at the lunar surface during an episode of hydrodynamic escape (Sossi et al. 2019; Tang & Young 2020), little isotopic fractionation could have been produced in the Moon.

More recently, Charnoz et al. (2021) made the case that they could explain the depletion and isotopic fractionation of MVEs by hydrodynamic escape from the LMO at 1600–1800 K. Tang & Young (2020) had prescribed a diffusive boundary layer in the vapor at the interface with the liquid, which limited escape efficiency and dampened evaporation-driven kinetic isotopic fractionation. Charnoz et al. (2021) argued that this assumption was unwarranted and made the case that tidally driven hydrodynamic escape could (i) generate sufficient undersaturation in the vapor to kinetically fractionate MVE isotopes and (ii) drive the depletions in K and Na observed in lunar rocks on a timescale of less 1000 yr, until a solid flotation crust formed at the top of the magma ocean (Elkins-Tanton et al. 2011; Tang & Young 2020).

Without going into the physics of hydrodynamic escape, we can still address the viability of this model by examining the rate of evaporation from the surface of the magma ocean in a vapor that is ~99% saturated (Figure 13(d)). Equation (28) still applies in this context. The main differences with the model of viscous drainage in the disk are the following:

(i) The volume-to-surface ratio $v/s$ of the Moon is $v/s = R/3 = 579$ km, which is larger than the value of 12 km calculated for the disk.

(ii) The temperature for evaporation from the LMO is only 1600–1800 K (we take 1700 K as the fiducial value). The equilibrium constant of the reaction $AO_{0.5,l} \leftrightarrow A_v + 0.25O_{2,v}$ has a strong temperature dependence (Figures 15(e), (f)), favoring higher equilibrium vapor pressure at higher temperature (Figures 15(a), (b)). The activity coefficients of $NaO_{0.5}$ and $KO_{0.5}$ also show some dependence on temperature (Figures 15(c), (d)), again favoring higher equilibrium vapor pressures at higher temperature (Figures 15(a), (b)).

(iii) The oxygen fugacity relevant to the LMO might have been $\Delta I W = -1.5$ (Karner et al. 2006), which has the following temperature dependence (Huebner 1971): $f_{O_2} = 10^{-27489/T+6.702+\Delta I W}$. While IW-1.5 and BSE (Visscher & Fegley 2013) (Equation (30)) $f_{O_2}$ differ by a factor of ~4 × 10⁵ at 1700 K, the effect of $f_{O_2}$ on evaporation timescale is relatively limited, because of the scaling $\propto f_{O_2}^{1/4}$ (each factor of 10 change in $f_{O_2}$ is a factor of 1.77 change in evaporation timescale; $(4 \times 10^5)^{1/4} \simeq 25$).

We calculated the evaporation timescales for vaporization from the LMO (and loss by hydrodynamic escape) using these parameters ($v/s = 579$ km, $T = 1700$ K, $f_{O_2} = \text{BSE}$) and an undersaturation of ~99%, and we find values of 12 and 18 Myr for Na and K, respectively (Figure 14). Using the same parameters but at $f_{O_2} = \text{IW-1.5}$, the evaporation timescales are 0.5 and 0.7 Myr for Na and K, respectively. Due to rapid cooling, a solid flotation crust of plagioclase would have formed at the surface of the Moon within ~1000 yr of its formation (Elkins-Tanton et al. 2011; Tang & Young 2020). This is much shorter than the depletion timescale calculated above, meaning that loss to space of Na and K from the LMO cannot explain the depletion of these elements from the Moon. The timescale for MVE depletion from the LMO is much
<figure><img src="image_18.png" /></figure>

**Figure 15.** Parameters influencing the kinetics of evaporation of Na and K in the protolunar disk and LMO. (a, b) Equilibrium vapor pressures of Na and K above a liquid of BSE composition as a function of temperature (see main text for “this study”; Visscher & Fegley 2013; Canup et al. 2015; Charnoz et al. 2021). Except otherwise noted (IW-1.5), all curves are given for an *f*<sub>O₂</sub> buffered by a liquid of BSE composition (Visscher & Fegley 2013). The lower Na and K equilibrium vapor pressures calculated here at BSE *f*<sub>O₂</sub> reflect the lower activity coefficients used (calculated from MELTS; Ghiorso et al. 2002). The higher equilibrium vapor pressures reported in Figure 14 of Charnoz et al. (2021) are inconsistent with the thermodynamic data that accompany that figure. (c, d) Activity coefficients Γ of NaO<sub>0.5</sub> and KO<sub>0.5</sub> in silicate liquid as a function of temperature (the results from MELTS are fitted as Γ<sub>NaO0.5</sub> = e<sup>-19201/T+2.914</sup>, Γ<sub>KO0.5</sub> = e<sup>-25614/T+3.7129</sup>, Ghiorso et al. 2002; Visscher & Fegley 2013; Canup et al. 2015; Charnoz et al. 2021) and comparison with existing experimental data in coal ash slags (Mueller et al. 2004). (e, f) Temperature dependence of the equilibrium constants of the reactions NaO<sub>0.5,l</sub> ↔ Na<sub>v</sub> + 0.25O<sub>2,v</sub> and KO<sub>0.5,l</sub> ↔ K<sub>v</sub> + 0.25O<sub>2,v</sub> (calculated from Equations (33) and (34) using the relationship lnK = -ΔG/RT; Lamoreaux and Hildenbrand 1984; Visscher & Fegley 2013; Canup et al. 2015; Charnoz et al. 2021). Increases in activity coefficient (panels (c) and (d)) and equilibrium constant (panels (e) and (f)) with increasing temperature lead to higher equilibrium vapor pressures and faster vaporization.

longer than depletion from the protolunar disk because of the
low temperature (low equilibrium vapor pressure and low
evaporation flux at a given undersaturation; Figure 14) and low
surface-to-volume ratio involved in the LMO setting.

Our conclusion that loss from the LMO is not a viable
explanation for the depletions in Na and K of the Moon
disagrees with the conclusion of Charnoz et al. (2021), who

calculated that such depletions could be achieved if hydro-
dynamic escape was sustained for 100–6000 yr. For compar-
ison, at BSE oxygen fugacity and a temperature of 1600–1800
K considered by Charnoz et al. (2021) (above 1800 K,
significant Li would be lost from the Moon, which is not
observed), we calculate depletion timescales of 3–61 Myr for
Na and 4–103 Myr for K (Figure 14). At a lower oxygen
<figure><img src="image_19.png" /></figure>

**Figure 16.** Evaporation timescales for reaching Moon-like depletions for Na in the model of Charnoz et al. (2021) of tidally driven hydrodynamic escape from the LMO. Here we focus on the evaporation step involving transfer of Na (panel (a)) and K (panel (b)) from the LMO to the lunar atmosphere, not transport in the hydrodynamic wind. While hydrodynamic escape may have been channeled along the line connecting the centers of Earth and the Moon (Charnoz et al. 2021), transfer to the atmosphere could have taken place over a wider surface area, and we conservatively calculate the evaporation flux over the whole lunar surface area (a solid angle of $4\pi$ sr). The oxygen fugacity is buffered by BSE liquid (Visscher & Fegley 2013; Canup et al. 2015; Equation (30)), the same assumption that was made by Charnoz et al. (2021). The blue solid curves are from this study, using the parameterization for $\gamma\Gamma$ from Zhang et al. (2021, Equations (31) and (32)) and a vapor saturation $S = 0.99$ (see text for details; also see Nie & Dauphas 2019). The black solid lines were calculated using the equilibrium vapor pressure from Canup et al. (2015) for K ($PK_{eq} = 10^{-4.7 \times 10^6/T^2 - 11935/T + 1.5825$), with $PK_{eq}$ in bar and $T$ in K) and Visscher & Fegley (2013) for Na ($PN_{eq} = 10^{(-4.5e6/T^2 - 11629/T + 3.0586)}$, $S = 0.99$, and evaporation coefficients $\gamma$ of 0.1 (Alexander 2002; Fedkin et al. 2006; Richter et al. 2011; Zhang et al. 2021). The dashed red lines were calculated by using the thermodynamic data provided by Charnoz et al. (2021) and assuming $S = 0.99$ and $\gamma = 0.1$. As shown, the calculated evaporation timescales at 1700 K (the temperature preferred by Charnoz et al. 2021) are all orders of magnitude longer than the lifetime of the exposed LMO (>0.8 Myr for K and >1.2 Myr for Na, compared to ~1000 yr for the LMO). The only way that we have found to decrease the evaporation timescales to values comparable to the lifetime of the exposed LMO is by taking $S = 0.8$ and $\gamma = 1$ and using the digitized K and Na equilibrium vapor pressures from Figure 14 of Charnoz et al. (2021). However, isotopic data call for evaporation in a near-saturated medium ($S = 0.99$; this study; Nie & Dauphas 2019; Charnoz et al. 2021), free evaporation experiments support an evaporation coefficient of $\gamma = 0.1$ (Alexander 2002; Fedkin et al. 2006; Richter et al. 2011; Zhang et al. 2021), and the K and Na equilibrium vapor pressures plotted in Figure 14 of Charnoz et al. (2021) are inconsistent with the thermodynamic data that accompany this figure.

fugacity of IW–1.5 (and the same temperature of 1600–1800 K), the evaporation timescales for Na and K decrease to 0.1–1.8 and 0.2–3.1 Myr, respectively. There are several possible compounding explanations for the discrepancy of several orders of magnitude between the present study and Charnoz et al. (2021) (Figure 16):

(i) The evaporation timescale is inversely proportional to the equilibrium vapor pressure. Visscher & Fegley (2013) calculated the Na and K equilibrium vapor pressure at a BSE-buffered $f_{O_2}$. The vapor pressure of K in equilibrium with BSE was later updated in Canup et al. (2015). We calculated the equilibrium vapor pressure at the same BSE-buffered $f_{O_2}$ and find values ~10 times lower than Visscher & Fegley (2013) for Na and ~4 times lower than Canup et al. (2015) for K at 1700 K. This is because (1) the activity coefficients that we use from MELTS are lower than those used by Visscher & Fegley (2013) and Canup et al. (2015) by factors of ~11.6 and ~7.3 for $NaO_{0.5}$ and $KO_{0.5}$, respectively, and (2) the equilibrium constants given by Equations (33) and (34) are higher than those of Visscher & Fegley (2013) and Canup et al. (2015) by factors of ~1.1 and ~1.8 for Na and K, respectively. Charnoz et al. (2021) used the activity coefficients for $NaO_{0.5}$ for the $Na_2O-SiO_2$ solid solution (Charles 1967; ~2 times lower than Visscher & Fegley 2013) and assumed that the activity coefficients for $KO_{0.5}$ were half of those for $NaO_{0.5}$ (~13.7 times higher than Canup et al. 2015). We used the thermodynamic data provided by Charnoz et al. (2021) to calculate the Na and K equilibrium vapor pressure at the BSE-buffered $f_{O_2}$ and find that the value for Na is a factor of ~1.8 lower than that calculated by Visscher & Fegley (2013), while the value for K is a factor of ~4.7 higher than that calculated by Canup et al. (2015). Differences in thermodynamic data can

only account for differences in evaporation timescales of factors of ~5 for Na and ~20 for K between Charnoz et al. (2021) and this study. We also digitized the Na and K equilibrium vapor pressure plotted in Figure 14 of Charnoz et al. (2021) and found that at 1700 K they were a factor of ~23 higher than the values that we calculate from the thermodynamic data that they provide for both Na and K, a factor of ~13 higher than the value given by Visscher & Fegley (2013) for Na, and a factor of ~109 higher than the values given by Canup et al. (2015) for K (Figure 16). If Charnoz et al. (2021) used the erroneous Na and K equilibrium vapor pressures from their Figure 14, this could shift the evaporation timescales by factors of ~135 and 460 for Na and K, respectively, relative to our estimates.

(ii) The evaporation timescale is inversely proportional to the evaporation coefficient ($\gamma_i$ in Equation (29)). For Na in silicate melts, the evaporation coefficients reported in the literature are 0.12 (Alexander 2002), 0.26 (Fedkin et al. 2006), 0.06–0.22 (Richter et al. 2011), 0.078 and 0.14 (Zhang et al. 2021). Charnoz et al. (2021) do not discuss the value of the evaporation coefficient of Na, but judging from their writing of the Hertz–Knudsen equation (their Equation (D13)), it is likely that they adopted a value of 1 for $\gamma$. If we conservatively take 0.1 for the evaporation coefficients of both Na and K as given in the literature, this could shift the evaporation timescale of Na and K by a factor of ~10.

(iii) The evaporation timescale is inversely proportional to $1 - S_i$, where $S_i$ is the vapor saturation. As shown by Nie & Dauphas (2019) and reiterated here, the heavy isotopic compositions of most MVEs require a vapor saturation of ~0.99. Although Charnoz et al. (2021) recognize that such a saturation is needed to explain the heavy isotope enrichments in K and Zn,
they do not use it explicitly as a constraint to calculate evaporation timescales. They calculate that at $T < 1800$ K, tidally driven hydrodynamic escape could sustain an undersaturation of 0.8 or more. Such a low saturation would, however, induce large kinetic isotopic fractionation, which is not seen. Using such an undersaturation could shift the evaporation timescale by a factor of $(1 - 0.8)/(1 - 0.99) = 20$.

The three factors highlighted above (possible erroneous equilibrium vapor pressures for Na and K, evaporation coefficients of 1, and 80% vapor undersaturation) would all contribute to lowering the Na and K evaporation timescale compared to what it should be. The overall shift could be as high as $23 \times 10 \times 20 = 4600$ for both Na and K, lowering the million-year evaporation timescale calculated here to the hundred-to-thousand-year timescale calculated by Charnoz et al. (2021). We conclude that evaporation kinetics is the main barrier to loss of MVEs from the LMO, and a flotation crust would form on the LMO before significant amounts of Na, K, and Rb can escape.

### 4.3. High-temperature Partial Condensation in Disk and Synestia Settings

We consider here a scenario involving partial condensation at high temperature (~3500 K), which would be like the models of Canup et al. (2015) and Lock et al. (2018). The vapor is assumed to have BSE isotopic composition. Kinetic isotopic fractionation will tend to enrich the partial condensate (i.e., the Moon) in the light isotopes, which is opposite to what is observed for all elements except Sn. In the context of the models of Canup et al. (2015) and Lock et al. (2018), we can therefore rule out condensation from a significantly supersaturated medium. We examine below whether the data can be explained by condensation under near-equilibrium conditions, for example, if the cooling timescale is much longer than the condensation timescale (Richter 2004; Hu et al. 2021). Isotopic fractionation during fractional condensation can be described using a Rayleigh distillation equation (each increment of liquid condensed is in equilibrium with the vapor, but once condensed, atoms are isolated),

$$ \delta_l \approx \delta_0 - \Delta_{v \to l - v} (1/f - 1) \ln(1 - f), \quad (36) $$

where f is the fraction of the element condensed (strictly speaking of the denominator isotope, e.g., ³⁹K) and Δ$_{v \to l - v}$ is the isotopic fractionation given by Equation (22), which will depend on temperature and the degree of supersaturation S. We set S = 1 because kinetic effects will enrich the condensate in the light isotopes, which is opposite to what is observed. The fraction of MVE condensed would be small, so using a batch equilibrium model where the bulk liquid stays in equilibrium with the vapor would not change the results significantly. We have used the known degrees of depletion of MVEs in the Moon (f; Table 2) and equilibrium fractionation factors (Δ$_{Eq}$; Section 3, Table 2) to calculate the isotopic composition of the condensate liquid for a model of fractional equilibrium condensation, and the results are plotted in Figure 17. As shown, the predicted isotopic fractionation in lunar rocks relative to terrestrial rocks is much smaller than what is measured. The large discrepancy is because equilibrium isotopic fractionation decreases with the square of the inverse of the temperature, and it is extremely small at 3500 K.

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Measured isotopic fractionations between the Moon and BSE in δ notation against predictions for fractional equilibrium condensation at 3500 K (Equation (35); Section 4.3). The model predictions correspond to what would be expected for the models of MVE depletion of Canup et al. (2015) and Lock et al. (2018), where MVE depletion happens very soon after the giant impact owing to incomplete condensation at high temperature and loss of MVEs to Earth.</figcaption></figure>

Furthermore, the condensation model does not provide any leverage to magnify instantaneous isotopic fractionation between liquid and vapor. Indeed, during a distillation the largest effects are produced in the residue, rather than the product reservoir. The largest fractionation is at the onset of condensation, and as condensation proceeds, the isotopic composition of the condensate approaches that of the starting material (silicate Earth). Our calculation thus indicates that the heavy K and Rb isotopic compositions of lunar rocks cannot be produced solely by partial condensation. Both the Canup et al. (2015) and Lock et al. (2018) models involve settings where melt parcels or moonlets could have experienced evaporative loss of alkali element and MVEs, so further work will be needed to evaluate whether the heavy isotopic enrichment of lunar rocks can be explained in the framework of those models. For example, in canonical impact models, evaporation under highly undersaturated conditions could have taken place while the disk underwent rapid collisional heating before it reached the quasi-steady-state conditions identified by Ward (2017). It is thus conceivable that the bulk of lunar MVE depletion was produced by partial condensation with little isotopic fractionation induced, and that the heavy isotopic compositions measured in lunar rocks came from a small fraction of material that experienced evaporation in a highly undersaturated medium, resulting in extensive isotopic fractionation.

## 5. Conclusion

Several giant impact models can explain the dynamical and physical properties of the Earth-Moon system, but these models cannot readily explain why the Moon is depleted in MVEs relative to Earth. Several recent studies have shown that lunar rocks were enriched in the heavy isotopes of several
MVEs, notably K and Rb. Because the equilibrium isotopic fractionation factors between condensate and vapor were unknown, the message that these heavy isotope enrichments conveyed about the conditions of volatile element depletion in the Moon was uncertain.

We have used an ab initio approach to calculate the reduced partition function ratios (β-factors) of various minerals and gas species and find that condensation or evaporation where each parcel of vapor is in equilibrium with the liquid cannot explain the heavy isotopic composition of the Moon relative to Earth. The heavy isotope enrichment of the Moon is more readily explained by partial vaporization in a gas at an undersaturation of ~99% for most MVEs, regardless of element volatility. Constant undersaturation during evaporation can be achieved if the vapor is advected away and elements are removed in proportion to their partial vapor pressures. A possible setting for such evaporation is viscous drainage onto Earth of the vapor layer in the protolunar disk, possibly powered by MRI. The timescale needed to deplete Na and K in the protolunar disk to levels encountered in lunar rocks is typically a month to a year, which is longer than the initial disk heating and vaporization phase that can last just a few days and is shorter than the lifetime of the disk of ~100 yr, meaning that it is a viable scenario. Loss of MVEs from the LMO by hydrodynamic escape is not favored because the timescale to achieve lunar-like depletion levels is ~0.1–103 Myr, which is much longer than the ~1000 yr lifetime of an exposed magma ocean free of flotation crust. The stark difference in timescales between protolunar disk and LMO settings stems from differences in temperature and evaporation geometries.

This work was supported by NASA grants NNX17AE86G, NNX17AE87G, 80NSSC17K0744, and 80NSSC20K0821 and NSF grant EAR-2001098 to N.D. DFT calculations were performed using the HPC resources from CALMIP (grant 2020-P1037) and the University of Chicago Research Computing Center. We thank two anonymous reviewers for their constructive comments and Edgard Rivera-Valentin for his handling of this manuscript.

## Appendix

### A.1. Bernoulli Expansion of β-factors for Diatomic Molecules

Supplementary Online Material: Mathematica code used to calculate the behavior of trace elements during LMO crystallization.

Bigeleisen (1958) carried out a Bernoulli expansion of the reduced partition function ratio, which for a diatomic molecule takes the form

$$ \ln \beta_{b/a} = \sum_{j=2}^{\infty} \frac{B_j}{j!} (u_a^j - u_b^j), \qquad (A1) $$

where $B_j$ are the Bernoulli numbers ($B_2 = 1/6$, $B_3 = 0$, $B_4 = -1/30$, $B_5 = 0$, $B_6 = 1/42$) and $u_a$ and $u_b$ are the reduced frequencies for the molecules substituted with isotopes a and b, respectively ($u = h\nu/kT$, with h the Planck constant, $\nu$ the vibrational frequency, k the Boltzmann constant, and T the temperature). For a diatomic molecule behaving as a harmonic oscillator made of two distinct atoms of masses m and M, we can relate the vibration frequency to the force constant of the

bond A through

$$ \nu = \frac{1}{2\pi} \sqrt{\frac{m + M}{mM}} A. \qquad (A2) $$

If we denote $m_a$ and $m_b$ the masses of isotopes a and b and inject Equation (A2) into Equation (A1), it follows that (with $\hbar = h/2\pi$)

$$ \ln \beta_{b/a} = \sum_{j=2}^{\infty} \frac{B_j}{j!j} \left( \frac{\hbar A^{1/2}}{kT} \right)^j \left[ \left( \frac{1}{M} + \frac{1}{m_a} \right)^{j/2} - \left( \frac{1}{M} + \frac{1}{m_b} \right)^{j/2} \right]. \qquad (A3) $$

We can write the last term as

$$ \left(\frac{1}{M} + \frac{1}{m_b}\right)^{j/2} = \left(\frac{1}{M} + \frac{1}{m_a} + \frac{1}{m_b} - \frac{1}{m_a}\right)^{j/2} \\ = \left(\frac{1}{M} + \frac{1}{m_a}\right)^{j/2} \left(1 + \frac{\frac{1}{m_b} - \frac{1}{m_a}}{\frac{1}{M} + \frac{1}{m_a}}\right)^{j/2}, \qquad (A4) $$

which can be approximated as

$$ \left(\frac{1}{M} + \frac{1}{m_b}\right)^{j/2} \approx \left(\frac{1}{M} + \frac{1}{m_a}\right)^{j/2} \left(1 + \frac{j \frac{1}{m_b} - \frac{1}{m_a}}{2 \frac{1}{M} + \frac{1}{m_a}}\right). \qquad (A5) $$

Injecting this formula in Equation (A3), we have

$$ 1000 \ln \beta_{b/a} \approx 1000 \left( \frac{1}{m_a} - \frac{1}{m_b} \right) \sum_{j=2}^{\infty} \frac{B_j}{j!2} \left( \frac{\hbar A^{1/2}}{kT} \right)^j \left( \frac{1}{M} + \frac{1}{m_a} \right)^{j/2-1}, \qquad (A6) $$

or using the reduced mass $\mu = 1/(1/M + 1/m) \approx 1/(1/M + 1/m_a)$,

$$ 1000 \ln \beta_{b/a} \approx 1000 \left( \frac{1}{m_a} - \frac{1}{m_b} \right) \sum_{j=2}^{\infty} \frac{B_j}{j!2} \left( \frac{\hbar A^{1/2}}{kT} \right)^j \mu^{1-j/2}. \qquad (A7) $$

We can also replace the force constant by the wavenumber $\tilde{\nu}$ (typically given in cm⁻¹) using the relationship $A = \mu(2\pi c\tilde{\nu})^2$,

$$ 1000 \ln \beta_{b/a} \approx 1000 \left( \frac{1}{m_a} - \frac{1}{m_b} \right) \mu \sum_{j=2}^{\infty} \frac{B_j}{j!2} \left( \frac{hc\tilde{\nu}}{kT} \right)^j. \qquad (A8) $$

For a diatomic molecule made of the same atoms (e.g., O₂ or N₂), the formulae would be

$$ 1000 \ln \beta_{b/a} \approx 1000 \left( \frac{1}{m_a} - \frac{1}{m_b} \right) \sum_{j=2}^{\infty} \frac{B_j}{j!2} \left( \frac{\hbar\sqrt{A}}{kT} \right)^j \left( \frac{2}{m_a} \right)^{j/2-1}, \qquad (A9) $$

$$ 1000 \ln \beta_{b/a} \approx 1000 \left( \frac{1}{m_a} - \frac{1}{m_b} \right) \frac{m_a}{2} \sum_{j=2}^{\infty} \frac{B_j}{j!} \left( \frac{hc\tilde{\nu}}{kT} \right)^j. \qquad (A10) $$
If we use the directionally averaged force constant $\langle F \rangle = A/3$, Equation (A7) takes the form

$$1000 \ln \beta_{b/a} \approx 1000 \left( \frac{1}{m_a} - \frac{1}{m_b} \right) \\ \sum_{j=2}^{\infty} \frac{B_j}{j!2} \left( \frac{\hbar \sqrt{3} \langle F \rangle}{kT} \right)^j \mu^{1-j/2}. \qquad (A11)$$

As an example, we apply Equations (A7), (A8), and (A11) to the $^{41}$K/$^{39}$K ratio in the diatomic molecule KO, which, truncated to the third order, yields

$$1000 \ln \beta \approx 1000 \left( \frac{1}{m^{39}_{\text{K}}} - \frac{1}{m^{41}_{\text{K}}} \right) \\ \times \left( \frac{\hbar^2 A}{24k^2T^2} - \frac{\hbar^4 A^2}{1440k^4T^4\mu} + \frac{\hbar^6 A^3}{60480k^6T^6\mu^2} \right), \qquad (A12)$$

$$1000 \ln \beta \approx 1000 \left( \frac{1}{m^{39}_{\text{K}}} - \frac{1}{m^{41}_{\text{K}}} \right) \\ \times \mu \left[ \frac{1}{24} \left( \frac{hc\tilde{\nu}}{kT} \right)^2 - \frac{1}{1440} \left( \frac{hc\tilde{\nu}}{kT} \right)^4 + \frac{1}{60480} \left( \frac{hc\tilde{\nu}}{kT} \right)^6 \right], \qquad (A13)$$

$$1000 \ln \beta \approx 1000 \left( \frac{1}{m^{39}_{\text{K}}} - \frac{1}{m^{41}_{\text{K}}} \right) \\ \times \left( \frac{\hbar^2 \langle F \rangle}{8k^2 T^2} - \frac{\hbar^4 \langle F \rangle^2}{160k^4 T^4 \mu} + \frac{\hbar^6 \langle F \rangle^3}{2240k^6 T^6 \mu^2} \right). \qquad (A14)$$

## A.2. Generic model of MVE depletion under constant undersaturation

The model of Nie and Dauphas (2019) of viscous drainage of vapor in the protolunar disk can explain the isotopic fractionations and depletions of all MVEs considered here, except Sn. Below we show how this model is part of a larger class of models involving advective removal of vapor species proportionally to their partial vapor pressures, which could be relevant to other settings such as formation of chondrules and refractory inclusions. According to the Hertz–Knudsen equation, the net flux for each element i during evaporation is given by

$$J_i = \sum_{j} \frac{\gamma_{j,i} (P_{eq,j,i} - P_{j,i})}{\sqrt{2\pi m_{j,i} RT}}, \qquad (A15)$$

where $j$ enumerates the vapor species (e.g., Na$_v$, NaO$_v$), $\gamma_{j,i}$ and $m_{j,i}$ are the condensation/evaporation coefficients and masses of the corresponding species, $R$ is the gas constant, $P_{j,i}$ and $P_{eq,j,i}$ are the partial and equilibrium vapor pressures, respectively. We consider a model where elements are transferred in a finite gas volume $V_g$, from which they are advected away in proportion to their total abundance in the gas volume ($J_{a,j,i}$ is the rate of advective loss in moles per unit time),

$$J_{a,j,i} = -a \frac{P_{j,i} V_g}{RT}, \qquad (A16)$$

with $a$ a constant. The change in the number of moles of the element considered in the gas is a balance between supply by evaporation from the solid of surface area $A$ (Equation (A15)) and removal by advection (Equation (A16)). For each species of the element considered we have

$$\frac{dn_{g,\mu,i}}{dt} = -a \frac{V_g}{RT} P_{\mu,i} + \frac{\gamma_{\mu,i} (P_{eq,\mu,i} - P_{\mu,i})}{\sqrt{2\pi m_{\mu,i} RT}} A \\ + \frac{V_g}{RT} \sum_{q, q' \neq \mu} (\kappa_{q,\mu,i} P_{q,i} - \kappa_{q,\mu,q,i} P_{\mu,i}), \qquad (A17)$$

where $n_{g,\mu,i}$ is the total number of moles of species $j$ of element $i$ in the gas volume $V_g$, $A$ is the surface area of the evaporating solid, and $\kappa$ are kinetic constants for the exchange reactions between the gas species ($\kappa_{q,\mu,i}$ is the rate of conversion of species $q$ into $\mu$ of element $i$; $\kappa_{q,\mu,q,i}$ is the rate of conversion of species $\mu$ into $q$ of element $i$). If we run a summation over all gas species, the kinetic terms cancel out and we have

$$\frac{dn_{g,i}}{dt} = -a \frac{V_g}{RT} P_i + \sum_{j} \frac{\gamma_{j,i} (P_{eq,j,i} - P_{j,i})}{\sqrt{2\pi m_{j,i} RT}} A. \qquad (A18)$$

In the following, we use the notations $S_i = P_i/P_{eq,i}$ and $S_{\mu,i} = P_{\mu,i}/P_{eq,\mu,i}$ to denote saturation. At steady state, $dn_{g,i}/dt = 0$ and it follows that

$$a \frac{V_g}{RT} S_i = \sum_{j} \frac{\gamma_{j,i} (1 - S_{\mu,i}) P_{eq,j,i}}{\sqrt{2\pi m_{\mu,i} RT} P_{eq,i}} A. \qquad (A19)$$

The ratio $x_{\mu,i} = P_{eq,\mu,i}/P_{eq,i}$ is known from equilibrium thermodynamics, and we have

$$a \frac{V_g}{RT} S_i = \sum_{j} \frac{\gamma_{j,i} (1 - S_{\mu,i}) x_{\mu,i}}{\sqrt{2\pi m_{\mu,i} RT}} A. \qquad (A20)$$

The different evaporated gas species can be out of equilibrium owing to differences in $\gamma_{\mu,i}$ and $m_{\mu,i}$, but chemical reactions in the gas would allow the different gas species to re-equilibrate partially or completely. Assuming that the gas residence time is sufficient to allow thermodynamic equilibration of the vapor species and that the proportions of the gas species correspond approximately to what they would be if the vapor was in equilibrium with the liquid, we write

$$P_{\mu,i} \approx \frac{P_{eq,\mu,i}}{P_{eq,i}} P_i. \qquad (A21)$$

With these assumptions the saturations of all species are approximately the same,

$$S_{\mu,i} \approx S_i. \qquad (A22)$$

Injecting this formula in Equation (A20), we have the following relationship for the overall saturation at steady state:

$$S_i = (1 - S_i) \sum_j \frac{A \gamma_{j,i} x_{j,i}}{a V_g} \sqrt{\frac{RT}{2\pi m_{j,i}}} . \qquad (A23)$$