**4.1 Modeling Parameters**

Our model considers downwelling driven by the dense IBC layer within a fully solidified Moon. Overturn driven by dense silicates in a partially solidified Moon45 is not considered here because it is unclear if this process took place after the isotopic closure of ferroan anorthosites (FANs), which is required for geochronology and the plutonic nature of the Mg-Suite8.

Our model Moon includes five layers: core, lower mantle (Layer 0), upper mantle (Layer 1), ilmenite-bearing cumulate layer (IBC, Layer 2), and crust (Layer 3), from bottom to top. The core is only energetically coupled with the lower mantle32,46. The IBC has lower viscosity and higher density than the mantle and is overlain by a less dense crust. We set the core-mantle boundary, the lower-upper mantle boundary, the IBC bottom, and the IBC-crust boundary at the nominal radii of 340 km, 1040 km, 1660 km, and 1710 km, respectively. Following previous work32, we treat the initial thickness of the IBC layer as a free parameter by modeling cases of 30, 50, 100, and 150km to explore the dynamic return flow patterns and timing of upwelling mantle from Layer 0.

The viscosity contrast between the IBC layer and the underlying mantle plays a key role in determining the dynamics of CMO36,47-49. Model viscosity is both temperature and compositionally dependent. Here we explore a range of IBC viscosities constrained by experiments36. The predicted viscosity of pure ilmenite is ~4 orders of magnitude lower than that of dry harzburgite (approximated by dry olivine). Mantle Layers 0 and Layer 1 therefore assume the rheology of dry harzburgite, ranging from 5x1020 - 1021 Pa s, and we vary the viscosity contrast of the IBC layer to be between 10-1 - 10-4 × that of dry harzburgite32.

Instability is induced via random distribution of chemical tracers50, meaning we assign no initial perturbation to the IBC-mantle interface. The evolution of the four silicate layers is solved with conservation of mass, momentum, and energy. The formulation and parameterization of melt generation have been described in32. Our mantle solidus is defined by38. The mantle thermal Rayleigh number32 is given to 6x105. We apply parameterized equations of51 to solve for the local production of melt using the mantle solidus from38. We first assume that melting in excess of 2% is mobilized and is detectable in the crust by remote techniques11, and later increase this threshold as a free parameter in our spatial analysis (detailed further in our results). Our model uses a 12×64×48×48 mesh based on CitcomS52, which gives an azimuthal resolution of 14 km and radial resolution of 22 km. The ilmenite fraction of the IBC layer is ~ 6-12 wt% over the range of IBC thicknesses32.

Following previous work46,53,54, our model assumes 50% of all heat producing elements (U, Th, and K) are present in the IBC layer7, while the remaining 50% is evenly distributed throughout the lunar mantle. The heat generation rate of these heat producing elements is calculated based on the bulk U and Th abundances of the Moon. The bulk U and Th abundances of the present day are taken as 25.7 and 102.8 ppb (Th/U = 4), respectively55. The Moon is highly depleted of the volatile element K56-58. We use the K/Th ratio of 2,50059. Additional details of our model are found in32.

**4.2 Natural Constraints**

Constraint 1 (Mg-suite volume) is defined by the estimated amount of Mg-suite material within the lunar crust. Although loosely constrained, a modest fraction of the lunar crust is estimated to be comprised of Mg-suite material, constituting as much as ~6-30 vol.% of the lunar crust19,31. We therefore consider the total melt volume generated by each model case in relation to the total volume of the lunar crust. We note that the total melt volumes reported here represent a conservative estimate as some Mg-suite melts may have assimilated crust in producing more Mg-suite material11,60. The reference volume of the lunar crust is estimated by assuming a spherical shell and using an average crustal thickness of 40 km61.

Constraint 2 (magmatic duration) is defined by the estimated duration of Mg-suite magmatism. Recent geochronological results for the Mg-suite tightly constrains magmatic events to occur within ± 9 Myrs8. The variance associated with the Mg-suite age of crystallization could reflect a real distribution of crystallization ages among each sample measured, or it could be attributed to analytical uncertainty. If the latter, then the real duration of Mg-suite crystallization would be even shorter than defined above. We therefore treat the ~18 My variance defined by8 as our duration constraint. We approximate the duration of mantle melting in our dynamical models using the full width at half maximum (FWHM) of melt production rates (Supplementary Fig. 1). Constraint 3 (formation interval) is defined by the maximum interval of time between the ages of FANs secondary Mg-suite samples (including the respective variance associated with each age). The most reliable crystallization ages stated in our introduction (Mg-suite: 4340 ± 9 Myrs, anorthositic crust: 4359 ± 9 Myrs) constrains the formation interval to ≤ 37 Myrs8. We note that this formation interval includes the time associated with IBC-layer formation following isotopic closure of the crust, magmatic emplacement, and cooling to the closure temperatures associated with the isotopic systems used to date the Mg-suite samples. Our model estimates the time of magmatic emplacement, prior to cooling and isotopic closure. The IBC layer is also fully formed at time zero in our model. In this study, we consider the interval between time zero of the model and the time step most closely associated with 50% of the cumulative melt volume for each case. The time to 50% cumulative melt volume therefore provides a relatively conservative estimate for the timing of Mg-suite formation compared to the onset of melting in each model. Models reaching 50% cumulative melt volume after 37 My are deemed unsuccessful based on Constraint 3 (Supplementary Fig. 2).

