# Hydrodynamical simulations of proto-Moon degassing

G. Madeira, L. Esteves, S. Charnoz, E. Lega, F. Moynier

► To cite this version:

G. Madeira, L. Esteves, S. Charnoz, E. Lega, F. Moynier. Hydrodynamical simulations of proto-Moon degassing. Earth and Planetary Science Letters, 2025, 651, 10.1016/j.epsl.2024.119163. insu-05080934

HAL Id: insu-05080934

https://insu.hal.science/insu-05080934v1

Submitted on 23 May 2025

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>

# Hydrodynamical simulations of proto-Moon degassing

G. Madeira<sup>a, ID,*</sup>, L. Esteves<sup>a,c, ID</sup>, S. Charnoz<sup>a</sup>, E. Lega<sup>b, ID</sup>, F. Moynier<sup>a</sup>

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

## ARTICLE INFO

Editor: O. Mousis

**Keywords:**

Moon, surface
Satellites, composition
Satellites, surfaces

## ABSTRACT

Modeling the isotopic and elemental abundance of the bulk silicate Moon represent major challenges. Similarities in the non-mass dependent isotopic composition of refractory elements with the bulk silicate Earth suggest that both the Earth and the Moon formed from the same material reservoir. On the other hand, the Moon's volatile depletion and isotopic composition of moderately volatile elements points to a global devolatilization processes, most likely during a magma ocean phase of the Moon. Here, we investigate the devolatilization of the molten Moon due to a tidally-assisted hydrodynamic escape, first proposed by Charnoz et al. (2021), with a focus on the dynamics of the evaporated gas. Unlike the 1D steady-state approach of Charnoz et al. (2021), we use 2D time-dependent hydrodynamic simulations carried out with the FARGOCA code modified to take into account the magma ocean as a gas source. Near the Earth's Roche limit, where the proto-Moon likely formed, evaporated gases from the lunar magma ocean form a circum-Earth disk of volatiles, with less than 30% of material being re-accreted by the Moon. We find that the measured depletion of K and Na on the Moon can be achieved if the lunar magma-ocean had a surface temperature of about 1800-2000 K. After about 1000 years, a thermal boundary layer or a flotation crust forms a lid that inhibits volatile escape. Mapping the volatile velocity field reveals varying trends in the longitudes of volatile reaccretion on the Moon's surface: material is predominantly re-accreted on the trailing side when the Moon-Earth distance exceeds 3.5 Earth radii. For $k_2/Q$ values of 0.0003 and 0.03, 60% and more than 99% of the volatile material, respectively, is re-accreted on the trailing side, suggesting a dichotomy in volatile abundances between the leading and trailing sides of the Moon. This dichotomy may provide insights on the tidal conditions of the early molten Earth. In conclusion, tidally-driven atmospheric escape effectively devolatilizes the Moon, matching the measured abundances of Na and K on timescales compatible with the formation of a thermal boundary layer or an anorthite flotation crust.

## 1. Introduction

The similarities in the isotopic composition of oxygen (Wiechert et al., 2001; Herwartz et al., 2014) and non-volatile elements such as Cr and Ti (Zhang et al., 2012; Mougel et al., 2018) between the bulk silicate Earth and the Moon indicate that both Earth's mantle and the bulk silicate Moon originated from material with a similar initial isotopic composition. The Moon's formation scenario is still debated: The 'canonical scenario', involving an impact with a Mars-sized object (Benz et al., 1989; Kokubo et al., 2000; Canup and Asphaug, 2001; Canup, 2008), successfully explains the anomalously high angular momentum of the Earth-Moon system and the Moon's bulk iron depletion. However, it reproduces the compositional similarities between the Earth and the

Moon only if the Earth and the impactor shared nearly identical compositions. This limitation arises because the simulations of the canonical scenario obtain a proto-lunar disk predominantly composed of impact material. Nevertheless, the formation of planet-sized objects with the required compositional similarities is highly improbable (Kaib and Cowan, 2015a,b).

High-energy, high-angular-momentum impacts (Ćuk and Stewart, 2012; Canup, 2012; Lock et al., 2018, 2020) produce a higher-temperature proto-lunar disk, where the homogenization of the proto-Earth and the impactor material is facilitated. However, the substantial vapor mass above the proto-lunar disk produced by these impacts can, in fact, inhibit the formation of a satellite as massive as the Moon (Nakajima et al., 2024), which has been shown to be a caveat of these high-energy models. Additionally, achieving the level of material mix-

<footnote><sup>a</sup> Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, Paris, F-75005, France</footnote>

<footnote><sup>b</sup> Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Lagrange (UMR7293), Nice cedex 4, F-06304, France</footnote>

<footnote><sup>c</sup> São Paulo State University, Grupo de Dinâmica Orbital e Planetologia, Guaratinguetá, SP, 12516-410, Brazil</footnote>

<footnote>* Corresponding author.
E-mail address: madeira@ipgp.fr (G. Madeira).</footnote>
ing necessary to explain the Earth-Moon system requires a significant excess of angular momentum (Meier et al., 2024).

Ćuk and Stewart (2012) propose that an ejection resonance with the Sun might have removed angular momentum from the post-impact system; however, it remains uncertain whether this mechanism alone can sufficiently account for the removal of all the angular momentum excess (Rufu and Canup, 2020; Ward et al., 2020). Regardless of the scenario, the physics governing the debris disk, the recondensation and reaccretion of evaporated silicates into condensates, and their coupled dynamics remain poorly understood due to the intrinsic complexity of the problem, leaving the Moon's formation as an open question.

In contrast to the similar abundance of refractory elements between the Earth and the Moon, the abundances of moderately volatile elements (e.g., Na, K, Zn, Rb, Cd) in lunar mare basalts exhibit a notable depletion when compared to terrestrial basalts (Ringwood, 1992; Lodders and Fegley Jr., 2015; Day and Moynier, 2014). Furthermore, it has been measured that volatile elements, including Zn (Paniello et al., 2012; Kato et al., 2015; Day et al., 2022), K (Wang and Jacobsen, 2016; Tian et al., 2020), Cl (Sharp et al., 2007; Boyce et al., 2015; Gargano et al., 2020), Ga (Kato and Moynier, 2017), and Rb (Pringle and Moynier, 2017), are enriched in the heavier isotopes in the Moon compared to the Earth, suggesting that the depletion of volatile elements in the Moon is attributable to global evaporation processes.

High-precision analysis of lunar magmatic rocks (Sossi et al., 2018) reveals an enrichment of lighter Cr isotopes compared to Earth's mantle, suggesting that Cr partitioning occurred near equilibrium in an oxygen-rich environment, at relatively low temperatures of around 1600–1800 K. Sossi et al. (2018) argue that these temperatures are similar to those expected from degassing at the surface of a magma ocean, indicating that volatile loss likely occurred during the magma ocean phase. However, we point out that it is also possible this process took place right after the giant impact, after the gas expanded and cooled to low temperatures. Nevertheless, due to the absence of hydrogen and the Earth's large mass, it seems unlikely that moderately volatile elements were lost through thermal escape (Nakajima and Stevenson, 2018). This is why several alternative scenarios have been proposed.

In the protolunar disk, volatile loss could happen due to incomplete condensation and liquid-vapor partitioning, caused by the viscous separation of the vapor from the liquid (Charnoz and Michaut, 2015) or by the gravitational torque exerted on the vapor by the proto-Moon (Canup et al., 2015). Additionally, Lock et al. (2018, 2020) propose that separation due to the gas drag of the growing droplets from the vapor may also play a role. However, the effectiveness of these scenarios depends heavily on the structure of the protolunar disk, which remains highly unknown. In fact, the existence of a protolunar disk is even debated. For example, Kegerreis et al. (2022) demonstrates through high-resolution simulations that the Moon may have formed immediately from a high-energy oblique impact of an object with the proto-Earth (also see Meier et al., 2024).

A plausible explanation is that the volatile loss occurred later in the system, after the Moon underwent wholesale melting following a high-energy accretion process and became covered by an ocean of magma (Wood et al., 1970). Isotopic evidence supporting this scenario comes from the variable Zn and Cl isotopic composition between lunar lithologies (Kato et al., 2015; Boyce et al., 2015). However, while isotopic data points to volatile loss during magma ocean degassing, the Moon's escape velocity is too large and would prevent the loss of moderately volatile elements.

To solve this conundrum Charnoz et al. (2021) developed a new model in which the volatile loss occurred immediately after the assembly of the proto-Moon, when the object was much closer to Earth. The influence of Earth's tides reduces the velocity needed for a gas to cross the Hill sphere and thermally escape from the proto-Moon. Exploring this scenario, Charnoz et al. (2021) find that between 3 and 6 Earth radii, the proto-Moon undergoes intense tidal devolatilization, and volatile material does not re-enter the satellite's Hill sphere when

under perturbative effects such as gas drag or radiation pressure. For lunar surface temperatures of approximately 1600-1800 K, Charnoz et al. (2021) demonstrated that lunar depletion in Na and K can be reproduced, provided the escape persists for at least 1000 years.

However, Charnoz et al. (2021) does not fully evaluate the problem by assuming a 1D steady-state approach. They solve the classical equations of hydrodynamic escape (Parker, 1963, 1965) only for the line connecting the Earth and the Moon, assessing whether the gas on the surface of the proto-Moon is hot enough to accelerate upwards to the Lagrange points, which are assumed to be situated on the Earth-Moon line. Consequently, their calculations of fluxes account for material released from only a limited region of the lunar surface. Moreover, they fail to systematically demonstrate that the material does not return to the Moon, relying instead on simulations with test particles that do not consider the hydrodynamic nature of the gas evolution. In fact, as we will demonstrate below, the presence of gas in the vicinity of the Moon acts to displace the Lagrange points away from the Earth-Moon line, challenging an important assumption of their model.

In this work, we take a step further than Charnoz et al. (2021) and investigate the tidally-driven atmospheric escape scenario using 2D time-dependent hydrodynamic simulations. This approach allows us to evaluate the evolution of gas released from the Moon's equator, quantifying both the amount of material that returns to the satellite and the locations on its surface where the gas is re-accreted. Additionally, we quantify the flow of material from the Moon to the Earth. In Section 2, we present our hydrodynamic set-up and our model for the Moon's magma ocean. In Section 3, we numerically simulate the evolution of gas released from the Moon and calculate the material fluxes. In Section 4, we examine the Moon's orbital evolution due to tides, calculating how the abundances of Na and K in the satellite change over time. Then, in Section 5, we explore how the formation of a stagnant lid can influence the total amount of material lost by the Moon. We discuss our results in Section 6 and draw our main conclusions in Section 7.

## 2. Methods

Our 2D numerical simulations were conducted using the FARGOCA code (Lega et al., 2014, FARGO with Co-latitude Added).<fnref n="1" /> FARGOCA is a parallelized, grid-based hydrodynamic code originally designed for studying protoplanetary disks. It is based on the FARGO code (Masset, 2000) and solves the fluid equations using a second-order upwind scheme within a time-explicit-implicit multi-step procedure.

We use polar coordinates $(r, \theta)$, where $r$ is the radial distance from the Earth (which is at the origin of the coordinate system) and $\theta$ is the azimuthal coordinate measured from the x-axis. We work in a coordinate system that rotates with the angular velocity of the Moon:

$$ \Omega_M = \sqrt{\frac{G(M_E + M_M)}{a_M^3}} \quad (1) $$

where $M_E$ is the mass of the Earth, $M_M$ is the current mass of the Moon, $G$ is the gravitational constant, and $a_M$ is the semi-major axis of the Moon. The Moon is assumed to be in a circular and planar spin-tidally locked orbit.

The dynamics of the gas are described by the integration of the Navier-Stokes equations, composed of a continuity equation and a set of two momentum equations:

### Continuity equation

The continuity equation is given by:

$$ \frac{\partial \Sigma}{\partial t} + \nabla \cdot (\Sigma \bar{v}) = 0 \quad (2) $$

<footnote n="1">¹ The simulations presented in this paper were obtained with a recently refactored version of the code available at https://gitlab.oca.eu/DISC/fargOCA.</footnote>
where Σ is the surface density of the gas and $\vec{v} = (v_r, v_\theta)$ is the velocity vector of the gas flow, with $v_\theta = r(\omega + \Omega_M)$ where $\omega$ is the azimuthal angular velocity in the rotating frame.

## Equations for the momenta

The Navier-Stokes equations for the radial momentum $J_r = \Sigma v_r$, and the angular momentum $J_\theta = \Sigma r v_\theta$ read:

