A South Pole-Aitken impact origin of the lunar compositional asymmetry.

The formation of the largest and most ancient lunar impact basin, South Pole-Aitken (SPA), was a defining event in the Moon's evolution. Using numerical simulations, we show that widespread mantle heating from the SPA impact can catalyze the formation of the long-lived nearside-farside lunar asymmetry in incompatible elements and surface volcanic deposits, which has remained unexplained since its discovery in the Apollo era. The impact-induced heat drives hemisphere-scale mantle convection, which would sequester Th- and Ti-rich lunar magma ocean cumulates in the nearside hemisphere within a few hundred million years if they remain immediately beneath the lunar crust at the time of the SPA impact. A warm initial upper mantle facilitates generation of a pronounced compositional asymmetry consistent with the observed lunar asymmetry.

INTRODUCTION

One of the fundamental mysteries of lunar evolution has been the origin of the prominent and ancient asymmetry between the Moon’s farside and nearside hemispheres (Fig. 1) (). The giant South Pole–Aitken (SPA) lunar impact basin overprints the farside hemisphere’s generally high albedo—diagnostic of the anorthositic lunar primary crust—and low abundances of elements such as Th and Ti. In contrast, Th and Ti abundances are elevated across the nearside, and swaths of low albedo covering a similar area denote lunar maria—dense basaltic lavas that preferentially resurfaced the nearside beginning ~4.3 billion years ago (Gya) (–). The nearside’s most extensive maria and highest Th concentrations encompass the SPA basin antipode and extend toward the southwest, marking the Procellarum KREEP Terrane (PKT) (), a province that hosts a unique geochemical component with elevated levels of potassium (K), rare earth elements (REE), phosphorous (P), and other trace and heat-producing elements such as Th (, ).

Fig. 1.

Maps showing the lunar surface compositional asymmetry.

(A) Lunar Reconnaissance Orbiter Wide Angle Camera mosaic (). (B) Surface Th concentration () [used as a proxy for KREEP distribution ()]. (C) Surface Ti concentration (). Left: Lunar nearside hemisphere. Right: Farside hemisphere. South Pole–Aitken (SPA) basin ellipse (dashed white line) and basin center antipode (cyan circle) () are marked. Procellarum KREEP Terrane (PKT) outline (pink) is denoted by the 3.5–parts per million Th contour (). Maps are plotted in Lambert equal-area azimuthal projection.

KREEP originated from the final magma ocean–derived residuum that crystallized beneath the lunar crust, alongside other late-stage cumulates bearing Ti-oxide minerals like ilmenite (–). The lunar compositional asymmetry is commonly hypothesized (, –) to have originated through the redistribution of a global layer of subcrustal late-stage cumulates that occurred as early as ~4.4 Gya (the end of lunar magma ocean crystallization (–)). However, the catalyst for such a redistribution has remained enigmatic. Here, we show that the impact event that formed the SPA basin by 4.3 Gya () was appropriately timed and imparted sufficient heating to initiate lunar mantle dynamics that entrain subcrustal KREEP- and Ti-rich late-stage cumulates in lateral flow toward the nearside hemisphere. Under warm upper mantle conditions favored by SPA basin formation constraints (), a pronounced compositional asymmetry forms that is consistent with surface observations of the ancient lunar asymmetry.

