A long-lived magma ocean on a young Moon

Fitting the isotopic composition of lunar rocks to a new thermal and crystallization model shows the Moon formed 4.40 to 4.45 Ga.


Comparison with existing models
To benchmark our model, we reproduced the simulations performed in (2) and (3). One of the most important features of the crystallization model is the crystallization temperature, noted to distinguish it from the crystallization temperature from our model, and which reads in ref. (2): where is in K, is the radius in km and ( ) = ( = )/ , with , the initial volume of the magma ocean. Using a surface temperature of 25 K, as done in (2) and (3), we reach a final volume of the magma ocean of 1% (which is their criteria for the end of solidification) in 28 Myr. This is more than the ≈10 Myr found in (2) but very similar to the 26 Myr obtained for the same set-up in (3), where this discrepancy was already pointed out.

Convergence tests for the convection simulations
3-D convection simulations are computationally extremely costly, and can only be performed using relatively coarse grids. In order to choose a resolution that both allows for a reasonable computation time and yields a solution that faithfully captures the dynamics of the system We also evaluated the discrepancy between 2-D and 3-D simulations since some tests (such as the thermo-chemical ones presented in Supplementary text section 2) were run on 2-D grids only.
Discrepancies in the heat flow and temperature between 2-D cylindrical and 3-D spherical simulations are well documented in steady-state e.g. (41), but less constrained for transient regimes such as those considered here. A common method to mimic a flow in 3-D spherical shells using a 2-D cylindrical shell is to re-scale the core radius in order to recover in 2-D the 3-D ratio of coremantle boundary surface to planet surface (42). This results in a decrease of the radius of the core, which is already quite small for the Moon. Using our fiducial set-up, we compared three cases: a 3-D simulation using a resolution of 132 radial shells (i.e. the geometry used for the main results of the study), a 2-D simulation with "real" core radius and 249 radial shells, which closely matches the simulation computed on a 2-D grid with 132 shells ( Figure S3), and another 2-D simulation using 249 shells, but with a re-scaled core. Figure S4 presents the time-series corresponding to these three cases. The 3-D simulation (green lines) yields the same onset time of heat piping as the 2-D simulation with non-scaled core (blue lines), but a slightly shorter magma ocean duration due to a lower heat-piping flux. The 2-D simulation with a re-scaled core (red lines) exhibits a later onset of convection, but a similar magma ocean lifetime as the non-scaled 2-D simulation due to the trade-off between late onset of heat-piping below thick crust and early onset below a thin crust, i.e. the mechanism described in the main text. This suggests that there is no need to rescale the core size in our 2-D runs, but to keep in mind that they likely slightly overestimate the duration of the magma ocean solidification.

Comparison between different effects included in the thermal model
The right-hand-side of equation [20] contains various terms. Some of them are necessary but it decreases rapidly as the crust thickens and the magma ocean cools down. The onset of convection in the cumulates, that results in heat piping, is identified in Figure S5A by the initial peak in the heat piping curve (red line). In Figure S5B, it corresponds to the reduction in the cooling rate of the magma ocean and in the growth of the plagioclase crust. Heat piping then slowly decreases as the cumulates cool down and less melt is produced. The internal heating diminishes as the magma ocean shrinks, although preferential partitioning of the heat-producing elements in the liquid is taken into account. Radioactive decay of the heat sources plays a minor role in the internal heating, but is still accounted for. The conductive heat flux from the cumulates is minor for most of the crystallization except the end, when the top of the cumulates solidifies following a steep temperature gradient associated with the crystallization of KREEP. This steep gradient results in increased conduction from the cumulates, and the end of the magma ocean crystallization is controlled by the balance between conductive heat fluxes from the cumulates and through the crust.

Effect of the various heat fluxes and sources
In Figure S6, for a case with = 2 W/m/K with heat diffusion in the cumulates (i.e. black lines corresponding to the two crystallization models, without internal heating. While the two lines are very close to each other, the blue and solid black lines corresponding to the two crystallization models, but with internal heating, depart significantly from one another.

Effect of a dense and weak, late crystallizing layer on heat piping
In order to assess the effects of compositional buoyancy of a late crystallized, dense and weak layer as introduced by (15) Figure S9 shows a closeup of the whole-rock ε 176 Hf radial distribution between radii 1450 km and 1500 km where this value reaches its maximum.  Figure S9B and Table   S5.

Hf isotopic composition of the lunar zircon
As discussed in the main text, our estimate of the age of the Moon is at odds with a recent estimate based on lunar zircons, according to which the age of the Moon should be 4.51±0.01 Ga (23). One potential problem with the lunar zircon Hf isotopic data is that they must be corrected for the effects of secondary neutron capture during cosmic ray exposure (CRE). CRE-effects on Hf isotopes in lunar samples can be quite significant (43,44), and accordingly, (23)

5))
. We further select the cases that also match ε However, bulk rock analyses of sample 14321 show no significant CRE effects on Hf isotopic composition (44,47), suggesting that the anomalous 178 Hf/ 177 Hf for 14321 zircons reported in an earlier study (45) are not due to CRE and that, therefore, these zircon data should not be corrected for CRE effects.  Figure S2.      In the thermo-chemical case, the viscosity in the dense layer is lowered by one order of magnitude, as per equation [11]. Both cases are ran on 2-D grids.   (24)) for all cases matching KREEP isotopic ratios (red and pink lines). The initial ε 176 Hf value for Kal 009 (25) corresponds to the grey shaded area (within standard deviations). (B) ε 176 Hf evolution of the highly fractionated reservoir of all cases selected by the KREEP fitting process. Those cases that also cross the uncertainty ellipsoid of Kal 009 (black point and error bars) are plotted in pink, while those which do not are plotted in red (as in Figure 4). (C) distribution, as plotted in Figure 4 of the main text.