$$ \begin{cases} \frac{\partial J_r}{\partial t} + \nabla \cdot (J_r \vec{v}) = \Sigma \left[ \frac{v_r^2}{r} - \frac{\partial \Phi}{\partial r} + \frac{1}{\Sigma} \left( f_r - \frac{\partial P_{2D}}{\partial r} \right) \right] \\ \frac{\partial J_\theta}{\partial t} + \nabla \cdot (J_\theta \vec{v}) = \Sigma r \left[ -\frac{\partial \Phi}{\partial \theta} + \frac{1}{\Sigma} \left( f_\theta - \frac{\partial P_{2D}}{\partial \theta} \right) \right] \end{cases} \quad (3) $$

where $P_{2D}$ is the vertically integrated pressure. The function $f = (f_r, f_\theta)$, which is the divergence of the stress tensor, is proportional to the disk's viscosity. In this study, we consider the fluid inviscid based on the results of previous one-dimensional calculations (Charnoz et al., 2021). Although an inviscid fluid is an appropriate approximation for low-density gases, as is the case for our vapor gas (Section 2.3), we point out that at the magma ocean interface, the gas density, and consequently its viscosity, may reach higher values, representing a caveat of our model.

The potential $\Phi$ acting on the disk consists of the contribution of the Earth $\Phi_E = -GM_E/r$ and the Moon $\Phi_M$, plus indirect terms that arise from the primary acceleration due to the Moon's and disk's gravity (Masset, 2002). The gravitational potential of the Moon on the disk is modeled as:

$$ \Phi_M = - \frac{GM_M}{\sqrt{r^2 + R_{sm}^2}} \quad (4) $$

where $R_{sm}$ is the smoothing length with $R_{sm} = 0.6H$, and $H = hr$ is the disk scale height, $h$ being the disk's aspect ratio (whose value is given in Section 2.1).

## Energy equation

To the Navier-Stokes system of equations, we add the following equation for the temporal variation of the 2D internal energy density:

$$ \frac{\partial e}{\partial t} + \nabla \cdot (e \vec{v}) = -P_{2D} \nabla \cdot \vec{v}, \quad (5) $$

where the term $-P_{2D} \nabla \cdot \vec{v}$ is the compressional heating, and the internal energy density is given by

$$ e = \Sigma c_v T, \quad (6) $$

with $T$ being the gas temperature and $c_v$ the specific heat at constant volume.

The vertically integrated pressure relates to the 2D internal energy density through the equation

$$ P_{2D} = (\gamma - 1)e, \quad (7) $$

where $\gamma$ is the adiabatic index.

In this work, our aim is to explore the dynamics of gas released from the surface of the Moon, deviating slightly from the scenarios for which the code was originally intended. Specifically, the FARGOCA code employs a point-mass approach for calculating the gravitational force of planets – or satellites, in our case – which is not appropriate for our simulations that require a physically resolved surface for the Moon. For this reason, we implemented modifications in the code, described below.

## 2.1. Implementations

We implemented the extended dimensions of the Moon in the code using boundary conditions. For each integration time-step, we calculate the distance $d_{i,M}$ from each grid cell $i$ to the Moon's center of mass. Cells within a given sphere contained by the Moon are assumed to have the

same properties as the Moon's center. Meanwhile, an evanescent boundary condition (de Val-Borro et al., 2006) is applied in the narrow shell between this sphere and the Moon's surface. This boundary condition is mathematically defined as:

$$ \Xi = \begin{cases} \Xi_M & \text{for } d_{i,M} \le R_{\text{shell}}, \\ \frac{\Xi + \lambda \Xi_M}{1 + \lambda} & \text{for } R_{\text{shell}} \le d_{i,M} \le R_M \end{cases} \quad (8) $$

where $\Xi$ represents any of the following quantities: velocity components, surface density, and energy. $\Xi_M$ is the quantity at the Moon's center, and $\lambda$, as defined by de Val-Borro et al. (2006), is the evanescent parameter responsible for controlling the rate of damping. The parameter relates to the integration time-step $T_{\text{step}}$ through the equation:

$$ \lambda = \frac{20T_{\text{step}}}{T_M} \left( \frac{R_M - d_{i,M}}{R_M - R_{\text{shell}}} \right)^2, \quad (9) $$

where $T_M$ is the Moon's orbital period.

The evanescent boundary condition applied to the Moon's surface is based on the approach originally used for the edges of the disk in the FARGOCA code. It is explicitly applied to the velocity, energy, and surface density of the gas. At the Moon's center, the radial and azimuthal components of the velocity in the rotating frame are set to $v_{r,M} = v_{\theta,M} = 0$, and the surface density is computed as (see Equations (6)-(7)):

$$ \Sigma_M = \frac{P_{2D}}{T_M} \frac{1}{c_v(\gamma - 1)}, \quad (10) $$

where $T_M$ is the temperature at the Moon's center. The gas energy at the Moon's center is computed using Equation (6). It should be noted that the condition $T = T_M$ at the Moon's center is applied implicitly in our code.

This approach ensures that all gas material within the sphere behaves as a single entity while promoting a smooth transition between the regions inside and outside the sphere. Simultaneously, if the surface density at the sphere, $\Sigma_M$, is higher than the disk's local surface density $\Sigma$, we naturally observe a flow of material from the Moon's surface to the disk. We interpret this as the vapor flux produced in the magma ocean.

## 2.2. Numerical setup-up

Our numerical simulations include the Earth as the central object, and the Moon, modeled following the implementation described in Section 2.1. Our calculations are performed in normalized units, assuming the Earth mass $M_E$ and radius $R_E$ as units of mass and distance, respectively ($G = M_E = R_E = 1$). We model the environment around the Earth using 1000 grid cells with arithmetic spacing in radius ($\delta r = 0.011$), distributed from 1 to 12 $R_E$, and 1000 grid cells in azimuth for the full circumference ($\delta \phi = 0.001$). This configuration results in a total of 10224 cells representing the Moon.

To ensure accuracy, we ran test simulations spanning 0.1 yr, measuring the fluxes of material ejected, re-accreted by the Moon, and lost at the disk boundaries. The outer disk boundary was increased from 9 $R_E$ in steps of 1 $R_E$, and grid resolution was refined from 500 to 1000 cells in both radial and azimuthal directions, in steps of 100. Fluxes converged to asymptotic values with increasing precision. We adopted a disk edge of 12 $R_E$ and 1000 grid cells, as these configurations resulted in relative flux differences below 10% compared to the immediately less precise cases. We also performed a set of test simulations varying the values of $R_{\text{shell}}$ to check the convergence of results with respect to this parameter, in the case of the Moon positioned at 3 $R_E$. These tests revealed that the results consistently converged for a minimum value of 0.95 $R_M$, which we adopted as the standard $R_{\text{shell}}$ in our study. The disk's inner and outer boundaries were treated as open.

Ideally, our simulations should begin with only the Earth and Moon present, to accurately analyze the flow of volatiles in the system. However, due to a code limitation, grid cells cannot be initialized as empty.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Partial pressure in the vapor above the lunar magma ocean as a function of lunar surface temperature, assuming a bulk silicate Earth (BSE) composition. Calculation was performed using the VAPOROCK code. Ptot stands for the total saturating vapor pressure.</figcaption></figure>

For simplicity, we start our simulations with an inviscid, adiabatic disk that has a constant surface density of $\Sigma = 1.6 \times 10^{-15}$ g/cm$^2$. This density is close to the minimum the code can handle and corresponds to a total mass of less than $< 10^{-22} M_E$. Moreover, this value is at least four orders of magnitude smaller than the $\Sigma_M$ values assumed in our simulations. The disk is characterized by a constant aspect ratio $h$, computed in normalized units as:

$$h = \sqrt{\frac{3RT_{BG}}{\mu}} \qquad (11)$$

where R is the gas universal constant, $\mu$ is the gas mean molecular mass, and $T_{BG}$ is the disk background temperature, assumed to be 300 K at $3 R_E$. Next, we detail our assumptions for the gas vapor released from the Moon's surface.

## 2.3. Magma ocean composition

We assume that the vapor above the lunar magma ocean consists of gas evaporated congruently by the oxides present in a liquid with a bulk silicate Earth (BSE) composition. The vapor composition is computed using the recently published VAPOROCK code (Wolf et al., 2023), that is an evolution of the MELT code. As in Visscher and Fegley (2013); Charnoz et al. (2021); van Buchem et al. (2023), the fraction of $O_2$ in the vapor is computed by assuming mass conservation of oxygen between the magma ocean and the vapor during the vaporization process.

Fig. 1 displays the resulting vapor composition in terms of partial pressures. Our bulk silicate Earth (BSE) vapor composition is very similar to that obtained with the LAVATMOS code (see Figure 4 of van Buchem et al., 2023) and the one published in Ito et al. (2015). However, significant differences appear when compared to the composition reported in Charnoz et al. (2021): our Na and K partial pressures computed with VAPOROCK are respectively about 5 and 10 times smaller at 2000 K than those in Charnoz et al. (2021), due to different models for activity coefficients. Note that, whereas our VAPOROCK calculation is in good agreement with LAVATMOS and Ito et al. (2015), there are still notable differences with the vapor composition obtained with the MAGMA code (see Figure 4 of van Buchem et al., 2023), which reports about an order of magnitude less Na and two orders of magnitude less K abundances at 2000 K than our calculations, probably due to different thermo-chemical data.

Table 1 shows the total pressure, the adiabatic index, the mean molecular mass, and atomic mass fractions of the gas above the magma

ocean as functions of the lunar surface temperature. Aluminum (Al) and Calcium (Ca) are also included in our gas composition; however, their mass fraction in the vapor are always less than 0.001 for all temperatures considered. In this work, we conducted a total of 65 numerical simulations, each with a timespan of 1 year, considering the magma ocean temperature $T_M$ and the Moon's orbital distance $a_M$ as free parameters. We assume $T_M = 1400, 1600, 1800, 2000$, and $2200$ K and $a_M$ ranging from 3 to 9 $R_E$, with a semi-major axis step of 0.5 $R_E$. As the Moon loses volatiles, the atomic concentrations in the magma ocean are expected to change over time, affecting the pressure, adiabatic index, and mean molecular mass of the vapor above the magma. For simplicity, however, we treat these parameters as constants throughout the simulations, determined solely by the gas temperature (Table 1), which is also kept constant.

## 3. Circum-Earth volatile disk and its dependency on lunar position and surface temperature

In this section, we explore the evolution and fate of the gas vapor released from the lunar surface. Initially, we focus on the case with the Moon at $3 R_E$, position at which lunar aggregates are expected to form (Ida et al., 2001; Salmon and Canup, 2012; Canup et al., 2015). Fig. 2 shows the outcome of the simulation with the Moon's surface temperature $T_M = 1800$ K. Panels 2a-2f show the surface density of the disk at different times, while panels 2g and 2h give the radial profile of surface density and temperature, respectively.

For a Moon located at $3 R_E$, we obtain that the gas vapor is energetic enough to cross the lunar Hill sphere and enter in orbit around the Earth, corroborating the scenario proposed by Charnoz et al. (2021). The gas exits the region gravitationally dominated by the Moon via spiral arms, resulting in the formation of a circum-Earth disk mainly composed of volatile species. The evolution of the disk is similar to that expected for a protoplanetary disk with a large planet immersed in it (e.g., Kley and Nelson, 2012; Morbidelli and Raymond, 2016; Armitage et al., 2019), including the formation of the gap at the lunar location.

The material released from the Moon expands throughout the entire disk, which eventually reaches an isothermal steady state. To understand this process, we must distinguish between the gas already present in the disk and the gas recently released from the Moon. The ejected gas, initially at the same temperature as the Moon, decompresses as it expands, resulting in a local drop in temperature. Simultaneously, this expanding gas heats the colder gas already present in the disk. This process repeats continuously, with the disk being progressively heated by the newly released gas until a steady state is reached, which occurs when the surface density of the disk approaches that of the Moon. At this point, the disk's temperature stabilizes at lunar temperature, and the gas only undergoes local compression or decompression in order to maintain a constant temperature.

In our model, we find that increases in the surface temperature lead to an increase in the amount of material released from the Moon. Nonetheless, the general dynamics of the circum-Earth disk are similar to those shown in Fig. 2. The main impact of lunar temperature on the structure of the circum-Earth disk is on the depth of the gap, which decreases with temperature. This result is somewhat expected since the disk height is proportional to the square root of temperature, and the proportional relation between the disk height and the depth of the gap is well-documented (Ziampras et al., 2020).

Fig. 3 illustrates the evolution of the gas material. Panel (a) shows the gas motion streamlines at a selected time, while panel (b) presents the longitudinal flux density at different distances from the Moon, as a function of lunar longitude.<fnref n="3" /> Panel (c) shows the fraction of re-accreted

