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

HAL
open science

# Small-Scale Overturn of High-Ti Cumulates Promoted by the Long Lifetime of the Lunar Magma Ocean

M. Maurice, N. Tosi, C. Hüttig

► To cite this version:

M. Maurice, N. Tosi, C. Hüttig. Small-Scale Overturn of High-Ti Cumulates Promoted by the Long Lifetime of the Lunar Magma Ocean. Journal of Geophysical Research. Planets, 2024, 129, 10.1029/2023JE008060. insu-04729905

HAL Id: insu-04729905

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

Submitted on 10 Oct 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Key Points:**

* Solid-state mantle convection during lunar magma ocean solidification enhances the overturn of high-Ti cumulates

* The overturn always occurs via small-scale downwellings and lasts at most a few hundred Myr

* Overturned material accumulates above the core-mantle boundary where it heats up and partially melts without undergoing any large-scale destabilization

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. Maurice,
maximemaurice@protonmail.com

**Citation:**

Maurice, M., Tosi, N., & Hüttig, C. (2024). Small-scale overturn of high-Ti cumulates promoted by the long lifetime of the lunar magma ocean. *Journal of Geophysical Research: Planets*, 129, e2023JE008060.
https://doi.org/10.1029/2023JE008060

Received 23 AUG 2023

Accepted 23 JAN 2024

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Small-Scale Overturn of High-Ti Cumulates Promoted by the Long Lifetime of the Lunar Magma Ocean

M. Maurice¹,² ID, N. Tosi² ID, and C. Hüttig²

¹Laboratoire de Météorologie Dynamique, Institut Pierre Simon Laplace, CNRS, Paris, France, ²German Aerospace Center, Planetary Science Institute, Berlin, Germany

**Abstract** The fractional crystallization of the lunar magma ocean (LMO) results in a gravitationally unstable layering, with dense Fe- and Ti-oxides overlying lighter mafic cumulates. Due to their high density, these are prone to overturn via Rayleigh-Taylor instability. However, this instability competes with the downward growth of the cold and stiff stagnant lid, which tends to trap the cumulates preventing their overturn. The fate of high-Ti cumulates (HTC) plays a major role in several aspects of the Moon's history, including mare volcanism, early magnetism, and the presence of a partially molten layer at the core-mantle boundary (CMB). To assess the extent of the overturn of HTC, we use 2D simulations of thermo-chemical mantle convection in the presence of a solidifying LMO. The long lifetime of the magma ocean caused by the insulating effect of the plagioclase crust delays the growth of the stagnant lid and promotes the onset of solid-state convection during magma ocean solidification. Both phenomena favor a high degree of mobilization of the HTC. Independent of the rheology of the mafic and HTC, the overturn is always characterized by small-scale instabilities and is completed within ~300–600 Myr. High-Ti material accumulates at the CMB where it undergoes partial melting until present-day, in agreement with the existence of a deep, partially molten layer inferred from geodetic data. Part of the overturned cumulates is entrained by mantle flow and can participate in secondary melting until ~1 Ga, in agreement with the age of high-Ti mare basalts.

**Plain Language Summary** The Moon likely underwent a hot accretion resulting in a deep magma ocean. An important consequence of the solidification of the lunar magma ocean is the formation of a shallow layer, rich in Titanium (Ti) and denser than the underlying mantle. This configuration is gravitationally unstable, with the Ti-rich layer that is prone to sink to the bottom of the mantle—an event referred to as the lunar mantle overturn. Yet, cooling of the lunar crust and shallow mantle can be efficient enough to stiffen this layer, preventing its mobilization and overturn. Understanding to what extent the Ti-rich layer sunk through the mantle is important because its fate influences many aspects of the lunar evolution, including volcanism, magnetism, and the present-day thermal state of the Moon. Here, we use numerical simulations to show that, if the Ti-rich layer started sinking while the magma ocean was still extant, the overturn would have been much more efficient than if it started only after the magma ocean solidified completely, as assumed in previous studies. Indeed, the hot magma ocean delays the cooling of the crust and shallow mantle and facilitates the mobilization of Ti-rich material as soon as it crystallizes.

## 1. Introduction

The idea that the Moon once hosted a deep magma ocean was an important outcome of the analysis of the Apollo samples. The crystallization of buoyant plagioclase crystals out of a large body of denser magma, and their accumulation at the surface to form a thick floating crust offer the best explanation for the ubiquity of anorthosite (anorthite-bearing rock, anorthite being the Ca end-member of plagioclase) in Apollo samples (e.g., Wood et al., 1970). A hot initial state for the Moon is furthermore a natural outcome of giant impact models of Moon formation (Lock et al., 2018; Pritchard & Stevenson, 2000; Sahijpal & Goyal, 2018), although recently Kegerreis et al. (2022) found that a rapidly forming Moon would only be superficially molten. The crystallization of this primordial lunar magma ocean (LMO) sets the structure from which the solid lunar mantle subsequently evolved. The stratification inherited from the crystallization of the LMO has thus a major influence on the long-term evolution of the Moon, and potentially on many of the processes that are associated with the dynamics of its mantle.

From a dynamical point of view, one of the most relevant events that characterize the mantle following LMO solidification is the Rayleigh-Taylor instability induced by the late crystallization at shallow depth of dense
cumulates. The LMO is thought to crystallize from the bottom up to the surface (Elkins-Tanton, 2012), and to become more and more enriched in incompatible elements (in particular Ti, Fe, and heat-producing elements). After 90%-95% of the LMO solidified, Ti-bearing phases started to crystallize (Charlier et al., 2018; Elkins-Tanton et al., 2011; Johnson et al., 2021; Lin et al., 2017; Rapp & Draper, 2018; Schwinger & Breuer, 2022; Snyder et al., 1992), forming a dense layer below the flotation crust, which led to a gravitationally unstable configuration (de Vries et al., 2009; Hess & Parmentier, 1995; Parmentier et al., 2002). The development of such a gravitational instability, however, is hindered by the growth of a cold and stiff stagnant lid, where the dense material may remain trapped before being mobilized (Elkins-Tanton et al., 2002; van Orman & Grove, 2000; Yu et al., 2019; Zhao et al., 2019). Depending on the rheology of the lunar mantle, late-stage high-Ti cumulates (HTC) may either remain locked where they initially formed, or overturn and sink partly or entirely to the core-mantle boundary (CMB). Here we use the term “High-Ti cumulates” rather than the widely used “Ilmenite-bearing cumulates” (IBC) because ilmenite is not a major Ti-bearing phase in our crystallization model, where most of the Ti is incorporated in clinopyroxene (recently, Schmidt and Kraettli (2022) found spinel to be the most important phase to incorporate Ti in the LMO crystallization sequence). As this discrepancy might be a caveat of our model, we base the range of our rheological parameters on results obtained for ilmenite (Dygert et al., 2016).

An event such as the giant impact that formed the South Pole Aitken (SPA) basin has also been suggested to trigger HTC overturn and lead to the lunar compositional asymmetry (Jones et al., 2022; Zhang et al., 2022). However, here we focus on the case where the overturn is purely driven by a convective instability, be it Rayleigh-Taylor or Rayleigh-Bénard in nature.

If the HTC overturn did occur, a dense layer must have formed above the CMB. Due to the high density of the cumulates, such a layer is likely stable against convection. However, HTC are also enriched in heat-producing elements. Therefore, it has been proposed that heating of the overturned layer over hundreds of millions to billions of years could increase the layer's thermal buoyancy sufficiently to offset its high intrinsic density, and destabilize it again (Hess & Parmentier, 1995; Stegman et al., 2003; Zhang et al., 2013, 2017; Zhong et al., 2000). The feasibility of this scenario depends on the degree of enrichment in heat-producing elements as well as on the density of the HTC. The time lag associated with the destabilization of the dense, overturned HTC layer has been suggested to be the cause for the late emplacement of the lunar maria between 4 and 3 Ga (Hess & Parmentier, 1995). Furthermore, Zhong et al. (2000) suggested that if it occurred via a degree-1 upwelling, this mobilization could explain the dichotomy in the lunar crustal thickness and volcanism. The increase in the CMB heat flow resulting from such an event (as hot, upwelling HTC are replaced with cold material at the CMB), led Stegman et al. (2003) to suggest that it could promote a core dynamo responsible for the early crustal magnetization of the Moon (Weiss & Tikoo, 2014). Finally, overturned HTC remaining at the CMB could lower the melting temperature, and help forming a partially molten layer inferred by geophysical and geodetic data (Briaud et al., 2023; Khan et al., 2007, 2014), although Walterová et al. (2023) recently questioned this conclusion.

The mineralogy of lunar samples has also been associated with HTC. Some ultramafic glasses among the Apollo samples are highly enriched in Ti (black glasses from Apollo 14 contain up to 18 wt.% Ti (Delano, 1986)), thus requiring a Ti-rich source, for which HTC (or a fraction thereof) is the best candidate. The assimilation of Ti by mafic melts from an olivine-orthopyroxene mantle reservoir percolating through a shallow HTC layer composed of ilmenite and clinopyroxene was initially invoked by Wagner and Grove (1997), but van Orman and Grove (2000) later showed that such a scenario failed to reproduce the composition of the ultramafic glasses. They suggested instead that if shallow HTC could remelt, the accompanying melt would be negatively buoyant and sink through the underlying mafic mantle, reacting with it to form a hybrid source for high-Ti mare basalts (Singletary & Grove, 2008). However, this scenario was also found unable to match the composition of the high-Ti glasses (Singletary & Grove, 2008). These studies considered a shallow (i.e., not overturned) emplacement for the HTC. Overturned high-Ti material, lying at the CMB, is an unlikely source for surficial high-Ti melt production, because HTC-derived melt could be negatively buoyant at depth (Mallik et al., 2019; van Kan Parker et al., 2012; Xu et al., 2022). Nevertheless, the entrainment of overturned HTC material in the bulk convecting mantle could lead to a redistribution of Ti-bearing materials and provide a mechanism for the generation of hybridized sources for the high-Ti mare basalts. The mineralogy of mantle excavated by the impact of the SPA basin has also been connected to HTC, and Moriarty et al. (2021) suggested that the overturn couldn't be complete at the formation time of SPA.

The fate of HTC is a topic that has received significant attention, and nearly three decades of modeling efforts have led to our current understanding of the process. Hess and Parmentier (1995) proposed a scenario in which the
initially thin layer of HTC that crystallizes near the end of the LMO progressively grows via small-scale Rayleigh-Taylor instabilities, ultimately forming a thicker layer that drives a large-scale overturn. Based on this idea, Parmentier et al. (2002) investigated the overturn of a thickened, diluted HTC layer using numerical simulations, and confirmed the possibility of a large-scale overturn. By contrast, Elkins-Tanton et al. (2002) and van Orman and Grove (2000) deemed HTC overturn unlikely because of the high viscosity associated with their low crystallization temperature. However, this conclusion was based on calculations considering a dry olivine rheology yielding a viscosity potentially several orders of magnitude higher than that of HTC (Dygert et al., 2016). Several studies have investigated the possibility of HTC overturn since then (de Vries et al., 2009; Li et al., 2019; Yu et al., 2019; Zhao et al., 2019), but none has reached an unequivocal conclusion concerning the fate of high-Ti material in the lunar mantle.

All these studies address the dynamics of HTC overturn in a completely solidified mantle, assuming that the crystallization of the LMO proceeds on a shorter timescale than that of the onset of mantle convection. However, recent studies have shown that if the magma ocean lifetime exceeds $10^5-10^6$ years, both Rayleigh-Bénard and Rayleigh-Taylor instabilities can trigger solid-state convection in the growing cumulates of Mars (Maurice et al., 2017), of the Earth (Ballmer et al., 2017), as well as of the Moon (Boukaré et al., 2018; Morison et al., 2019). The existence of melting/freezing boundaries between the magma ocean and its cumulates can further help the onset of convection (Agrusta et al., 2019; Morison et al., 2019). The long lifetime of the LMO (~200 Myr) recently inferred by Maurice et al. (2020) would strongly favor this scenario. Noticeably, Zhao et al. (2019) introduced in their model setup a hot (and thus weak) layer in order to mimic the presence of a residual magma ocean. They found that the presence of such a weak layer enhances the amount of dense cumulates that can be mobilized. However, a self-consistent approach is still needed to capture the influence of the magma ocean on the dynamics of the overturn as well as to reveal the feedbacks between the mantle flow and the thermal evolution of the LMO (Maurice et al., 2020).