We examine lunar evolution under two end-member scenarios for the state of the interior at the time of the SPA impact (Fig. 2, first row). All numerical convection simulations in this study (see Supplementary Text) begin with a global layer of ilmenite-bearing late-stage cumulates overlying the lunar mantle at the time of the SPA impact, in accordance with recent remote sensing analyses (). Reference scenario 1 (RS1) (Fig. 2A) considers an SPA impact into the stratified mantle predicted to exist immediately following solidification of the primordial lunar magma ocean (). On the basis of experimental constraints on the rheology of ilmenite-bearing cumulates (), the unmixed mantle of RS1 allows for late-stage cumulates that are weakened (i.e., lower viscosity) by a moderate ilmenite content (see Materials and Methods). Alternatively, for RS2 (Fig. 2, B and C), we consider an SPA impact into a thermally homogenized post-magma ocean cumulate mantle, in line with constraints for the lunar mantle thermal structure at the time of impact [e.g., ()]. To account for the possibility of some ilmenite depletion by mixing of the underlying mantle (), the subcrustal late-stage cumulate layer is less ilmenite-rich than in RS1 and has no compositional weakening effect. However, despite the difference in compositional effects, the temperature dependence of viscosity in our simulations means that the warm upper mantle of RS2, including the late-stage cumulate layer, flows more easily than the cooler upper mantle of RS1.

Fig. 2.

Cross sections of 3D lunar interior evolution simulations showing the effects of the SPA impact.

All cross sections are oriented with the farside SPA basin at the north pole. The lunar core (inner white circle) is shown to scale, while the crust (dark gray outer layer) is illustrated in postprocessing with exaggerated thickness. (A) Reference scenario 1 (RS1) (stratified post-magma ocean mantle) shown immediately before (top row) and after (second row) the SPA impact as well as 300 Ma (third row) and 600 Ma (bottom row) after impact. Left: Composition, distinguishing for simplicity between the material that is primarily late-stage cumulates (pink) and primarily other lunar mantle (gray). The dashed circle denotes ~500-km depth, the approximate maximum source depth of Ti-rich igneous rocks taken from the lunar nearside surface (). Right: Temperature. (B) Same as (A) for RS2 (mixed mantle after magma ocean solidification). (C) Same as (A) for RS2B (mixed mantle, more energetic impact). (D) Schematic cross-sectional illustrations of SPA-induced shock heating and simplified SPA-induced convection pattern. Shock heating at the time of the SPA impact (second row) is represented with an arbitrary isotherm (dashed orange line), radially symmetric from a point at depth below the impact.

RS1 and RS2 consider best-fit SPA impact parameters from vertical basin formation simulations (85-km radius projectile, 10-km/s impact) (, ). The more likely formation of the elliptical SPA basin by an oblique impact () allows for a more energetic collision [e.g., ()] with more extensive heating, because basin radius is a function of the vertical component of impact velocity (). To further highlight the range of possible outcomes, RS2B allows for an approximately order of magnitude more energetic SPA impact with the same initial mantle conditions as RS2 (Fig. 2).

RESULTS AND DISCUSSION

At the start of our simulations, we determine impact shock–induced heating using an analytic model that scales with the energy of the impact (, ). We find, for the range of plausible mantle parameters, that the SPA impact emplaces a mantle thermal anomaly (Fig. 2, second row) that exerts control over lunar interior evolution for hundreds of millions of years. The SPA-induced thermal anomaly drives a hemisphere-scale upwelling beneath the SPA basin on the lunar farside (Fig. 2, third and bottom rows), and the larger thermal anomaly of RS2B leads to a similarly sized downwelling beneath the impact antipode within the nearside PKT region (location marked in Fig. 1).

In every scenario examined, the SPA-induced hemisphere-scale convection pattern brings about lateral upper mantle flow toward the SPA antipode that entrains late-stage cumulates and generates a hemispheric compositional asymmetry, which is especially pronounced given warmer initial upper mantle conditions (Fig. 3). The development of the asymmetry in KREEP- and ilmenite-bearing late-stage cumulates as a result of the SPA impact agrees with the timing and distribution of multiple observed features of the lunar asymmetry.

Fig. 3.

Plots of the simulated and observed lunar compositional asymmetry.