<footnote n="2">² Animations of Figs. 2 and 4 can be found at gusmadeira.github.io/data.</footnote>

<footnote n="3">³ We define zero longitude on the line connecting the Earth and the Moon, on the near side. The west corresponds to the leading side, and the east to the trailing side.</footnote>
<caption>**Table 1**</caption>

Total pressure, adiabatic index, mean molecular mass, and atomic mass fraction of the lunar magma ocean, initially assumed to have the same composition as the bulk silicate Earth.

<table><thead><tr><th></th><th>Symbol</th><th>1400 K</th><th>1600 K</th><th>1800 K</th><th>2000 K</th><th>2200 K</th></tr></thead><tbody><tr><td>Total pressure (bar)</td><td>P</td><td>3.2×10<sup>-8</sup></td><td>1.2×10<sup>-6</sup></td><td>2.0×10<sup>-5</sup></td><td>1.9×10<sup>-4</sup></td><td>1.2×10<sup>-3</sup></td></tr><tr><td>adiabatic index</td><td>γ</td><td>1.3552</td><td>1.2806</td><td>1.2278</td><td>1.1938</td><td>1.1725</td></tr><tr><td>mean molecular mass (g/mol)</td><td>μ</td><td>25.3</td><td>25.4</td><td>25.5</td><td>25.9</td><td>26.6</td></tr><tr><td>atomic mass fraction</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>O</td><td>μ<sub>O</sub></td><td>0.253</td><td>0.252</td><td>0.252</td><td>0.254</td><td>0.264</td></tr><tr><td>Si</td><td>μ<sub>Si</sub></td><td>&lt;0.001</td><td>&lt;0.001</td><td>0.001</td><td>0.006</td><td>0.024</td></tr><tr><td>Mg</td><td>μ<sub>Mg</sub></td><td>&lt;0.001</td><td>&lt;0.001</td><td>&lt;0.001</td><td>0.001</td><td>0.004</td></tr><tr><td>Fe</td><td>μ<sub>Fe</sub></td><td>0.001</td><td>0.004</td><td>0.013</td><td>0.030</td><td>0.056</td></tr><tr><td>Na</td><td>μ<sub>Na</sub></td><td>0.698</td><td>0.689</td><td>0.675</td><td>0.648</td><td>0.593</td></tr><tr><td>K</td><td>μ<sub>K</sub></td><td>0.048</td><td>0.054</td><td>0.059</td><td>0.061</td><td>0.059</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Outcome of the simulation with the molten Moon located at 3 $R_E$ and a surface temperature of 1800 K. Panels (a)-(f) show the surface density of the circum-Earth disk, with green circles indicating the lunar physical radius and the Hill sphere. Panels (g) and (h) show the radial surface density and radial temperature of the disk, respectively. The system is presented in the rotating frame with the Moon.</figcaption></figure>

material as a function of lunar longitude, and panel (d) illustrates the amount of material released from the Moon, re-accreted by the satellite, and lost at the inner and outer boundaries of the disk.

Material is expected to enter orbit around the Earth when crossing the lunar Hill sphere (solid green curve in Fig. 3a), which is theoretically computed using the simplified relation $R_H = a_M(M_M/(3M_E))^{1/3} \approx 0.6a_M/R_E$ in Moon radii. Nonetheless, when the Moon is close to Earth, planetary tides significantly distort the lunar Roche lobe (Kopal, 1972; Murray and Dermott, 2000), meaning that a perfectly symmetrical Hill sphere does not accurately represent the region where the Moon's gravity is dominant. This can be seen by the yellow and pink

pure-gravitational equipotential curves in Fig. 3a, related to the $L_1$ and $L_2$ points of the system,<fnref n="4" /> respectively. To fully understand the circulation of material around the Moon, Fig. 3a shows the longitudinal flux density not only at the circular Hill radius but also at 1.2 Moon radii (dotted green curve), approximately corresponding to the closest distance of the Roche lobe to the Moon.

<footnote n="4">⁴ The equipotential curves are computed without accounting for the pressure effect of the circum-Earth disk. Pressure effects are known to only slightly change the position of $L_1$ and $L_2$ (Fig. 10, Casoli and Masset, 2009).</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** (a) Surface density and motion streamlines in the vicinity of the Moon, where the solid black line indicates the physical radius of the Moon, and the dashed and solid green lines correspond to circumferences with radii equal to the Hill radius $R_H$ and 1.2 Moon radii, respectively. The yellow and pink curves represent pure-gravitational equipotential curves related to the $L_1$ and $L_2$ points, respectively. (b) Longitudinal flux density at 1.2 Moon radii (dotted lines) and at the Hill radius (solid lines) at different times, as a function of lunar longitude. (c) Fraction of re-accreted material at different longitudes of the Moon. (d) The amount of material released from the Moon, accreted by Earth, re-accreted by the Moon, and lost at the outer boundary of the disk. The Moon is located at $3 R_E$ and has a surface temperature of 1800 K.</figcaption></figure>

The pressure effect of the gas disk is observed to displace the system's Lagrange points (Casoli and Masset, 2009), causing $L_1$ and $L_2$ to librate around 30° W and 150° E, respectively. Although volatile gas is released from the entire cross-section of the Moon, we find that the gas flows out of the Moon's Roche lobe via spiral arms passing through $L_1$ and $L_2$. These spiral arm regions correspond to areas of positive flux in Fig. 3b, with higher flux observed on the near side of the Moon. Volatile re-enter the Roche lobe in regions other than the spiral arms, with a fraction of this material flowing toward the Moon and being re-accreted at roughly the same longitudes where it first penetrates the lobe. The circum-Earth disk material re-enters Moon's Roche lobe mostly after orbiting the Earth for a timescale of a few days.

Over the course of one year of simulation, a total mass of $1.1 \times 10^{-6} M_M$ is released and escapes the Moon's Roche lobe, with approximately 18% of the material eventually re-entering this region and being re-accreted by the satellite. Consequently, the net loss of volatile is reduced to $9 \times 10^{-7} M_M$. Of the material not re-accreted, about 18% is lost at the inner boundary of the disk – that is, it is accreted by Earth – and 72% is lost at the outer boundary of the disk. The remainder remains in the circum-Earth disk until the end of the 1-year simulation.

Figs. 4 and 5 show the same panels as Figs. 2 and 3, but for the Moon initially at $8 R_E$. In this case, the theoretical Hill radius corresponds to almost 5 times the Moon's physical radius, an obvious consequence of the new lunar position. The total amount of material released from the Moon is greater than that obtained for $a_M = 3 R_E$ ($4.1 \times 10^{-6} M_M$), however, only 5% of the material is actually put into orbit around Earth. Because of this, a much lighter disk of volatile is formed when the Moon is further away from the planet.

For $a_M = 3 R_E$, the surface density of the disk reaches peak values of $0.2 \text{ g/cm}^2$, while for $a_M = 8 R_E$, it is always lower than $0.05 \text{ g/cm}^2$. For the latter case, it is observed that 40% of the disk material penetrates the Moon's Roche lobe after few dozen lunar orbits and is re-accreted by the satellite, which implies a small net loss of volatiles in the Moon. The vast majority of the material released from the Moon (95% of the total) never crosses the lunar Roche lobe and remains in the lunar extended atmosphere, cycling around the satellite for a few hours before being re-accreted. This material tends to move away from the satellite on its leading side and approach it on the trailing side. In fact, in all our simulations with $a_M > 3.5 R_E$, volatiles are preferentially re-accreted on the trailing side of the Moon (Fig. 5c), likely due to the satellite's orbital motion. In the next section, we analyze the fluxes of gas in the system.

### 3.1. Fluxes of ejection and accretion of volatile

Fig. 6 gives the fluxes and fate of material ejected from the Moon in our 1-year simulations, plotted as functions of Earth-Moon distance and surface temperature. As seen in Table 1, the pressure is highly sensitive to surface temperature. Consequently, the lunar surface density – and therefore the ejection flux – increases with temperature (dashed lines in Fig. 6a), despite the inverse relation between $\Sigma_M$ and $T_M$ in Equation (10). For $T_M = 1400$ K, the ejection fluxes are of the order of $10^{-9} M_M$/year, while for $T_M = 2200$ K, they increase to $10^{-4} M_M$/year. It is noteworthy that although the surface density does not depend on the Moon's semi-major axis, the height of the background disk does, which influences the ejection fluxes. This explains the increase in ejection flux with distance.

Nevertheless, the ejection flux is not representative of the material actually lost by the Moon, since some of it returns to the satellite. As discussed earlier, an increasing fraction of material is reaccreted by the Moon without ever crossing the Roche lobe as the Moon's semi-major axis expands. In addition, some of the material actually put in orbit around the Earth is re-accreted by the satellite after a few orbital periods. The re-accretion fraction (Fig. 6b), combines the re-accretion in these both scenarios, and shows a weak dependence on the surface temperature. For $a_M = 3 R_E$, ~15-20% of the ejected material is reaccreted
<figure><img src="image_7.png" /><figcaption>Fig. 4. Outcome of the simulation with the molten Moon located at 8 R_E with a surface temperature of 1800 K. The surface density of the circum-Earth disk is shown in panels (a) to (f) for different times. The system is in the rotating frame with Moon. The green circles represent the lunar physical radius and the Hill sphere. Panels (g) and (h) show the radial surface density and disk temperature, respectively.</figcaption></figure>

by the Moon, while for $a_M = 9 R_E$, only about 1% of the ejected material does not return to the satellite.

The net loss flux from the Moon is given by the solid lines in Fig. 6. A direct comparison of the fluxes obtained in this study with those of the dry model of Charnoz et al. (2021) (their Fig. 7) is not feasible due to the different assumptions about the vapor gas. Specifically, Charnoz et al. (2021) assumes constant molar mass (35 g/mol) and adiabatic index (≈1.5) while assuming the pressure given by the adiabatic equation of state. Our model, on the other hand, assumes these quantities as functions of the lunar temperature. As a result, we obtain a much smoother decrease in the net loss flux with the distance, compared to the estimates of Charnoz et al. (2021). This suggests that the tidally-assisted hydrodynamic escape may be more efficient than previously envisioned, especially considering that our results take into account the material re-accreted by the Moon – a factor not considered in Charnoz et al. (2021).

The Fig. 6c presents the fate of the material not re-accreted by the Moon, showing that the lunar temperature has a small impact on the fraction of material lost at the edges of the disk. In general, about 10% of the material remains in the disk at the end of the simulation, while most of it is lost at the outer edge, which we interpret as a tendency of the material to be released into the interplanetary environment. Some of the gas is also accreted by the Earth, indicating that the tidally assisted hydrodynamic escape not only devolatilizes the lunar mantle, but also contributes to the enrichment of the Earth with volatile materials. The

fraction accreted by the Earth increases with distance, reaching approximately 10% of the material not recreated by the Moon for $a_M = 3 R_E$, and reaching ~35% for $a_M = 9 R_E$. However, it is important to note that in absolute terms, the amount of material accreted by Earth when the Moon is at 3 $R_E$ is greater than it is at larger distances, due to the greater amount of gas injected in the circum-Earth disk.

## 4. Moon's orbital expansion and composition

In this section, we take into account the migration of the Moon due to tidal torques and calculate its composition over time. The dynamic evolution of the early Moon is a puzzle and various scenarios have been explored over the years (Webb, 1982; Čuk and Stewart, 2012; Čuk et al., 2016, 2021; Tian et al., 2017). In particular, the uncertainty about the tidal effects on a molten Earth and Moon turns the post-formation evolution of the satellite into an intricate question (Chen and Nimmo, 2016). For simplicity, we adopt a simplified constant-Q tidal model for the evolution of the semi-major axis, described by (Goldreich and Soter, 1966):

$$ \frac{da_M}{dt} = 3 \frac{k_2}{Q} \frac{M_M}{M_E} \sqrt{\frac{GM_E}{a_M}} \left(\frac{R_E}{a_M}\right)^5 \quad (12) $$

where $k_2$ and Q are the Love number and quality factor of Earth, respectively. Currently, these values are $k_2 = 0.3$ and $Q \approx 10$ ($k_2/Q \approx 0.03$, Williams and Boggs, 2015). Whilst the quality factor of a molten Earth
<figure><img src="image_8.png" /></figure>

