Density structure of the Moon’s crust beneath the polar regions as revealed by gravity data.

The cumulative effects of impact cratering on the Moon have given rise to a surficial fine-grained impact gardened regolith that is 10s of meters thick, along with deeper layers of thick impact ejecta and in situ fractured crust. Collectively, the thick ejecta deposits and underlying fractured crust are referred to as the megaregolith, but its physical properties and thickness remain poorly understood. Impact generated porosity can affect the bulk density, seismic velocities, and thermal conductivity of the crust. Knowledge of how porosity varies with depth and with location is thus important for understanding the Moon’s impact history, geologic evolution and thermal evolution.

The gravity field of a planet is directly related to its internal density structure, and the high-resolution measurements made by the Gravity Recovery and Interior Laboratory (GRAIL) mission allow to investigate how density (and by inference porosity) varies both laterally and with depth in the crust. The work of Besserer et al. (2014) made use of an early GRAIL-derived gravity model to invert for the depth dependence of density in thecrust using assumed linear and exponential density profiles. Since then, substantially improved gravity models have been developed using the complete set of observations collected by the GRAIL mission. Recent models reach spherical harmonic degrees as high as 1200 (Goossens et al., 2020), 1500 (Konopliv et al., 2014), and even 1800 (Park et al., 2025). These new gravity models, combined with improved analysis techniques, allows us to revisit the subsurface density structure of the Moon.

In this study, we employed a localized spectral analysis method (Wieczorek & Simons, 2005; 2007) to compute the effective density spectrum which was then used to constrain the depth dependence of density in the lunar crust. We restricted our analyses to the polar regions where the quality of the gravity field was best. In contrast to previous studies that imposed a density profile of a specific form, we assumed that the density structure could be approximated by Nconstant density layers, and we investigated models with 2 to 5 independent layers. As a result of the large number of inversion parameters, we used a Bayesian inversion approach with a Markov Chain Monte Carlo sampling of the density profiles. For our localized analyses we used a spherical cap size of 15° with a spectral bandwidth of 43 and employed 14 windows with concentration factors greater than 0.99. When computing the misfit between the model and observations, we used spherical harmonic degrees between 250 and 780, in contrast to 250 to 550 as used in the earlier study by Besserer et al. (2014). For the model results presented here, we only use density profiles increasing with depth.

Our results reveal a distinct two-layer density structure near both poles (see Figure 1 for an example centered on the South Pole). The models are characterized by an upper megaregolith layer that is about 1-2 km thick with porosities of ~15–20% and a less porous underlying layer with porosities of 7-13%. These two layers are interpreted as being composed of allochthonous impact ejecta and autochthonous in situ fractured crust, respectively. The depth of the lower layer extends to at least 11 km, beyond which our inversion approach loses sensitivity to density. We show these porosities substantially reduce crustal bulk seismic velocities and thermal conductivity.

Figure 1. Posterior distributions of density as a function of depth for layered models centered over the South Pole, where density is constrained to increase with depth. Panels (a) to (d) show the results for 2-, 3-, 4-, and 5-layer models, respectively. The yellow line in each subplot indicates the median density as a function of depth; the solid and dashed red lines represent the maximum posterior density and its uncertainty in the high- and low-density groups, and the solid and dashed horizontal green lines represent our estimate of the transition depth and its uncertainty. Density and depth were sampled in the ranges of 1000–4000 kg/m³ and 0–50 km, respectively.

Reference:

Besserer, J., Nimmo, F., Wieczorek, M. A., Weber, R. C., Kiefer, W. S., McGovern, P. J., Andrews‐Hanna, J. C., Smith, D. E., & Zuber, M. T. (2014). GRAIL gravity constraints on the vertical and lateral density structure of the lunar crust. Geophysical Research Letters, 41(16), 5771–5777.

Goossens, S., Sabaka, T. J., Wieczorek, M. A., Neumann, G. A., Mazarico, E., Lemoine, F. G., Nicholas, J. B., Smith, D. E., & Zuber, M. T. (2020). High‐Resolution Gravity Field Models from GRAIL Data and Implications for Models of the Density Structure of the Moon’s Crust. Journal of Geophysical Research: Planets, 125(2), e2019JE006086.

Konopliv, A. S., Park, R. S., Yuan, D., Asmar, S. W., Watkins, M. M., Williams, J. G., Fahnestock, E., Kruizinga, G., Paik, M., Strekalov, D., Harvey, N., Smith, D. E., & Zuber, M. T. (2014). High‐resolution lunar gravity fields from the GRAIL Primary and Extended Missions. Geophysical Research Letters, 41(5), 1452–1458.

Park, R. S., Berne, A., Konopliv, A. S., Keane, J. T., Matsuyama, I., Nimmo, F., Rovira-Navarro, M., Panning, M. P., Simons, M., Stevenson, D. J., & Weber, R. C. (2025). Thermal asymmetry in the Moon’s mantle inferred from monthly tidal response. Nature, 641, 1188–1192.

Wieczorek, M. A., & Simons, F. J. (2005). Localized spectral analysis on the sphere. Geophysical Journal International, 162(3), 655–675.

Wieczorek, M. A., & Simons, F. J. (2007). Minimum-variance multitaper spectral estimation on the sphere. Journal of Fourier Analysis and Applications, 13(6), 665–692.

Wieczorek, M. A. (2024). Lunar shape models (LDEM_shape_pa) [Dataset]. SHTOOLS. https://shtools.github.io/SHTOOLS/python-datasets.html.