Vapor Drainage in the Protolunar Disk as the Cause for the Depletion in Volatile Elements of the Moon

The prevailing view of the origin of the Moon is that it formed by a collision between the proto-Earth and a large impactor that has been named Theia, which produced a disk of liquid and vapor from which the Moon grew (Hartmann & Davis 1975; Cameron & Ward 1976). Simulations of this giant impact can explain the angular momentum of the Earth–Moon system and the small mass of the lunar core, but other chemical and isotopic constraints are more challenging to reproduce (Canup & Asphaug 2001; Pahlevan & Stevenson 2007; Canup 2012; Cuk & Stewart 2012; Reufer et al. 2012; Nakajima & Stevenson 2015; Lock et al. 2018). Among these is the lunar depletion in moderately volatile elements (MVEs) such as Rb, K, and Zn (Ringwood & Kesson 1977; O'Neill 1991; Albarède et al. 2015; Righter 2019). The settings considered to explain these depletions are (i) incomplete condensation of a vapor cloud surrounding the Earth (synestia) in the aftermath of a very energetic impact event (Lock et al. 2018), (ii) accretion of the Moon from a partially condensed, volatile-depleted silicate melt carried across the Roche limit by viscous spreading (Canup et al. 2015), (iii) removal of volatile elements from a stratified protolunar disk by accretion of the vapor layer onto Earth (Charnoz & Michaut 2015; Wang et al. 2019), and (iv) loss to space from the lunar magma ocean (Day & Moynier 2014; Saxena et al. 2017). One difficulty in assessing these scenarios is the lack of thermodynamic data constraining the condensation/evaporation behaviors of MVEs at the very high temperatures (presumably >3500 K; Canup et al. 2015; Lock et al. 2018) involved in the formation of the Moon. Isotopic analyses of several MVEs have revealed heavy isotope enrichments of lunar rocks compared to Earth (Herzog et al. 2009; Paniello et al. 2012; Kato et al. 2015; Wang & Jacobsen 2016; Kato & Moynier 2017; Pringle & Moynier 2017; Wimpenny et al. 2019), but the extents of those fractionations are uncertain, and a quantitative understanding of how they relate to the conditions of the Moon formation is still missing.

We report here new and high-precision Rb isotopic analyses of lunar and terrestrial samples that we use to understand why the Moon is depleted in MVEs relative to Earth. In Section 2, we report the results of Rb isotopic analyses of Apollo lunar samples (primarily basalts) indicating that the Moon is enriched in the heavy isotopes of Rb relative to Earth. In Section 3, we calculate the saturation of Rb (∼99%) in the vapor medium needed to explain the heavy Rb isotope enrichment of lunar rocks by evaporation, and show that the same degree of saturation can account for previously documented heavy isotope enrichments in other MVEs K, Ga, Cu, and Zn. In Section 4, we show that this saturation can be explained in the context of viscous drainage of vapor in the protolunar disk onto Earth if the α-viscosity of the vapor layer was 10⁻³–10⁻², which is consistent with viscosities predicted if the magnetorotational instability (MRI) was active.

We have measured the isotopic compositions of Rb in six Apollo lunar samples comprising five mare basalts and one norite (Table 1) to better constrain the bulk Rb isotopic composition of the Moon, which is expressed in δ⁸⁷Rb notation (the permil departure of the ⁸⁷Rb/⁸⁵Rb ratio relative to reference material NIST SRM984). To compare to the Earth, we have also measured the Rb isotopic compositions of terrestrial basalts and granites (Table 1). Pringle & Moynier (2017) previously measured the Rb isotopic compositions of lunar rocks and found a barely resolvable difference in δ⁸⁷Rb of +0.17 ± 0.13‰ between the Moon and Earth using a protocol that presumably did not fully separate Rb from K (see the Appendix). The motivations for our new Rb isotopic analyses were therefore to improve on the precision of these analyses and test their accuracy using a protocol that fully separated K from Rb (see the Appendix for details). Basalt samples have been prioritized because the chemical and isotopic compositions of other lithologies could have been more severely affected by lunar magma ocean differentiation. Indeed, Rb is incompatible during magmatic processes but it can partition in feldspar. It is thus conceivable that it would have been fractionated isotopically when an anorthositic flotation crust formed from the lunar magma ocean. For example, the cataclastic norite sample 77215 from Apollo 17 measured in this study, which is a white-colored fragment containing abundant plagioclase, shows a fractionated Rb isotopic composition (−0.15 ± 0.03‰) compared to mare basalts (average +0.03 ± 0.03‰). All basalts have similar δ⁸⁷Rb and δ⁴¹K values (Figure 1; Tables 1 and 2), so it is unlikely that their isotopic compositions were affected by Rb and K volatilization during eruption at the lunar surface.

In the following, we compare Rb with K as these are two alkali elements with similar geochemical behaviors (O'Neill 1991). To mitigate the effect of magmatic fractionation, we plot the δ⁸⁷Rb (this study and Pringle & Moynier 2017) and δ⁴¹K (Wang & Jacobsen 2016) values of lunar and terrestrial rocks against the La/U ratios. The motivation for doing so is that La and U have a similar geochemical behavior to Rb and K during magmatic processes (O'Neill 1991), but unlike the latter, La and U are two refractory elements that were immune to volatile-loss processes and should therefore be in chondritic proportion in the bulk Moon and Earth. Samples with near-chondritic La/U ratios (∼27.6) are therefore less likely than others to have had their Rb and K isotopic compositions fractionated by magmatic differentiation. Focusing on samples with La/U ratios within 0.5× to 1.5× the chondritic value, we estimate that the Moon and Earth differ in their δ⁸⁷Rb and δ⁴¹K values by +0.16 ± 0.06 and +0.40 ± 0.05‰, respectively.

An alternative approach for estimating the Rb and K isotopic compositions of the Earth and Moon is to interpolate the trends between δ⁸⁷Rb, δ⁴¹K, and La/U ratios to the chondritic La/U ratio (Figure 1). Variations in the La/U ratio most likely reflect magmatic differentiation processes (including mixing with a KREEP component characterized by (La/U)/(La/U)_CI ≃ 0.65; Warren 1989), and we take the δ⁸⁷Rb and δ⁴¹K values at (La/U)/(La/U)_CI = 1 (unfractionated La/U ratio) as representative of the bulk Moon. Accordingly, we estimate that the δ⁸⁷Rb and δ⁴¹K values of the Moon are heavier than those of the Earth by +0.16 ± 0.04‰ and +0.41 ± 0.07‰, respectively, which agrees with the values given above by taking the average of all samples with near-chondritic La/U ratios. The slightly positive slope of lunar δ⁸⁷Rb values against La/U ratios is controlled by the two heaviest basalts, but excluding them from the regression does not change the estimated bulk lunar δ⁸⁷Rb value.

Zoom In Zoom Out Reset image size

In some models of lunar formation (scenarios (i) and (ii) above in Section 1), the depletions in Rb and K are interpreted to reflect incomplete condensation (Canup et al. 2015; Lock et al. 2018). However, this process cannot explain the heavy isotopic compositions of Rb and K in the Moon because at the temperature of ∼3500 K considered in these models, equilibrium isotopic fractionation between vapor (as Rb and K atoms) and condensate is too small to account for the observed values (Dauphas et al. 2018b), and kinetic isotopic fractionation would enrich the partial condensate (the Moon) in the light isotopes (Richter 2004), which is opposite to what is observed. Most likely, the elevated δ⁸⁷Rb and δ⁴¹K values of the Moon relative to Earth reflect partial evaporation. A model of evaporation from the lunar magma ocean (scenario (iv) above) indicates that for Na, whose volatility is similar to K, the cumulative loss over the lifetime of the magma ocean is small (5%–20%; Saxena et al. 2017), which is insufficient to account for the six-fold depletion of K and Rb in the Moon relative to Earth. More work, however, is needed to fully assess this scenario, especially for the more volatile element Zn (Day & Moynier 2014; Charnoz et al. 2019; Young & Tang 2019).

Loss of MVEs could have happened through drainage and loss of vapor from the protolunar disk onto the proto-Earth (scenario (iii) above; Charnoz & Michaut 2015). In this model, it is assumed that the protolunar disk is stratified, comprising a liquid layer at the midplane surrounded above and below by a vapor layer. MVEs are evaporated from the liquid layer into the vapor layer. The fate of MVEs and their depletions in the Moon depend sensitively on the viscosity of the vapor layer. The reason is that viscosity mediates transport of mass and angular momentum across the disk. Little loss of liquid from the disk to the Earth could have taken place over much of the cooling history of the disk, because a gravitationally stable low-viscosity region would be present within ∼1.7 R_⊕ that would have acted as a barrier preventing liquid from being accreted by the Earth. If the viscosity of the vapor was low, it would have remained stagnant and no volatile loss would have occurred. If the viscosity of the vapor was high due to MRI (Carballido et al. 2016; Gammie et al. 2016), the vapor layer carrying MVEs would have been efficiently accreted by the Earth, leaving a disk depleted in volatile elements.

Drainage of the vapor layer onto Earth is a complex process that would have involved some evaporation from the liquid to the vapor layer (Charnoz & Michaut 2015). The kinetics and isotopic fractionation associated with such evaporation depend on the mode of transport inside (advective versus diffusive) and across (equilibrium versus kinetic) the liquid and vapor layers (Craig & Gordon 1965; Richter et al. 2002, 2007; Richter 2004; Dauphas et al. 2015). Given the rate of heat loss experienced by the disk through radiation to space, much of the vapor layer would have convected vigorously (Thompson & Stevenson 1988; Ward 2011). Furthermore, given the difference in orbital velocities between the vapor and liquid layers, it is likely that Kelvin–Helmholtz instabilities developed (Thompson & Stevenson 1988; Charnoz & Michaut 2015) that prevented the development of a viscous layer at the vapor–liquid interface. In this context, where transport within the liquid and vapor layers is entirely advective, the extent of isotopic fractionation of element i during evaporation depends (through the Hertz–Knudsen equation) on its vapor saturation (Richter et al. 2002, 2007; Dauphas et al. 2015), expressed as ${S}_{i}={P}_{i}/{P}_{i,\mathrm{eq}}$, where P_i is the partial vapor pressure of element i and P_i,eq is its equilibrium vapor pressure at the relevant temperature.

If S_i = 1, there is no net evaporative flux, the vapor is in thermodynamic equilibrium with the liquid, and the isotopic fractionation between vapor and liquid is entirely equilibrium in nature (noted as ${{\rm{\Delta }}}_{\mathrm{eq},i}^{v-l}$ hereafter). Ab initio calculations give the equilibrium fractionation between monoatomic K vapor (the dominant gas species in conditions relevant to lunar formation) and K-feldspar (taken as a proxy for K in silicate melt) to be ${{\rm{\Delta }}}_{{\rm{e}}{\rm{q}},{\rm{K}}}^{v-l}\simeq -0.012$ to −0.024‰ at 3500 K and −0.055 to −0.107‰ at 1650 K (Dauphas et al. 2018b; Zeng et al. 2018; Li et al. 2019; Zeng et al. 2019). Note that in silicate melts with (K, Na)/Al ratios of <1 relevant to magmas of bulk silicate Earth or Moon compositions, K and Rb charge balance Al³⁺ in tetrahedral coordination, so to first order feldspar is a reasonable model structure for those melts (Seifert et al. 1982; Toplis et al. 1997). Rubidium forms chemical bonds of similar strengths to K (Zeng et al. 2019), and we estimate that the vapor–liquid equilibrium fractionation for Rb should be ${{\rm{\Delta }}}_{{\rm{e}}{\rm{q}},{\rm{R}}{\rm{b}}}^{v-l}\simeq -0.003$ to −0.007‰ at 3500 K and −0.016 to −0.031‰ at 1650 K $({{\rm{\Delta }}}_{{\rm{e}}{\rm{q}},{\rm{R}}{\rm{b}}}^{v-l}/{{\rm{\Delta }}}_{{\rm{e}}{\rm{q}},{\rm{K}}}^{v-l}\simeq 0.3)$.

If S_i = 0, the flux is unidirectional (from the liquid to the vapor) and the isotopic fractionation is mostly kinetic in nature (noted as ${{\rm{\Delta }}}_{\mathrm{kin},i}^{v-l}$ hereafter). Experiments of K evaporation in vacuum give a kinetic isotopic fractionation factor ${{\rm{\Delta }}}_{\mathrm{kin},{\rm{K}}}^{v-l}\,=[{(39/41)}^{0.43}-1]\times 1000=-22\unicode{x02030}$ (Richter et al. 2011). No data are available for Rb evaporation, but we can reasonably assume that it follows the same mass dependence', which gives ${{\rm{\Delta }}}_{\mathrm{kin},\mathrm{Rb}}^{v-l}=[{(85/87)}^{0.43}-1]\times 1000=-10\unicode{x02030}$.

If evaporation took place in a medium that was partially saturated (0 < S_i < 1), the instantaneous isotopic fractionation factor between vapor and liquid (${{\rm{\Delta }}}_{i}^{v-l}$) would be intermediate between the equilibrium and kinetic values given above (Dauphas et al. 2015):

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}^{v-l}={{\rm{\Delta }}}_{\mathrm{eq},i}^{v-l}+\left(1-{S}_{i}\right){{\rm{\Delta }}}_{\mathrm{kin},i}^{v-l}.\end{eqnarray} \tag{ 1 }$$

