This new understanding of the basin has important implications for the upcoming robotic and human exploration of the lunar south pole, a key goal of which is to sample the ejecta of the SPA basin41. Proposed south polar Artemis landing sites are now seen to be situated within the Th-rich ejecta blanket at the downrange end of the basin. Thus, the rocks sampled by Artemis may constrain not only the age of the basin and history of lunar bombardment but also the composition of the late-stage magma ocean and timing of its solidification.
In this study, we identify the basin outline using the discrete transition in Bouguer gravity (the gravity anomaly corrected for the effects of topography) and topography at the outer edge of the central zone in which the crust has been thinned by the impact. The outlines of the comparison basins were identified using the radial distance at which the contrast from the basin centre to exterior (in Bouguer gravity or topography) reached 50% of the maximum value. This choice makes the radius relatively insensitive to data resolution and also less sensitive to the natural variability in the data than identifying points closer to the geological rim of the basin at a greater radial distance at which the quantity of interest is closer to the surrounding value. As a result, basin dimensions are smaller than previously published values in which the rim is typically defined using the high point in topography outside the steep transition leading into the basin centre. For multiring basins, our basin dimensions are much smaller than the outer ring diameters but are comparable with the diameters of the positive Bouguer anomalies. For basins that are largely isostatically compensated, topography, Bouguer gravity and crustal thickness all yield similar results. However, topography is more strongly affected by external processes such as impacts and erosion, and crustal thickness models typically require strong filtering to stabilize the solution, thereby decreasing their resolution. For these reasons, Bouguer gravity is taken as the preferred dataset. For the buried Utopia impact basin, the proximity of the basin to the dichotomy boundary complicates the selection of the outline and we manually select the Bouguer gravity value to use as the threshold at a value closer to the basin floor to mitigate this effect. For all basins, we mask out areas that are affected by other later processes (for example, impact craters superposed on the Crisium and Smythii rims, volcanic features and erosion on the rim of Hellas and the downrange tail region in the south and tectonic modification on the northern edge of Sputnik).
For the SPA basin, the above criterion yields a very irregular outline. The topography of the basin has been heavily modified by subsequent cratering, making the topographically defined basin shape unreliable. The Bouguer gravity and crustal thickness signature of the basin provides a better outline, but the transition between the basin floor and surrounding highland crust is slightly more gradational and irregular, particularly in the west, leading to large outward excursions in the rim location defined using a simple threshold. As discussed in the text, this may reflect the effect of the subsurface magma ocean on the collapse of the transient cavity in the western half of the basin, as predicted by models. The western outline in Bouguer gravity is similarly affected by the large linear gravity anomalies related to intrusive bodies. Moreover, the superposition of the basin on the global crustal and Bouguer gravity asymmetry has the potential to introduce errors in a simple threshold-based approach.
Instead, for tracing the outline of the SPA basin, we use gravity gradients, which we calculate as the maximum amplitude eigenvalue of the tensor of second horizontal derivatives of the Bouguer potential. The short-wavelength gravity gradients are strongly sensitive to long-wavelength variations in surface topography owing to the strong attenuation/amplification of the gravity anomalies with upward/downward continuation, resulting in much stronger values in areas of high topography when evaluated at a constant radius. This effect is removed by calculating the gravity and gravity gradients at the local surface of a smooth (spherical harmonic degree 30) topography model, so that the gravity gradients are everywhere evaluated at the local surface. In the gravity gradients, sets of discontinuous arcuate features circumferential to the basin define the basin outline. We identify the gravity gradient ring most closely matching the point at which the Bouguer anomaly reaches half of its maximum contrast in the eastern portions of the basin and trace this feature around the basin. In the northeastern quadrant of the basin, the gravity gradient ring most closely matching the Bouguer threshold approach takes a step inward towards the basin centre, leading to a more irregular outline, whereas in the northwest, the Bouguer threshold deviates out away from the basin rim. Our preferred outline tracing is based on a compromise between proximity of the gravity gradient feature to the 50% Bouguer threshold, continuity of the gravity gradient arcs and a smoothly varying basin rim.
As well as the gravity-based outline above and the topographic outline based on discontinuous massifs and scarps discussed in the main text, we also evaluate two different choices of basin outline for the SPA basin using gravity gradients and Bouguer gravity (Extended Data Fig. 1). The second outline is designed to more closely follow the 50% Bouguer gravity threshold, resulting in a more irregular basin rim with a distinct shift in radius in the eastern basin rim. The third outline selects the outermost prominent gravity gradient ring, although this seems to highlight an outer ring in the eastern basin compared with the western basin. We also test one extra topographic outline tracing for the SPA basin that stays further from the basin centre in the southern part of the basin. Ultimately, the precise choice of outline for the SPA basin is uncertain, but all choices taper southward, as shown in the next section.
With a discontinuous and irregular basin outline, it is not possible to simply measure the tapering of the basin outside some model of a fit to the shape of the outline. To fit the basin outlines, we first identified the best-fit centre by fitting the radial distance of the basin outline as a function of azimuthal angle as defined from the approximate basin centre with a Fourier transform and using the first Fourier mode to determine the offset between the true basin centre and the initial choice. The second Fourier mode reveals the long-axis orientation. We then use a Metropolis-Hastings Markov chain Monte Carlo algorithm to iteratively adjust the outline fit based on the misfit to the data.
To fit the data to a simple tapered ellipse, we use a shape defined by:
with the basin centred at the origin, in which |k| < 1/a. For the case of k = 0, this equation defines a simple ellipse with semimajor axis a and semiminor axis b. Positive and negative values of k define tapers in opposite directions and thus provide a means of evaluating the probability of a northward taper and northward projectile direction for the SPA basin. In the Markov chain Monte Carlo algorithm, a, b and k are varied, as well as the basin centre point location. This approach allows for a simple elliptical basin with no preference in tapering or impact direction. The basin tapering factor f is calculated by taking the difference between the basin width at the points midway between the basin centre and the uprange/downrange ends and dividing by the along-track distance between these points. Thus, this tapering represents the gradient in basin width of the middle 50th percentile of the basin. Fits were obtained for all basins discussed in this work (Extended Data Fig. 1). The mean and 1σ ranges in parameters were taken from the posterior distributions of the model (Extended Data Table 1) after 25,000 iterations with a burn-in period of 2,000 iterations. As well as the tapered ellipse, we also fit the basin outlines using a geometric model based on the physics of the formation of a tapered elongated impact basin from an oblique impact (see Supplementary Discussion).
We analyse the surface composition using elemental data derived from the low-altitude, high-resolution (0.5°) Lunar Prospector Gamma Ray Spectrometer (GRS) data. All data are reprojected in a basin-centred polar projection. Regions of interest are defined on the basin floor, in the western ejecta blanket, in the eastern ejecta blanket and in the highlands north of the basin using sectors in an axisymmetric reference frame (Extended Data Fig. 2). In each region of interest, we calculate the area-weighted means and standard deviations of Th, Ti and FeO. The standard deviations are converted to standard errors (the uncertainties on the means) by scaling by n, in which n is the effective number of independent observations in a given area assuming a dataset resolution of 5° of arc (we use a lower resolution than the actual grid resolution to be conservative, as the data are derived from a range of measurement heights). For Ti, the data are highly pixellated, so we smooth the data using a spherical harmonic filter with a low-pass cosine transition between degrees 50 and 70 (half-wavelength resolution of about 90 km) to highlight regional trends.
Further to the west (lower left in the figures), outside the probable extent of the SPA basin continuous ejecta blanket, the Th concentration remains high within the region surrounding Mare Australe. The highest Th concentrations are found within the mare units themselves and associated cryptomare deposits and thus are unrelated to the SPA basin ejecta deposit. These areas are excluded from the analysis. By contrast, increased Th concentrations within the broad SPA basin ejecta blanket are not associated with known mare or cryptomare deposits and remain at similar higher levels through the cratered highlands in that area. To the north-northwest of the SPA basin, two small (roughly 200 km) areas of strongly enriched Th concentrations are found, one of which is associated with a known cryptomare deposit. These localized features are not consistent with the expected broad distribution of a well-mixed ejecta blanket of relatively uniform composition and are more likely related to Th-rich features found globally, including the Compton-Belkovich anomaly, and to the near-side Th-rich red spots, some of which are associated with local magmatic and volcanic processes. Excluding these spots, Th concentrations north of the basin are indistinguishable from the far-side highlands, as discussed in the text. A similar Th-rich spot east of the basin is also excluded from the analysis.
Although there are many unknowns about early lunar evolution and the development of the present-day asymmetry, a simple model can account for the compositional asymmetry based on known or well-accepted aspects of lunar structure and evolution as described in the text. The Moon has an asymmetric crust that is much thicker on the far side than on the near side and that crust is largely isostatically supported apart from an important but small-amplitude fossil figure. The crust of the Moon crystallized from a magma ocean, which caused the progressive enrichment of incompatible elements including KREEP in that magma ocean. From these very simple aspects of lunar structure and evolution, the concentration of the final dregs of the magma ocean on the near side can be seen to be a natural consequence of the crustal asymmetry.
For a crust and mantle in vertical hydrostatic equilibrium, the remnant magma ocean trapped between the crust and cumulate mantle at any point will be thicker where the crust is thinner. This simple process can account for the near-side KREEP terrane without invoking any further geodynamic complexities, although it does not explain the formation of the crustal asymmetry to begin with nor the later concentration of Ti-rich cumulates on the near side, which required an extra process. The fundamental assumption made by this model is that the crustal asymmetry was largely established by the late stages of magma ocean crystallization, but it does not otherwise depend on the mechanism for generating the asymmetry.
We model the lateral distribution and thickness of the late-stage magma ocean assuming that the crust and mantle are in vertical hydrostatic equilibrium. We assume that the long-wavelength crustal thickness patterns outside the main impact basins were established during crustal formation. The modelled crustal thickness within the PKT is affected by the presence of the maria, so we begin with a model of the lunar crust before the maria, which has a thinner crust within the PKT than models that do not account for the presence of the maria. We interpolate the crustal thickness across the 38 largest impact basins, including the SPA basin itself, using an inverse-distance-squared algorithm. For all basins smaller than Serenitatis, there is commonly a pronounced thickening of the crust by ejecta outside the central positive Bouguer anomaly and we scale the positive Bouguer anomaly diameter by a factor of two before basin removal. For the SPA basin, we scale the diameter by 1.2 to avoid any rim and ejecta effects. After the interpolation, the base of the crust is filtered in the spherical harmonic domain with a low-pass cosine taper between degrees 45 and 60 to remove short-wavelength structures that are probably of later origin and to smooth out the interpolation boundaries. We assume that the base of the magma ocean is a spherical interface until it intersects the base of the crust, ignoring the small effects of rotational-tidal distortion that would affect the top of the mantle and base of the crust similarly.
We then map out the magma ocean thickness for a range of basal interfaces in 0.1-km increments beginning at the point at which the thickest crust first comes into contact with the underlying cumulates and compute the corresponding magma ocean volume and PCS (Extended Data Fig. 3). At these late stages of crystallization, the bulk cumulate mixture (including mafic and feldspathic minerals) is negatively buoyant. It is unknown whether these minerals segregate by density and the floatation crust continues to thicken or whether the bulk crystal mixture sinks. We consider both scenarios, assuming that plagioclase makes up 20% of the column of crystallizing minerals in each increment for the former case. For models in which the crust continues to thicken during these late stages, a model beginning with the present-day crustal thickness will end with a crust approximately 7 km too thick where the crust is thinnest. We iterate the model, adjusting the starting thickness to end up with the desired final pre-SPA basin crustal thickness model. The nominal model assumes a magma ocean depth of 1,000 km (bottom radius at 737 km) but we also consider magma oceans extending to the core-mantle boundary at 400 km (ref. ). The assumed initial thickness of the magma ocean only affects the calculation of the PCS at each step in the model. The nominal model assumes that the crustal thickness generated at each incremental step of crystallization is 20% of the crystallized magma ocean, but we also consider the case in which no further crust is generated because the mixed cumulates together are negatively buoyant.
We compare the predicted distribution of the magma ocean remnants to the Th anomalies in the PKT and the SPA basin, defined using the 3-ppm Th concentration contour. For the SPA basin, the distribution of Th on the basin floor as a function of azimuth is similar to that in the ejecta blanket, with the highest Th concentrations in the western half of the basin. For the SPA basin, the early stages of magma ocean crystallization will always have magma ocean beneath the Th anomalies on the basin floor, but the impact occurred at a transitional time in which the magma ocean was pinching out to zero thickness in the northeast of the basin. Many areas outside the basin probably did have remnant magma ocean at the time of impact and so cannot be used to evaluate the models. We compute the areal fraction of the high-Th zone within the basin underlain by magma ocean and identify the model with 50% coverage as the optimal model and the middle 50th percentile as the acceptable range (Extended Data Table 2). A qualitative comparison of the distribution of magma ocean relative to the distribution of Th-rich ejecta outside the basin yields similar results. For the PKT at the latest stage of magma ocean crystallization, the desired solution will have magma ocean below the PKT but not beneath the surroundings and thus models are penalized for both the lack of magma ocean within the PKT and the presence of magma ocean outside it. The optimal model has the best areal fit to the observed PKT and we identify the range of models that differ from the best fit by less than half of the total area of the PKT as the acceptable range. In the discussion in the text, we focus on the optimal values from the nominal model with a 1,000-km-deep magma ocean and continued crustal formation, and the range from all models, yielding a PCS for the SPA basin of 98.8% (98.4-99.3%) and for the PKT of 99.90% (99.82-99.98%).
The results of this study are not strongly sensitive to the details of the pre-SPA basin crustal thickness model. The topography and crustal thickness of the Moon outside the main impact basins is dominated by global spherical harmonic degrees 1 and 2 patterns. The location of the SPA basin southwest of the peak far-side crustal thickness and the orientation of its long axis towards the northwest dictates that the crust will be thicker to the northeast and thinner to the southwest, as in our reconstructed pre-SPA basin crust. The lack of any strong crustal thickness anomalies on the scale of the SPA basin or smaller, apart from identified impact basins, indicates that the presence of any strong departures from the global trends in the pre-SPA basin crust is unlikely. The preservation state of the SPA basin indicates that the crustal thickness distribution has not changed appreciably since that time, apart from later formed impact basins. If we instead take only the spherical harmonic degrees 1 and 2 components of the pre-SPA basin crustal model, the late-stage magma ocean liquids still concentrate first in the southwest half of the SPA basin and later in the PKT (Extended Data Fig. 4). Other studies have attempted to reproduce the pre-SPA basin crustal thickness distribution under different sets of assumptions but still produce the northeast-southwest gradient in crustal thickness across the basin that is responsible for the predicted excavation of late-stage magma ocean liquids in our work. Any model of the pre-SPA basin crust with a gradual decrease in crustal thickness away from the peak value in the far-side highlands will yield similar results.
The estimated concentration of Th in the reservoir depends on the assumptions made about the mixing of materials in the ejecta. As a lower-bound estimate on the Th concentration in the excavated reservoir at the SPA basin, we take the observed Th concentration in the west ejecta of 1.8 ppm assuming that the deepest materials are exposed directly at the surface in the overturned ejecta layer. More likely, the west ejecta blanket is a mixture of excavated crust and magma ocean liquids. Assuming that the late-stage magma ocean liquids had an FeO abundance of 25 wt% (ref. ), the observed FeO excess of the west ejecta relative to the far-side highlands of 2.2 wt% is compatible with a mixture of roughly 9% magma ocean liquids with crustal material, equivalent to excavating a 4-km layer beneath the 40-km-thick crust. For that mixing ratio, the Th excess of the west ejecta relative to the far-side highlands of 0.8 ppm would correspond to a Th concentration in the underlying reservoir of 9 ppm, which we consider to be a more realistic estimate. This value is comparable with the peak Th concentration of 5 ppm on the basin floor, although the basin floor is probably dominated by the differentiated melt sheet. Alternatively, if the ejecta includes material excavated from depths up to 80 km (ref. ), including a residual magma ocean thickness of 5 km only in the southwestern half of the basin, the Th excess relative to the east ejecta (which also excavated mantle material) of 0.6 ppm corresponds to a reservoir concentration of 9.6 ppm, whereas the Th excess relative to the far-side highlands would correspond to a reservoir concentration of 12.8 ppm. Together, these considerations support a range of 1.8-12.8 ppm, with a preferred value of 9 ppm.
Later, after the establishment of the near-side PKT, the concentration of Th in the Imbrium ejecta reveals the concentration in the final solidified KREEP reservoir. Assuming that the peak Th concentrations at the surface correspond to pure KREEP places a lower bound of 13 ppm, whereas accounting for the mixing of the excavated KREEP-free crust and KREEP-rich material provides a more realistic estimate of 36 ppm and further assuming that the ejecta also mixes with target rock by ballistic sedimentation raises this estimate to 75 ppm (ref. ). Similarly, scaling the Th excess of the peak concentrations in the Imbrium ejecta relative to the far-side highlands by the ratio of the FeO concentration in the late-stage magma ocean to the FeO excess in the ejecta yields a Th concentration of 28 ppm. These considerations support a range in Th concentrations of 13-75 ppm, with a preferred value of 28-36 ppm.
The evolving concentration of Th in the magma ocean is dependent on the partitioning of Th between the melt and solid phases on crystallization. However, estimates of the partition coefficient in the range 10 to 10 (refs. ) are effectively indistinguishable from perfect exclusion of Th from the solid phase. Instead, the partitioning of Th between the growing cumulate mantle and diminishing magma ocean is dominated by the presence of small quantities of melt trapped in the interstices of the cumulate mantle. Extraction of melt from the cumulate mantle would have been limited by percolation, which becomes inefficient at melt fractions less than about 2%. Previous studies have typically assumed a trapped melt fraction of approximately 5%, with values as high as 20% deemed unlikely. Because trapped melt leads to incorporation of incompatible elements in the cumulate mantle, these melt fractions correspond to an effective Th partition coefficient in the range 10 to 10, exceeding the partition coefficient for simple crystallization. We model the evolving Th concentration in the magma ocean as a function of PCS for melt fractions of 0, 5 and 10% and assuming an average Th partition coefficient of 5 × 10 (ref. ). The initial Th concentration is assumed to be the bulk silicate Earth value of 79.5 ppb (ref. ). The true bulk lunar Th concentration may differ from bulk silicate Earth depending on the relative contributions of the Earth's crust, Earth's mantle and projectile material to the debris disk that formed the Moon.
In discussing the timing of the SPA basin-forming impact, we use the constraint that the formation of the final KREEP reservoir in the PKT was largely complete by 4.34-4.37 Ga, as supported by several studies (see review in ref. ). However, geochronological estimates based on unshocked zircons indicate the survival of some residual magma ocean beneath the near side until about 3.90 Ga (refs. ). An extended duration for some fraction of the magma ocean on the near side is expected owing to the high concentrations of heat producing elements ultimately concentrated in the PKT, but this does not affect the interpretations or conclusions of this work.