**Fig. 5.** (a) Surface density and motion streamlines in the vicinity of the Moon, initially with $a_M = 8 R_E$ and $T_M = 1800$ K. The solid black line indicates the physical radius of the Moon, and the dashed and solid green lines correspond to circumferences with radii equal to the Hill radius and 3.1 Moon radii, respectively. The yellow and pink curves represent equipotential curves related to the $L_1$ and $L_2$ points, respectively. (b) Longitudinal flux density at 3.1 Moon radii (dotted lines) and at the Hill radius (solid lines) at different times. (c) Fraction of the total re-accreted mass as a function of lunar longitude. (d) Amount of material released from the Moon, accreted by Earth, re-accreted by the Moon, and lost at the outer boundary of the disk.

is unknown, it is believed to be higher at that time due to lower internal dissipation (Chen and Nimmo, 2016; Charnoz et al., 2021). Simultaneously, $k_2$ is expected to have been larger due to the absence of a shear modulus in the planet (e.g., see Raevskiy et al., 2015). Seen this, we investigate the early Moon's evolution by assuming different val-

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

**Fig. 6.** Panel (a) shows the flux of gas vapor ejected from the Moon (dashed lines) and the flux of material effectively lost by the satellite (solid lines). Panel (b) presents the fraction of the ejected material that is re-accreted by the satellite. Panel (c) shows the fraction of material effectively lost from the Moon that is accreted by Earth (solid lines), lost at the outer boundary (dashed lines), and that remains in the circum-Earth disk (dotted lines) at the end of the simulation. All plots are presented as functions of the Earth-Moon distance.

ues of the ratio $k_2/Q$, specifically $k_2/Q = 0.3, 0.03, 0.003$, and $0.0003$. In addition, we explore the scenario proposed by Ćuk and Stewart (2012), in which the Moon is temporarily captured in ejection resonance. Such a resonance occurs when the precession period of lunar perigee is equal to the orbital period of the Earth. For this scenario, we adopt the time evolution given in Figure 3 of Ćuk and Stewart (2012), for which $k_2/Q \approx 0.003$. The time evolution of the Moon's semi-major axis is shown in Fig. 7.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Time evolution of the Earth-Moon distance. In the case of the colored lines, the Moon is migrating only due to the Earth's tides, with different values of $k_2/Q$. The black line, on the other hand, corresponds to the scenario in which the Moon is also involved in an ejection resonance, as shown in Figure 3 of Ćuk and Stewart (2012).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Mass fractions of Na (black solid line) and K (red dashed line) in the Moon, in relation to BSE. The lunar surface temperature is 1800 K and the Earth's tidal parameter is 0.003. The time evolution is given on the lower x-axis, while the respective semi-major axis of the Moon is shown on the upper x-axis. The black and red regions outline, respectively, the ranges of Na and K mass fractions estimated in Visscher and Fegley (2013). The markers delimit the interval of time and semi-major axis in which the mass fraction of both elements simultaneously match the estimates.</figcaption></figure>

In the several tide migration scenarios considered, we account for the material lost from the Moon over time by linearly interpolating the net loss fluxes shown in Fig. 6. Our particular interest lies in the abundances of Na and K, which are evaluated based on the fractions of atoms shown in Table 1 (see Charnoz et al., 2021). The analysis of the abundance of moderately volatile elements is interesting because they are measurable in limited amount of materials (Wolf and Anders, 1980), allowing a direct comparison between our results and laboratory data.

In Fig. 8, we present the fractions of Na and K in relation to BSE as a function of time, for a scenario in which $T_M = 1800$ K and $k_2/Q = 0.003$. Initially, the Moon's composition reflects that of the BSE and the mass fractions are 1. Over time, the Moon migrates outwards and loses material through hydrodynamic escape. The abundance of sodium, initially greater than that of potassium, decreases more rapidly. Charnoz et al. (2021) obtains that the loss of volatile via hydrodynamic escape ceases before the Moon reaches a distance of 6 $R_E$, leaving a content of moderately volatile elements in the Moon potentially similar to the estimates. According to Visscher and Fegley (2013), the current abundances of Na and K relative to BSE in the Moon vary within the ranges of 0.11-0.63 and 0.11-0.8, respectively.

In our analysis, the net loss flux decreases more slowly with distance than observed in Charnoz et al. (2021) and we obtain in the simulation shown in Fig. 8, the Moon being totally depleted of Na and K in less than 25,000 years. This suggests that mechanisms other than lunar migration must be at work in the system to prevent the total loss of volatiles. As will be discussed in Section 5, the formation of a stagnant lid could

<figure><img src="image_12.png" /><figcaption>Fig. 9. (a) Time interval for the end of volatile loss and respective (b) lunar positions for our entire set of simulations, with the Moon starting at 3 $R_E$. The x-axis shows the temperature of the lunar surface, while the different colors correspond to different Earth tidal parameters. Our calculations were carried out up to $10^6$ years.</figcaption></figure>

be one of these mechanisms. Fig. 9 gives, for each of our simulations, the time interval and the corresponding lunar positions at which the loss of volatile must be stopped, in order to explain the observed lunar abundances of Na and K.

This time interval corresponds to the period during which the mass fractions of Na and K simultaneously match the estimates of Visscher and Fegley (2013). The midpoint of this range is calculated as the average between the minimum and maximum values. Our analysis extends up to $10^6$ years. As expected, higher temperatures on the Moon's surface result in faster volatile loss, requiring an earlier halt to this process. At the same time, lower $k_2/Q$ leads to slower migration, which means that the satellite spends a longer period closer to Earth. Therefore, in scenarios with lower $k_2/Q$, the interruption of volatile loss also needs to occur earlier. For $k_2/Q = 0.003$, the time interval for the loss of volatile remains similar, regardless of whether the Moon is in ejection resonance or not. This phenomenon occurs because the Moon loses most of its volatile before being captured in the resonance. For a surface temperature of 1400 K, it is observed that the Moon never loses enough of its volatile content.

In the scenario of the Moon forming from a proto-lunar disk, although the first aggregates coalesce at about 3 $R_E$, the Moon fully forms only at around 5 $R_E$ (Salmon and Canup, 2012). Seen this, we performed the same analysis as in Fig. 9, but starting with the Moon at 5 $R_E$. To account for the direct formation scenario, we also considered a case where the Moon begins at 7 $R_E$, a location consistent with those obtained in Kegerreis et al. (2022); Meier et al. (2024). The average times for the cessation of volatile loss in these two cases are shown in Fig. 10. For the Moon at 5 $R_E$, the cessation times shows only a twofold increase, maintaining the same order of magnitude as those represented in Fig. 9. For a Moon starting at 7 $R_E$, the volatile loss flux is significantly lower, and in cases of rapid tidal migration, the Moon does not lose sufficient
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Average time for the end of volatile loss for cases where the Moon starts at 5 $R_E$ (solid lines with circles) and 7 $R_E$ (dotted lines with triangles). The lunar surface temperature is shown on the x-axis, and the different colors represent different Earth tidal parameters.</figcaption></figure>

volatile material. The Na and K mass fractions match the estimates only for $k_2/Q \le 0.003$.

However, we draw the reader's attention to the fact that these results apply to a Moon in a circular orbit. For the disk formation scenario, this approximation is reasonable, as the Moon is expected to form with $e \sim 0.01$ (Salmon and Canup, 2012). For the direct formation scenario, in contrast, this assumption may not hold, as the Moon is expected to form in a highly eccentric orbit ($e \sim 0.3$, Meier et al., 2024). Assuming a semi-major axis of 7 $R_E$ and $e = 0.3$, we performed a simple estimation of the mass loss over one orbital period by interpolating the flux at different orbital locations. This analysis indicates an approximately 30% increase in mass loss compared to the circular case, suggesting that the cessation time of a directly formed Moon might be shorter than that shown in Fig. 10, although we do not expect the values to differ by orders of magnitude. Nonetheless, this should be further investigated in future publications. In the next section, we calculate the timescales for the crystallization of the stagnant lid.

**5. Formation of an early stagnant lid on the Moon**

The formation of a conductive lid over a vigorously convecting magma ocean on the Moon may result both from the buoyancy of suspended crystals within the magma during mantle solidification (Elkins-Tanton et al., 2011; Elkins-Tanton, 2012) or the establishment of a stagnant lid due to rapid surface cooling into a magma with strongly temperature-dependent viscosity (e.g. Davaille and Jaupart, 1993). This conductive lid can effectively extend magma ocean cooling duration (e.g. Watson et al., 2022; Michaut and Neufeld, 2022; Lichtenberg et al., 2023) and has the potential to reduce or even prevent volatile degassing.

Elkins-Tanton et al. (2011) estimate that the formation of an anorthite flotation crust (plagioclase) on the lunar magma ocean should happen on a timescale of ~10³ years, assuming instantaneous crystal flotation. Anorthite floats due to its lower density compared to the interstitial liquid. Nonetheless, anorthite's contribution is negligible for the formation of a conductive lid until solidification of the lunar magma ocean reaches ~80-90% or a depth of ~100 km. A stagnant lid that emerges early enough could limit heat loss from the magma ocean and therefore delay the formation of the anorthite crust (Perera et al., 2018; Michaut and Neufeld, 2022). Here, to estimate the onset of the conductive lid, we calculate the thermal evolution of the lunar magma ocean alongside the growth of a stagnant lid. Details about our model for stagnant lid formation are provided in Appendix A.

By integrating Equation (A.5), we obtain the thermal evolution of the lunar magma ocean over 10⁵ years, assuming logarithmically spaced timesteps. At each interval, we assess the criteria expressed in Equation (A.8) for stagnant lid formation. Upon matching this criterion, a stagnant lid forms and the heat flux extracted from the magma ocean due to convection decreases. Our preliminary calculations show that the

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Temporal evolution of (a) temperatures, (b) stagnant lid thickness, and (c) depth of the lunar magma ocean. The color-coding represents numerical calculations with lid emergence of sizes of 1 m, 10 m, 100 m, and without a lid. The purple dashed line denotes the start of anorthite solidification when magma ocean depth drops below 100 km.</figcaption></figure>

criterion is reached quite early, around 10 years, when the magma ocean has a potential temperature of ≈ 2800 K. However, the thickness of the stagnant lid developed during such a few years is on the order of 10⁻⁷ m, which in practical terms is negligible.

To address this limitation, we propose that the stagnant lid becomes effective in inhibiting degassing and reducing magma ocean convection with the surface once it reaches a certain size, $h_0$. Determining this critical size falls beyond the scope of this paper; hence, we consider $h_0$ to span three orders of magnitude from 1 m to 100 m. Therefore, for every timestep in which Equation (A.8) holds true, we calculate the size of the stagnant lid $h$ via Equation (A.12). Also, when $h \ge h_0$, we alter the convective heat flux in the magma ocean from $F(\Delta T)$ to $F(T_\mu)$.

Fig. 11 illustrates the evolution of the lunar magma ocean with an initial potential temperature of 3000 K for a case without stagnant lid formation and for three different values of $h_0$. The blue lines in Fig. 11(a) denote the evolution of potential temperature (solid line) and surface temperature (dash-dotted line) for the case with no lid formation, in which the temperatures show an evolution similar to those obtained for planets without an atmosphere (e.g. Lebrun et al., 2013; Nikolaou et al., 2019). Fig. 11(c) displays the evolution of magma ocean depth. To calculate the depth, we first obtain the radial temperature profile $T(r,t)$ given by the adiabat with $T_p(t)$, which also provides the solid fraction in the magma ocean. Then, we assume that all formed solids sink and instantly settle in the liquid magma, adding to the magma ocean's bottom. We can estimate the onset of anorthite crust formation when the magma ocean depth reaches 100 km (e.g. Elkins-Tanton et al., 2011; Tian et al., 2017). In Fig. 11(c), the purple dashed line indicates the onset of anorthite crust according to this criterion.

We observe that the anorthite crust formation occurs in approximately 1800 years for the magma ocean evolving without stagnant lid formation, a finding consistent with multi-composition fractional crys-
tallization models (Elkins-Tanton et al., 2011). For $h_0 = 1$ m, lid emergence occurs around 501 years. At the time of lid emergence, Fig. 11(a) shows the surface temperature drops due to the decrease in convective heat flux from $F(\Delta T)$ to $F(T_\mu)$ (see Appendix A for details) and the magma ocean's potential temperature remains slightly higher than in the case with no lid. Lid emergence for $h_0$ equal to 10 m and 100 m occurs at 698 years and 1012 years, respectively. In Fig. 11(c), we note that the onset of the stagnant lid precedes the beginning of anorthite solidification for all $h_0$ values considered. Yet, the early emergence of the stagnant lid delays plagioclase formation to the order of $10^5$ years.

For the sake of completeness, we have computed the case with $h_0 = 1$ km, obtaining that the time to form the lid is roughly the same as the time to form the flotation crust. For different initial potential temperatures, we find that the variation in the timescale for the emergence of the lid is marginal. This is explained by the fact that convection, and therefore the heat flux escaping from the magma ocean, is proportional to the mantle's potential temperature.