In this study, we investigate the impact of a solidifying LMO on the efficiency of HTC overturn. For the first time, we couple thermo-chemical convection simulations with a thermal evolution model of the LMO crystallization to self-consistently address the mechanical interplay between the crystallization of the LMO and the fate of HTC. We introduce the model setup in Section 2 (in particular the additions and differences with respect to the model of Maurice et al. (2020)). The results concerning the efficiency of the HTC overturn are described in Section 4, and Section 5 elaborates on more far-reaching implications for the lunar evolution.

## 2. Layered Structure of the LMO Cumulates

### 2.1. LMO Crystallization

We use the same crystallization sequence of the LMO as Maurice et al. (2020), whose calculation is described in Schwinger and Breuer (2022). It results from an initially 1,000-km-deep magma ocean undergoing pure fractional crystallization, producing a 44-km-thick crust by flotation of anorthosites. The incompatibility of Fe and Ti also results in their enrichment in the late crystallized (shallow) layers, where they enter clinopyroxene and Fe-Ti-oxides, increasing the cumulates density (Figure 1a) and forming the HTC. A melt fraction of 0.5% is retained in the solid mantle, providing a reservoir for heat-producing elements, which are otherwise considered perfectly incompatible. This results in a typical fractionation abundance profile, with an effective mineral-melt partition coefficient of 0.05 (Figure 1b). The temperature at which each cumulate layer is formed is also an output of the crystallization model and provides the thermal initial condition for our mantle convection simulations (Figure 1c). Below 1,000 km depth, we consider an accretion-like temperature profile in the unmolten lower mantle. This is followed by the crystallization temperature of the cumulates up to ~110 km. Above this depth, corresponding to the initial depth of the LMO for our simulations, we consider an isothermal profile (vertical orange line in Figure 1c) characterizing the LMO. The LMO evolution is calculated self-consistently with our model (see Section 3.3 for details). The LMO shrinks over time because of the downward growth of the plagioclase crust and the upward rise of the solidification front (horizontal dotted black line in Figure 1c), which follows the crystallization temperature (blue line in Figure 1c from 110 km depth up to the surface).

### 2.2. 3-Layer Lunar Mantle Model

The density and rheology of the HTC have a first-order influence on both their overturn and on the possibility that an overturned layer is re-mobilized and participates in mantle convection over a longer timescale. Such a
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Initial conditions obtained from a model of lunar magma ocean fractional crystallization. Blue lines denote density (a), heat production (b), and crystallization temperature (c) profiles resulting from this model (see text and Maurice et al. (2020) for details). The top 44 km correspond to the plagioclase crust. The orange dashed lines in panels (a) and (b) indicate volume-averaged values used in the simulations and listed in Table 1. In panel (c), the orange line is the initial temperature profile characterized by an accretion-like shape in the unmolten lower mantle below 1,000 km depth, by the crystallization temperature of the cumulates up to 110 km, which is the depth at which the plagioclase crust starts to grow and marks the initial position of the crystallization front (horizontal dashed line), and an isothermal magma ocean up to the surface (see text and Section 3.3 for details).</figcaption></figure>

re-mobilization, being driven by the cumulates internal heat production rate, also constitutes an essential quantity to monitor. We thus track the presence of HTC in the convecting mantle and scale density, viscosity and heat production accordingly. Similar to previous studies (Li et al., 2019; Yu et al., 2019; Zhao et al., 2019), we approximate the sequence resulting from the fractional crystallization of the LMO (blue lines in Figure 1) with a three-layer profile for the density and heat production. We adopt the same model used by Yu et al. (2019) (to which we refer for a detailed description) consisting of mantle, HTC, and crust as represented by the orange dashed line in Figures 1a and 1b. The transition between these three layers occurs at the largest density gradients. The values of density and internal heat production for each layer are provided in Table 1.

Our crystallization model does not constrain the viscosity of the cumulates (which is controlled not only by the composition, but also by the water content, the grain size, the creep regime etc.). Therefore, to account for the weakening induced by the presence of HTC, following Yu et al. (2019), we multiply the intrinsic viscosity of the material initially present in this layer by a factor between $10^{-1}$ and $10^{-3}$. For comparison, pure ilmenite has been shown to be up to four orders of magnitude weaker than olivine (Dygert et al., 2016). The Fe-Ti oxides content in the HTC predicted by our crystallization model does not exceed a few tens of wt.%. Therefore, we do not expect the viscosity of this layer to be as low as that of pure ilmenite.

# 3. Numerical Implementation

## 3.1. Mantle Convection

We solve the non-dimensional conservation equations of mass (1), linear momentum (2), thermal energy (3), and composition (4) in a 2-D cylindrical grid using the finite-volume mantle convection code GAIA (Hüttig et al., 2013):

<caption>**Table 1**</caption>

*Initial Depth, Average Density and Heat Production of the Three Layers Used to Approximate the Lunar Magma Ocean Crystallization Sequence*

<table><thead><tr><th>Layer</th><th>Top and bottom depth (km)</th><th>Density ρ<sub>i</sub> (kg/m<sup>3</sup>)</th><th>Internal heating h<sub>i</sub> (W/kg)</th></tr></thead><tbody><tr><td>Crust (i = 3)</td><td>0–44</td><td>2,715</td><td>1.01 × 10<sup>-10</sup></td></tr><tr><td>HTC (i = 2)</td><td>44–70</td><td>3,775</td><td>1.27 × 10<sup>-10</sup></td></tr><tr><td>Mantle (i = 1)</td><td>70–1,350</td><td>3,206</td><td>1.43 × 10<sup>-11</sup></td></tr></tbody></table>
$$ \vec{\nabla} \cdot \vec{v}' = 0, \qquad (1) $$

$$ \vec{\nabla} \cdot (\eta'(\vec{\nabla}\vec{v}' + \vec{\nabla}\vec{v}^T)) - \vec{\nabla}p' = (\mathrm{Ra}T' - \mathrm{Ra}_c \sum_{i=1}^{3} \Delta\rho'_i c_i) g'\vec{e}_r, \qquad (2) $$

$$ \frac{\partial T'}{\partial t'} + \vec{v}' \cdot \vec{\nabla} T' - \nabla^2 T' = \frac{\mathrm{Ra}_H}{\mathrm{Ra}} \sum_{i=1}^{3} h'_i c_i, \qquad (3) $$

$$ \frac{\partial c_i}{\partial t} + \vec{v}' \cdot \vec{\nabla} c_i = 0, \quad i = 1, 2, 3, \qquad (4) $$

where $\vec{v}'$ is the velocity vector, $\eta'$ is the viscosity, $p'$ is the dynamic pressure, $T'$ is the temperature, $Ra$, $Ra_c$, and $Ra_H$ are the thermal, compositional and internal heating Rayleigh numbers, respectively, $g'$ is the non-dimensional radial profile of gravity (see below), $\Delta\rho'_i$ and $h'_i$ ($i = 1, 2, 3$) corresponding to the different layers given in Table 1) are non-dimensional density anomalies and internal heating rates, respectively, and $c_i$ are the volume fractions of each of the three components, which always add up to one. This constraint is enforced in the initial condition, and the divergence-free nature of the flow ensures that it is satisfied over time. All quantities in Equations 1–4 are non-dimensional. All dimensional variables and parameters are reported in Table 2. The scalings used to adimensionalize them are listed in Table 3. The three Rayleigh numbers are defined as follows:

$$ \mathrm{Ra} = \frac{\alpha \rho_1 g_0 \Delta T D^3}{\eta_{\mathrm{ref}} \kappa}, \qquad (5) $$

$$ \mathrm{Ra}_c = \frac{(\rho_2 - \rho_1) g_0 D^3}{\eta_{\mathrm{ref}} \kappa}, \qquad (6) $$

$$ \mathrm{Ra}_H = \frac{\rho_1^2 g_0 \alpha h_1 D^5}{\eta_{\mathrm{ref}} \kappa k}, \qquad (7) $$

where $\alpha$ is the thermal expansivity, $g_0$ is the gravitational acceleration at the surface, $\kappa$ is the thermal diffusivity, $k$ is the thermal conductivity, $D$ is the thickness of the lunar mantle, $\eta_{\mathrm{ref}}$ is the reference viscosity, and $\Delta T$ is the reference temperature contrast.

The viscosity is calculated according to a pressure- and temperature-dependent Ahrrenius law, including a composition-dependent term:

$$ \eta(P, T, c_2) = \eta_{\mathrm{ref}} \exp\left(\frac{E^* + PV^*}{RT} - \frac{E^* + P_{\mathrm{ref}}V^*}{RT_{\mathrm{ref}}}\right) (\Delta\eta_{\mathrm{HTC}})^{c_2}, \qquad (8) $$

where $\eta_{\mathrm{ref}}$ is the reference viscosity, $P$ is the lithostatic pressure, $E^*$ and $V^*$ are the activation energy and volume of diffusion creep, respectively, $R$ is the gas constant, $T_{\mathrm{ref}}$ and $P_{\mathrm{ref}}$ are the reference temperature and pressure at which $\eta_{\mathrm{ref}}$ is attained, $\Delta\eta_{\mathrm{HTC}}$ is a non-dimensional weakening factor of HTC, and $c_2$ is the HTC volume fraction.

In order to account for the non-constant gravity through the lunar mantle (because of the small lunar core), following Maurice et al. (2020), we calculate the radial profile of gravity assuming a homogenous mantle density equal to $\rho_1$:

$$ g(r) = \frac{4/3\pi G}{r^3} (\rho_c R_c^3 + \rho_1 (r^3 - R_c^3)), \qquad (9) $$

where $r$ is the radius, $G$ is the universal gravitation constant, and $R_c$ and $\rho_c$ are the core radius and density, respectively.

The top boundary of the computational domain where Equations 1–4 are solved corresponds to the bottom of the LMO while this is extant, and to the actual lunar surface afterward (the crust being added to the computational
<table><caption>**Table 2**
Model Variables and Parameters</caption><thead><tr><th>Symbol</th><th>Quantity</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>v̄</td><td>Velocity</td><td>Variable</td><td>m/s</td></tr><tr><td>η</td><td>Viscosity</td><td>Equation 8</td><td>Pa s</td></tr><tr><td>p</td><td>Dynamic pressure</td><td>Variable</td><td>Pa</td></tr><tr><td>P</td><td>Lithostatic pressure</td><td>Variable</td><td>Pa</td></tr><tr><td>T</td><td>Temperature</td><td>Variable</td><td>K</td></tr><tr><td>t</td><td>Time</td><td>Variable</td><td>s</td></tr><tr><td>c<sub>i</sub></td><td>Volume fraction of <i>i</i> material</td><td>Variable</td><td>—</td></tr><tr><td>ρ<sub>i</sub></td><td>Density of <i>i</i> material</td><td>Table 1</td><td>kg/m<sup>3</sup></td></tr><tr><td>h<sub>i</sub></td><td>Internal heating of <i>i</i> material</td><td>Table 1</td><td>W/kg</td></tr><tr><td>Δη<sub>HTC</sub></td><td>HTC weakening factor</td><td>10<sup>0</sup>–10<sup>−3</sup> (study parameter)</td><td>—</td></tr><tr><td>Ra</td><td>Rayleigh number</td><td>Equation 5</td><td>—</td></tr><tr><td>Ra<sub>c</sub></td><td>Compositional Rayleigh number</td><td>Equation 6</td><td>—</td></tr><tr><td>Ra<sub>H</sub></td><td>Internal heating Rayleigh number</td><td>Equation 7</td><td>—</td></tr><tr><td>α</td><td>Thermal expansivity</td><td>3 × 10<sup>−5</sup></td><td>1/K</td></tr><tr><td>ΔT</td><td>Temperature scale</td><td>1,930</td><td>K</td></tr><tr><td>D</td><td>Mantle thickness</td><td>1,350</td><td>km</td></tr><tr><td>η<sub>ref</sub></td><td>Reference viscosity</td><td>10<sup>20</sup>–10<sup>21</sup> (study parameter)</td><td>Pa s</td></tr><tr><td>κ</td><td>Thermal diffusivity</td><td>10<sup>−6</sup></td><td>m<sup>2</sup>/s</td></tr><tr><td>g<sub>0</sub></td><td>Surface gravity</td><td>1.62</td><td>m/s<sup>2</sup></td></tr><tr><td>k</td><td>Thermal conductivity</td><td>4</td><td>W/m/K</td></tr><tr><td>E*</td><td>Activation energy</td><td>150, 335 (study parameter)</td><td>kJ/mol</td></tr><tr><td>V*</td><td>Activation volume</td><td>4 × 10<sup>−6</sup></td><td>m<sup>3</sup>/mol</td></tr><tr><td>R</td><td>Gas constant</td><td>8.314</td><td>J/mol/K</td></tr><tr><td>P<sub>ref</sub></td><td>Reference pressure</td><td>3</td><td>GPa</td></tr><tr><td>T<sub>ref</sub></td><td>Reference temperature</td><td>1,600</td><td>K</td></tr><tr><td>ρ<sub>c</sub></td><td>Core density</td><td>7,200</td><td>kg/m<sup>3</sup></td></tr><tr><td>c<sub>P,c</sub></td><td>Core heat capacity</td><td>840</td><td>J/K/kg</td></tr><tr><td>R<sub>c</sub></td><td>Core radius</td><td>390</td><td>km</td></tr></tbody></table>

domain after full crystallization of the LMO). At the top boundary, we set free-slip boundary conditions and a prescribed temperature, corresponding to the LMO temperature while this is extant, and to 250 K afterward. For the bottom boundary (i.e., the CMB), we also use free-slip conditions and prescribe the core temperature, which we calculate at each timestep considering the secular cooling of an isothermal core obtained from the heat flux at the CMB, that is:

