Large impact cratering during lunar magma ocean solidification

The lunar cratering record is used to constrain the bombardment history of both the Earth and the Moon. However, it is suggested from different perspectives, including impact crater dating, asteroid dynamics, lunar samples, impact basin-forming simulations, and lunar evolution modelling, that the Moon could be missing evidence of its earliest cratering record. Here we report that impact basins formed during the lunar magma ocean solidification should have produced different crater morphologies in comparison to later epochs. A low viscosity layer, mimicking a melt layer, between the crust and mantle could cause the entire impact basin size range to be susceptible to immediate and extreme crustal relaxation forming almost unidentifiable topographic and crustal thickness signatures. Lunar basins formed while the lunar magma ocean was still solidifying may escape detection, which is agreeing with studies that suggest a higher impact flux than previously thought in the earliest epoch of Earth-Moon evolution.


Supplementary Note 1: Numerical impact modelling setup
The numerical impact simulations assumed a simplified three-layered target that was assumed flat for projectiles smaller than 90 km in diameter and curved for projectiles larger or equal to 90 km in diameter. The layers consisted of (a) a crust that was 10, 25 or 50 km thick, (b) an extremely low viscosity layer with a thickness of 10, 25 or 50 km (from here on referred to as the melt layer), and (c) a solid mantle. The justification of the crustal thickness range comes from the GRAILderived crustal thickness map of the Moon 2 (Supplementary Table 1) and estimates of the thickness of the flotation crust from LMO solidification calculations 3,4 . The presence of even the smallest melt fraction in rocks is known to significantly decrease their viscosity 5,6 . We do not address the chemistry of the melt, which determines its solidus temperature, because the melt composition is uncertain and changes as the magma ocean solidifies.
Two temperature profiles were employed for the majority of our simulations. For the first (Initial Moon A in Supplementary Fig. 1), the temperature gradient at the surface was set to 50 K/km, and the temperature was limited to maximum value of 1400 K at 25 km depth to ensure that the mantle was everywhere below the solidus temperature. Below 25 km depth, the mantle temperature profile followed an adiabat. Though this model has a surface temperature of about 80 K (which is somewhat lower than the true average value), however, the near-surface temperature rises to 280 K within 4 km depth due to the steep gradient applied. We note that iSALE does not provide much freedom for customizing temperature profiles when working with curved surfaces. This temperature profile was used for simulating the largest basins, for which size the curved surface plays a role in basin formation 7,8 . Regardless, tests using a Cartesian geometry showed that this lower surface temperature had no significant impact on the final basin morphology when compared with the other tested temperature profiles ( Supplementary Fig. 3a). For our second temperature profile, we used the profile H-0LB of Laneuville et al. 9 at 4.5 Ga which is the initial condition for their post-magma ocean thermal evolution simulations. This temperature profile (Initial Moon B in Supplementary Fig. 1) is cooler in the crust and hotter in the mantle than the first temperature profile. This thermal profile was used for small to mid-size basins, where the target was assumed flat.
Although the two initial temperature profiles described above were used for the majority of our simulations, to further test how the choice of the temperature profile affects the final basin morphology, we also tested a much hotter temperature profile (Initial Moon C in Supplementary  Fig. 1). This temperature profile is representative of the interior of the Procellarum KREEP Terrane at 4 Ga 9,10 . Though the temperature of the crust is between our two previous profiles, the uppermost mantle is considerably hotter and exceeds the mantle solidus. Tests show that even when the hottest thermal profile was used, it was the presence or absence of a low viscosity melt layer that played the most dominant role in the final crater morphology. Supplementary Fig. 3 shows that our final crustal thickness profiles are nearly the same for all three temperature profiles. We further investigated how much impact melt was generated using each of these temperature profiles, and found that there was little difference in the size and extent of the melt pool (see, Supplementary Fig. 4).
Supplementary Fig. 1: Temperature as a function of depth below the surface. Initial Moon A and B are two thermal profiles used in this work that are representative of the period immediately following magma ocean crystallization. Neither profile has temperatures above the material's solidus. Initial Moon C is the hottest tested temperature profile, which is representative of the interior of the nearside Procellarum KREEP Terrane at 4 Ga. Solidus and liquidus appropriate for dunite 9 are shown in black and gray, respectively.
Supplementary Fig. 2 shows our computed yield strength as a function of depth for the three temperature profiles A, B and C. The panel on the left shows the yield strength profile with depth when there is no melt layer, for the two thermal profiles used in this work (Initial Moon A and B).
Considering that these thermal profiles do not cross the set solidus, the yield strength is non-zero. This is changed when we introduce a melt layer, which causes the yield strength to drop to zero where the melt layer is located ( Supplementary Fig. 2  The melt layer was treated as a cohesionless fluid with constant low viscosity. Because of the melt being treated as a cohesionless fluid, its yield strength is zero in the computations. Solomatov (2007) 10 suggests that a typical viscosity of near-liquidus materials in a magma ocean is 0.01 Pa s with a factor of 10 uncertainty. Our tests showed that the chosen viscosity of the melt layer had little influence on the final basin morphology for all values less than about 10 10 Pa s, which is significantly lower than the viscosity of a solid rock. In this work, we used a constant viscosity of 100 Pa s, simply to avoid using (near) zero values in calculation.
Supplementary Fig. 2: Yield strength profiles as a function of depth below the surface. (left) Yield strength predicted using the two temperature profiles that do not include melt in the mantle (Initial Moon A and B). (right) Yield strength profiles after introducing a 25 km thick melt layer beneath the crust for the three different temperature profiles (A, B, and C). In this plot, the crustal thickness is assumed to be 25 km, with the exception of temperature profile B that also considers a 50 km crust.
We note that our simulations with a melt layer are not entirely self-consistent with the temperature profiles in Supplementary Fig. 1. Our approach was to investigate the consequence of adding a melt layer at the base of the crust while keeping all other variables constant. Given that the composition of the melt layer is uncertain, estimating its liquidus and solidus temperatures would also be uncertain. Furthermore, the composition of the melt layer changes as the magma ocean continues to crystallize. To investigate how a melt layer would affect the basin morphology, we thus simply changed the viscosity of the material to a low, non-zero, value (100 Pa s), that is appropriate for molten materials within magma oceans 11 while leaving the temperature of the melt unchanged. Table 2. These models are very similar to the models used extensively in previous lunar basin modelling 8,10,[12][13][14][15][16][17][18] . Our previous work 10 investigated the effects of acoustic fluidisation on large crater formation and found that its inclusion had little to no influence on the final basin morphology. In that study, the basin depth was largely unchanged when using acoustic fluidization with warm and hot temperature profiles in the target. However, it was shown that it was necessary to include acoustic fluidisation in order to match the observed depth of the basin floor when using colder temperature profiles. Acoustic fluidisation is included in this work because temperature profile B is comparable to some of the cold profiles used in Miljkovic et al. 10 . Acoustic fluidisation was applied in the same way on targets both with and without a melt layer. Furthermore, for the simulations that included a melt layer, test simulations were made with acoustic fluidisation turned on and off, which demonstrated that the final crater morphology was insignificantly different for the two cases. This further establishes that the melt layer is the dominant factor that controls the final crater morphology in our simulations ( Supplementary Fig. 3c).

The list of input parameters for our iSALE simulations is shown in Supplementary
In our simulations, we used both 17 km/s and 10 km/s for the bolide impact speed. The two speeds did not show significant differences in basin morphologies when the outcomes are analyzed in terms of the kinetic energies of the impactor. We chose to focus on the 17 km/s impact speed in order to limit the number of free parameters in our simulations. Simulations were run until the crater modification stage was completed, which was confirmed by the crater profile (both at the surface and the interface between the crust and the mantle) reaching a stable position and not moving more than a couple of cells over a significant time interval. This equilibrium was reached typically within ~3 h following impact depending on the basin size, which is similar to previous works 10, [13][14][15][16][17]25 .
In our work, we compare the relief of the surface and crust-mantle interface with the known surface topography of the Moon and GRAIL-derived crustal thickness models (Model 1 of ref. 2). An alternative approach could have been instead to compare directly the gravity field predicted by iSALE with GRAIL observations. However, to do so, it would have been necessary to specify the density and porosity of the crust, which are both uncertain. Any undamped numerical oscillations in the vertical direction at the end of the iSALE simulations would also have affected the gravity signal. We note that isostatic adjustment of the basin (which is largely vertical in nature) would modify the predicted gravity field but would not affect significantly the crustal thickness. For these reasons, we chose to the analyse the crustal thickness profiles instead of the predicted gravity, as they are less affected by the details of how the gravity was computed.
In this work, we focus on the final crustal thickness variations, and the formation of basin rings and their spacing. Though the final surface relief is also an important outcome of basin formation, this was not easy to interpret for many of our iSALE simulations. For example, for the same sized impact basin, the cell dimension was 2.5 km by 2.5 km when basin formation was modelled until completion, but we used 500 by 500 m cells when simulating the formation of faults (rings). Many of our highest resolution simulations that focused on faulting did not run to completion because of computational resource considerations. In simulations where the melt layer was present, even after 3 h following basin formation, there were vertical oscillations (equal to a couple of cells moving up-down) that affected equally the surface and crust-mantle interface, but not the crustal thickness that is the difference of the two. For this reason, we did not interpret the final surface relief predicted by our models, but concentrated instead almost exclusively on the resulting crustal thickness profile. We note that in this study, we made between 150 and 200 simulations all of which required day to weeks-long runtimes.
Supplementary  The impactor diameter was 60 km impactor, the impact speed was 17 km/s, and the crustal thickness was 25 km. Heavy black contours delimit the crust, initial melt layer and mantle boundaries, whereas the thin connected lines denote initially horizontal tracer particles that were separated by a 20-km depth interval. The impact melt region is approximated by temperatures greater than 1500 K. Tracers were removed from the simulations when they become too distorted, and the central region that lacks tracer lines closely corresponds to the region of impact melted materials. The impact melted region is roughly bowl-shaped, is no more than 200 km deep in the basin centre, and extends up to 300 km radially at the surface level. Most of the impact melted materials are derived from the mantle. Though the impact melt pool is slightly larger in the hot case, the total impact melt volume and geometry are similar for the two temperature profiles. The basin morphology is only weakly dependent on the thickness of the melt layer, when the melt layer is thicker than 10 km. For the case where there is no melt layer, a peak ring and main rim form at radial distances of ~220 and ~440 km, respectively. For the thinnest melt layer thickness of 10 km, it is possible that a single ring-like fault forms at the surface at a radial distance of ~300 km. However, any ring-like structure is ambiguous for thicker melt layers (panels C and D). This result is compatible with the observation that there is usually only one possible observable ring for the oldest pre-Nectarian basins (Supplementary Table 1).

Supplementary
Supplementary Fig. 7: Basin morphology dependence on melt layer thickness. From top to bottom, the melt layer is 0, 10, 25, and 50 km thick. The projectile diameter was 60 km and impacted into a 25 km thick crust at 17 km/s. White, light grey and dark grey shades denote the crust, melt layer, and mantle, respectively. Massless tracer particles (denoted by connected lines) were placed initially at every 5 km depth in the pre-impact target.

and melt layer thickness
Supplementary Fig. 9 shows crustal thickness signatures as a function of basin size and melt layer thickness. Each panel corresponds to a single projectile size (15,30, 60, 90, 120 and 160 km diameter), and the crustal profiles correspond to different melt layer thicknesses, ranging from 0 (purple), 10 km (blue), 25 km (grey), to 50 km (red). When no melt layer is present, all basins show a prominent thinning of the crust in the basin interior. When the melt layer is >10 km, all basins show a very similar profile, where the crustal thickness is larger in the basin interior and grades more smoothly to the ambient pre-impact value than for the case where no melt layer is present. The South Pole-Aitken basin is thought to be the largest and oldest impact basin on the Moon 7 . Our best-fit simulation for this basin was for a 160 km diameter projectile, similar to what was used in previous studies 7 . When a melt layer is present (of any thickness larger than 10 km), we find extensive crustal inflow into the basin centre giving rise to a crustal cap covering the entire impact melt sheet. When no melt layer is present, similar to previous investigations, basin formation ends with an exposed mantle-rich melt pool at the surface 7 . Although, the large size of the South Pole-Aitken basin in comparison to the lunar lithospheric thickness may explain the lack of prominent inner and outer ring structures, the final crustal signature when including a melt layer is the most similar to the GRAIL-derived crustal thickness profile ( Supplementary Fig. 10) 2 . Our simulations suggest that the crustal inflow is composed mostly of broken rafts flowing back into the basin centre that are composed primarily of overturned and jumbled crust originating from lower and mid-level crust levels. f