This expression gives the instantaneous isotopic fractionation between vapor and liquid. In the model of Charnoz & Michaut (2015) (scenario (iii) above), the vapor is continuously removed by accretion onto Earth through viscous spreading, leaving behind an MVE-depleted liquid layer. This is a complex process that to first order can be modeled using a Rayleigh distillation:

$$\begin{eqnarray}&&{\delta }_{i,l}-{\delta }_{i,0}\simeq \left[{{\rm{\Delta }}}_{\mathrm{eq},i}^{v-l}+\left(1-{S}_{i}\right){{\rm{\Delta }}}_{\mathrm{kin},i}^{v-l}\right]\mathrm{ln}{f}_{i},\end{eqnarray} \tag{ 2 }$$

where ${\delta }_{i,l}$ is the isotopic composition (in ‰) of the residual liquid after evaporation, δ_i,0 is the initial isotopic composition, and f_i is the fraction of the element remaining after vapor loss. In the context of lunar formation, the present isotopic composition of the Earth would be the starting composition (because in the Earth–Moon system, most of Rb and K resides in Earth), the lunar composition would be the residual liquid, and the fraction of either Rb or K remaining would be ∼0.17 (corresponding to a six-fold depletion; Ringwood & Kesson 1977; O'Neill 1991; Albarède et al. 2015). Wang et al. (2019) argued that the heavy K isotopic composition of the Moon relative to Earth could be explained through evaporative loss where each step in the distillation involves vapor–liquid equilibrium. Using ab initio data (Dauphas et al. 2018b; Zeng et al. 2018; Li et al. 2019; Zeng et al. 2019), we calculate that the δ⁴¹K value of the residue remaining after 83% loss of K by distillation under equilibrium conditions (S_i = 1 in Equation (2)) should be +0.021 to +0.042‰ at 3500 K (Canup et al. 2015; Lock et al. 2018) and +0.097 to +0.190‰ at 1650 K (Wang et al. 2019). These values are smaller than the observed Moon–Earth difference of +0.41 ± 0.07‰. We therefore disagree with the assessment of Wang et al. (2019) that the heavy K isotopic composition of the Moon (and by extension Rb and possibly Zn) can be explained by evaporative loss under equilibrium conditions. The data instead call for loss in an undersaturated medium involving kinetic isotopic fractionation.

Using Equation (2), we calculate the vapor saturation needed to explain the measured isotopic fractionations between lunar and terrestrial rocks for Rb and K. We find values of S_Rb = 0.991 and S_K = 0.990. If evaporation had occurred in a more undersaturated medium, kinetic isotope effects would have left the Moon more isotopically fractionated in Rb and K than measured. Conversely, if the evaporation had occurred in a medium closer to saturation, the instantaneous fractionation would approach the equilibrium value, leaving the Moon less isotopically fractionated than it is.

In Figure 2, we test if this calculated level of saturation could also explain the isotopic fractionations documented for other MVEs. Given the very high temperature involved, we can neglect the equilibrium term in Equation (2) and to a good approximation, we have

$$\begin{eqnarray}&&{\delta }_{i,l}-{\delta }_{i,0}\simeq \left(1-{S}_{i}\right){{\rm{\Delta }}}_{\mathrm{kin},i}^{v-l}\mathrm{ln}{f}_{i}.\end{eqnarray} \tag{ 3 }$$

If all MVEs experienced evaporation under the same saturation, we would expect to find a linear relationship between δ_Moon–δ_Earth and ${{\rm{\Delta }}}_{\mathrm{kin},i}^{v-l}\mathrm{ln}{f}_{i}$ whose slope is $1-S$. The MVEs Rb, Ga, Cu, K, and Zn indeed define a straight line corresponding to S = 0.989 ± 0.002 (Figure 2; weighted linear regression). Refractory elements such as Ca (Simon & DePaolo 2010) and Ti (Millet et al. 2016) are not shown in Figure 2, but they would plot at the origin (like Li, Mg, and Si) and would therefore follow the trend defined by MVEs. The only exception is Sn, which seems to be enriched in the lighter isotopes in lunar basalts compared to terrestrial rocks (Wang et al. 2019). This could be due to a large vapor–liquid equilibrium isotopic fractionation and a large evaporation coefficient for Sn (Section 4), such that the equilibrium term would dominate over the kinetic term in Equation (2) $\left(\left|{{\rm{\Delta }}}_{\mathrm{eq},i}^{v-l}\right|\gt \left|\left(1-{S}_{i}\right){{\rm{\Delta }}}_{\mathrm{kin},i}^{v-l}\right|\right)$. It could also reflect the fact that the lunar and terrestrial Sn isotopic compositions are not well known, as this element can exist in several oxidation states in planets (Sn⁰, Sn²⁺, and Sn⁴⁺), which can drive large isotopic fractionation (Badullovich et al. 2017; Dauphas et al. 2018a; Wang et al. 2018; Roskosz et al. 2019).

Zoom In Zoom Out Reset image size

The fact that elements showing great diversity in their volatilities (e.g., Zn is ∼200-fold depleted while K and Rb are ∼6-fold depleted; Ringwood & Kesson 1977; O'Neill 1991; Albarède et al. 2015) experienced evaporation under similar saturation S ≃ 0.99 is a fundamental observation and critical test for scenarios of lunar volatile depletion. As demonstrated below, this is a natural outcome of the model of volatile loss by viscous drainage of the vapor layer of the protolunar disk onto Earth. Still, further work is needed to evaluate whether other scenarios of Moon formation (Pahlevan & Stevenson 2007; Canup et al. 2015; Nakajima & Stevenson 2015; Saxena et al. 2017; Lock et al. 2018; Charnoz et al. 2019; Righter 2019; Young & Tang 2019) can similarly account for MVE loss by evaporation under a saturation of 99%.

The fact that the vapor was undersaturated means that there was a net evaporative flux from the liquid to the vapor layer, which must have been balanced by vapor removal. We consider a scenario (Figure 3) where the disk is composed of a liquid layer at the midplane of the protolunar disk, overlain by a vapor layer that is drained efficiently to the Earth due to its high viscosity, presumably powered by MRI (Charnoz & Michaut 2015; Carballido et al. 2016; Gammie et al. 2016). The timescale for producing a vapor layer with a hydrostatic profile for element i at saturation S_i is obtained by dividing the total mass of an element in the vapor layer by the flux across the liquid/vapor interface given by the Hertz–Knudsen equation (see the Appendix)

$$\begin{eqnarray}&&{t}_{S,i}=\displaystyle \frac{\pi {S}_{i}}{{\gamma }_{i}\left(1-{S}_{i}\right)}\sqrt{\displaystyle \frac{{m}_{i}{R}^{3}}{\overline{m}{{GM}}_{\oplus }}},\end{eqnarray} \tag{ 4 }$$

where γ_i is a dimensionless evaporation coefficient (available experimental results give values between 0.017 and 0.13 for evaporation of K from silicate melt; Fedkin et al. 2006; Richter et al. 2011); m_i and $\overline{m}$ are the mass of the vapor species for element i (e.g., ∼39 g mol⁻¹ for K because monoatomic K dominates the vapor) and the mean molecular mass of the bulk vapor (∼30 g mol⁻¹; Ward 2011; Charnoz & Michaut 2015), respectively; R is the distance to Earth's center (∼1.7 R_⊕ = 11 × 10⁶ m if we take the outer edge of the gravitationally stable liquid disk; Charnoz & Michaut 2015); and M_⊕ is Earth's mass. Importantly, this timescale is independent of the volatility (P_i,eq) of the element considered (and the uncertainties attached to it at the elevated temperature considered) because t_S,i is the ratio of the mass of an element in the vapor layer divided by the evaporation flux across the liquid/vapor boundary, which both scale linearly with P_i,eq. In our preferred scenario (Figure 3), vapor removal would occur by accretion onto the Earth through viscous spreading, which takes place on a timescale

$$\begin{eqnarray}&&{t}_{a}=\displaystyle \frac{2{R}^{2}}{3\nu },\end{eqnarray} \tag{ 5 }$$

where ν is the kinematic viscosity (in m² s⁻¹).

Zoom In Zoom Out Reset image size

Both timescales (t_S,i and t_a) are independent of volatility (P_i,eq), meaning that to first order, the level of saturation in the vapor layer should be similar for all elements. More volatile elements would still be more efficiently partitioned into the vapor layer and more efficiently lost than less volatile elements, explaining the different extents of depletion of MVEs measured in lunar rocks. Some variations may still exist in the saturation level due to differences in atomic/molecular mass (m_i) and more importantly evaporation coefficients (γ_i), so further experimental work is needed to better characterize the latter.

If t_a ≪ t_S,i, vapor is removed more rapidly than it can be replenished by vaporization from the liquid such that the saturation would be lower than S ≃ 0.99, and the isotopic fractionation would be larger than that measured in lunar rocks. Conversely, if t_a ≫ t_S,i, the saturation would be higher than calculated, and the isotopic fractionation would be smaller than measured in lunar rocks. Equating these two timescales, we can calculate the viscosity needed to maintain the vapor at the saturation level dictated by isotopes

$$\begin{eqnarray}&&\nu =\displaystyle \frac{2{\gamma }_{i}\left(1-{S}_{i}\right)}{3\pi {S}_{i}}\sqrt{\displaystyle \frac{\overline{m}{{RGM}}_{\oplus }}{{m}_{i}}},\end{eqnarray} \tag{ 6 }$$

or alternatively the saturation set for each element by a given viscosity (we use the approximation S_i ∼ 1)

$$\begin{eqnarray}&&{S}_{i}\simeq 1-\displaystyle \frac{3\pi \nu }{2{\gamma }_{i}}\sqrt{\displaystyle \frac{{m}_{i}/\overline{m}}{{{RGM}}_{\oplus }}}.\end{eqnarray} \tag{ 7 }$$

As shown in Figure 2, most MVEs record the same saturation level S ≃ 0.99. In the context of our model, the only element-specific controls on the saturation level are the mass of the vapor species ($\sqrt{{m}_{i}/\overline{m}}$) and the evaporation coefficient (γ_i). The $\sqrt{{m}_{i}/\overline{m}}$ term varies little from one element to another (1.14 for K, 1.46 for Cu, 1.48 for Zn, 1.52 for Ga, 1.69 for Rb) but the similarity in saturation also requires that the evaporation coefficients γ_i be similar (Equation (7)). As discussed in Section 3, the isotopic composition of Sn in lunar rocks seems to support evaporation under equilibrium conditions (S = 1; Wang et al. 2019). A possible explanation for this peculiar behavior is that the evaporation coefficient for Sn is much higher than for other MVEs (Equation (7)), which can be tested in the future by performing evaporation experiments.

To assess whether the calculated kinematic viscosity of the vapor layer (Equation (6)) is realistic or not, we put it in the context of an α-disk, where the viscosity is expressed as ν = αc²/Ω, with α a dimensionless number, $c=\sqrt{{kT}/\overline{m}}$ the sound speed (k is the Boltzmann constant), and ${\rm{\Omega }}=\sqrt{{{GM}}_{\oplus }/{R}^{3}}$ the Keplerian angular velocity. We calculate that α must have been on the order of 10⁻³–10⁻² in order to explain the inferred saturation of ∼0.99. Carballido et al. (2016) and Gammie et al. (2016) estimated that in the protolunar disk, where a fraction of alkalis could have been ionized, MRI could have sustained an α value on the order of 10⁻². Our estimate of α based on the isotopic compositions of Rb, K, and other volatile elements is thus entirely consistent with theoretical expectations, suggesting that MRI-powered viscous drainage is a viable mechanism for explaining the depletion of MVEs in the Moon. At the present time, we cannot tell whether such viscous drainage of vapor would have been a steady process occurring throughout the lifetime of the protolunar disk, or whether it was episodic.

Discussions with F. M. Richter, F. Ciesla, D. J. Stevenson, R. M. Canup, C. W. Visscher, S. Charnoz, S. J. Lock, J. Hu, H. Zeng, M. Meheut, M. Blanchard, J. Z. Zhang, S. M. Aarons, and A. W. Heard were greatly appreciated. We thank CAPTEM and R. A. Zeigler for providing the Apollo lunar samples analyzed in this study. This work was supported by a NASA NESSF fellowship (NNX15AQ97H) to N.X.N., NASA grants NNX17AE86G (LARS), NNX17AE87G (Emerging Worlds), and 80NSSC17K0744 (Habitable Worlds) to N.D.

A.1. Saturation Timescale

In this section, we calculate the timescale for replenishing a hydrostatic atmosphere above the disk midplane such that the surface pressure is at a set saturation S. The vertical forces acting on a parcel of gas are the pressure forces and the vertical component of the gravitational pull exerted by the Earth

$$\begin{eqnarray}&&\tfrac{1}{{\rho }_{g}}\tfrac{\partial P}{\partial z}=-\tfrac{{{GM}}_{\oplus }z}{{\left({R}^{2}+{z}^{2}\right)}^{3/2}},\end{eqnarray} \tag{ 8 }$$

where P is the gas pressure, z is the vertical distance from the midplane, ρ_g is the gas density, G is the gravitational constant, M_⊕ is the mass of the Earth, and R is the distance from Earth's center. The gas density ρ_g is related to the partial pressure P through the ideal gas law and we have

$$\begin{eqnarray}&&\tfrac{1}{P}\tfrac{\partial P}{\partial z}=-\tfrac{\overline{m}{{GM}}_{\oplus }z}{{{kT}\left({R}^{2}+{z}^{2}\right)}^{3/2}},\end{eqnarray} \tag{ 9 }$$

where $\overline{m}$ is the mean molecular mass of the vapor (∼30 g mol⁻¹, equivalent to ∼5 × 10⁻²⁶ kg molecule⁻¹; Ward 2011), k is Boltzmann's constant (1.38 × 10⁻²³ m² kg s⁻² K⁻¹), and T is the temperature (in K). Integration of this differential equation with $P\left(z=0\right)={P}_{0}$ yields

$$\begin{eqnarray}&&P={P}_{0}{e}^{-\tfrac{\overline{m}{{GM}}_{\oplus }}{{kTR}}\left(1-\tfrac{1}{\sqrt{1+{z}^{2}/{R}^{2}}}\right)}.\end{eqnarray} \tag{ 10 }$$

With z ≪ R (the scale height is much smaller than the distance to the Earth), we have

$$\begin{eqnarray}&&P={P}_{0}{e}^{-\tfrac{\overline{m}{{GM}}_{\oplus }}{2{{kTR}}^{3}}{z}^{2}},\end{eqnarray} \tag{ 11 }$$

which can be rewritten as

$$\begin{eqnarray}&&P={P}_{0}{e}^{-\tfrac{{z}^{2}}{2{h}^{2}}},\end{eqnarray} \tag{ 12 }$$

with $h=\sqrt{\tfrac{{{kTR}}^{3}}{\overline{m}{{GM}}_{\oplus }}}$ the scale height of the vapor layer. Assuming that the vapor layer is well mixed, the scale height of all elements is the same and the partial pressure of an element i as a function of altitude can be written as

$$\begin{eqnarray}&&{P}_{i}={P}_{i,0}{e}^{-\tfrac{{z}^{2}}{2{h}^{2}}}.\end{eqnarray} \tag{ 13 }$$

Introducing the sound speed $c=\sqrt{{kT}/\overline{m}}$ and the angular Keplerian velocity ${\rm{\Omega }}=\sqrt{{{GM}}_{\oplus }/{R}^{3}}$, the scale height can be rewritten as h = c/Ω. We now evaluate what the timescale t_S,i is for an element to be maintained at a given level of saturation S_i. The total surface density for an element (for the half-space; in kg m⁻²) is given by the integral

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{i}={\int }_{0}^{+\infty }\tfrac{{m}_{i}{P}_{i}\left(z\right)}{{kT}}{dz}=\tfrac{{m}_{i}{P}_{i,0}}{{kT}}{\int }_{0}^{+\infty }{e}^{-\tfrac{{z}^{2}}{2{h}^{2}}}{dz}=\sqrt{\tfrac{\pi }{2}}\tfrac{{{hm}}_{i}{P}_{i,0}}{{kT}},\end{eqnarray} \tag{ 14 }$$

where m_i is the mass of the vapor species for element i (in kg). The net mass flux of an element (in kg m⁻² s⁻¹) across the liquid/vapor interface is given by the Hertz–Knudsen equation (note that the more commonly used molar flux would have m_i at the denominator)

$$\begin{eqnarray}&&{\phi }_{i}={\gamma }_{i}\sqrt{\tfrac{{m}_{i}}{2\pi {kT}}}\left({P}_{i,\mathrm{eq}}-{P}_{i,0}\right),\end{eqnarray} \tag{ 15 }$$

where γ_i is an evaporation coefficient. The saturation timescale t_S,i is the total mass surface density of an element in the vapor layer at saturation S_i (Equation (14) with ${P}_{i,0}={S}_{i}{P}_{i,\mathrm{eq}}$) divided by the mass flux across the liquid/vapor layer at the same saturation given by the Hertz–Knudsen equation (Equation (15), again with ${P}_{i,0}={S}_{i}{P}_{i,\mathrm{eq}}$)

$$\begin{eqnarray}&&{t}_{S,i}=\tfrac{{{\rm{\Sigma }}}_{i}}{{\phi }_{i}}=\tfrac{\pi {S}_{i}}{{\gamma }_{i}\left(1-{S}_{i}\right)}\sqrt{\tfrac{{m}_{i}{R}^{3}}{\overline{m}{{GM}}_{\oplus }}}.\end{eqnarray} \tag{ 16 }$$

This saturation timescale is akin to a residence time.

A.2. Chemical Purification of Rubidium

A new Rb purification procedure was developed to separate Rb from rock matrices for isotopic analysis. The Rb blank of the procedure (digestion and column chemistry) is ∼0.14 ng, which accounts for less than 0.5% of total Rb from a typical sample (≥40 ng). The total Rb purification yields exceed 95% for all the samples.

The separation of Rb from matrix elements is difficult because it is a trace element and it behaves very similarly to K, which is usually in much higher concentration (the K/Rb weight ratio of terrestrial rocks is ∼325). Our procedure uses (i) AG50W-X8 cation exchange resin and HNO₃ to first separate Rb from matrix elements other than K, and (ii) Eichrom crown ether Sr resin (4,4'(5')-di-t-butylcyclohexano 18-crown-6 in 1-octanol) and HNO₃ to then separate Rb from K. Previous Rb isotope studies (Nebel et al. 2011; Pringle & Moynier 2017) used mainly cation exchange resin and HCl to separate Rb from K and matrix elements. The second step in our procedure is similar to the chemistry used by Zhang et al. (2018), who used Sr resin but in a shorter column to separate Rb from K.

Samples of about 100 mg or less were digested in three steps using mixtures of concentrated HF–HNO₃–HCl–HClO₄ acids:

• (i)  
  4 ml 28 M HF + 2 ml 15 M HNO₃ + 1 ml 10 M HClO₄ was added to each sample in a 30 ml fluoropolymer beaker that was closed and left on a hotplate at 130 °C for 24 hr.
• (ii)  
  The solution was dried down at 130 °C on a hotplate, taken up in 4.5 ml 11 M HCl + 1.5 ml 15 M HNO₃, and left on a hotplate at 130 °C for 24 hr with the lid closed. This step was repeated twice.
• (iii)  
  After evaporation to dryness, the sample was redissolved in 4 ml 15 M HNO₃, and left on a hotplate for 24 hr with the lid closed. If a sample was not fully digested, the residue was transferred to a fluoropolymer beaker to be further digested in a Parr bomb at 175–180°C for at least 3 days and it was then combined with the previously digested solution.

The combined solution was dried and redissolved in 1 M HNO₃ for column chemistry, which consists of five columns:

• (i)  
  The first column (Figure 4(a)) uses 16 ml cation resin (AG50W-X8, 200–400 mesh resin in 20 ml Bio-Rad Econo-Pac columns of 15 mm diameter and 140 mm length). The samples were loaded in 4 ml 1 M HNO₃, and 160 ml of 1 M HNO₃ was then passed through the column to prepurify Rb from most matrix elements (Strelow 1960; Strelow et al. 1965; Dybczynski 1972). This step removes elements with a higher distribution coefficient than Rb that are retained on the column (e.g., Fe, Mg, Al, Ca, Cr, and Mn) while Rb is eluted. All the matrix elements with distribution coefficients lower than, or similar to, Rb are collected together with Rb (1–160 ml; light-blue shaded area in Figure 4(a)). The elution scheme was designed to process relatively large sample masses (typically 100 mg). We collect all elution before the Rb peak to avoid Rb loss as excessive matrix element concentrations could potentially cause Rb to be eluted early by competing for ion exchange sites on the resin. This step removes a significant portion of matrix elements and all Rb is recovered.
• (ii)  
  The second column (Figure 4(b)) uses the same type of resin (16 ml AG50W-X8 200–400 mesh in 20 ml Bio-Rad Econo-Pac columns of 15 mm diameter and 140 mm length) but a different HNO₃ molarity of 0.5 M. The reason for using this lower acid molarity is that the elements are better separated from each other (Strelow 1960; Strelow et al. 1965), the downside being a larger elution volume compared to the first column (the Rb peaks are at 280 versus 120 ml, respectively; Figure 4). Samples were loaded in 4 ml of 0.5 M HNO₃, followed by 360 ml 0.5 M HNO₃ for elution. Rubidium is collected in 130–360 ml elution volume (light-blue shaded area in Figure 4(b)). Elements that have a higher or lower partition coefficient than Rb are eliminated. Most Ti is removed in this step, while most K is kept in the Rb elution fraction.
• (iii)  
  After step (ii), a small amount of Ti may remain with Rb. The quantitative removal of remaining Ti was achieved by using a small anion column (1 ml resin of AG1-X8 200–400 mesh in a Bio-Rad Poly-Prep 10 ml column of 8 mm diameter and 90 mm length). The samples were loaded on the resin in 0.5 ml of 2 M HF and then eluted with 10 ml of the same acid. In 2 M HF, Ti sticks to the anion exchange resin while Rb is eluted (Faris 1960; Nelson et al. 1960).
• (iv)  
  This step separates Rb from K. We use a custom-made fluoropolymer column (40 cm in length and 0.4 cm in diameter) filled with Eichrom Sr resin (50–100 μm) to separate Rb from K. Element partitioning data for the Eichrom Sr resin can be found in Philip Horwitz et al. (1992). Samples were loaded and eluted in 3 M HNO₃. As shown in Figure 4(c), complete separation of K and Rb was achieved. This step is similar to the one reported by Zhang et al. (2018) but we use a much longer column (40 cm versus 13 cm) to better separate Rb from K.
• (v)  
  The solutions from the previous Sr resin step contained pure Rb but were relatively viscous, likely due to dissolution of small amounts of organics from the Sr resin. An additional small clean-up column (1 ml AG50W-X8 resin in a Bio-Rad Poly-Prep 10 ml column of 90 mm height and 8 mm diameter) was therefore used. Samples were loaded onto the column in 1 ml of 0.5 M HNO₃, followed by another 1 ml 0.5 M HNO₃, and Rb was then collected in 10 ml 6 M HNO₃. This step is not necessary for all Rb samples, especially for samples with high Rb contents. Our tests showed that solutions before and after this column gave the same Rb isotopic compositions when measured by MC-ICPMS.

Zoom In Zoom Out Reset image size

Pringle & Moynier (2017) used 0.5 M HCl and a 1 ml homemade column (0.6 cm inner diameter and 3.5 cm length) filled with AG50W-X8 resin to separate Rb from K. They claimed that this column could separate K from Rb but did not show any elution curve. We tested the method, and as shown in Figure 5, the elution peaks of K and Rb largely overlap. In contrast, our method (Figure 4(c)) allows a complete separation of Rb from K.

Zoom In Zoom Out Reset image size

All purified Rb solutions were checked for the presence of potential matrix elements (the elements shown in Figure 4) before measuring their Rb isotopic compositions. After the procedure, K was reduced to about the same concentration level as Rb or lower, and all other matrix elements had lower concentrations than Rb. The remaining K did not affect the Rb isotopic analyses. Our tests with K-doped Rb solutions showed that the measured isotopic composition of Rb was shifted only when K/Rb concentration ratios (ppm/ppm) were >50 (Figure 6).

Zoom In Zoom Out Reset image size

Even small quantities of Sr in the Rb solution can produce a significant isobaric interference of ⁸⁷Sr⁺ on ⁸⁷Rb⁺. The contribution of ⁸⁷Sr was corrected for by monitoring ⁸⁸Sr and assuming a constant ⁸⁷Sr/⁸⁸Sr ratio of 0.085. Our tests with Sr-doped Rb solutions showed that solutions with ⁸⁸Sr/⁸⁵Rb intensity ratios (V/V) of ≤0.001 yielded accurate δ⁸⁷Rb isotopic compositions within ±0.03‰ (Figure 7). The ⁸⁸Sr/⁸⁵Rb intensity ratios (V/V) of our purified samples were all <0.0005 (for most samples the ratio was much lower), which is 10 times lower than the ratio after purification in Pringle & Moynier (2017) (⁸⁸Sr/⁸⁵Rb < 0.005).

Zoom In Zoom Out Reset image size

A.3. Mass Spectrometry

The Rb isotopic analyses were performed using the Thermo Scientific Neptune multi-collector inductively coupled plasma mass spectrometer (MC-ICPMS) at the University of Chicago. A standard-sample bracketing technique was used. Rubidium sample and standard solutions of ∼15–25 ppb (∼1–1.5 V for ⁸⁵Rb in low resolution) in 0.3 M HNO₃ medium were introduced into the MC-ICPMS at a flow rate of 100 μl minute⁻¹ via a dual cyclonic-Scott-type quartz spray chamber. The samples and standards were matched for Rb concentrations within 1% for MC-ICPMS analyses. Normal Ni sampler and skimmer cones were used, and all measurements were performed in low-resolution mode. The MC-ICPMS at the University of Chicago is equipped with nine Faraday collectors. Rubidium-85 and -87 were measured on L2 and axial (A) Faraday collectors, respectively, and the isobaric interference of Sr was monitored at mass ⁸⁸Sr on H1. All three collectors were equipped with the 10¹¹ Ω amplifiers. Isobaric interference ⁸⁷Sr was corrected for by assuming a constant ⁸⁷Sr/⁸⁸Sr ratio of 0.085, which is the terrestrial Sr ratio. For our purified samples, the correction of the ⁸⁷Sr interference on ⁸⁷Rb is small enough that uncertainties in the ⁸⁷Sr/⁸⁸Sr ratio are largely inconsequential (Figure 7). For example, adopting two extreme ⁸⁷Sr/⁸⁸Sr values covering documented variations among chondrites, terrestrial, and lunar rocks (0.0835 and 0.0885, which correspond to ⁸⁷Sr/⁸⁶Sr ratios of 0.70 and 0.74, respectively) would shift δ⁸⁷Rb values of the most Sr-rich purified solutions (⁸⁸Sr/⁸⁵Rb intensity ratio of ∼0.0005) by 0.008‰ at most.

Data were collected as a single block of 25 cycles of 4.194 s integration time, with a take-up time of 90 s and a washout time (0.45 M HNO₃ for wash) of 60 s. The clean acid solution (0.3 M HNO₃) that was used for diluting Rb for isotopic analysis was measured before and after each sample and standard under conditions identical to the sample or the standard (with the same integration, take-up, and washout times). For each sample and standard, the average intensity of the two bracketing blank acids (0.3 M HNO₃) was subtracted. The Rb background was normally around 0.001–0.003 V for ⁸⁵Rb and increased by 0.001 V after a measurement session (typically ∼15 hr or less).

The sample solutions were measured 5–12 times, depending on the total Rb amount, and the average δ⁸⁷Rb values were calculated. The uncertainty for a sample was calculated using the formula of 2 × σ/$\sqrt{n}$, with n the number of replicates for the sample and σ the standard deviation of the δ⁸⁷Rb values of the standards calculated by treating the standards as if they were samples and calculating the δ⁸⁷Rb values by considering the two nearest standards bracketing each standard (Dauphas et al. 2009). All Rb isotopic compositions are expressed in δ⁸⁷Rb notation, which is the departure in ‰ of the ⁸⁷Rb/⁸⁵Rb ratio of a sample from that of the reference material NIST SRM984,

$$\begin{eqnarray*}\begin{array}{rcl}{\delta }^{87}\mathrm{Rb}(&permil; ) & = & [{\left({}^{87}\mathrm{Rb}{/}^{85}\mathrm{Rb}\right)}_{\mathrm{sample}}/{\left({}^{87}\mathrm{Rb}{/}^{85}\mathrm{Rb}\right)}_{\mathrm{SRM}984}-1]\\ & & \times \,1000.\end{array}\end{eqnarray*}$$

A.4. Test Results

The yield and accuracy of the procedure were checked with the reference standard NIST SRM984 treated as a sample, and with various synthetic and natural samples. Synthetic peridotite samples were made by mixing geostandard (peridotite) powders containing very little Rb with the Rb reference standard SRM984. The matrix elements in these synthetic samples are found in natural peridotites while their Rb isotopic compositions should be identical or very similar to SRM984 (i.e., δ⁸⁷Rb = 0). Two synthetic samples were made: DTS-2b (dunite) + SRM984 and PCC-1 (peridotite) + SRM984. As shown in Figure 8, the two synthetic Rb-doped peridotite samples, and SRM984 treated as sample, gave δ⁸⁷Rb values of zero within error.

Zoom In Zoom Out Reset image size

We also measured several geostandards, including basalts and granites (BHVO-2, BCR-2, BE-N, W-2, AGV-2, GSR-1, GS-N, G-A, and G-3). The majority of samples were digested and measured more than once and yielded reproducible results, which are in good agreement with previously published results (Pringle & Moynier 2017; Zhang et al. 2018; Figure 8). The Allende carbonaceous chondrite (powder from the Smithsonian Institution) was also measured to test the method on a low-Rb sample. Compared with the geostandards, Allende has much lower Rb concentration (∼1 ppm versus tens of ppm or more in the geostandards). Its Rb concentration is, however, similar to mare basalts (∼1 ppm). It was digested and measured three times and all the three measurements yielded reproducible results (Figure 8). Our measured Rb isotopic compositions of lunar samples (Table 1) in general agree with the data reported in Pringle & Moynier (2017) (Table 2).