We also compare the convective timescales of the magma ocean with the timescales during which gas parcels remain cycling in the lunar extended atmosphere, which primarily occurs for $a_M \ge 3.2 R_M$. The convective velocities of the magma ocean during most of the convecting-only phase is on the order of 1 m/s, resulting in overturn timescales of about 20 days. We calculate the timescale for a gas parcel to be accreted by Moon to be about 0.4 day, depending on the assumed lunar surface temperature. This implies that ~50 cycles of gas in atmosphere occur before the magma ocean undergoes a complete overturn and news its surface.

Given the simplicity of the model used, the strong conclusion of this section is that, regardless of the stability of the stagnant lid over the lunar magma ocean, the different scenarios indicate that volatile loss should be inhibited in a timescale of $10^3$ years. This timescale is consistent with our analysis of the end of volatile loss (Fig. 9), provided that the lunar temperature is around 1800 K-2000 K. This result is particularly valid for $k_2/Q$ values lower than that of the present-day Moon, while for equal or higher values, the best match occurs at temperatures around 2000 K. For the direct formation scenario, in which the Moon forms at $7 R_E$, $k_2/Q \approx 0.0003$ is required to explain the concentrations of Na and K in the satellite, with the best matches occurring for temperatures around 2000 K.

It should be noted that the various temperatures employed are derived from distinct problems and assumptions. The lunar temperature of previous sections is a mean temperature of the magma ocean and is fixed, while both potential and surface temperatures in this section vary in time. Nonetheless, by comparing the time-weighted average potential temperature of the magma ocean ($\bar{T}_p$), as discussed herein, with the lunar temperature ranging from 1800 K to 2000 K, we find good correspondence. For the smallest $h_0$, $\bar{T}_p$ is equal to 1990.77 K and when $h_0$ is set to 100 m, $\bar{T}_p$ gives 1844.11 K.

# 6. Discussion

In this paper, we put to trial the hypothesis of the tidally-driven atmospheric escape scenario. Charnoz et al. (2021) assume a 1D steady-state approach and apply the hydrodynamic escape formalism to the Earth-Moon line to calculate the material crossing the lunar Hill sphere. Here, we evaluate the problem in a more general way by performing 2D time-dependent hydrodynamic simulations, which allows us to track the dynamic evolution of the volatile gas. As predicted by Charnoz et al. (2021), we find that when the Moon is closer to the Earth, most of the gas material that crosses the Moon's Roche lobe does not return to the satellite. The fraction of re-accretion by the Moon is less than 30% when the satellite is located up to $3.5 R_E$. Our results therefore contradict those of Dauphas et al. (2022), which state that hydrodynamic escape from the Moon is insufficient to explain the depletion of moderate volatile elements on the satellite. However, the physics we consider is more realistic as we discuss below.

## 6.1. On the work of Dauphas et al. (2022)

Dauphas et al. (2022) reported a critical comparison of different models of volatile loss. By using the Hertz-Knudsen (HK) evaporation theory – which describes the kinetic evaporation of a gas from a free surface –, they link the measured isotopic fractionation to the evaporating flux (see, e.g. Richter et al., 2002, for a detailed model). In fact, the evaporation flux ($J_i$) in the kinetic theory of gases is proportional to the difference between the gas pressure (P) and the equilibrium saturating vapor pressure of the liquid ($P_{sat}$), so for a species i, $J_i = \gamma_s(P - P_{sat})/(2\pi\mu kT)$, where $\gamma_s$ is a sticking efficiency coefficient – measured experimentally and between 0 and 1 – and $\mu$ is the mass of the molecule (see, e.g. Hirth and Pound, 1963; Richter et al., 2002). Different isotopes with slightly different values of $\mu$ will have slightly different evaporating fluxes and different vapor pressures. These different evaporation fluxes induce a net isotopic fractionation of the material that remains in the liquid, leading to progressive enrichment in the heavier isotopes. Based on this argument, Dauphas et al. (2022) proposed that the measured isotopic variations of K in lunar rocks imply that the effective vapor pressure above the lunar magma must be $P \sim 0.99 P_{sat}$ (with $T \sim 1700$ K). Due to the close-to-total saturation and low temperature, it is concluded that about 10 Myrs would be necessary to explain the measured isotopic composition, a time by which a stagnant lid would have formed (~1000 years).

Whereas we do not dispute the simple calculations presented in Dauphas et al. (2022), the approach developed in this paper and in Charnoz et al. (2021) relies on completely different physics: the Navier-Stokes equation – conservation of mass, energy, and momentum – is solved and the net flux at the Hill sphere of the Moon is computed, not simply at the magma ocean surface as in the HK model. Many additional physical complexities arise that make the HK calculation too simplistic: for example, the escaping flux is highly anisotropic, potentially leading to non-isotropic fractionation and mixing in the lunar atmosphere before escape. Also, molecular and turbulent isotopic diffusion in the lunar magma ocean and the atmosphere could alter isotopic fractionation before escape (see Fig. 5). Finally, Hertz-Knudsen equation is derived from the kinetic theory of gases and is adapted to evaporation from a free liquid surface with a vapor with zero net velocity (Hirth and Pound, 1963). To account for overly simplistic physics, the HK theory uses several multiplicative factors to correct for physics not included in the equation (e.g. the sticking coefficient $\gamma_s$). In the problem considered in Charnoz et al. (2021), the gas flow is forced by a strong pressure gradient – this is a Parker Wind – with non-zero net velocity at the magma ocean surface, which is a priori outside the domain of validity of the HK theory. Therefore, the kinetic fractionation of HK developed in the Dauphas et al. (2022) model does not capture all the complexity of the escape process presented in Charnoz et al. (2021).

This timescale discrepancy between dynamical processes and atomic fractionation processes is not limited to the tidally assisted escape model. In Dauphas et al. (2022), it is shown that viscous drainage from a disk, using the MRI model of Charnoz and Michaut (2015), could lead to an isotopic fractionation of K compatible with lunar samples in less than 400 days. However, dynamical studies of lunar disk cooling and spreading show that spreading and cooling timescales are about 50-100 years (see Thompson and Stevenson, 1988; Ward, 2012), leading again to an apparent paradox where dynamical and isotopic fractionation timescales are very different.

In addition, the viscous drainage model requires an ionized environment (Charnoz et al., 2021) with a non-zero local magnetic field to foster the magnetorotational turbulence. Isotopic fractionation processes in an ionized environment are not well understood and can produce mass-independent fractionation processes (Kuga et al., 2017; Robert et al., 2021), which can further complicate the problem. We conclude that it is still very difficult to reconcile the fractionation timescale deduced from the mere application of the HK theory with dynamical models of
hydrodynamical escape or lunar disk evolution. Improvements could be made on both sides. Thus, for now, this question remains open.

## 6.2. Comparison with Charnoz et al. (2021) and general remarks

In spite of the differences in methodology and in the properties of the lunar magma ocean between our work and that of Charnoz et al. (2021), we obtain escape fluxes that are fairly consistent. This is especially evident, at least in order of magnitude, with those reported by Charnoz et al. (2021) when the Moon is at about 3 $R_E$. The Earth's tides reduce the energy needed for thermal escape and, closer to the Earth, it turns out that all the gas content formed by the evaporation of the magma ocean is energetic enough to cross Moon's Roche lobe and be placed in orbit around the Earth. Only a fraction of this gas is re-accreted by the Moon, resulting in an efficient net loss of volatiles.

As the Moon migrates outwards due to tidal torques with the Earth, the energy required for thermal escape increases and less material is able to cross Moon's Roche lobe. However, we find a much smoother drop in the loss flux with distance than Charnoz et al. (2021). While the authors conclude that the loss of volatiles from the Moon is halted up to 6 $R_E$, we obtain that, at this position, the net loss flux is reduced by only one order of magnitude (Fig. 6). This implies that tidally-driven atmospheric escape is even more efficient than previously estimated. Regardless of whether the loss of volatiles starts already in the lunar building blocks, in which the loss flux should be even more efficient due to their smaller masses, or starts only with the fully formed Moon at around 5 $R_E$, the loss is expected to be strong enough to completely deplete the objects in volatiles in Na and K, if the surface temperature is above 1600 K.

In an isothermal expanding gas cloud, a fraction of the volatiles is expected to condense due to pressure changes as the gas flows toward the lunar Roche lobe, an effect not considered in this study. Nonetheless, condensate droplets small enough ($\lesssim$ µm) are expected to remain well-coupled to the gas, being carried outward by inertia. Furthermore, if the early Earth's photosphere had temperatures $\gtrsim$ 2000 K, micrometer-sized condensates would have been subjected to radiation pressure, possibly being blown out of the system (Hyodo et al., 2018). Therefore, we expect that condensation would increase the loss flux, making hydrodynamic escape even more efficient. Indeed, Charnoz et al. (2021) showed that this effect generates a 40% increase in loss flux, and we leave the question of how condensation quantitatively impacts the net loss flux for future publications.

Additionally, Earth tides would distort the Moon, giving it a nonspherical shape (Chandrasekhar, 1989). For example, assuming the satellite is at 3 $R_E$ and approximating it as a perfect fluid, we calculate that its projection in the XY plane corresponds to an ellipse with semi-axes of 0.95 and 1.15 Moon radii (see Madeira et al., 2023). In general, planetary tides increase the satellite's surface area, which leads to larger ejection fluxes. Test simulations accounting for the Moon's distortion show that the net volatile flux can increase by approximately 30% at 3 $R_E$. Therefore, Moon's distortion is also expected to make tidally assisted hydrodynamic escape more efficient, even though it is not expected to change the order of magnitude of our results.

These results highlight the need for another mechanism, in addition to the Moon's orbital evolution, to explain the abundances of moderate volatile elements measured in lunar mare basalts. We suggest that a conductive lid could develop over the magma ocean, stopping the loss of volatiles and surface recycling, either due to a high temperature contrast – thermal stagnant lid – or the accumulation of buoyant solids – anorthite flotation crust. Considering a simple model of thermal evolution and stagnant lid growth, we find that a stagnant lid ranging in size from 1 m to 100 m should form on top of the convective layer within an interval of about 500-1000 years. If this thermal stagnant lid is stable, it would block magma ocean degassing within 1000 years. However, if the stagnant lid is not stable, the lunar magma ocean is expected to lose heat more rapidly, initiating anorthite solidification and the formation of a flotation crust around 1800 years. Regardless of the scenario envisioned

in our model, volatile loss should stop on a timescale of $10^3$ years. To achieve such timing, we constrain the temperature of the early Moon's magma ocean between 1800-2000 K., which is in line with models based on the Cr isotopic composition of the Moon (Sossi et al., 2018).

For the lunar positions $\gtrsim$3.2 $R_E$, we obtain that the material is predominantly re-accreted on the trailing side of the Moon, which, in principle, could result in a dichotomy of volatile element concentrations between the east and west sides of the satellite. This is only valid if the satellite lost most of its volatile content when fully formed, since the impacts between its building blocks would erase any possible dichotomy they might have had. In addition, other mechanisms might act over time to reduce this dichotomy, such as the redistribution of surface material due to meteorite impacts and the reorientation of lunar longitudes due to migration or other events – the Moon's current zero longitude meridian does not necessarily correspond to the same one from billions of years ago.

However, if measured, the dichotomy in the content of volatiles in the Moon could give clues about the migration of the satellite. For the same temperature, the time needed to obtain the current abundances of Na and K in the satellite are of the same order, regardless of tidal migration (Fig. 9). However, the location of the Moon at this time depends on the migration parameter. For the Moon's current $k_2/Q$ and temperature of 2000 K, gas loss is predominantly away from the Earth, and more than 99% of the material is re-accreted on the trailing side of the Moon, implying a strong dichotomy. Now, for a slow migration with $k_2/Q = 0.0003$, for example, the loss of volatiles occurs within 3.2 $R_E$, and only ~60% of the material is re-accreted on the trailing side.

## 6.3. Caveats

The early Moon's degassing is an intricate problem, with many potentially relevant effects being not accounted for in our model. These effects, however, should be explored in future studies. The first of these, discussed in the previous section, concerns the possible condensation of volatiles – an effect that cannot be currently computed given that the FARGOCA code, in its present version, is not a multi-phase code. Also mentioned, the Moon distortion by Earth tides is another of these effects. We assume the Moon to be spherical in our simulations, which is a reasonable approximation only if the early Moon had a high internal cohesion (see Madeira et al., 2023). However, the lack of knowledge of physical parameters of the early Moon and early Earth force us to always rely on assumptions regarding tidal distortion, and investigating the Moon's shape effect in the hydrodynamic escape process is beyond the scope of this work.