<table><caption>**Table 3**
Non-Dimensional Quantities and Corresponding Scaling</caption><thead><tr><th>Quantity</th><th>Scaling</th></tr></thead><tbody><tr><td>v̄'</td><td>D/κ × v̄</td></tr><tr><td>T'</td><td>(T - T<sub>s</sub>)/ΔT</td></tr><tr><td>p'</td><td>D<sup>2</sup>/(κη<sub>ref</sub>) × p</td></tr><tr><td>η'</td><td>1/η<sub>ref</sub> × η</td></tr><tr><td>t'</td><td>κ/D<sup>2</sup> × t</td></tr><tr><td>Δρ<sub>i</sub>'</td><td>(ρ<sub>i</sub> - ρ<sub>1</sub>)/(ρ<sub>2</sub> - ρ<sub>1</sub>)</td></tr><tr><td>h<sub>i</sub>'</td><td>h<sub>i</sub>/h<sub>1</sub></td></tr><tr><td>g'</td><td>g/g<sub>0</sub></td></tr></tbody></table>

$$ \rho_c c_{P,c} \frac{4}{3} \pi R_c^3 \frac{dT_c}{dt} = -4\pi R_c^2 k \left. \frac{\partial T}{\partial r} \right|_{r=R_c}, \quad (10) $$

where $c_{P,c}$ and $T_c$ are the heat capacity and temperature of the core, respectively.

All simulations are run on a 2-D cylindrical grid, which has been shown to reproduce well results for a 3-D sphere in the case of a Moon-like core size (Maurice et al., 2020), with 251 radial shells each containing 1,234 cells, corresponding to a radial resolution of 5.4 km and a lateral resolution ranging from 2.0 km at the CMB to 8.9 km at the surface.
## 3.2. Advection of Composition

Equation 4 controlling the advection of the composition is solved using the particle-in-cell method (Plesa et al., 2012), suitably adapted to track multiple phases. Each particle corresponds to one of the three phases, which are identified according to their initial position. As convection mixes the mantle, the volume fraction of the ith component in a given cell is calculated as the number of particles of type i divided by the total number of particles in that cell. This implementation makes the scheme potentially dependent on both the grid resolution and on the initial particles distribution density. We ran convergence tests on both these quantities on a reference case and found the effect to be negligible.

Clustering of particles due to the divergence of the interpolated velocity field can also be an issue (Pusok et al., 2016; Wang et al., 2015). The velocity field obtained by GAIA is divergence-free insofar that the sum of the fluxes computed at the cell walls of control volumes vanishes. However, the field obtained at arbitrary locations by interpolating velocities from grid points is, in general, not divergence-free, making a correction necessary. We therefore add a corrective term to a standard barycentric interpolation scheme, to make sure that the interpolated velocity field seen by particles is actually divergence-free.

## 3.3. LMO Thermal Evolution

The evolution of the LMO determines the position of the boundary as well as the value of the thermal boundary condition at the top of the solidifying mantle. It is calculated self-consistently following the procedure described in Maurice et al. (2020) to which we refer for all the details. Here we provide a summary.

The thermal evolution of the LMO is treated by ensuring heat conservation of the shrinking magma ocean. Before the formation of the plagioclase flotation crust, the LMO radiates heat to space and solidification is fast compared to the timescale for the onset of solid-state mantle convection. When plagioclase starts to crystallize, its floatation leads to the formation of a solid lid that insulates the interior and delays LMO cooling. For resolution reasons, we begin our simulations with a 5-km-thick plagioclase crust. At this stage, the LMO has a depth of ~110 km (Figure 1c). We balance the cooling flux by thermal conduction flux through the growing flotation crust (which we denote $f_{crust}$) with the conductive heat flux from the top of the cumulates pile, beneath the LMO ($f_{cond}$), and the heat piping flux generated upon the extraction of hot melts produced by pressure-release melting in the convecting mantle (Maurice et al., 2020) ($f_{HP}$), as well as the internal heat production due to the decay of long-lived radioactive isotopes $^{235}$U, $^{238}$U, $^{232}$Th, and $^{40}$K (expressed as an equivalent flux at the surface, $f_{rad}$).

Figure 2a schematically shows the setup of our model. The LMO (yellow region in Figure 2a) is treated as an isothermal shell whose internal radius corresponds to the top radius of the solid cumulates (light brown layer in Figure 2a). The outer radius of the LMO corresponds to the bottom radius of the crust (gray layer in Figure 2a). Below the LMO, we solve numerically Equations 1–4 in the upward-growing solid cumulates as described in Section 3.1. The conductive flux through the crust is obtained by solving the heat conduction equation in a downward-growing 3D spherical shell with a fixed surface temperature of 250 K and a bottom temperature equal to the LMO temperature. Notice that the assumption of spherical symmetry allows us to treat the crust in 3-D without significant time cost. We account for internal heat production both in the crust and in the solid cumulates (with heating rates given in Table 1). The cooling flux for the LMO is calculated as $f_{crust} = k_{crust}dT/dr$ ($r = R_{crust}$), where $k_{crust}$ is the thermal conductivity of the crust, and the temperature gradient is evaluated at the bottom of the crust.

The conductive heat flux from the underlying cumulates into the LMO is $f_{cond} = -k_{cum}dT/dr$ ($r = R_{MO}$), where $k_{cum}$ is the heat conductivity of the cumulates (i.e., of the mantle), and the temperature gradient is calculated at the top of the cumulates, from the temperature field calculated by the mantle convection model. The heat piping flux is obtained by calculating the energy necessary to buffer the temperature at the solidus at each timestep, assuming that this energy is consumed by latent heat during melting. Finally, the internal heating rate accounts for the partitioning of heat producing elements in the melt as the LMO solidifies, considering that they are perfectly incompatible and only enter the solid cumulates due to trapped interstitial melt.

## 3.4. Parameter Space

We vary the reference viscosity $\eta_{ref}$ between $10^{20}$ and $10^{21}$ Pa s to account for the uncertainty on the viscosity of the early lunar mantle. The thermal conductivity of the plagioclase crust, found by Maurice et al. (2020) to play a
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) High-Ti cumulates (HTC) overturn in the presence of a solidifying lunar magma ocean (LMO) or (b) after complete LMO solidification. The second is the traditional paradigm for HTC overturn (adopted e.g., by Li et al. (2019), Yu et al. (2019), and Zhao et al. (2019)), while the first is the new setup introduced in this study. On both panels, time advances clockwise, from the top quadrant representing the initial conditions, through the right one corresponding to the onset of the Rayleigh-Taylor instability due to the presence of shallow, dense HTC, and to the bottom quadrant, corresponding to a stage in which the cold stagnant lid has grown and engulfed—partly or fully—the HTC. Note that layer thicknesses are not drawn to scale.</figcaption></figure>

central role for the lifetime of the LMO, is held constant at 2 W/m/K, the fiducial value adopted in Maurice et al. (2020). The crust keeps its low thermal conductivity also after the end of LMO solidification, thereby affecting the long-term thermal evolution of the solid mantle (Ziethe et al., 2009).

Similar to Yu et al. (2019), the distribution of heat-producing elements as well as the density profile are obtained from the crystallization model described in Section 2.2. We thus keep them constant, and focus on the role of the rheology. As in Yu et al. (2019), we vary the value of the activation energy (E*) using either 335 kJ/mol (the value for diffusion creep previously adopted), or 150 kJ/mol, a decreased value aimed at mimicking the effect of dislocation creep (see e.g., Schulz et al. (2020)).

We investigate the role of HTC rheology by varying the viscosity contrast Δη<sub>HTC</sub> between 1 and 10<sup>-3</sup>, the lowest value approaching the four orders of magnitude reduction found by Dygert et al. (2016) for pure ilmenite rheology compared to dry olivine, which is a lower bound for HTC where Fe-Ti-oxides are diluted.

To assess the influence of the presence of the magma ocean on the overturn of HTC, we run simulations either starting from a completely solidified lunar mantle (the traditional scenario as e.g., in Li et al. (2019), Yu et al. (2019) and Zhao et al. (2019)), or coupled with the LMO thermal evolution model. In the first case, the initial temperature follows the crystallization temperature (solid blue line in Figure 1c), that is, we assume that the cumulates retain their solidification temperature during the whole solidification of the magma ocean. In the second case, each newly crystallized cumulates shell starts from its solidification temperature, and the LMO temperature (i.e., the top thermal boundary condition) is set according to its depth.

**4. HTC Overturn During LMO Crystallization**

