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

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.


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 coremantle boundary (CMB).Here we use the term "High-Ti cumulates" rather than the widely used "Ilmenitebearing 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., 2013Zhang et al., , 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., 2007Khan et al., , 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.

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 enrichement in the late crystallized (shallow) layers, where they enter clinopyroxene and Fe-Tioxides, 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).

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 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.

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):  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).
where ⃗v′ is the velocity vector, η′ 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 nondimensional radial profile of gravity (see below), Δρ′ 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: where α is the thermal expansivity, g 0 is the gravitational acceleration at the surface, κ is the thermal diffusivity, k is the thermal conductivity, D is the thickness of the lunar mantle, η ref is the reference viscosity, and ΔT is the reference temperature contrast.
The viscosity is calculated according to a pressure-and temperature-dependent Ahrrenius law, including a composition-dependent term: where η 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 ref and P ref are the reference temperature and pressure at which η ref is attained, Δη 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 ρ 1 : where r is the radius, G is the universal gravitation constant, and R c and ρ 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 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: 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.

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.

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 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 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.

Parameter Space
We vary the reference viscosity η 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 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 Δη HTC between 1 and 10 3 , 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  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.

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 (selfconsistently 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 (Ψ HTC ), which is defined as follows (Yu et al., 2019):  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.
where Ψ HTC is expressed in percent, V HTC is the total volume of HTC, and V′ 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 Ψ HTC .The final value, Ψ HTC,∞ (reached at the end of the timeevolution), is of particular interest to quantify the overturn efficiency.

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 η ref = 10 21 Pa s.When using this value (Figure 3a), the amount of overturned HTC is negligible (Ψ HTC,∞ < 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 Δη HTC also favor the overturn.For E* = 150 kJ/mol, Ψ HTC,∞ exceeds ∼50% for all values of Δη 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 Δη HTC = 10 3 could not be run past LMO crystallization, and we consider Ψ HTC at the end of LMO crystallization for Ψ HTC,∞ .This could underestimate Ψ HTC,∞ , making our observed trend conservative.
For the low value of the reference viscosity (η ref = 10 20 Pa s, Figure 3b), Ψ HTC,∞ 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 (Ψ HTC,∞ > 50%), while Ψ HTC,∞ ∼ 0% when LMO is absent, for all values of Δη HTC (except for Δη HTC = 10 3 , for which the case without LMO could not be run).Notice that the case E* = 335 kJ/mol and Δη HTC = 10 3 could not be run past the LMO crystallization.However its Ψ 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 Ψ HTC,∞ across the whole parameter space for cases that consider a crystallizing LMO (Figure 4a) and cases that do not (Figure 4b).All values of Ψ HTC,∞ are reported in Tables A1 and A2.

Morphology of the Overturn
When the crystallization of the LMO is considered, for η 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 η ref = 10 21 Pa s, E* = 150 kJ/mol, and Δη 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 Δη HTC (i.e., as HTC are made weaker).For these cases, the onset time of convection (which in the following we will indicate with τ conv ) decreases with Δη HTC (see Table A1).
Different from the case with E* = 150 kJ/mol, for E* = 335 kJ/mol (and η 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 (Δη 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 η ref = 10 21 Pa s, E* = 335 kJ/mol and Δη 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.
For the lower value of the reference viscosity (η ref = 10 20 Pa s), thermal convection starts early, prior to HTC crystallization, regardless of the value of the activation energy (and τ conv is accordingly independent of Δη HTC , see Table A1).The corresponding wavelength of the flow is similar to the cases with η ref = 10 21 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 shortwavelength 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.

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 (T ref = 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 T ref , 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 3 .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 21 Pa s, Ψ HTC,∞ > 50% when Δη HTC = 10 3 .By contrast, the overturn is nearly negligible for Δη HTC ≥ 10 2 .For the cases with η ref = 10 20 Pa s, Ψ HTC,∞ consistently increases with decreasing Δη HTC .

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(Maurice et al., , 2020)).Except for cases with η ref = 10 21 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 21 Pa s, Δη HTC = 10 3 , 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   HTC content and temperature colormaps are clipped as in Figure 5.
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  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.
flow and moderate (although not negligible) heat-piping, while the largerscale 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 η 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).

Influence of the HTC Overturn on the Long-Term Evolution of 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 deepseated, partially molten layer (Section 5.5).

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(Khan et al., , 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).

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., 2013Zhang et al., , 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 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 (Δρ = 574 kg/m 3 in our case compared to 300 kg/m 3 in Zhong et al. (2000), <165 kg/m 3 in Stegman et al. (2003), and 90 kg/m 3 in Zhang et al. (2013), where Δρ is the density difference between HTC and mantle).This difference (partly due to the asumption 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 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).
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 3 is "largely hypothetic."This corresponds to Δρ = 294 kg/m 3 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. (2013Zhang et al. ( , 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 η ref = 10 21 Pa s and ∼500 Myr for η ref = 10 20 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 η ref = 10 20 Pa s and Δη HTC = 10 1 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 presentday (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 produce.

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 postprocessing 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 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 Tibearing 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 η 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 (Δη 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 factor Δη 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.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 Δη 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 η 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 Δη 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 Δη HTC and the cumulative quantity of Ti-bearing melt produced or the duration of its production, although it is noteworthy that in the case Δη 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 η 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.

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 = α c T c k c g(R c )/c P,c , where α 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 η 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 η ref = 10 20 Pa s, E* = 335 kJ/mol and Δη HTC = 1 and 10 1 (which do not feature HTC overturn), and E* = 150 kJ/mol and Δη 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.
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.

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. (2007Khan et al. ( , 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 η ref = 10 21 Pa s, E* = 150 kJ/mol and Δη HTC = 10 1 , Figure 13 shows the evolution of the laterally averaged 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.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 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.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 η ref = 10 21 Pa s, E* = 335 kJ/mol, and Δη HTC = 10 3 .Results for all cases are presented in Tables A1 and A2.

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.

Figure 1 .
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 Table1.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).

Figure 2 .
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.

Figure 3 .
Figure 3. Time-series of Ψ HTC for η ref = 10 21 Pa s (a) and η ref = 10 20 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 Δη HTC from 1 (i.e., no viscosity contrast between high-Ti cumulates (HTC) and mantle) to 10 3 .Solid and dashed lines correspond to activation energies of 150 and 335 kJ/mol, respectively.The cases with LMO 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 , respectively, did not run past LMO crystallization, and we prolong their final value of Ψ HTC .The cases without LMO with η ref = 10 20 Pa s, Δη HTC = 10 3 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.

Figure 4 .
Figure 4. Final value of Ψ HTC (Ψ HTC,∞ ) for the whole parameter space investigated in this study.Panel (a) correspond 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 = 10 20 Pa s and Δη HTC = 10 3 could not be run.

Figure 5 .
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 η ref = 10 21 Pa s, E* = 150 kJ/mol, and Δη 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.

Figure 6 .
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 η ref = 10 21 Pa s, E* = 335 kJ/mol, and Δη 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.

Figure 7 .
Figure7.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 η ref = 10 21 Pa s and Δη HTC = 10 3 ; the case on the top row has E* = 335 kJ/mol (same case as in Figure6), while the case on the bottom row has E* = 150 kJ/mol (same case as in Figure5).Downward spikes in f HP in panel (d) are artifacts corresponding to simulation restarts that do not affect the evolution.

Figure 8 .
Figure 8. Lifetime of the lunar magma ocean (LMO) over the complete (η ref E* Δη HTC ) parameter space investigated in this study.

Figure 9 .
Figure 9. Present-day thermal state of the Moon for η 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 byKhan et al. (2007) (K07) andKhan et al. (2014), respectively (K14).

Figure 10 .
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 η ref = 10 21 Pa s, E* = 150 kJ/mol and Δη 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).

Figure 11 .
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), η ref = 10 21 Pa s, while in panels (b) and (d), η ref = 10 20 Pa s.Line colors and styles refer to different values of the HTC weakeningfactor Δη 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.

Figure 12 .
Figure 12.Time-series of the core-mantle boundary (CMB) heat flux for cases with η ref = 10 21 Pa s (a) and η 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 Figure1c).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.

Figure 13 .
Figure 13.Time evolution of the laterally averaged profiles of (a) temperature and (b) high-Ti cumulates (HTC) fraction for the case η ref = 10 21 Pa s, E* = 150 kJ/mol and Δη 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 Figure5, the HTC content colormap is not saturated.

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

Table 2
Model Variables and Parameters