Another limitation of our model is related to the 2D approximation. Due to the high temperature of the material ejected from the Moon, a thick circum-Earth disk forms, similar to the expected to the protolunar disk (Nakajima and Stevenson, 2014; Ward, 2011). At the end of our simulations, the disk height at 3 $R_E$ reaches about 0.5 Moon radii, and at 8 $R_E$, it exceeds 2 Moon radii. While our 2D approximation may be marginally appropriate at 3 $R_E$, where hydrodynamic escape is strongest, we acknowledge that we may be overlooking some physical effects when the Moon is at greater distances from the Earth. The extended atmosphere of a Moon embedded in a thick disk is expected to exchange material with the circum-Earth disk not only through spiral arms (a 2D effect), but also via meridional gas circulation near the gap (see details in Teague et al., 2019; Schulik et al., 2020). These meridional flows deliver material to the satellite and consequently, the re-accretion flux onto the Moon might be higher than obtained by us. A more realistic model of the disk would require highly computationally expensive 3D hydrodynamical simulations, and our computations were already very demanding with the current setup.

We also mention the uncertainty surrounding the initial composition of the Moon. While it is well-established that the Moon and BSE formed with the same concentration of refractory elements (Ringwood, 1992), our assumption that the Moon has a BSE-like concentration of volatile
elements is more speculative. To better understand how this assumption might impact the system, we recomputed the vapor composition by varying the mass fractions of Na and K in the magma ocean. Our findings show that the lunar surface density – which determines the ejection flux – is proportional to the concentration of these elements. For instance, when the mass fractions of Na and K are doubled, the lunar surface density increases by a factor of 1.4.

The lunar surface density impacts the ejection flux (Fig. 6a), but it is not expected to significantly affect the fraction of material that is reaccreted by the Moon (Fig. 6b), implying that the net loss flux will also be proportional to the initial volatile content. Therefore, it is possible that the overall timescale for the end of volatile loss we obtained, at least in terms of order of magnitude, remains unchanged for different volatile contents in the magma ocean. This is because any initial excess or deficit of volatiles relative to the BSE might be balanced by corresponding changes in depletion rates. We will leave the exploration of this scenario for a future publication.

To conclude, we address a comment on the gas adiabatic index. The VAPOROCK code provides the actual adiabatic index of the vapor gas, which may differ from the value expected in the vertically integrated (2D) scenario simulated here. This difference arises because, in the latter, there is an additional degree of freedom in the system due to vertical disk expansion (Goldreich et al., 1986), and the 2D adiabatic index should be lower than the actual value used in our simulations. To quantify the impact of this parameter, we calculated the net loss flux for $T_M = 1800$ K and $a_M = 3 R_E$, varying the adiabatic index. A 10% reduction in this parameter increases the flux by an order of magnitude, whereas reducing the flux by the same amount requires a 60% increase in the adiabatic index. Within the range $0.9\gamma - 1.6\gamma$ (1.105–1.964 for $T_M = 1800$ K and $a_M = 3 R_E$), the volatile loss timescale is expected to remain of the same order of magnitude as the results presented here, leaving our conclusions unaffected.

## 7. Conclusion

We conclude that tidally driven atmospheric escape is a solid mechanism to explain the depletion of moderate volatile elements in the Moon, as it does not rely on any specific formation model. Our model is able to explain the current observed depletion of Na and K on the satellite, without requiring any prior loss of volatiles in the proto-Moon disk, which is uncertain to have existed (Kegerreis et al., 2022; Sossi et al., 2024). Our results show that the hydrodynamic escape of the early Moon is more intense than previously predicted by Tang and Young (2020); Charnoz et al. (2021); Dauphas et al. (2022), with the formation of a stagnant lid being the key mechanism preventing the complete depletion of Na and K in the satellite. To reproduce the current abundances of these elements, the temperature of the Moon’s magma ocean must be around 1800-2000 K.

### CRediT authorship contribution statement

**G. Madeira:** Writing – review & editing, Writing – original draft, Investigation, Formal analysis, Data curation, Conceptualization. **L. Esteves:** Writing – review & editing, Writing – original draft, Investigation. **S. Charnoz:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Conceptualization. **E. Lega:** Writing – review & editing, Writing – original draft, Software, Methodology, Conceptualization. **F. Moynier:** Writing – review & editing, Writing – original draft, Resources, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