We run simulations of the dynamical evolution of the lunar mantle during and after LMO crystallization, focusing on the fate of HTC. We compare cases where the LMO crystallization is virtually instantaneous—as previously assumed by for example, Yu et al. (2019) and Zhao et al. (2019)—with cases where the LMO is long-lived (self-consistently computed as in Maurice et al. (2020)), and solid-state convection starts in the cumulates underlying the solidifying LMO. We track the evolution of the volume fraction of overturned HTC (Ψ<sub>HTC</sub>), which is defined as follows (Yu et al., 2019):
$$ \Psi_{\text{HTC}} = \frac{V_{\text{HTC}} - V'_{\text{HTC}}}{V_{\text{HTC}}} \times 100, \qquad (11) $$

where $\Psi_{\text{HTC}}$ is expressed in percent, $V_{\text{HTC}}$ is the total volume of HTC, and $V'_{\text{HTC}}$ is the volume of cumulates that has not yet been mobilized. In this section, we discuss the influence of each parameter on the efficiency of the HTC overturn in light of the time-series of $\Psi_{\text{HTC}}$. The final value, $\Psi_{\text{HTC},\infty}$ (reached at the end of the time-evolution), is of particular interest to quantify the overturn efficiency.

## 4.1. Influence of the LMO on the Overturn

Zhao et al. (2019) suggested that the presence of a weak layer (i.e., the LMO) at the top of the solid cumulates pile would enhance the overturn by mechanically decoupling the convecting mantle from the cold and stiff crust. They accounted for this effect by using a temperature-dependent rheology and introducing a hot layer with low viscosity overlying the HTC. Here we perform the first simulations where this interface is modeled self-consistently by using a free-slip boundary condition at the top of the cumulates pile and a continuous temperature profile, as we generally do not expect the LMO to be hotter than the underlying cumulates. As shown in Figure 3, our results largely support the idea that the overturn is facilitated by the presence of the LMO, in particular for a high mantle reference viscosity $\eta_{\text{ref}} = 10^{21}$ Pa s. When using this value (Figure 3a), the amount of overturned HTC is negligible ($\Psi_{\text{HTC},\infty} < 1\%$) for all cases that do not account for the presence of the LMO but one (thick green line), while it becomes significant for the corresponding cases considering a crystallizing LMO. Decreasing the activation energy and, to a lower extent, decreasing the HTC weakening factor $\Delta\eta_{\text{HTC}}$ also favor the overturn. For $E^* = 150$ kJ/mol, $\Psi_{\text{HTC},\infty}$ exceeds ~50% for all values of $\Delta\eta_{\text{HTC}}$ when the presence of a crystallizing LMO is accounted for, while for $E^* = 335$ kJ/mol, only the weakest HTC yield a significant overturn. The case $E^* = 150$ kJ/mol and $\Delta\eta_{\text{HTC}} = 10^{-3}$ could not be run past LMO crystallization, and we consider $\Psi_{\text{HTC}}$ at the end of LMO crystallization for $\Psi_{\text{HTC},\infty}$. This could underestimate $\Psi_{\text{HTC},\infty}$, making our observed trend conservative.

For the low value of the reference viscosity ($\eta_{\text{ref}} = 10^{20}$ Pa s, Figure 3b), $\Psi_{\text{HTC},\infty}$ approaches 100% regardless of the presence of the LMO for $E^* = 150$ kJ/mol. For $E^* = 335$ kJ/mol, the LMO plays a crucial role, allowing for extensive overturn ($\Psi_{\text{HTC},\infty} > 50\%$), while $\Psi_{\text{HTC},\infty} \sim 0\%$ when LMO is absent, for all values of $\Delta\eta_{\text{HTC}}$ (except for $\Delta\eta_{\text{HTC}} = 10^{-3}$, for which the case without LMO could not be run). Notice that the case $E^* = 335$ kJ/mol and $\Delta\eta_{\text{HTC}} = 10^{-3}$ could not be run past the LMO crystallization. However its $\Psi_{\text{HTC}}$ already reached ~100% by the end of LMO crystallization, so its subsequent evolution is not relevant.

Figure 4 summarizes the above results showing the value of $\Psi_{\text{HTC},\infty}$ across the whole parameter space for cases that consider a crystallizing LMO (Figure 4a) and cases that do not (Figure 4b). All values of $\Psi_{\text{HTC},\infty}$ are reported in Tables A1 and A2.

## 4.2. Morphology of the Overturn

When the crystallization of the LMO is considered, for $\eta_{\text{ref}} = 10^{21}$ Pa s, two different mechanisms operate depending on the value of the activation energy. For $E^* = 150$ kJ/mol, solid-state convection is first triggered by a Rayleigh-Taylor instability that grows when HTC start to crystallize after ~100 Myr. This situation is illustrated in Figure 5 for the case $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 150$ kJ/mol, and $\Delta\eta_{\text{HTC}} = 10^{-3}$. The overturn always begins with small-scale diapirs, which join to form five downwellings as convection becomes also thermally driven. The flow wavelength then decreases as new HTC diapirs begin to sink, particularly for low values of $\Delta\eta_{\text{HTC}}$ (i.e., as HTC are made weaker). For these cases, the onset time of convection (which in the following we will indicate with $\tau_{\text{conv}}$) decreases with $\Delta\eta_{\text{HTC}}$ (see Table A1).

Different from the case with $E^* = 150$ kJ/mol, for $E^* = 335$ kJ/mol (and $\eta_{\text{ref}} = 10^{21}$ Pa s), thermally driven convection starts as early as after 31 Myr (see Table A1), long before HTC begin crystallizing. Therefore, the Rayleigh-Taylor instability accompanying HTC crystallization occurs in an already established flow, whose wavelength is set by Rayleigh-Bénard convection. Actually in this case, only for the weakest HTC ($\Delta\eta_{\text{HTC}} = 10^{-3}$) does an overturn occur. The overturn then induces a transient reorganization of the flow, which is interrupted by the end of LMO crystallization. This case is shown in Figure 6 for $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 335$ kJ/mol and $\Delta\eta_{\text{HTC}} = 10^{-3}$. The later evolution of the overturn is heavily hampered by the presence of the crust, which effectively sets a strong shear stress at the top of the HTC layer, which prevents its efficient mobilization.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time-series of Ψ<sub>HTC</sub> for η<sub>ref</sub> = 10<sup>21</sup> Pa s (a) and η<sub>ref</sub> = 10<sup>20</sup> Pa s (b), comparing cases including a crystallizing lunar magma ocean (LMO) (thin lines) and cases starting after LMO solidification (thick lines). Circles represent the end of LMO solidification, corresponding to the starting time of the simulations that do not account for LMO crystallization. Colors denote different viscosity contrasts Δη<sub>HTC</sub> from 1 (i.e., no viscosity contrast between high-Ti cumulates (HTC) and mantle) to 10<sup>-3</sup>. Solid and dashed lines correspond to activation energies of 150 and 335 kJ/mol, respectively. The cases with LMO with η<sub>ref</sub> = 10<sup>21</sup> Pa s, E* = 150 kJ/mol and Δη<sub>HTC</sub> = 10<sup>-3</sup>, and η<sub>ref</sub> = 10<sup>20</sup> Pa s, E* = 335 kJ/mol and Δη<sub>HTC</sub> = 10<sup>-3</sup>, respectively, did not run past LMO crystallization, and we prolong their final value of Ψ<sub>HTC</sub>. The cases without LMO with η<sub>ref</sub> = 10<sup>20</sup> Pa s, Δη<sub>HTC</sub> = 10<sup>-3</sup> and E* = 150 and 335 kJ/mol could not be run. Extreme viscosity gradients was the reason why these various cases could not be run.</figcaption></figure>

For the lower value of the reference viscosity (η<sub>ref</sub> = 10<sup>20</sup> Pa s), thermal convection starts early, prior to HTC crystallization, regardless of the value of the activation energy (and τ<sub>conv</sub> is accordingly independent of Δη<sub>HTC</sub>, see Table A1). The corresponding wavelength of the flow is similar to the cases with η<sub>ref</sub> = 10<sup>21</sup> Pa s (six downwellings for E* = 150 kJ/mol, eight for E* = 335 kJ/mol). The HTC overturn, which is significant for all cases, is initiated along pre-existing cold downwellings, but rapidly disrupts the pattern established by the thermal instability, and makes the flow highly time-dependent for E* = 150 kJ/mol, with a tendency to promote short-wavelength features upon decreasing the strength of the HTC.

In general, while the overturn of HTC alters the flow pattern, the pre-existence of thermally driven convection structures in most cases affects the growth of the Rayleigh-Taylor instability associated with the crystallization of HTC. Therefore, calculations of overturn wavelength based on classical scaling of Rayleigh-Taylor instability, which have often been invoked in the literature (Evans & Tikoo, 2022; Hess & Parmentier, 1995), may fail to capture the complexity of the overturn. Even though the flow pattern is affected by the 2-D geometry, a similar mechanism will likely operate also in 3-D.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Final value of ΨHTC (ΨHTC,∞) for the whole parameter space investigated in this study. Panel (a) corresponds to cases where lunar magma ocean crystallization is accounted for and panel (b) to those where it is ignored. Lighter tones correspond to more extensive high-Ti cumulates overturn. The two cases in panel (b) with ηref = 1020 Pa s and ΔηHTC = 10−3 could not be run.</figcaption></figure>

**4.3. Influence of Mantle Rheology**

The mantle rheology is controlled by the reference viscosity and activation energy, the pressure span in the lunar mantle being too small to significantly affect the viscosity. To isolate the effect of these rheological parameters, for cases accounting for convection during the LMO crystallization, we turn again to Figure 4, which shows the values of ΨHTC,∞ across the parameter space, and in particular as a function of ηref, E*, and ΔηHTC. The effect of the reference viscosity is straightforward: an overall weaker mantle facilitates the overturn. The effect of the activation energy (i.e., the temperature dependence of the viscosity) is also clear: decreasing E* increases the overturn efficiency. This is due to the fact that the reference temperature (Tref = 1,600 K, at which the effective viscosity coincides with the reference viscosity) is higher than the temperature at which HTC crystallize. Increasing the temperature-dependence of the viscosity increases the stiffness of the mantle wherever it is colder than Tref, in particular that of the crystallizing HTC, thereby preventing their mobilization.

In order to investigate the impact of a weak HTC rheology on the overturn (Dygert et al., 2016), we also vary the HTC weakening factor (ΔηHTC) between 1 and 10⁻³. This is illustrated by the vertical axis of Figure 4. The trend is clear and confirms the expectations: the weaker the HTC are, the more prominent is their overturn. For the cases with ηref = 10²¹ Pa s, ΨHTC,∞ > 50% when ΔηHTC = 10⁻³. By contrast, the overturn is nearly negligible for ΔηHTC ≥ 10⁻². For the cases with ηref = 10²⁰ Pa s, ΨHTC,∞ consistently increases with decreasing ΔηHTC.

**4.4. Influence of the Overturn on the LMO Lifetime**

In Maurice et al. (2020) we investigated the influence of convection and secondary melting on the LMO lifetime, yet considering a purely thermal model. We found that the heat-piping flux—that is, heating of the LMO due to partial melting associated with mantle convection—could significantly prolong the lifetime of the LMO. Here we additionally consider the effect of compositional heterogeneities (i.e., the role of Rayleigh-Taylor instability). As these can have a significant impact on the mantle flow, we expect the LMO cooling chronology to be also impacted. As long as HTC have not started to crystallize, convection in the solid cumulates, if ongoing, is purely thermal (in fact, we assume the mantle beneath the HTC to have a homogeneous composition and hence a uniform density). Nevertheless, as also shown above, solid cumulates may start convecting before the crystallization of HTC (Figure 6a). In such cases, the onset of solid-state convection in the cumulates is thermally driven, caused by the unstable thermal profile of the LMO crystallization temperature (Maurice et al., 2017, 2020). Except for cases with ηref = 10²¹ Pa s and E* = 335 kJ/mol, the overturn of HTC, if it occurs, starts after the onset of thermal convection, and is associated with a later increase in the intensity of mantle flow.

Figure 7 shows the time-series of the depth of the LMO and of the bottom of the crust, as well as the evolution of the heat budget of the LMO for the two cases illustrated in Figures 5 and 6: ηref = 10²¹ Pa s, ΔηHTC = 10⁻³, and E* = 335 kJ/mol (top row) and E* = 150 kJ/mol (bottom row). In the first case, the effect of the onset of mantle
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Snapshots of the temperature field (left column) and high-Ti cumulates (HTC) content (right column) at three successive times during the lunar magma ocean crystallization (a, b after 100 Myr; c, d after 150 Myr; e, f after 300 Myr), for the case $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 150$ kJ/mol, and $\Delta\eta_{\text{HTC}} = 10^{-3}$. Here, mantle convection is triggered by the Rayleigh-Taylor instability of the HTC layer (a, b). For better visibility, the colormap of the temperature field is clipped at 1,000 K, and the colormap of the HTC content is saturated at 0.1.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Snapshots of the temperature field (left column) and high-Ti cumulates (HTC) content (right column) at three successive times during the lunar magma ocean crystallization (a, b after 50 Myr; c, d after 200 Myr; e, f after 250 Myr), for the case $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 335$ kJ/mol, and $\Delta\eta_{\text{HTC}} = 10^{-3}$, that is, as in Figure 5 but for a higher activation energy. Here, mantle convection is triggered by Rayleigh-Bénard instability (a, b), long before the HTC layer starts to crystallize (c, d). HTC content and temperature colormaps are clipped as in Figure 5.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Time-series of the depth of the bottom of the lunar magma ocean (blue curve) and of the crust (orange curve) (a, c), and of the heat fluxes (b, d): conductive cooling through the crust (blue curve, $f_{crust}$), conductive heating from the top of the cumulates pile (orange curve, $f_{cond}$), internal heating by radioactive decay (green curve, $f_{rad}$), and heat piping flux (red curve, $f_{HP}$). Both cases have $\eta_{ref} = 10^{21}$ Pa s and $\Delta\eta_{HTC} = 10^{-3}$; the case on the top row has $E^* = 335$ kJ/mol (same case as in Figure 6), while the case on the bottom row has $E^* = 150$ kJ/mol (same case as in Figure 5). Downward spikes in $f_{HP}$ in panel (d) are artifacts corresponding to simulation restarts that do not affect the evolution.</figcaption></figure>

convection on the evolution of the LMO becomes relevant early on, after less than 20 Myr. As already pointed out by Maurice et al. (2020) for a similar case (also prior to HTC crystallization), it slows down the LMO solidification. Here, HTC start crystallizing and overturning via Rayleigh-Taylor instability after ~180 Myr. The return flow associated with the overturn causes mantle melting and is accompanied by a heat-piping flux (red curve in Figure 7b), which heats the LMO causing it to slightly thicken (see blue and orange curves in Figure 7a at ~200 Myr).

For this combination of parameters, the late increase in the heat piping flux prolongs the LMO lifetime by about 100 Myr compared to purely thermal cases that neglect the effect of HTC overturn. In other cases, the role of heat-piping can be still more pronounced. In the case shown in Figures 7c and 7d, mantle convection is triggered by the Rayleigh-Taylor instability associated with HTC crystallization and overturn starting after ~100 Myr. Since the mantle is stiffer, the HTC layer needs to grow thicker before it can be destabilized. Instead of a progressive overturn of the HTC layer as it crystallizes, the overturn takes place as one intense episode, occurring when the LMO is insulated below an already thickened crust (see Figure 7c and videos during the LMO solidification in Maurice et al. (2023b)). The LMO lifetime is mainly controlled by the balance between heat-piping flux and conduction through the crust (red and blue curves in Figures 7b and 7d, respectively), which in turn depends on the thickness of the crust. The late onset of heat-piping (Figures 7c and 7d) induces a delay in heat extraction through the crust (compare with Figures 7a and 7b where $f_{crust}$ overshoots the peak in $f_{HP}$), leading to a significant re-thickening of the LMO and eventually a much longer lifetime (~400 Myr).

Figure 8 shows the lifetime of the magma ocean over the whole parameter space investigated. In general, we find that cases with a low reference viscosity feature efficient, small-scale overturns, that trigger only limited return
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Lifetime of the lunar magma ocean (LMO) over the complete ($\eta_{\text{ref}} - E^* - \Delta\eta_{\text{HTC}}$) parameter space investigated in this study.</figcaption></figure>

flow and moderate (although not negligible) heat-piping, while the larger-scale overturn (if any) associated with higher reference viscosity provides much more secondary melt that re-heat the LMO. However, the influence of the various parameters on the LMO lifetime is not systematic and is made quite complex by the significant time-dependence of the flow that accompanies the HTC overturn. We find that the longest LMO lifetime (close to 400 Myr) is reached for the stiffest cases with $\eta_{\text{ref}} = 10^{21}$ Pa s and $E^* = 150$ kJ/mol.

Using numerical simulations, two recent studies have suggested that the formation of the SPA basin could have induced a large-scale reorganization of mantle flow, ultimately responsible for the formation of the lunar compositional asymmetry (Jones et al., 2022; Zhang et al., 2022). These models, however, assume that the LMO crystallization was completed at the time of SPA formation (~4.3 Ga according to Evans et al. (2018)). On the one hand, according to our results, at this time, if the LMO solidification was completed, a significant part of the HTC may have already overturned. On the other hand, if the LMO was still solidifying, as predicted by most of our simulations, the SPA-forming impact would have occurred on a Moon that was still hosting an

internal magma ocean, a process whose consequences remain to be investigated in detail with dedicated simulations (Miljković et al., 2021).

# 5. Influence of the HTC overturn on the Long-Term Evolution of the Moon

The fate of the HTC layer is key to address several issues related to the long-term evolution of the lunar interior. In this section, we discuss the influence of the HTC overturn on the present-day thermal state of the Moon, comparing the predictions of our simulations with the temperature distribution inferred from geophysical and geodetic data (Section 5.1). We describe the limited mobilization of the overturned HTC (Section 5.2) and the extent to which their entrainment by mantle flow results in partial melting of mafic cumulates mixed with HTC that may be related to the generation of Ti-rich mare basalts (Section 5.3). We show the influence of the HTC overturn on the CMB heat flux that contributes to drive convection and dynamo generation in the core (Section 5.4). Finally, we discuss the formation of melt atop the CMB in light of the inferred existence of a deep-seated, partially molten layer (Section 5.5).

## 5.1. Present-Day Thermal State

Constraints of the present-day thermal state of the Moon have been obtained based on seismic and gravity data (Khan et al., 2007, 2014). Figure 9 compares these estimates with the results of our simulations after 4.5 Gyr, both with and without accounting for LMO crystallization and its influence on the HTC overturn. As already observed, a weaker rheology allows for a more extensive overturn. Since it also allows for a faster cooling (compare left and right columns of Figure 9), one could expect a correlation between overturn and present-day thermal state. However, the important role of the HTC weakening makes this picture more complicated. We observe that cases where the HTC overturn is extensive are characterized by a hotter lower mantle, and a colder upper mantle. Furthermore, cases that account for the solidification of the LMO generally have a slightly warmer upper mantle (compare top and bottom rows of Figure 9). Our results, while on the warm side, compare well with the estimates of Khan et al. (2007), and overlap with estimates from Khan et al. (2014) in the mid mantle. The range of CMB temperatures reached by our simulations is also in agreement with recent alternative estimates by Briaud et al. (2023) (~1,500–1,850 K), lower than those by Khan et al. (2014).

## 5.2. Lack of Large-Scale Mobilization of the Overturned HTC

Whenever the overturn occurs, we do not observe a later global re-mobilization of the sunken HTC layer, which was reported by some previous modeling studies (Stegman et al., 2003; Zhang et al., 2013, 2017; Zhong et al., 2000). These studies did not simulate the overturn, but prescribed a post-overturn density stratification as initial condition. They focused on the evolution of overturned HTC while making assumptions on their thickness, density, heat-sources content, and rheology. Here, we self-consistently simulate both the HTC overturn and its
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Present-day thermal state of the Moon for $\eta_{\text{ref}} = 10^{21}$ Pa s (resp. $10^{20}$ Pa s) on the left column (resp. right column). In the top row, lines correspond to the present-day temperature profiles of cases that do account for the initial presence of a crystallizing lunar magma ocean (LMO) while the blue-shaded area corresponds to the envelope defined by the cases that do not account for the initial presence of the LMO. Conversely, in the bottom row, lines correspond to cases that do not account for the initial presence of the LMO while the pink shaded area indicates the envelope defined by the cases that do account for it. On all panels, the orange-shaded and hatched areas correspond to the range of present-day temperature of the lunar mantle inferred from seismic and gravity data by Khan et al. (2007) (K07) and Khan et al. (2014), respectively (K14).</figcaption></figure>

subsequent evolution, employing initial conditions based on a LMO crystallization model and thus avoiding the need to make the above assumptions. With respect to these previous works, we generally find a significantly denser overturned HTC layer as a result of inefficient mixing with the background mantle during the overturn, as well as of an intrinsically higher HTC density predicted by our crystallization model ($\Delta\rho = 574 \text{ kg/m}^3$ in our case compared to $300 \text{ kg/m}^3$ in Zhong et al. (2000), $<165 \text{ kg/m}^3$ in Stegman et al. (2003), and $90 \text{ kg/m}^3$ in Zhang et al. (2013), where $\Delta\rho$ is the density difference between HTC and mantle). This difference (partly due to the assumption originally made by Hess and Parmentier (1995) that dense cumulates dilute into underlying mafic cumulates, which we do not observe in our simulations) explains the overall lack of large-scale mobilization of the overturned HTC layer that we systematically observe. The density contrast obtained from our crystallization model is in agreement with the recent experiments of Schmidt and Kraettli (2022) if perfect separation between
plagioclase and mafic crystals is assumed. Yet, Schmidt and Kraettli (2022) also question the validity of this assumption and suggest that a cumulate density above 3,500 kg/m³ is “largely hypothetic.” This corresponds to Δρ = 294 kg/m³ and still on the upper limit of the range previously investigated.

The main driver of re-mobilization in the works of Stegman et al. (2003), Zhang et al. (2013, 2017), and Zhong et al. (2000) is internal heating by heat-producing elements concentrated in the HTC, which we also take into account (in the simulations of Stegman et al. (2003), the core can also play a role). While we do observe heating of the lower mantle, the induced positive thermal buoyancy never overcomes the effect of the negative compositional buoyancy of the overturned HTC (see Maurice et al., 2023b). Cases with E* = 150 kJ/mol exhibit a quiescent period of 1.5–2 Gyr for η<sub>ref</sub> = 10<sup>21</sup> Pa s and ~500 Myr for η<sub>ref</sub> = 10<sup>20</sup> Pa s, during which the lower mantle heats up (see the low rms velocity between 4 and 3 Ga for the corresponding cases in Figure S1 in Supporting Information S1), followed by the formation of long-lived convecting piles at the top of a static dense HTC layer (Figures 10e and 10f, and Maurice et al., 2023b). For E* = 335 kJ/mol, cold downwellings largely characterize the convection during the early stages following the LMO solidification, before basal heating by the overturned layer takes over as the main convection driver. The resulting structure is nevertheless similar, although the erosion of the overturned layer is more pronounced (yet still far from a large-scale re-mobilization). The case η<sub>ref</sub> = 10<sup>20</sup> Pa s and Δη<sub>HTC</sub> = 10<sup>-1</sup> is the only one where the overturn occurs mainly after the end of the LMO phase, and is too limited to produce a global layer at the CMB. Instead, two antipodal piles persist until present-day (Maurice et al., 2023b). Notice that our discretization of the crystallization sequence (Figure 1) is more likely to overestimate the heat production in HTC than their density, thus making our model conservative, and the lack of HTC re-mobilization a robust result.

Another effect that can enhance the buoyancy of the overturned HTC layer is its partial melting. U et al. (2023) carried out simulations of the evolution of the lunar mantle assuming the presence above the CMB of a dense HTC layer enriched in heat-producing elements. As the layer heats up and undergoes partial melting, its density decreases because, for a similar composition, liquid HTC, although denser than the mafic lunar mantle, are less dense than solid HTC. Thus partial melting of the HTC layer can have a net positive effect on its buoyancy. U et al. (2023) found that this effect, combined with melt migration, results in some degree of HTC mobilization toward the shallower mantle. Yet, the compressibility of HTC melts—an important parameter controlling this effect—remains to be constrained. We also notice that the initial HTC distribution chosen by U et al. (2023) make the basal layer more prone to be mobilized, with a concentration increasing exponentially with depth and reaching high values only very close to the CMB, unlike what our simulation self-consistently produces.

5.3. Re-Melting of Entrained HTC and Production of High-Ti Mare Basalts

As pointed out in the previous section, despite heating due to the enrichment in long-lived radioactive isotopes, the overturned HTC layer is never fully destabilized. Instead, its erosion causes the transport of Ti-rich materials in the bulk convecting mantle (Figures 10c and 10d) that can then undergo secondary (pressure-release) melting and thus contribute to the production of high-Ti mare basalts. From our simulations we estimate an upper bound for the mass of Ti-bearing melt produced by comparing the temperature in each grid cell that contains HTC with melting curves from Wyatt (1977) for HTC-analog material (29% ilmenite—71% clinopyroxene). Assuming a linear dependence of the HTC melt fraction on temperature between the HTC solidus (HTC melt fraction equal to 0) and HTC liquidus (HTC melt fraction equal to 1), we can estimate how much HTC melt is produced. However, it is unclear at which depth Ti-bearing melt can be buoyant enough to be extracted to the lunar surface and form high-Ti mare basalts (Mallik et al., 2019; van Kan Parker et al., 2012; Xu et al., 2022). HTC melt in the lower mantle is likely negatively buoyant (Mallik et al., 2019; Xu et al., 2022), otherwise the survival of a supposed partially molten layer atop the CMB until today is hard to explain (see Section 5.5). Therefore, we focus here on HTC melt produced within the pressure window of 1–2.5 GPa inferred from multiple saturation points of high-Ti mare basalt samples (Elkins-Tanton et al., 2002). Our estimate yields an upper bound because we do not account for mantle depletion in HTC material due to melt extraction. As we only calculate HTC melt production as a post-processing step, HTC parcels that have already experienced melting can be counted multiple times. While this simple approach potentially yields a large overestimate, it is still insightful as it provides an upper bound on how long and how much Ti-bearing melt can be produced.

In general, we observe that secondary melting of HTC-bearing materials post-LMO solidification is always fed by
the convecting piles at the top of the overturned layer, rather than HTC downwellings, which do not have time to
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Snapshots of the temperature (left column) and high-Ti cumulates (HTC) content (right column) at the end of lunar magma ocean crystallization (top row), after 1.1 Gyr (middle row), and 4.5 Gyr (bottom row), for the case $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 150$ kJ/mol and $\Delta\eta_{\text{HTC}} = 10^{-1}$. Overlay contours correspond to areas where the temperature is high enough for HTC to melt (solid white lines in panels b and d), and to the pressure window of inferred production of the parent magmas of high-Ti mare basalts (1–2.5 GPa, dashed orange lines in panels b, d and f). Notice that the downwelling plumes originating from the shallow HTC layer in panels (d) and (f) are trapped in the stagnant lid and are no longer involved in the overturn (which is long finished).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** High-Ti melt produced by secondary partial melting of high-Ti cumulates (HTC) that are entrained in the mantle after the end of the lunar magma ocean crystallization and overturn, expressed as thickness of a global equivalent layer (GEL) (panels a and c), and melt production rate expressed as GEL thickness per Myr (panels b and d). In panels (a) and (b), $\eta_{ref} = 10^{21}$ Pa s, while in panels (b) and (d), $\eta_{ref} = 10^{20}$ Pa s. Line colors and styles refer to different values of the HTC weakening factor $\Delta\eta_{HTC}$ and of the activation energy, respectively. Lines that do not appear in some panels (e.g., dashed lines other than green in panels a and b) plot below the range of the y-axis, corresponding to cases that produce no Ti-bearing melt.</figcaption></figure>

heat up to their melting temperature before they reach the CMB. Furthermore, at the end of the overturn, downwellings break up, leaving down-pointing plumes that get trapped in the cold stagnant lid, and no longer experience significant heating in spite of their heat-producing elements content.

Figure 11 shows the time-series of the cumulative (left column) and instantaneous (right column) amount of Ti-bearing melt produced between 1 and 2.5 GPa, represented as a global equivalent layer (GEL) thickness. The left column tracks the total quantity of Ti-bearing melts produced (again, providing an upper bound), while the right column tracks the lifetime of Ti-bearing melt production. The trends are very different depending on the mantle reference viscosity. For $\eta_{ref} = 10^{21}$ Pa s (top row), as already pointed out, a high activation energy generally correlates with the absence of HTC overturn, which prevents Ti-bearing material from participating in mantle convection and undergoing secondary melting, except for the most extreme HTC weakening ($\Delta\eta_{HTC} = 10^{-3}$). In this case, due to the higher activation energy, the convecting part of the mantle—which is hotter than $T_{ref}$—is weaker than for cases with $E^* = 150$ kJ/mol (see Equation 8). This results in an intense convection and better
mixing of entrained HTC material. The calculated cumulative Ti-bearing melt production in this case is several tens of km, exceeding the maximum of ~25 km corresponding to the complete extraction of HTC. This limit is exceeded because we neglect extraction of Ti-bearing melts. However, a more intense convection is also associated with faster cooling of the mantle, which acts to shorten the secondary melting phase, resulting in a similar duration of the high-Ti melt production period with respect to cases with the lower activation energy, extending it until 1.1 to 1.5 Ga. For cases with $E^* = 150$ kJ/mol, there is no clear trend between the HTC weakening factor $\Delta\eta_{HTC}$ and either the cumulative amount of Ti-bearing melt produced or the duration of Ti-bearing secondary melt production, due to complex trade-offs between HTC mixing in the mantle and its thermal evolution. The melt production time-series follow the evolution of the intensity of convection (see Figure S1 in Supporting Information S1), with a quiescent period around 3 Ga before heating of the lower mantle drives convective piles (see Maurice et al., 2023b).

For $\eta_{ref} = 10^{20}$ Pa s, the situation is somewhat different: here, cases with $E^* = 150$ kJ/mol cool down rapidly and fail to produce Ti-bearing secondary melt after LMO crystallization, except for $\Delta\eta_{HTC} = 10^{-3}$ where a small amount of Ti-bearing melt is produced just after the end of the LMO solidification. For $E^* = 335$ kJ/mol, there is still no trend between $\Delta\eta_{HTC}$ and the cumulative quantity of Ti-bearing melt produced or the duration of its production, although it is noteworthy that in the case $\Delta\eta_{HTC} = 10^{-2}$, the Ti-bearing melt production is almost 1 Gyr longer than in cases with stiffer HTC. The cumulative amount of melt produced is slightly higher than for the higher reference viscosity cases (a few km), and the duration of melt production is significantly shorter, reflecting the faster cooling of the mantle (see Figure S1 in Supporting Information S1). For both values of $\eta_{ref}$, all cases can sustain melt production of Ti-rich melts between 1 and 2.5 GPa until ~3 Ga, which is in good agreement with the age of high-Ti mare basalts (e.g., Hiesinger et al., 2003). Ti-bearing melt production for all cases is presented in Tables A1 and A2.

### 5.4. Influence of the HTC Overturn on the Lunar Core Dynamo

The early magnetization of the Moon is arguably one of its most intriguing features (Weiss & Tikoo, 2014). Among the hypotheses for its origin, a thermally driven core dynamo has been suggested. Such a mechanism requires a sufficiently large CMB heat flux, in excess of the heat flux conducted along the core adiabat ($f_{c,ad} = \alpha_c T_c k_c g(R_c)/c_{P,c}$, where $\alpha_c = 10^{-4}$ K$^{-1}$ is the thermal expansivity of the core (Stegman et al., 2003), $T_c$ its temperature, and other quantities are given in Table 2). Stegman et al. (2003) proposed that the late destabilization of overturned HTC can result in such a heat flux, in which case the destabilization time lag would account for the delayed magnetization of the lunar crust (approximately 1 Gyr after the Moon formation).

As pointed out earlier, our simulations do not support such a large-scale destabilization of the overturned HTC. Figure 12 shows time-series of the CMB heat flux from simulations that account for LMO crystallization. Interestingly, for cases that feature HTC overturn, for $\eta_{ref} = 10^{21}$ Pa s (respectively $10^{20}$ Pa s), we find a period extending over 2 Gyr (respectively 1 Gyr) during which the CMB heat flux drops below 0, centered at 3 Ga (respectively 3.5 Ga). It corresponds to the period of quiescent convection already that we alluded to in the context of HTC mobilization (or lack thereof), illustrated in Figure S1 in Supporting Information S1. During this period, the core (initially colder than the lower mantle due to the accretion-like initial temperature profile (Figure 1c)) is heated by the overturned HTC.

Except for the cases with $\eta_{ref} = 10^{20}$ Pa s, $E^* = 335$ kJ/mol and $\Delta\eta_{HTC} = 1$ and $10^{-1}$ (which do not feature HTC overturn), and $E^* = 150$ kJ/mol and $\Delta\eta_{HTC} = 10^{-3}$ (which does), all the time during which the CMB heat flux is high enough to drive a core dynamo is spent during LMO crystallization, and this phase occurs too early to match the magnetization record of the lunar crust. Furthermore, for the two exceptions mentioned, the case where HTC overturn occurs only produces a potential core dynamo very late (~0.5 Ga), also in disagreement with the lunar magnetic chronology. Thus, our results do not support a thermal core dynamo origin for the lunar crust magnetization.

Laneuville et al. (2013) also did not observe a sufficiently high CMB heat flow to drive a long-lived core dynamo in the case of pure thermal convection. They suggested instead that inner core growth can provide an additional buoyancy source to power a core dynamo driven by compositional convection, in line with the recent study of Briaud et al. (2023) that favors the existence of a solid inner lunar core.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Time-series of the core-mantle boundary (CMB) heat flux for cases with $n_{\text{ref}} = 10^{21}$ Pa s (a) and $n_{\text{ref}} = 10^{20}$ Pa s (b). The initial negative CMB heat flux is due to the accretion-like temperature profile that we prescribe in the unmolten lower mantle (see Figure 1c). The blue shaded area corresponds to the CMB heat flux exceeding the adiabatic heat flux out of the core, thus potentially driving a thermal dynamo. While the minimum value is different for each case, its variation between the cases is not visible with the adopted y-axis resolution, and we adopt the value of 6.9 mW/m$^2$, which is the average over all cases.</figcaption></figure>

Recently, Evans and Tikoo (2022) proposed an alternative scenario where the required CMB heat flux is provided by latent heat consumption associated with melting of the overturned HTC. This requires a protracted overturn encompassing ~1 Gyr. While the authors showed, based on scaling arguments, that this timescale is feasible under the assumption that the overturn occurs via sinking of small (~40 km in diameter) diapirs, our simulations show that HTC overturn is completed within at most a few hundreds of Myr from Moon formation, lending no support to this scenario. Nevertheless, the consumption of latent heat associated with re-melting of HTC is a feature that will need to be included in our model in order to address more thoroughly the issue of a lunar core dynamo driven by mantle flow.

## 5.5. Formation of a Partially Molten Layer in the Deep Lunar Mantle

The present-day existence of a partially molten layer atop the lunar CMB has been inferred from seismic, geodetic, and geophysical data (see Briaud et al. (2023) and Khan et al. (2007, 2014) and references therein). While alternative rheology models have been showed to also comply with the data without requiring such a layer (Walterová et al., 2023), its existence would place a strong constraint on the present-day thermal state of the Moon, with which thermal evolution models such as the one presented in this study should comply. The presence of an overturned layer of HTC helps satisfying this constraint as the melting temperature of HTC is significantly lower than that of the lunar mantle (Mallik et al., 2019; Wyatt, 1977). In order to test our model against this idea, we compare our laterally averaged temperature profiles with the HTC melting curves from Wyatt (1977). For the case $n_{\text{ref}} = 10^{21}$ Pa s, $E^* = 150$ kJ/mol and $\Delta n_{\text{HTC}} = 10^{-1}$, Figure 13 shows the evolution of the laterally averaged
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Time evolution of the laterally averaged profiles of (a) temperature and (b) high-Ti cumulates (HTC) fraction for the case $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 150$ kJ/mol and $\Delta\eta_{\text{HTC}} = 10^{-1}$. The contour lines in both panels corresponds to 5% and 30% melt-fraction of HTC (the actual melt fraction may be 0 if the local HTC content is 0 and the temperature is below the mafic mantle solidus). Notice that unlike in Figure 5, the HTC content colormap is not saturated.</figcaption></figure>

temperature (top panel) and of the laterally averaged HTC content (bottom panel), including a contour plot of the melt fraction of HTC at local *P*-*T* conditions. Because of the small pressure range of the lunar mantle, the contour plot follows almost exactly the isotherms of Figure 13a. In order to infer the actual melt fraction from the HTC melt fraction, one needs to multiply by the local HTC fraction (assuming that the *P*-*T* conditions are below the mantle solidus). However, the layer atop the CMB is nearly pure HTC, and the melt fraction at the CMB can be directly read from the contour plot. In the case presented in Figure 13, a ~30-km-thick partially molten layer (with melt fractions up to 12%) survives until present day above the CMB.

Among all our simulations, those which keep a warm CMB due to limited mantle cooling (i.e., stiffer cases that exhibit HTC overturn) achieve a higher CMB melt fraction (and a thicker partially molten layer). Mallik et al. (2019) recently obtained melting temperatures for HTC-analog material lower than Wyatt (1977), with respect to which they obtained a 100 K offset on account of the lower Mg# in their starting material, which can be relevant for HTC concentrations that we find. Using the melting curves for the HTC composition of Mallik et al. (2019) rather than those from Wyatt (1977) leads to significantly more cases having a partially molten layer
at the CMB. In all our cases accounting for LMO solidification, the largest CMB melt fraction (50%) and the thickest partially molten layer (187 km) are obtained for $\eta_{\text{ref}} = 10^{21}$ Pa s, $E^* = 335$ kJ/mol, and $\Delta\eta_{\text{HTC}} = 10^{-3}$. Results for all cases are presented in Tables A1 and A2.

## 6. Conclusions

By performing the first lunar overturn simulations accounting for a solidifying magma ocean, we have shown the following:

* the mechanical decoupling between the stiff plagioclase crust and the convecting mantle in the presence of a magma ocean significantly enhances the overturn of dense, late-stage HTC;

* the mantle and HTC rheology exert a strong influence on the overturn mechanism, whereby weak mantle and HTC promote a high time-dependence of the flow associated with small-scale dripping of HTC as soon as they crystallize, while stiffer mantle and HTC yield a steadier and larger-scale overturn after significant growth of the HTC layer;

* in cases where the overturn occurs while the LMO is still extant, it can significantly alter the LMO lifetime through an increase of the heat-piping flux associated with the return flow caused by the overturn. In most cases, the LMO lifetime is between 200 and 300 Myr, with the most extreme case peaking at 400 Myr.

Using our setup to investigate the influence of the overturn on the long-term evolution of the Moon, we have reached the following conclusions:

* despite internal heating, large-scale re-mobilization of the overturned HTC layer is unlikely to occur due to its high density. Instead, the overturned layer forms small piles that tend to be eroded with time, injecting high-Ti material in the convecting lunar mantle;

* Ti-bearing material entrained from the overturned HTC is prone to pressure-release melting, and under the assumption that the generated melt is positively buoyant within the 1–2.5 GPa pressure range, Ti-bearing melt can be produced until ~1 Ga, in agreement with estimates of the age of mare basalt production;

* the overturned HTC heat up the CMB, yielding a protracted time period (between 3.8 and 3 to 2.3 Ga) of negative CMB heat flow, during which a dynamo powered by thermal convection in the core is unlikely;

* the heat provided by the overturned HTC layer as well as its low melting temperature allow sustaining a partially molten layer above the CMB until present day, in agreement with geophysical constraints.

## Appendix A: Results Tables

Tables A1 and A2 present quantitative results over the complete parameter space investigated, for cases starting during LMO crystallization and after it, respectively. $\tau_{\text{LMO}}$ is the LMO lifetime, $\delta_{\text{W77}}$ and $\delta_{\text{M19}}$ indicate the thickness of the partially molten layer at the CMB (defined as having a melt fraction larger than 1%) calculated using the melting curves from Mallik et al. (2019) and Wyatt (1977), respectively, $\phi_{\text{W77}}$ and $\phi_{\text{M19}}$ indicate the melt fraction (in %) at the CMB calculated from Mallik et al. (2019) and Wyatt (1977), respectively, and $\delta_{\text{HTC}}$ is the thickness of the Ti-bearing mare basalt volcanism in GEL, all calculated at the present-day.
<table><caption>**Table A1**

Cases Accounting for Lunar Magma Ocean (LMO) Crystallization</caption><thead><tr><th>η<sub>ref</sub> (Pa s)</th><th>E* (kJ/mol)</th><th>Δη<sub>HTC</sub></th><th>τ<sub>conv</sub> (Myr)</th><th>τ<sub>LMO</sub> (Myr)</th><th>Ψ<sub>HTC,∞</sub> (%)</th><th>δ<sub>W77</sub> (km)</th><th>δ<sub>M19</sub> (km)</th><th>φ<sub>W77</sub> (%)</th><th>φ<sub>M19</sub> (%)</th><th>δ<sub>HTC</sub> (m)</th></tr></thead><tbody><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>0</sup></td><td>187</td><td>378</td><td>51</td><td>24</td><td>84</td><td>10</td><td>41</td><td>2,361</td></tr><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>-1</sup></td><td>138</td><td>400</td><td>75</td><td>52</td><td>133</td><td>18</td><td>43</td><td>596</td></tr><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>-2</sup></td><td>120</td><td>410</td><td>82</td><td>79</td><td>182</td><td>25</td><td>46</td><td>423</td></tr><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>-3</sup></td><td>107</td><td>402</td><td>84</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>0</sup></td><td>31</td><td>220</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>-1</sup></td><td>31</td><td>222</td><td>4</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>-2</sup></td><td>31</td><td>224</td><td>6</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>-3</sup></td><td>31</td><td>296</td><td>59</td><td>79</td><td>187</td><td>38</td><td>50</td><td>53,488</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>0</sup></td><td>15</td><td>312</td><td>86</td><td>0</td><td>133</td><td>0</td><td>31</td><td>0</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>-1</sup></td><td>15</td><td>322</td><td>95</td><td>0</td><td>106</td><td>0</td><td>26</td><td>0</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>-2</sup></td><td>15</td><td>330</td><td>97</td><td>0</td><td>84</td><td>0</td><td>23</td><td>0</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>-3</sup></td><td>15</td><td>330</td><td>98</td><td>0</td><td>57</td><td>0</td><td>13</td><td>558</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>0</sup></td><td>3</td><td>212</td><td>16</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2,020</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>-1</sup></td><td>3</td><td>252</td><td>56</td><td>0</td><td>46</td><td>0</td><td>19</td><td>3,276</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>-2</sup></td><td>3</td><td>282</td><td>64</td><td>57</td><td>176</td><td>23</td><td>45</td><td>1,894</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>-3</sup></td><td>3</td><td>308</td><td>79</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note. $τ_{LMO}$ is the lifetime of the LMO, $δ_{W77}$ and $δ_{M19}$ are the thickness of the partially molten layer computed with the HTC melting curves from Mallik et al. (2019) and Wyatt (1977), respectively, $φ_{W77}$ and $φ_{W77}$ the corresponding melt fractions at the CMB, and $δ_{HTC}$ the GEL thickness of HTC-bearing produced over the whole lunar evolution (computed with the melting curves from Wyatt (1977)). Notice that the cases with $η_{ref} = 10^{21}$ Pa s, $E* = 150$ kJ/mol and $Δη_{HTC} = 10^{-3}$, and $η_{ref} = 10^{20}$ Pa s, $E* = 335$ kJ/mol and $Δη_{HTC} = 10^{-3}$ could not run until 4,500 Myr due to computational difficulties, so we did not calculate their present-day state (the indicated value of $Ψ_{HTC,∞}$ for these cases correspond to the value at the end of LMO crystallization).

<table><caption>**Table A2**

Cases Not Accounting for Lunar Magma Ocean Crystallization</caption><thead><tr><th>η<sub>ref</sub> (Pa s)</th><th>E* (kJ/mol)</th><th>Δη<sub>HTC</sub></th><th>Ψ<sub>HTC,∞</sub> (%)</th><th>δ<sub>W77</sub> (km)</th><th>δ<sub>M19</sub> (km)</th><th>φ<sub>W77</sub> (%)</th><th>φ<sub>M19</sub> (%)</th><th>δ<sub>HTC</sub> (m)</th></tr></thead><tbody><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>0</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>-2</sup></td><td>10</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1,482</td></tr><tr><td>10<sup>21</sup></td><td>150</td><td>10<sup>-3</sup></td><td>72</td><td>89</td><td>198</td><td>55</td><td>57</td><td>18,548</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>0</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>-2</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>21</sup></td><td>335</td><td>10<sup>-3</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>0</sup></td><td>80</td><td>0</td><td>0</td><td>0</td><td>0</td><td>37,968</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>-1</sup></td><td>93</td><td>0</td><td>89</td><td>0</td><td>33</td><td>38,298</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>-2</sup></td><td>98</td><td>24</td><td>149</td><td>13</td><td>41</td><td>18,506</td></tr><tr><td>10<sup>20</sup></td><td>150</td><td>10<sup>-3</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>0</sup></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>-1</sup></td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>-2</sup></td><td>8</td><td>84</td><td>285</td><td>34</td><td>49</td><td>3,681</td></tr><tr><td>10<sup>20</sup></td><td>335</td><td>10<sup>-3</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note. The cases with $η_{ref} = 10^{20}$ Pa s and $Δη_{HTC} = 10^{-3}$ (for both values of $E*$) could not run until 4,500 Myr due to computational difficulties, so we did not calculate their present-day state.
## Data Availability Statement

All data, along with a python notebook to reproduce Figures 1 and 3–12 are available at Maurice et al. (2023a). Animations of the evolution of the temperature and HTC content of all simulations can be retrieved at Maurice et al. (2023b). Plots of the evolution of the laterally averaged temperature and HTC can be retrieved at Maurice et al. (2023c). GAIA is a proprietary code of DLR. Readers interested in using it should contact Nicola Tosi (nicola.tosi@dlr.de) and Christian Hüttig (christian.huettig@dlr.de).

## Acknowledgments

We thank two anonymous reviewers for their constructive comments that helped improve a previous version of this work. We also thank Sabrina Schwinger for her help with the LMO crystallization model. Part of this work was supported by the Helmholtz Association (Grant VH-NG-1017). Maxime Maurice received support from NASA Grant 80NSSC18K0828 during the preparation of the manuscript. Open Access funding enabled and organized by Projekt DEAL.

## References

Agrusta, R., Morison, A., Labrosse, S., Deguen, R., Alboussière, T., Tackley, P. J., & Dubuffet, F. (2019). Mantle convection interacting with magma oceans. *Geophysical Journal International*, 220(3), 1878-1892. https://doi.org/10.1093/gji/ggz549

Ballmer, M. D., Lourenço, D. L., Hirose, K., Caracas, R., & Nomura, R. (2017). Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle. *Geochemistry, Geophysics, Geosystems*, 18(7), 2785-2806. https://doi.org/10.1002/2017GC006917

Boukaré, C.-E., Parmentier, E. M., & Parman, S. W. (2018). Timing of mantle overturn during magma ocean solidification. *Earth and Planetary Science Letters*, 491, 216-225. https://doi.org/10.1016/j.epsl.2018.03.037

Briaud, A., Ganino, C., Fienga, A., Mémin, A., & Rambaux, N. (2023). The lunar solid inner core and the mantle overturn. *Nature*, 617(7962), 743-746. https://doi.org/10.1038/s41586-023-05935-7

Charlier, B., Grove, T. L., Namur, O., & Holtz, F. (2018). Crystallization of the lunar magma ocean and the primordial mantle-crust differentiation of the Moon. *Geochimica et Cosmochimica Acta*, 234, 50-69. https://doi.org/10.1016/j.gca.2018.05.006

Delano, J. W. (1986). Pristine lunar glasses: Criteria, data, and implications. *Journal of Geophysical Research*, 91(B4), D201-D213. https://doi.org/10.1029/JB091Ib04p0D201

de Vries, J., van den Berg, A., & van Westrenen, W. (2009). Formation and evolution of a lunar core from ilmenite-rich magma ocean cumulates. *Earth and Planetary Science Letters*, 292(1-2), 139-147. https://doi.org/10.1016/j.epsl.2010.01.029

Dygert, N., Hirth, G., & Liang, Y. (2016). A flow law for ilmenite in dislocation creep: Implications for lunar cumulate mantle overturn. *Geophysical Research Letters*, 43(2), 532-540. https://doi.org/10.1002/2015gl0066546

Elkins-Tanton, L. T. (2012). Magma oceans in the inner solar system. *Annual Review of Earth and Planetary Sciences*, 40(1), 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 and Planetary Science Letters*, 304(3-4), 326-336. https://doi.org/10.1016/j.epsl.2011.02.004

Elkins-Tanton, L. T., van Orman, J. A., Hager, B. H., & Grove, T. L. (2002). Re-examination of the lunar magma ocean cumulate overturn hypothesis: Melting or mixing is required. *Earth and Planetary Science Letters*, 196(3-4), 239-249. https://doi.org/10.1016/s0012-821x(01)00613-6

Evans, A. J., Andrews-Hanna, J. C., Head, J. W., III., Soderblom, J. M., Solomon, S. C., & Zuber, M. T. (2018). Reexamination of early lunar chronology with grail data: Terranes, basins, and impact fluxes. *Journal of Geophysical Research: Planets*, 123(7), 1596-1617. https://doi.org/10.1029/2017je005421

Evans, A. J., & Tikoo, S. M. (2022). An episodic high-intensity lunar core dynamo. *Nature Astronomy*, 6(3), 325-330. https://doi.org/10.1038/s41550-021-01574-y

Hess, P. C., & Parmentier, E. M. (1995). A model for the thermal and chemical evolution of the Moon's interior: Implications for the onset of mare volcanism. *Earth and Planetary Science Letters*, 134(3-4), 501-514. https://doi.org/10.1016/0012-821x(95)00138-3

Hiesinger, H., Head, J., III, Wolf, U., Jaumann, R., & Neukum, G. (2003). Ages and stratigraphy of mare basalts in Oceanus Procellarum, Mare Nubium, Mare Cognitum, and Mare Insularum. *Journal of Geophysical Research*, 108(E7), 5065. https://doi.org/10.1029/2002je001985

Hüttig, C., Tosi, N., & Moore, W. B. (2013). An improved formulation of the incompressible Navier-Stokes equations with variable viscosity. *Physics of the Earth and Planetary Interiors*, 220, 11-18. https://doi.org/10.1016/j.pepi.2013.04.002

Johnson, T., Morrissey, L., Nemchin, A., Gardiner, N., & Snape, J. (2021). The phases of the Moon: Modelling crystallisation of the lunar magma ocean through equilibrium thermodynamics. *Earth and Planetary Science Letters*, 556, 116721. https://doi.org/10.1016/j.epsl.2020.116721

Jones, M. J., Evans, A. J., Johnson, B. C., Weller, M. B., Andrews-Hanna, J. C., Tikoo, S. M., & Keane, J. T. (2022). A south pole-Aitken impact origin of the lunar compositional asymmetry. *Science Advances*, 8(14), eabm8475. https://doi.org/10.1126/sciadv.abm8475

Kegerreis, J. A., Ruiz-Bonilla, S., Eke, V. R., Massey, R. J., Sandness, T. D., & Theodorο, L. F. A. (2022). Immediate origin of the Moon as a post-impact satellite. *The Astrophysical Journal Letters*, 937(2), L40. https://doi.org/10.3847/2041-8213/ac8d96

Khan, A., Connolly, J. A. D., Maclennan, J., & Mosegaard, K. (2007). Joint inversion of seismic and gravity data for lunar composition and thermal state. *Geophysical Journal International*, 168(1), 243-258. https://doi.org/10.1111/j.1365-246X.2006.03200.x

Khan, A., Connolly, J. A. D., Pommier, A., & Noir, J. (2014). Geophysical evidence for melt in the deep lunar interior and implications for lunar evolution. *Journal of Geophysical Research: Planets*, 119(10), 2197-2221. https://doi.org/10.1002/2014JE004661

Laneuville, M., Wieczorek, M. A., Breuer, D., & Tosi, N. (2013). Asymmetric thermal evolution of the Moon. *Journal of Geophysical Research: Planets*, 118(7), 1435-1452. https://doi.org/10.1002/jgre.20103

Li, H., Liang, N. Z. Y., Wu, B., Dygert, N. J., Huang, J., & Parmentier, E. M. (2019). Lunar cumulate mantle overturn: A model constrained by ilmenite rheology. *Journal of Geophysical Research: Planets*, 124(5), 1357-1378. https://doi.org/10.1029/2018je005905

Lin, Y., Tronche, E., Steenstra, E., & van Westrenen, W. (2017). Evidence for an early wet Moon from experimental crystallization of the lunar magma ocean. *Nature Geoscience*, 10(14), 14-18. https://doi.org/10.1038/ngeo2845

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. *Journal of Geophysical Research: Planets*, 123(4), 910-951. https://doi.org/10.1002/2017JE005333

Mallik, A., Ejaz, T., Shcheka, S., & Garapic, G. (2019). A petrologic study on the effect of mantle overturn: Implications for evolution of the lunar interior. *Geochimica et Cosmochimica Acta*, 250, 238-250. https://doi.org/10.1016/j.gca.2019.02.014

Maurice, M., Tosi, N., & Hüttig, C. (2023a). Figures scripts. https://doi.org/10.6084/m9.fishshare.24018063.v1

Maurice, M., Tosi, N., & Hüttig, C. (2023b). HTC overturn movies. https://doi.org/10.6084/m9.fishshare.23985891.v1

Maurice, M., Tosi, N., & Hüttig, C. (2023c). HTC overturn plots. https://doi.org/10.6084/m9.fishshare.23985882.v1

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. *Journal of Geophysical Research: Planets*, 122(3), 577-598. https://doi.org/10.1002/2016JE005250
Maurice, M., Tosi, N., Schwinger, S., Breuer, D., & Kleine, T. (2020). A long-lived magma ocean on a young Moon. *Science Advances*, **6**(28), eaba8949. https://doi.org/10.1126/sciadv.aba8949

Miljković, K., Wieczorek, M. A., Laneuville, M., Nemchin, A., Bland, P. A., & Zuber, M. T. (2021). Large impact cratering during lunar magma ocean solidification. *Nature Communications*, **12**(1), 5433. https://doi.org/10.1038/s41467-021-25818-7

Moriarty, D. P. I., Dygert, N., Valencia, S. N., Watkins, R. N., & Petro, N. E. (2021). The search for lunar mantle rocks exposed on the surface of the Moon. *Nature Communications*, **12**(1), 4659. https://doi.org/10.1038/s41467-021-24626-3

Morison, A., Labrosse, S., Deguen, R., & Alboussière, T. (2019). Timescale of overturn in a magma ocean cumulate. *Earth and Planetary Science Letters*, **516**, 25–36. https://doi.org/10.1016/j.epsl.2019.03.037

Parmentier, E. M., Zhong, S., & Zuber, M. T. (2002). Gravitational differentiation due to initial chemical stratification: Origin of lunar asymmetry by the creep of dense KREEP? *Earth and Planetary Science Letters*, **201**(3–4), 473–480. https://doi.org/10.1016/S0012-821X(02)00726-4

Plesa, A.-C., Tosi, N., & Hüttig, C. (2012). Thermochemical convection in planetary mantles: Advection methods and magma ocean overturn simulations. In C.-P. Rückemann (Ed.), *Integrated information and computing systems for natural, spatial, and social sciences* (pp. 302–323). IGI Global. https://doi.org/10.4018/978-1-4666-2190-9.ch015

Pritchard, M. E., & Stevenson, D. J. (2000). Thermal aspects of a lunar origin by a giant impact. In R. M. Canup & K. Righter (Eds.), *Origin of the Earth and Moon* (pp. 179–196). The University of Arizona Press.

Pusok, A. E., Kaus, B. J. P., & Popov, A. A. (2016). On the quality of velocity interpolation schemes for marker-in-cell method and staggered grids. *Pure and Applied Geophysics*, **174**(3), 1071–1089. https://doi.org/10.1007/s00024-016-1431-8

Rapp, J. F., & Draper, D. S. (2018). Fractional crystallization of the lunar magma ocean: Updating the dominant paradigm. *Meteoritics & Planetary Sciences*, **53**(7), 1423–1455. https://doi.org/10.1111/maps.13086

Sahijpal, S., & Goyal, V. (2018). Thermal evolution of the early Moon. *Meteoritics & Planetary Sciences*, **53**(10), 2193–2211. https://doi.org/10.1111/maps.13119

Schmidt, M. W., & Kraettli, G. (2022). Experimental crystallization of the lunar magma ocean, initial selenotherm and density stratification, and implications for crust formation, overturn and the bulk silicate Moon composition. *Journal of Geophysical Research: Planets*, **127**(5), e2022JE007187. https://doi.org/10.1029/2022JE007187

Schulz, F., Tosi, N., Tosi, A.-C., & Breuer, D. (2020). Stagnant-lid convection with diffusion and dislocation creep rheology: Influence of a non-evolving grain size. *Geophysical Journal International*, **220**(1), 18–36. https://doi.org/10.1093/gji/ggz417

Schwinger, S., & Breuer, D. (2022). Employing magma ocean crystallization models to constrain structure and composition of the lunar interior. *Physics of the Earth and Planetary Interiors*, **322**, 106831. https://doi.org/10.1016/j.pepi.2021.106831

Singletary, S., & Grove, T. L. (2008). Origin of lunar high-titanium ultramafic glasses: A hybridized source? *Earth and Planetary Science Letters*, **268**(1–2), 182–189. https://doi.org/10.1016/j.epsl.2008.01.019

Snyder, G. A., Taylor, L. A., & Neal, C. R. (1992). A chemical model for generating the sources of mare basalts: Combined equilibrium and fractional crystallization of the lunar magmasphere. *Geochimica et Cosmochimica Acta*, **56**(10), 3809–3823. https://doi.org/10.1016/0016-7037(92)90172-F

Stegman, D. R., Jellinek, A. M., Zatman, S. A., Baumgardner, J. R., & Richards, M. A. (2003). An early lunar core dynamo driven by thermochemical mantle convection. *Nature*, **421**(6919), 143–146. https://doi.org/10.1038/nature01267

U, K., Kameyama, M., & Ogawa, M. (2023). The volcanic and radial expansion/contraction history of the Moon simulated by numerical models of magmatism in the convective mantle. *Journal of Geophysical Research: Planets*, **128**(9), e2023JE007845. https://doi.org/10.1029/2023JE007845

van Kan Parker, M., Sanloup, C., Sator, N., Guillot, B., Tronche, E. J., Perillat, J. P., et al. (2012). Neutral buoyancy of titanium-rich melts in the deep lunar interior. *Nature Geoscience*, **5**(3), 186–189. https://doi.org/10.1038/ngeo1402

van Orman, J. A., & Grove, T. L. (2000). Origin of lunar high-titanium ultramafic glasses: Constraints from phase relations and dissolution kinetics of clinopyroxene-ilmenite cumulates. *Meteoritics & Planetary Sciences*, **35**(4), 783–794. https://doi.org/10.1111/j.1945-5100.2000.tb01462.x

Wagner, T. P., & Grove, T. L. (1997). Experimental constraints on the origin of lunar high-Ti ultramafic glasses. *Geochimica et Cosmochimica Acta*, **61**(6), 1315–1327. https://doi.org/10.1016/s0016-7037(96)00387-0

Walterová, M., Běhounková, M., & Efroimsky, M. (2023). Is there a semi-molten layer at the base of the lunar mantle? *Journal of Geophysical Research: Planets*, **128**(7), e2022JE007652. https://doi.org/10.1029/2022JE007652

Wang, H., Agrusta, R., & van Hunen, J. (2015). Advantages of a conservative velocity interpolation (CVI) scheme for particle-in-cell methods with application in geodynamic modeling. *Geochemistry, Geophysics, Geosystems*, **16**(6), 2015–2023. https://doi.org/10.1002/2015GC005824

Weiss, B. P., & Tikoo, S. M. (2014). The lunar dynamo. *Science*, **346**(6214), 1246753. https://doi.org/10.1126/science.1246753

Wood, J. A., Dickey, J. S., Jr., Marvin, U. B., & Powell, B. N. (1970). Lunar anorthosites and a geophysical model of the Moon. In *Proceedings of the Apollo 11 lunar science conference* (Vol. 1).

Wyatt, B. A. (1977). The melting and crystallisation behaviour of a natural clinopyroxene-ilmenite intergrowth. *Contributions to Mineralogy and Petrology*, **61**, 1–9. https://doi.org/10.1007/BF00375941

Xu, M., Jing, Z., van Orman, J. A., Yu, T., & Wang, Y. (2022). Experimental evidence supporting an overturned iron-titanium-rich melt layer in the deep lunar interior. *Geophysical Research Letters*, **49**(13), e2022GL099066. https://doi.org/10.1029/2022GL099066

Yu, S., Tosi, N., Schwinger, S., Maurice, M., Breuer, D., & Xiao, L. (2019). Overturn of ilmenite-bearing cumulates in a rheologically weak lunar mantle. *Journal of Geophysical Research: Planets*, **124**(2), 418–436. https://doi.org/10.1029/2018JE005739

Zhang, N., Ding, M., Zhu, M.-H., Li, H., & Yue, Z. (2022). Lunar compositional asymmetry explained by mantle overturn following the South Pole–Aitken impact. *Nature Geoscience*, **15**(1), 37–41. https://doi.org/10.1038/s41561-021-00872-4

Zhang, N., Dygert, N., Liang, Y., & Parmentier, E. H. (2017). The effect of ilmenite viscosity on the dynamics and evolution of an overturned lunar cumulate mantle. *Geophysical Research Letters*, **44**(13), 6543–6552. https://doi.org/10.1002/2017GL073702

Zhang, N., Parmentier, E. M., & Liang, Y. (2013). A 3-D numerical study of the thermal evolution of the Moon after cumulate mantle overturn: The importance of rheology and core solidification. *Journal of Geophysical Research: Planets*, **118**(9), 1789–1804. https://doi.org/10.1029/jgre.20121

Zhao, Y., de Vries, J., van den Berg, A. P., Jacobs, M. H. G., & van Westrenen, W. (2019). The participation of ilmenite-bearing cumulates in lunar mantle overturn. *Earth and Planetary Science Letters*, **511**, 1–11. https://doi.org/10.1016/j.epsl.2019.01.022

Zhong, S., Parmentier, E. M., & Zuber, M. T. (2000). A dynamic origin for the global asymmetry of lunar mare basalts. *Earth and Planetary Science Letters*, **177**(3–4), 131–140. https://doi.org/10.1016/S0012-821X(00)00041-8

Ziethe, R., Seiferlin, K., & Hiesinger, H. (2009). Duration and extent of lunar volcanism: Comparison of 3D convection models to mare basalt ages. *Planetary and Space Science*, **57**(7), 784–796. https://doi.org/10.1016/j.pss.2009.02.002