(Left axis) Simulation results for surface-equivalent thickness of KREEP- and Ti-bearing lunar late-stage cumulates (pink) at 600 Ma after the SPA impact, azimuthally averaged with respect to distance from the SPA basin. (Right axis) Mean (black line) and 1σ variation (gray shaded region) of present-day lunar surface Th concentration () azimuthally averaged with respect to distance from the SPA basin. Simulated late-stage cumulate thicknesses are normalized to the initial layer thickness and are calculated from the total amount of late-stage cumulates above 500-km depth, the approximate maximum source depth of nearside Ti-rich igneous rocks (). The left and right vertical axes are independently scaled to show the full range of each dataset. Extent of the SPA basin () and approximate extent of the Th-rich PKT are denoted as blue shaded regions. The Imbrium (IMB) basin center position () is marked by a vertical dashed line at the shortest great circle distance from the center of SPA.

An early, global, subcrustal late-stage cumulate layer is supported by evidence for a global lunar magma ocean and the confirmed presence of Ti-rich basalts and KREEP on the lunar nearside (–) as well as recent remote sensing analyses, which indicate KREEPy late-stage cumulate lithologies that are associated with Th anomalies in the SPA basin (). The expected lateral and vertical distributions of emplaced SPA ejecta [e.g., ()] are proposed by Moriarty et al. () to agree with the distributions of observed lithologies and surface Th on the farside after processes including postimpact basin modification and inward slumping of Th-rich ejecta, thus supporting the presence of late-stage cumulates that were intact beneath the lunar crust and excavated during formation of the SPA basin.

In addition, crater retention ages suggest that the PKT is younger than the SPA basin (), and the most ancient nearside mare basalts erupted as early as 4.3 Gya (, ) [up to ~200 million years (Ma) after the proposed age of the SPA impact ()] followed by intense episodes of nearside mare volcanism beginning ~3.8 Gya () [~500 Ma to 700 Ma after the SPA impact ()]. The timing and location of the intense nearside mare volcanism are well explained by an early sequestration of heat-producing KREEP in the lunar nearside (). Our findings show that within 300 Ma (Fig. 2, third row) and until at least 600 Ma after the SPA impact (Figs. 2, bottom row, and 3), late-stage cumulates are concentrated in the nearside, which would lead to a KREEP and thermal asymmetry that could be linked with the formation of the PKT and nearside mare basalts.

Toward the end of our simulations, at ~3.9 Gya (400 Ma to 600 Ma after the SPA impact) (), the Imbrium impact would have excavated the locally concentrated KREEP material on the nearside (Fig. 3). Localization of KREEP-rich material to the lunar nearside followed by its excavation during the Imbrium impact has previously been suggested to have generated the Th signature across much of the lunar surface ().

The simulations shown in Fig. 2 represent outcomes along a continuum where the degree of nearside KREEP localization catalyzed by the SPA basin–forming impact is dependent on the thermal state and rheology of the lunar mantle as well as the energy of the SPA impact. For the relatively cooler and more rigid near-surface conditions of the stratified mantle in RS1, lateral transport of late-stage cumulates is moderate (farside thickness is reduced to about half and nearside thickness increases slightly), and the layer remains intact beneath the crust (Figs. 2 and 3). However, for warmer and less rigid upper mantle conditions predicted at the time of the SPA impact [e.g., ()], as in RS2 and RS2B, nearly all late-stage cumulates are transported away from the farside and accumulate concentrically around the SPA antipode. Constraints on the origin of Th-rich ejecta on the lunar surface suggest that basin-forming impacts other than Imbrium (i.e., those outside of the PKT) excavated little to no KREEP (), and the majority of extensive maria are limited to the area within and proximal to the PKT (Fig. 1); these observations point to the pronounced late-stage cumulate localization demonstrated by RS2 and RS2B.