GM and FM thank the European Research Council (#101001282, METAL). LE acknowledges the financial support from FAPESP through grants 2021/00628-6 and 2023/09307-3. HPC resources from “Mesocentre SIGAMM”, hosted by Observatoire de la Côte d’Azur were used. We thank Alain Miniussi for maintaining and developing the code FAR-GOCA and Alessandro Morbidelli, Aurelien Crida, Angela Limare, and Maylis Landeau for helpful discussions. We thank the anonymous reviewers for their comments, which greatly helped us to improve the article.

## Appendix A. Magma ocean model

In this work, we adapted the model of mixed interior convection below a growing conductive stagnant lid from Watson et al. (2022) to start from a convecting-only lunar magma ocean without a lid. Table A.2 describes the general parameters used in the model. As in Watson et al. (2022), we adopt solidus and liquidus temperatures linearized with respect to pressure $p$ (Maurice et al., 2017), while the viscosity for the liquid and solid phases of the magma follows the Arrhenius law, where the viscosity contrast between solid and liquid is on the order of $10^{22}$ Pa·s. In the mixing phase, the viscosity remains continuous following the function $f(\phi)$ of Costa et al. (2009). The top panel of Fig. A.12 shows the viscosity function for $p = 0$ GPa. Empirical constants $\delta, \gamma$, and the critical solid fraction $\phi^*$ are chosen as in Table A.2, and $\xi$ is used as a free parameter to ensure viscosity continuity at the limit $\phi \to 1$.

Given that the magma ocean is convecting, the temperature profile follows an adiabat as in Lebrun et al. (2013), but with a pressure range of 0 to ~8 GPa for the Moon. The variation of temperature with respect to pressure is given by:

$$ \frac{dT}{dp} = \frac{\alpha T}{\rho C_p}, \qquad (A.1) $$

where $T(p=0)$ is the potential temperature $T_p$, $\alpha$ is the thermal expansion, and $C_p$ is the specific heat capacity. By calculating the adiabatic temperature profiles for various $T_p$, it shows that even in the relatively low-pressure environment of the Moon’s mantle, the adiabats first intersect the liquidus (and solidus) at the highest pressures. This indicates bottom-up solidification.

Integrating the energy conservation equation over the volume of the magma ocean ($V_m$), we obtain the evolution of the potential temperature in time (Abe, 1997):

$$ \int_m dV_m \rho C_p \frac{dT_p}{dt} = \int_m (-F + \rho L \frac{d\phi}{dt}) dV_m, \qquad (A.2) $$

with $L$ being the latent heat of solidification. The heat flux, $F$, follows the Rayleigh-Bénard convection equation and is given as a function of the temperature difference across the unstable boundary $\Delta T$ and the Rayleigh number $Ra$:

$$ F(\Delta T) = C \frac{k \Delta T}{d} Ra(\Delta T)^{1/3}, \qquad (A.3) $$

where

$$ Ra(\Delta T) = \frac{\rho \alpha g d^3}{\kappa \mu} \Delta T. \qquad (A.4) $$

This leads to:

$$ F(\Delta T) = Ck \left( \frac{\rho \alpha g}{\kappa \mu} \right)^{1/3} \Delta T^{4/3}. \qquad (A.5) $$

Here, $F(\Delta T)$ is independent of the height of the convective layer $d$. $k = \kappa \rho C_p$ is the thermal conductivity, $\kappa$ is the thermal diffusivity and $C$ is an experimental constant of proportionality. While the magma ocean evolves without a stagnant lid, the heat flux is measured immediately
<table><caption>**Table A.2**
Stagnant lid evolution model quantities.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Description</th><th>Reference</th></tr></thead><tbody><tr><td colspan="4"><strong>MOON PARAMETERS</strong></td></tr><tr><td>R</td><td>1737</td><td>Moon radius (km)</td><td></td></tr><tr><td>R<sub>c</sub></td><td>380</td><td>Moon core radius (km)</td><td></td></tr><tr><td>ρ</td><td>3328</td><td>mantle mean density (kg m<sup>-3</sup>)</td><td></td></tr><tr><td>ρ<sub>c</sub></td><td>5171</td><td>core mean density (kg m<sup>-3</sup>)</td><td></td></tr><tr><td>g(r)</td><td><sup>G(M<sub>mol</sub> + M<sub>molk</sub>)(r)</sup>/<sub>r<sup>2</sup></sub></td><td>Moon gravity a radius r (m s<sup>-1</sup>)</td><td></td></tr><tr><td colspan="4"><strong>VISCOSITY PARAMETERS</strong></td></tr><tr><td>T<sub>liq</sub></td><td>1400 + 149.5p</td><td>silicate liquidus and solidus. p (GPa), T (K)</td><td></td></tr><tr><td>T<sub>sol</sub></td><td>1977 + 64.1p</td><td></td><td></td></tr><tr><td>μ<sub>l</sub></td><td>10<sup>21</sup> exp[<sup>E<sub>l</sub></sup>/<sub>R<sub>l</sub></sub>(<sup>1</sup>/<sub>T<sub>l</sub></sub> - <sup>1</sup>/<sub>T<sub>liq</sub></sub>)]</td><td>solid phase viscosity (Pa s)</td><td>Maurice et al. (2017)</td></tr><tr><td>μ<sub>s</sub></td><td>0.1 exp[<sup>E<sub>s</sub></sup>/<sub>R<sub>s</sub></sub>(<sup>1</sup>/<sub>T<sub>s</sub></sub> - <sup>1</sup>/<sub>T<sub>sol</sub></sub>)]</td><td>liquid phase viscosity (Pa s)</td><td></td></tr><tr><td>μ<sub>l+s</sub></td><td>μ<sub>l</sub>f(φ)</td><td>partially molten phase viscosity (Pa s)</td><td></td></tr><tr><td>f(φ)</td><td>See Costa et al. (2009)</td><td>empirical function for relative viscosity.</td><td>Watson et al. (2022) Costa et al. (2009)</td></tr><tr><td>φ*</td><td>0.5</td><td>critical solid fraction</td><td></td></tr><tr><td>δ<sub>s</sub>, γ</td><td>7, 4</td><td>empirical constants for f(φ)</td><td></td></tr><tr><td>B<sub>E</sub></td><td>2.5</td><td>Einstein coefficient</td><td></td></tr><tr><td colspan="4"><strong>MAGMA OCEAN PARAMETERS</strong></td></tr><tr><td>α<sub>0</sub></td><td>3 × 10<sup>-5</sup></td><td>thermal expansivity at p = 0 (K<sup>-1</sup>)</td><td></td></tr><tr><td>α</td><td>α<sub>0</sub>(0.02p + 1)<sup>-3/4</sup></td><td>thermal expansivity (K<sup>-1</sup>).</td><td></td></tr><tr><td>C<sub>p</sub></td><td>10<sup>3</sup></td><td>specific heat capacity (K<sup>-1</sup>)</td><td></td></tr><tr><td>κ</td><td>8.8 × 10<sup>-7</sup></td><td>thermal diffusivity (m<sup>2</sup> s<sup>-1</sup>)</td><td></td></tr><tr><td>k</td><td>κρC<sub>p</sub></td><td>thermal conductivity (W m<sup>-1</sup> K<sup>-1</sup>)</td><td>Lebrun et al. (2013) Nikolaou et al. (2019)</td></tr><tr><td>L</td><td>4 × 10<sup>5</sup></td><td>latent heat for silicate solidification (J kg<sup>-1</sup>)</td><td></td></tr><tr><td>ε</td><td>1</td><td>black body emissivity</td><td></td></tr><tr><td>a</td><td>0.3</td><td>Earth's albedo</td><td></td></tr><tr><td>σ</td><td>5.670374 × 10<sup>-8</sup></td><td>Stefan-Boltzmann constant (W m<sup>-2</sup> K<sup>-4</sup>)</td><td></td></tr><tr><td>S<sub>0</sub></td><td>1361</td><td>solar flux (W m<sup>-2</sup>)</td><td></td></tr></tbody></table>

at the surface with $\Delta T = T_p - T_{surf}$ representing the temperature difference at the top of the convective layer and C assumes the value of 0.089 (Lebrun et al., 2013). Note that the dynamical viscosity $\mu$ is both temperature and pressure dependent due to our assumed formulation for viscosity, while g and $\alpha$ can be considered as $g(r=0)$ and $\alpha_0$, respectively.

We calculate the surface temperature by balancing the convective flux with the irradiating flux at the surface:

$$F(\Delta T) = F_{\text{rad}} = \sigma\epsilon(T_{\text{surf}}^4 - T_{\text{eq}}^4) \qquad (A.6)$$

where the equilibrium temperature $T_{eq}$ is given in terms of albedo $a$, Stefan-Boltzmann constant $\sigma$ and Solar flux $S_0$ (see Table A.2):

$$T_{\text{eq}}^4 = \frac{(1-a)}{4\sigma} S_0. \qquad (A.7)$$

At each timestep, we compute the condition for the existence of the stagnant lid due to the viscosity contrast at the top of the convection layer, given by Davaille and Jaupart (1993)

$$T_{\mu} = -\frac{\mu(T, p)}{d\mu(T, p)/dT} < \frac{(T_p - T_{\text{surf}})}{B}. \qquad (A.8)$$

The above condition depends on the factor B that can be constrained by laboratory experiments (Davaille and Jaupart, 1993) or numerical models (Thiriet et al., 2019). We use $B = 2.54$ provided by Thiriet et al. (2019). If the condition is satisfied, a stagnant lid forms at the top of the magma ocean reducing the heat flux from magma ocean convection and the surface temperature. In the model, the heat flux drop is represented by the fact that $F(\Delta T)$ is then calculated with $\Delta T = T_\mu$ which is smaller than $T_p - T_{surf}$ and the constant C becomes 0.47 as measured in the experiments of Davaille and Jaupart (1993). Another consequence is a new boundary interface between the stagnant lid and the magma ocean at the depth of the stagnant lid $h$ as in Watson et al. (2022). By assuming a temperature drop at the boundary between the magma ocean and the stagnant lid given by $T_\mu B$, we can obtain the relation for the temperature at the bottom of the stagnant lid, $T_{bsl}$, as

$$T_{\mu}B = T_{p} - T_{bsl}. \qquad (A.9)$$

The convective flux from the magma ocean $F(T_\mu)$ is now balanced with the heat conduction flux through the stagnant lid $F_{\text{cond}}$, and the growth of the stagnant lid is proportional to the difference $F_{\text{cond}} - F(T_\mu)$. The conductive flux can be calculated as

$$F_{\text{cond}} = -k \left. \frac{\partial T}{\partial z} \right|_{z = R-h} \qquad (A.10)$$

As we are only interested in the time period during which the stagnant lid emerges and has very small size relative to the convection layer, we consider that the timescale for heat diffusion through the lid is very short so that the temperature profile through the lid is approximately linear. This implies that

$$F_{\text{cond}} = -k \frac{(T_{\text{surf}} - T_{\text{bsl}})}{h} \approx F(T_{\mu}) \qquad (A.11)$$

Then,

$$h = -k \frac{(T_{\text{surf}} - T_{\text{bsl}})}{F(T_{\mu})}. \qquad (A.12)$$

This approximation does not produce significant differences in the calculation of the stagnant lid thickness $h$, and it is much faster compared to integrating a heat conduction equation for such a thin stagnant lid (< 1 km). Note that due to this approximation, $F(T_\mu)$ remains the heat flux reaching the surface and therefore is used to determine $T_{surf}$.

In the Fig. A.12, we plot the viscosity $\mu(T, p)$, rheological temperature scale $T_\mu$ and heat flux $F_{\text{conv}}(T_\mu)$ of a magma ocean below a stagnant lid. As the potential temperature decreases due to the cooling of the magma ocean, the viscosity increases, transitioning from values on the order of $10^{-1}$ Pa·s to $10^{21}$ Pa·s. With a potential temperature ($T_p$) of 1977 K, the temperature of the magma ocean crosses the liquidus, initiating crystallization at a pressure of $p = 0$ GPa. From this point, the function $f(\phi)$ drives the viscosity from $\mu_l$ to $\mu_s$ during the increase of the fraction of solids in the magma (see Costa et al., 2009). $T_\mu$ depends on the derivative of viscosity that is very sensitive to slope variations along the viscosity curve. Consequently, both $T_\mu$ and $F_{\text{conv}}(T_\mu)$ experience abrupt changes upon crossing the liquidus and solidus temperatures. In
<figure><img src="image_15.png" /><figcaption>Fig. A.12. Viscosity, rheological temperature scale and convective heat flux in function of the potential temperature.</figcaption></figure>

the calculations addressed in this work, the stagnant lid remains quite thin (< 1 km), so the thermal boundary between the lid and the magma ocean remains close to the surface, with *p* ≈ 0 GPa.

## Data availability

Data will be made available on request.

## References

* Abe, Y., 1997. Thermal and chemical evolution of the terrestrial magma ocean. Phys. Earth Planet. Inter. 100, 27-39. https://doi.org/10.1016/S0031-9201(96)03229-3.

* Armitage, P.J., Kley, W., Kley, W., 2019. Planet Formation and Disk-Planet Interactions. From Protoplanetary Disks to Planet Formation. Sas-a-Fee Advanced Course, vol. 45. Swiss Society for Astrophysics and Astronomy, pp. 151-260.

* Benz, W., Cameron, A., Melosh, H., 1989. The origin of the Moon and the single-impact hypothesis iii. Icarus 81, 113-131.

* Boyce, J.W., Treiman, A.H., Guan, Y., Ma, C., Eiler, J.M., Gross, J., Greenwood, J.P., Stolper, E.M., 2015. The chlorine isotope fingerprint of the lunar magma ocean. Sci. Adv. 1, e1500380.

* Canup, R.M., 2008. Lunar-forming collisions with pre-impact rotation. Icarus 196, 518-538.

* Canup, R.M., 2012. Forming a Moon with an Earth-like composition via a giant impact. Science 338, 1052-1055.

* Canup, R.M., Asphaug, E., 2001. Origin of the Moon in a giant impact near the end of the Earth's formation. Nature 412, 708-712.

* Canup, R.M., Visscher, C., Salmon, J., Begley Jr, B., 2015. Lunar volatile depletion due to incomplete accretion within an impact-generated disk. Nat. Geosci. 8, 918-921.

* Casoli, J., Masset, F.S., 2009. On the horseshoe drag of a low-mass planet. I. Migration in isothermal disks. Astrophys. J. 703, 845-856. https://doi.org/10.1088/0004-637X/703/1/845. arXiv:0907.4677.

* Chandrasekhar, S., 1989. The equilibrium and the stability of the Roche ellipsoids. Selected Papers, Volume 4: Plasma Physics, Hydrodynamic and Hydromagnetic Stability, and Applications of the Tensor-Virial Theorem 4, 430.

* Charnoz, S., Michaut, C., 2015. Evolution of the protolunar disk: dynamics, cooling timescale and implantation of volatiles onto the Earth. Icarus 260, 440-463.

* Charnoz, S., Sossi, P.A., Lee, Y.N., Siebert, J., Hyodo, R., Allibert, L., Pignatatale, F.C., Landeau, M., Oza, A.V., Moynier, F., 2021. Tidal pull of the Earth strips the proto-moon of its volatiles. Icarus 364, 114451.

* Chen, E.M., Nimmo, F., 2016. Tidal dissipation in the lunar magma ocean and its effect on the early evolution of the Earth-moon system. Icarus 275, 132-142.

* Costa, A., Caricchi, L., Bagdassarov, N., 2009. A model for the rheology of particle-bearing suspensions and partially Molten rocks. Geochem. Geophys. Geosyst. 10, Q03010. https://doi.org/10.1029/2008GC002138.

* Ćuk, M., Stewart, S.T., 2012. Making the Moon from a fast-spinning Earth: a giant impact followed by resonant despinning. Science 338, 1047-1052.

* Ćuk, M., Hamilton, D.P., Lock, S.J., Stewart, S.T., 2016. Tidal evolution of the Moon from a high-obliquity, high-angular-momentum Earth. Nature 539, 402-406.

* Ćuk, M., Lock, S.J., Stewart, S.T., Hamilton, D.P., 2021. Tidal evolution of the Earth-moon system with a high initial obliquity. Planet. Sci. J. 2, 147.

* Dauphas, N., Nie, N.X., Blanchard, M., Zhang, Z.J., Zeng, H., Hu, J.Y., Meheut, M., Visscher, C., Canup, R., Hopp, T., 2022. The extent, nature, and origin of k and rb depletions and isotopic fractionations in Earth, the Moon, and other planetary bodies. Planet. Sci. J. 3, 29.

* Davaille, A., Jaupart, C., 1993. Transient high-Rayleigh-number thermal convection with large viscosity variations. J. Fluid Mech. 253, 141-166. https://doi.org/10.1017/S0022112093001740.

* Day, J., Moynier, F., Ishizuka, O., 2022. A partial melting control on the Zn isotope composition of basalts. Geochem. Perspect. Lett. 23, 11-16.

* Day, J.M., Moynier, F., 2014. Evaporative fractionation of volatile stable isotopes and their bearing on the origin of the Moon. Philos. Trans. R. Soc. A, Math. Phys. Eng. Sci. 372, 20130259.

* de Val-Borro, M., Edgar, R., Artymowicz, P., Ciecielag, P.e., Cresswell, P., d'Angelo, G., Delgado-Donate, E., Dirksen, G., Fromang, S., Gawryszczak, A., et al., 2006. A comparative study of disc-planet interaction. Mon. Not. R. Astron. Soc. 370, 529-558.

* Elkins-Tanton, L.T., 2012. Magma oceans in the inner solar system. Annu. Rev. Earth Planet. Sci. 40, 113-139. https://doi.org/10.1146/annurev-earth-042711-105503.

* Elkins-Tanton, L.T., Burgess, S., Yin, Q.Z., 2011. The lunar magma ocean: reconciling the solidification process with lunar petrology and geochronology. Earth Planet. Sci. Lett. 304, 326-336. https://doi.org/10.1016/j.epsl.2011.02.004.

* Garcia, R.F., Gagnepain-Beyneix, J., Chevrot, S., Lognonné, P., 2011. Very preliminary reference Moon model. Phys. Earth Planet. Inter. 188, 96-113. https://doi.org/10.1016/j.pepi.2011.06.015.

* Gargano, A., Sharp, Z., Shearer, C., Simon, J.I., Halliday, A., Buckley, W., 2020. The cl isotope composition and halogen contents of Apollo-return samples. Proc. Natl. Acad. Sci. 117, 23418-23425.

* Goldreich, P., Soter, S., 1966. Q in the solar system. Icarus 5, 375-389.

* Goldreich, P., Goodman, J., Narayan, R., 1986. The stability of accretion tori-I. Long-wavelength modes of slender tori. Mon. Not. R. Astron. Soc. 221, 339-364.

* Herwartz, D., Pack, A., Friedrichs, B., Bischoff, A., 2014. Identification of the giant impactor Theia in lunar rocks. Science 344, 1146-1150.

* Hirth, J., Pound, G., 1963. Condensation and evaporation: nucleation and growth kinetics. Prog. Mater. Sci., Macmillan. https://books.google.fr/books?id=joFFAQAAIAAJ.

* Hyodo, R., Genda, H., Charnoz, S., Pignatatale, F.C., Rosenblatt, P., 2018. On the impact origin of Phobos and Deimos. IV. Volatile depletion. Astrophys. J. 860, 150.

* Ida, S., Kokubo, E., Takeda, T., 2001. N-body simulations of Moon accretion. In: Collisional Processes in the Solar System, pp. 203-222.

* Ito, Y., Ikoma, M., Kawahara, H., Nagahara, H., Kawashima, Y., Nakamoto, T., 2015. Theoretical emission spectra of atmospheres of hot rocky super-earths. Astrophys. J. 801, 144. https://doi.org/10.1088/0004-637X/801/2/144. arXiv:1501.05393.

* Kaib, N.A., Cowan, N.B., 2015a. Brief follow-up on recent studies of Theia's accretion. Icarus 258, 14-17.

* Kaib, N.A., Cowan, N.B., 2015b. The feeding zones of terrestrial planets and insights into Moon formation. Icarus 252, 161-174.

* Kato, C., Moynier, F., 2017. Gallium isotopic evidence for extensive volatile loss from the Moon during its formation. Sci. Adv. 3, e1700571.
Kato, C., Moynier, F., Valdes, M.C., Dhaliwal, J.K., Day, J.M., 2015. Extensive volatile loss during formation and differentiation of the Moon. Nat. Commun. 6, 7617.

Kegerreis, J.A., Ruiz-Bonilla, S., Eke, V.R., Massey, R.J., Sandnes, T.D., Teodoro, L.F., 2022. Immediate origin of the Moon as a post-impact satellite. Astrophys. J. Lett. 937, L40.

Kley, W., Nelson, R., 2012. Planet-disk interaction and orbital evolution. Annu. Rev. Astron. Astrophys. 50, 211–249.

Kokubo, E., Canup, R., Ida, S., Righter, K., 2000. Lunar accretion from an impact-generated disk. In: Origin of the Earth and Moon, pp. 145–163.

Kopal, Z., 1972. The Roche model and its applications to close binary systems. Adv. Astron. Astrophys. 9, 1–65. https://doi.org/10.1016/B978-0-12-003209-9.50006-9.

Kuga, M., Cernogora, G., Marrocchi, Y., Tissandier, L., Marty, B., 2017. Processes of noble gas elemental and isotopic fractionations in plasma-produced organic solids: cosmochemical implications. Geochim. Cosmochim. Acta 217, 219–230. https://doi.org/10.1016/j.gca.2017.08.031.

Lebrun, T., Massol, H., ChassefièRe, E., Davaille, A., Marcq, E., Sarda, P., Leblanc, F., Brandeis, G., 2013. Thermal evolution of an early magma ocean in interaction with the atmosphere. J. Geophys. Res., Planets 118, 1155–1176. https://doi.org/10.1002/jgre.20068.

Lega, E., Crida, A., Bitsch, B., Morbidelli, A., 2014. Migration of Earth-sized planets in 3d radiative discs. Mon. Not. R. Astron. Soc. 440, 683–695.

Lichtenberg, T., Schaefer, L.K., Nakajima, M., Fischer, R.A., 2023. Geophysical evolution during rocky planet formation. In: Inutsuka, S., Aikawa, Y., Muto, T., Tomida, K., Tamura, M. (Eds.), Protostars and Planets VII, p. 907. arXiv:2203.10023.

Lock, S.J., Stewart, S.T., Petaev, M.I., Leinhardt, Z., Mace, M.T., Jacobsen, S.B., Cuk, M., 2018. The origin of the Moon within a terrestrial synestia. J. Geophys. Res., Planets 123, 910–951.

Lock, S.J., Stewart, S.T., Ćuk, M., 2020. The energy budget and figure of Earth during recovery from the Moon-forming giant impact. Earth Planet. Sci. Lett. 530, 115885.

Lodders, K., Fegley Jr, B., 2015. Chemistry of the solar system. R. Soc. Chem.

Madeira, G., Charnoz, S., Zhang, Y., Hyodo, R., Michel, P., Genda, H., Winter, S.G., 2023. Exploring the recycling model of Phobos formation: rubble-pile satellites. Astron. J. 165, 161.

Masset, F., 2000. Fargo: a fast Eulerian transport algorithm for differentially rotating disks. Astron. Astrophys. Suppl. Ser. 141, 165–173.

Masset, F.S., 2002. The co-orbital corotation torque in a viscous disk: numerical simulations. Astron. Astrophys. Suppl. Ser. 387, 605–623. https://doi.org/10.1051/0004-6361:200202040. arXiv:astro-ph/0205211.

Maurice, M., Tosi, N., Samuel, H., Plesa, A.C., Hüttig, C., Breuer, D., 2017. Onset of solid-state mantle convection and mixing during magma ocean solidification. J. Geophys. Res., Planets 122, 577–598. https://doi.org/10.1002/2016JE005250.

Meier, T., Reinhardt, C., Timpe, M., Stadel, J., Moore, B., 2024. A systematic survey of Moon-forming giant impacts. II. Rotating bodies. arXiv preprint. arXiv:2409.02746.

Michaut, C., Neufeld, J.A., 2022. Formation of the lunar primary crust from a long-lived slushy magma ocean. Geophys. Res. Lett. 49, e95408. https://doi.org/10.1029/2021GL095408.

Morbidelli, A., Raymond, S.N., 2016. Challenges in planet formation. J. Geophys. Res., Planets 121, 1962–1980.

Mougel, B., Moynier, F., Göpel, C., 2018. Chromium isotopic homogeneity between the Moon, the Earth, and enstatite chondrites. Earth Planet. Sci. Lett. 481, 1–8.

Murray, C.D., Dermott, S.F., 2000. Solar System Dynamics. Cambridge University Press.

Nakajima, M., Stevenson, D.J., 2014. Investigation of the initial state of the Moon-forming disk: bridging sphs simulations and hydrostatic models. Icarus 233, 259–267.

Nakajima, M., Stevenson, D.J., 2018. Inefficient volatile loss from the Moon-forming disk: reconciling the giant impact hypothesis and a wet Moon. Earth Planet. Sci. Lett. 487, 117–126. https://doi.org/10.1016/j.epsl.2018.01.026. arXiv:1812.10502.

Nakajima, M., Atkins, J., Simon, J.B., Quillen, A.C., 2024. The limited role of the streaming instability during Moon and exomoon formation. Planet. Sci. J. 5, 145.

Nikolaou, A., Katyal, N., Tosi, N., Godolt, M., Grenfell, J.L., Rauer, H., 2019. What factors affect the duration and outgassing of the terrestrial magma ocean? Astrophys. J. 875, 11. https://doi.org/10.3847/1538-4357/ab08ed. arXiv:1903.07436.

Paniello, R.C., Day, J.M., Moynier, F., 2012. Zinc isotopic evidence for the origin of the Moon. Nature 490, 376–379.

Parker, E., 1963. Kinematical hydromagnetic theory and its application to the low solar photosphere. Astrophys. J. 138, 552.

Parker, E., 1965. Dynamical theory of the solar wind. Space Sci. Rev. 4, 666–708.

Perera, V., Jackson, A.P., Elkins-Tanton, L.T., Asphaug, E., 2018. Effect of reimplementing debris on the solidification of the lunar magma ocean. J. Geophys. Res., Planets 123, 1168–1191. https://doi.org/10.1029/2017JE00551210.1002/essoar.976474f02c7aac3ca1320398ca3470b1.

Pringle, E.A., Moynier, F., 2017. Rubidium isotopic composition of the Earth, meteorites, and the Moon: evidence for the origin of volatile loss during planetary accretion. Earth Planet. Sci. Lett. 473, 62–70.

Raevskiy, S., Gudkova, T., Kuskov, O., Kronrod, V., 2015. On reconciling the models of the interior structure of the Moon with gravity data. Izv. Phys. Solid Earth 51, 134–142.

Richter, F.M., Davis, A.M., Ebel, D.S., Hashimoto, A., 2002. Elemental and isotopic fractionation of Type B calcium-, aluminum-rich inclusions: experiments, theoretical considerations, and constraints on their thermal evolution. Geochim. Cosmochim. Acta 66, 521–540. https://doi.org/10.1016/S0016-7037(01)00782-7.

Ringwood, A., 1992. Volatile and siderophile element geochemistry of the Moon: a reappraisal. Earth Planet. Sci. Lett. 111, 537–555.

Robert, F., Chaussidon, M., Gonzalez-Cano, A., Mostefaoui, S., 2021. Oxygen and magnesium mass-independent isotopic fractionation induced by chemical reactions in plasma. Proc. Natl. Acad. Sci. 118, e2114221118. https://doi.org/10.1073/pnas.2114221118.

Rufu, R., Canup, R.M., 2020. Tidal evolution of the ejection resonance/quasi-resonance and the angular momentum of the Earth-moon system. J. Geophys. Res., Planets 125, e2019JE006312.

Salmon, J., Canup, R.M., 2012. Lunar accretion from a Roche-interior fluid disk. Astrophys. J. 760, 83.

Schulik, M., Johansen, A., Bitsch, B., Lega, E., Lambrechts, M., 2020. On the structure and mass delivery towards circumplanetary discs. Astron. Astrophys. 642, A187.

Sharp, Z., Barnes, J., Brearley, A., Chaussidon, M., Fischer, T., Kamenetsky, V., 2007. Chlorine isotope homogeneity of the mantle, crust and carbonaceous chondrites. Nature 446, 1062–1065.

Sossi, P.A., Moynier, F., Van Zuilen, K., 2018. Volatile loss following cooling and accretion of the Moon revealed by chromium isotopes. Proc. Natl. Acad. Sci. 115, 10920–10925.

Sossi, P.A., Nakajima, M., Khan, A., 2024. Composition, structure and origin of the Moon. arXiv e-prints arXiv:2408.16840. https://doi.org/10.48550/arXiv.2408.16840.

Tang, H., Young, E., 2020. Evaporation from the lunar magma ocean was not the mechanism for fractionation of the Moon's moderately volatile elements. Planet. Sci. J. 1, 49.

Teague, R., Bae, J., Bergin, E.A., 2019. Meridional flows in the disk around a young star. Nature 574, 378–381.

Thiriet, M., Breuer, D., Michaut, C., Plesa, A.C., 2019. Scaling laws of convection for cooling planets in a stagnant lid regime. Phys. Earth Planet. Inter. 286, 138–153. https://doi.org/10.1016/j.pepi.2018.11.003.

Thompson, C., Stevenson, D.J., 1988. Gravitational instability in two-phase disks and the origin of the Moon. Astrophys. J. 333, 452. https://doi.org/10.1086/166760.

Tian, Z., Wisdom, J., Elkins-Tanton, L., 2017. Coupled orbital-thermal evolution of the early Earth-moon system with a fast-spinning Earth. Icarus 281, 90–102.

Tian, Z., Jolliff, B.L., Korotev, R.L., Fegley Jr, B., Lodders, K., Day, J.M., Chen, H., Wang, K., 2020. Potassium isotopic composition of the Moon. Geochim. Cosmochim. Acta 280, 263–280.

van Buchem, C.P.A., Miguel, Y., Zilinskas, M., van Westrenen, W., 2023. LavAtmos: an open-source chemical equilibrium vaporization code for lava worlds. Meteorit. Planet. Sci. 58, 1149–1161. https://doi.org/10.1111/maps.13994. arXiv:2210.10463.

Visscher, C., Fegley, B., 2013. Chemistry of impact-generated silicate melt-vapor debris disks. Astrophys. J. Lett. 767, L12.

Wang, K., Jacobsen, S.B., 2016. Potassium isotopic evidence for a high-energy giant impact origin of the Moon. Nature 538, 487–490.

Ward, W.R., 2011. On the vertical structure of the protolunar disk. Astrophys. J. 744, 140.

Ward, W.R., 2012. On the vertical structure of the protolunar disk. Astrophys. J. 744, 140. https://doi.org/10.1088/0004-637X/744/2/140.

Ward, W.R., Canup, R.M., Rufu, R., 2020. Analytical model for the tidal evolution of the ejection resonance and the timing of resonance escape. J. Geophys. Res., Planets 125, e2019JE006266.

Watson, C., Neufeld, J.A., Michaut, C., 2022. Early asymmetric growth of planetary stagnant lids. J. Fluid Mech. 952, A3. https://doi.org/10.1017/jfm.2022.864.

Webb, D., 1982. Tides and the evolution of the Earth—moon system. Geophys. J. Int. 70, 261–271.

Wiechert, U., Halliday, A.N., Lee, D.C., Snyder, G.A., Taylor, L.A., Rumble, D., 2001. Oxygen isotopes and the Moon-forming giant impact. Science 294, 345–348.

Williams, J.G., Boggs, D.H., 2015. Tides on the Moon: theory and determination of dissociation. J. Geophys. Res., Planets 120, 689–724.

Wolf, A.S., Jäggi, N., Sossi, P.A., Bower, D.J., 2023. VapoRock: thermodynamics of vaporized silicate melts for modeling volcanic outgassing and magma ocean atmospheres. Astrophys. J. 947, 64. https://doi.org/10.3847/1538-4357/acbcc7. arXiv:2208.09582.

Wolf, R., Anders, E., 1980. Moon and Earth: compositional differences inferred from siderophiles, volatiles, and alkalis in basalts. Geochim. Cosmochim. Acta 44, 2111–2124.

Wood, J.A., Dickey Jr, J.S., Marvin, U.B., Powell, B., 1970. Lunar anorthosites and a geophysical model of the Moon. In: Levinson, A.A. (Ed.), Geochimica et Cosmochimica Acta Supplement, Volume 1. In: Proceedings of the Apollo 11 Lunar Science Conference Held 5-8 January, 1970 in Houston, TX. Volume 1: Mineralogy and Petrology. Pergamon Press, New York, p. 965.

Zhang, J., Dauphas, N., Davis, A.M., Leya, I., Fedkin, A., 2012. The proto-Earth as a significant source of lunar material. Nat. Geosci. 5, 251–255.

Ziampras, A., Kley, W., Dullemond, C.P., 2020. Importance of radiative effects in gap opening by planets in protoplanetary disks. Astron. Astrophys. 637, A50.