Constraints 4 and 5 (exposure proportion and distance, respectively) account for the spatial distribution of Mg-Suite rocks. We take the spatial distribution of exposed Mg-Suite rocks from21. This study examines the mineralogy exposed in 164 craters across the global surface (33-199 km in diameter). Our selection of their dataset is based on their extensive search of craters and basins and the use of a common approach to mineral identification. Here, craters exhibiting spectral signatures of olivine, orthopyroxene, and spinel, with no evidence for clinopyroxene (i.e., a marker for mare basalts) are considered positive Mg-suite detections (Fig. 1). The Mg-suite exposure proportion based on the dataset of21 is 0.52 detections per crater examined.

Although the complete distribution of subsurface Mg-suite is unknown, the observed spatial distribution of Mg-suite detections can be quantified by measuring the current distance between each detection and its farthest neighbor (Constraint 5). If Mg-suite detections are confined to a small region of the Moon, the farthest neighbor distance for each detection will be relatively short compared to if Mg-suite is more globally distributed. We calculate the average farthest neighbor of observed Mg-suite rocks to be 5103 ± 243 km, which is nearly half the circumference of the Moon (~5460 km, or the maximum farthest neighbor distance achievable) indicating a broad distribution of Mg-suite exposures across the global surface (Fig. 1).

**References**

53. Boukaré, C.-E., Parmentier, E.M., Parman, S. Timing of mantle overturn during magma ocean solidification. Earth Planet. Sci. Lett., **491**, 216-225 (2018).

54. Zhang, N., Parmentier, E.M., Liang, Y. A 3D numerical study of the thermal evolution for the Moon after cumulate mantle overturn: The importance of rheology and core solidification. J. Geophys. Res. Planets, **118**, 1789-1804 (2013a).

55. Elkins-Tanton, L.T., Van Orman, J.A., Hager, B.H., and Grove, T.L. Reexamination of the lunar magma ocean cumulate overturn hypothesis: melting or mixing is required. Earth Planet. Sci. Lett., **196**, 249-259 (2002).

56. Parmentier, E.M., Zhong, S., Zuber, M.T. Gravitational differentiation due to initial chemical stratitfication: Origin of lunar symmetry by the creep of dense KREEP. Earth Planet. Sci. Lett., **201**, 473-480 (2002).

57. Zhao, Y., de Vries, J., van den Berg, V., Jacobs, M., van Westrenen, W. The participation of ilmenite-bearing cumulates in lunar mantle overturn. Earth Planet. Sci. Lett., **511**, 1-11 (2019).

58. Scheinberg, A., Elkins-Tanton, L., Zhong, S. Timescale and morphology of Martian mantle overturn. J. Geophys. Res. Planets. **119**, 454-467 (2014).

59. Katz, R. Spiegelman, M. Langmuir, C. A new parameterization of hydrous mantle melting. Geochem. Geophys. Geosys.,** 4**, 1073 (2003).

60. Zhong, S., et al. A benchmark study on mantle convection in a 3-D spherical shell using CitcomS. Geochem. Geophys. Geosys. **9**, Q10017 (2008).

61. Zhong, S. Parmentier, M., Zuber, M. A dynamic origin for the global asymmetry of lunar mare basalts. Earth Planet. Sci. Lett., **177**, 131-140 (2000).

62. Zhang, N., Parmentier, E.M., Liang, Y. Effects of lunar cumulate mantle overturn and megaregolith on the expansion and contraction history of the Moon. Geophys. Res. Lett., **40**, 5019-5023 (2013b).

63. Taylor, R.S. *Planetary Science: A lunar perspective*, (p. 481). Houston, Texas: Lunar Planet. Inst (1982).

64. Jones, J.H. Palme, H. Geochemical constraints on the origin of the Earth and Moon. In R. Canup, and K. Righter (Eds.), *Origin of the Earth and Moon* (pp. 197-216). Tuscon: Univ. of Ariz. Press (2000).

65. Albarede, F., Albalat, E. Lee, C.-T. An intrinsic volatility scale relevant to the Earth and Moon and the status of water in the Moon. Meteorit. Planet. Sci. **50**, 568-577 (2014).

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

67. Laneuville, M., et al. Asymmetric thermal evolution of the Moon. J. Geophys. Res. Planets. **118**, 1435-1452 (2013).

68. Warren, P.H. Anorthosite assimilation and the origin of the Mg/Fe-related bi-modality of pristine Moon rocks: support for the magmasphere hypothesis. J. Geophys. Res. **91**, D331-D343 (1986).

69. Wieczorek, M.A. et al. The crust of the Moon as seen by GRAIL. Science **339**, 671-675 (2013).