With a warm upper mantle as well as a more energetic SPA impact than constrained by axisymmetric models (, ), late-stage cumulates may have been later driven deep into the lunar mantle. An SPA-catalyzed downwelling of late-stage cumulates as highlighted by RS2B (Fig. 2C) would be consistent with several lines of evidence that suggest that dense, subcrustal, KREEP- and ilmenite-bearing lithologies were transferred deeper into the lunar interior. Volcanic glasses found on the lunar nearside—with melt source depths constrained to ~500-km—assimilated high-Ti and KREEP components, necessitating the sinking of late-stage cumulates (). Reanalysis of lunar seismic data indicates a low-viscosity layer surrounding the lunar core (), which is proposed to be consistent with a stable layer of ilmenite-bearing lunar magma ocean cumulates that foundered to the core-mantle boundary (CMB) (, ). Furthermore, constraints from crater retention ages and remote sensing analyses suggest that the SPA impact occurred before the removal of early late-stage cumulates from beneath the crust (, ), allowing for an SPA-induced downwelling of late-stage cumulates.

A nearside downwelling of dense ilmenite-bearing cumulates, due to passive gravitational instability, has previously been put forth to explain the origin of the lunar asymmetry (), although the development of such a hemisphere-scale downwelling depends greatly on the rheology of late-stage cumulates (, ). A passive density-driven mechanism for the downwelling of late-stage cumulates has also been demonstrated as unlikely to entrain the KREEP-rich magma ocean residuum (). Our results show that hemisphere-scale convection in the aftermath of the SPA impact (a known basin-forming event) can initiate downwelling of late-stage cumulates even if the layer is not weaker or denser than the mantle. SPA-induced convection therefore offers a viable stand-alone catalyst that could be augmented by passive gravitational instability.

Our simulations of interior evolution for plausible early lunar mantle states demonstrate that heating from the SPA basin–forming impact initiates long-lived hemisphere-scale mantle convection, which, especially under warm mantle conditions favored by SPA basin formation constraints (), is capable of generating a compositional asymmetry consistent with noted observations. Although the major dynamical influence of interior heating by large impacts has already been examined extensively for Mars [e.g., (, )], near-surface processes discussed here may provide insight to the giant impact origin hypothesis of Mars’s north-south hemispheric dichotomy [e.g., (, )]. Furthermore, neither the reorganization of lunar interior dynamics nor the nearside concentration of late-stage cumulates caused by SPA impact heating have been widely considered in investigations of lunar evolution. Our results presented here indicate that SPA-induced convection is a fundamental consideration for lunar history, and the associated hemisphere-scale processes could have many implications beyond what is discussed here. For example, the nearside localization and possible downwelling of dense late-stage cumulates may influence the offset between the Moon’s center of figure and center of mass.

MATERIALS AND METHODS

Thermochemical evolution

Our lunar interior evolution simulations use a modified version of CitcomS, a three-dimensional, finite-element, thermochemical evolution code with Lagrangian tracer–based composition tracking (–). We consider a lunar mantle that can be described using the Boussinesq approximation and infinite Prandtl number. The governing equations of conservation of mass, momentum, energy, and composition take the following nondimensional forms (, )where u is velocity, P is dynamic pressure, η is viscosity, Ra is the Rayleigh number, T is temperature, δ is the Kronecker delta (where subscript r denotes the radial direction), t is time, Q is the internal heat production, and C is composition. Subscripts i and j are spatial indices, and a subscript comma indicates the partial derivative with respect to the following index or variable. Reference values used in nondimensionalization are listed in table S1.

In CitcomS, the basal heating Rayleigh number Ra is defined as ()where ΔT is the temperature drop from the CMB to the surface (in units of K), R0 is the planetary radius, and subscript 0 denotes the dimensional reference value of a quantity used in nondimensionalization. Relative to typical definitions of Ra, which normalize with the convective layer thickness, the CitcomS definition using the full radius R0 leads to a value of Ra that is approximately 20% larger for the Moon given our selected lunar parameters (table S1). As detailed in the Lunar interior models section, our simulations exclude the lunar crust, so ΔT and R0 are defined so that the outer surface of the model is the boundary between the crust and the underlying mantle (including late-stage cumulates). Hence, the model domain is the extent of the lunar mantle, from the CMB to the crust-mantle boundary. Boundary conditions at both surfaces are isothermal and free slip. The effects of rigid body rotation are removed.

Our simulations assume Newtonian rheology where temperature-, depth-, and composition-dependent viscosity is computed according to the dimensionless equationwhere A is a composition-dependent viscosity prefactor, is the dimensionless reference viscosity (with a magnitude of unity under the nondimensionalization scheme of CitcomS), Ea is the activation energy, Va is the activation volume, is the dimensionless temperature at the outer surface of the model, Rc is the core radius, and all other symbols are as defined above. Dimensionless activation energy Ea and volume Va are defined aswhere Ea and Va are the dimensional reference values for activation energy and activation volume, respectively, R is the universal gas constant, and all other symbols are as defined above.

The domain of our models consists of a spherical shell divided into 12 spherical caps, each with a grid resolution of 21 by 21 by 39 nodes in the directions of latitude, longitude, and radius, respectively. The global mesh is composed of (20 × 20 × 38) × 12 = 182,400 brick elements with radial length of ~36.1 km and lateral length ranging from ~18.3 km at the CMB to ~94.4 km at the outer surface of the model.

Composition is tracked using advective tracers, which allow for composition-dependent viscosity (). Our simulations use two tracer compositions—one for late-stage cumulates and one for the rest of the lunar mantle. Simulations are initialized with 50 tracers per element. At each time step of a simulation, the ratio of each tracer composition is used to compute the viscosity reduction coefficient A in each brick element of the mesh as a geometric mean between the late-stage cumulate compositional viscosity prefactor ALSC and unity.

Time step size is computed such that the size of each step is 70% of the dynamically computed maximum stable advection time step. For our model parameters, this yields time steps on the order of ~10−1 Ma to 100 Ma.

Impact-induced shock heating

Large impacts like that which formed the SPA basin impart substantial instantaneous heating to planetary interiors, and the resultant thermal anomaly can significantly influence subsequent mantle dynamics (, ). We model the emplacement of the SPA impact–induced thermal anomaly (impact parameters listed in table S1) by prescribing a temperature increase throughout the lunar mantle according to the Hugoniot release method of Watters et al. ()whereand where shock-induced temperature increase ΔTs is a function of shock-induced pressure increase Pδ, c is the specific heat capacity, ρ0 is the reference density, and S and C are the slope and intercept of the linear Hugoniot shock equation of state (EOS), respectively (C is also the acoustic velocity in the target material).

Instantaneous shock pressure Pδ (hence ΔTs) is approximately uniform within a spherical region referred to as the isobaric core, outside of which Pδ and ΔTs decay exponentially with radial distance (, ). The isobaric core is located at depth dc and has radius rc (). For the impact parameters used in this study (table S1)where values for the radius of the impactor Rimp = 85 km and velocity of the impact vimp = 10 km/s are from best-fit constraints of SPA basin formation simulations (, ). The pressure within the isobaric core, Pc, iswhere uc is the shock-induced particle velocity within the isobaric core under the assumption that the target and projectile materials are the same, e.g., dunite (), and other symbols are as defined above. Outside this region, the peak shock pressure, Ps, decays according to the inverse power lawwhere r is the radial distance from the origin of the isobaric core, and other symbols are as defined above.

Pδ at a given point is defined as the difference between the peak shock pressure Ps (or Pc when r ≤ rc) at that point and the preimpact pressure P0 at that point, i.e., Pδ = Ps − P0.

Our modified version of CitcomS calculates shock-induced temperature increase at each node of the finite element mesh at an impact time defined relative to the simulation start time (our simulations use timpact = 5 Ma after initialization). Once the start time (units of Ma) of a step exceeds the specified impact time and the energy equation has been solved for that step, shock heating is calculated and added to the temperature at each node of the mesh. Following the example of Roberts and Arkani-Hamed (), we assume that any melt produced in the impact returns to the mantle solidus within one time step because the time scale of crystallization of a body of impact-generated melt at the surface is much quicker than the time scale of solid-state convection (). That is, the maximum temperature Ts,max calculated during our shock heating routine at any node is limited by the mantle solidus.

From the equations and assumptions defined above, the magnitude and extent of the lunar mantle thermal anomaly generated by the SPA impact depend in part on the lunar mantle shock EOS (from S and C) and the energy of the SPA impact (from Rimp and vimp). Across a range of materials potentially relevant to the lunar mantle, differences in constrained values of S and C lead to little enough variation in the SPA-induced thermal anomaly (fig. S1) that we expect similar outcomes with any reasonable choice of shock EOS. To produce representative scenarios (Figs. 2 and 3), we therefore select values that lead to an intermediate thermal anomaly within the range of potentially relevant materials [S = 1.5 and C = 5.2 km/s ()].

Lunar interior models

To examine the first-order effects of the SPA-induced thermal anomaly on lunar interior dynamics, focusing especially on the evolution of an initially global, subcrustal layer (–, ) of KREEP- and Ti-bearing late-stage cumulates, we consider a two-layer lunar mantle (fig. S2). At the start of our models, a 60-km-thick layer of late-stage cumulates overlies a compositionally uniform silicate mantle that extends to the lunar CMB. Because of the temperature dependence of viscosity in our models and the expected low temperature of the lunar crust (, ), the crust is taken to be much more viscous than the mantle and is excluded from our models. This also reduces the overall viscosity contrast in the simulations and ensures that CitcomS is able to maintain numerical stability. Compositional density is uniform throughout the model. A dense late-stage cumulate layer could alter our simulations by promoting Rayleigh-Taylor instability of the layer [e.g., (, )]; however, we expect our conclusions to be largely unaltered as short-wavelength instabilities would likely be entrained in the stronger SPA-induced lateral flow field and long-wavelength instabilities may augment the compositional redistribution that we observe and facilitate a nearside downwelling of late-stage cumulates.

Simulations are initialized with one of two end-member thermal structures (fig. S2) that represent the range of plausible lunar mantle thermal states at the time of the SPA impact. Mantle model 1 allows for an SPA impact that occurs immediately after solidification of the lunar magma ocean and formation of the late-stage cumulate layer beneath the crust. The thermal structure of model 1 is adapted from model results for the stratified lunar mantle state at the end of magma ocean solidification (). Mantle model 2 considers a thermally homogenized cumulate portion of the post-magma ocean lunar mantle, in line with some considerations of the lunar thermal state at the time of the SPA impact [e.g., ()]. The late-stage cumulate layer is initially intact beneath the lunar crust in both models on the basis of recent remote sensing analyses that suggest SPA excavated a subcrustal KREEP-rich layer (). The deep, early lunar thermal structure is not well constrained; our models allow for a warm, basal primordial mantle layer that remains unmixed with the overlying magma ocean–derived cumulate mantle at the time of the SPA impact, based on models of magma ocean solidification (), but thermal constraints relevant for the upper mantle () simply assume an approximately adiabatic profile all the way to the CMB.

Since the temperature of the outer surface of each mantle model is different, ΔT also differs for each (table S2). The basal heating Rayleigh number Ra and nondimensional heating number Q thus differ between each, whereand where H is internal heating in terms of mass. For mantle model 1, Ra = 1.3 × 106 and Q = 13.72. For mantle model 2, Ra = 4.4 × 105 and Q = 41.16. The influence of internal heating and the convective driving force of the SPA-induced thermal anomaly mean that Ra as reported here underrepresents the vigor of convection.

We additionally consider compositional viscosity reductions on the late-stage cumulate layer. Recent rheological experiments indicate that ilmenite fractions of approximately 5 to 15%, as expected for lunar late-stage magma ocean cumulates, are capable of reducing the viscosity of silicate rock by up to ~2 to 3 orders of magnitude due to the low viscosity of ilmenite (). To broadly characterize the relationship between relative ilmenite content and the outcome of SPA-induced convection, we examine model scenarios with zero-, one-, and two-order magnitude viscosity reductions on the late-stage cumulate layer.