# High-Resolution Gravity Field Models from GRAIL Data and Implications for Models of the Density Structure of the Moon's Crust

S. Goossens, T. Sabaka, M. Wieczorek, G. Neumann, E. Mazarico, F.

Lemoine, J. Nicholas, D. Smith, M. Zuber

► To cite this version:

S. Goossens, T. Sabaka, M. Wieczorek, G. Neumann, E. Mazarico, et al. 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, 2020, 125 (2), 10.1029/2019JE006086. hal-02980457

HAL Id: hal-02980457

https://hal.science/hal-02980457

Submitted on 27 Jan 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

JGR Planets

RESEARCH ARTICLE

10.1029/2019JE006086

**Key Points:**

* We apply a special constraint based on topography to determine a high-resolution lunar gravity field model

* Our new model has high correlations with topography over its entire degree range

* We derive maps for the vertical and lateral density structure of the Moon's crust at high resolution

**Correspondence to:**

S. Goossens,
sander.j.goossens@nasa.gov

**Citation:**

Goossens, S., Sabaka, T. J., Wieczorek, M. A., Neumann, G. A., Mazarico, E., Lemoine, F. G., et al. (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, e2019JE006086. https://doi.org/10.1029/2019JE006086

Received 12 JUN 2019
Accepted 15 OCT 2019
Accepted article online 6 NOV 2019

# High-Resolution Gravity Field Models from GRAIL Data and Implications for Models of the Density Structure of the Moon's Crust

S. Goossens¹,² ID, T. J. Sabaka² ID, M. A. Wieczorek³ ID, G. A. Neumann² ID, E. Mazarico² ID, F. G. Lemoine² ID, J. B. Nicholas²,⁴, D. E. Smith⁵, and M. T. Zuber⁵ ID

¹CRESST, University of Maryland, Baltimore County, Baltimore, MD, USA, ²NASA Goddard Space Flight Center, Greenbelt, MD, USA, ³Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France, ⁴Emergent Space Technologies, Greenbelt, MD, USA, ⁵Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA

**Abstract** We present our latest high-resolution lunar gravity field model of degree and order 1200 in spherical harmonics using Gravity Recovery and Interior Laboratory (GRAIL) data. In addition to a model with the standard spectral Kaula regularization constraint, we determine models by applying a constraint based on topography called rank-minus-one (RM1). The new models using this RM1 constraint have high correlations with topography over the entire degree range by design. The RM1 models allow the determination of apparent crustal densities at all spatial scales (called effective density) covered by the model, whereas the Kaula-constrained model can only be used globally up to spherical harmonic degree 700. We find that the effective density spectrum has a smaller slope for the high degrees when compared to the medium degrees. We interpret this as indicative of a global average surface density, as opposed to an ever-decreasing effective density as one approaches the surface. We use the RM1 models to derive maps of lateral and vertical density variations in the lunar crust. These models allow us to increase the resolution of this analysis compared to previous studies, by increasing the degree range over which to fit theoretical models of vertical density variations, and by decreasing the size of the spherical caps used in a localized analysis. Several regions on the Moon, such as South Pole-Aitken and Mare Orientale, are distinct from their surroundings in terms of surface densities. The RM1 models are especially valuable in (localized) spectral studies of the structure of the lunar crust.

**Plain Language Summary** The Gravity Recovery and Interior Laboratory (GRAIL) mission was designed to investigate the Moon's interior structure and advance the knowledge of its history of heat flow by mapping the Moon's variations in its gravity field to high precision. We present new models of the Moon's gravity field (with a resolution of 4.5 km by 4.5 km at the Moon's equator), where we have used information of the Moon's topography. With such models, we can investigate the structure of the crust at all spatial scales, whereas standard models can only be used globally up to lower resolution (on average, 8 km by 8 km at the most). We find indications for a global average surface density, and we present maps of the lateral and vertical density structure of the Moon. Several regions, such as South Pole-Aitken and Mare Orientale, stand out as having different surface densities from their surroundings. Owing to their high resolution and improved correlations with topography, our models will be especially useful in local studies of the structure of the Moon's crust.

## 1. Introduction

The Gravity Recovery and Interior Laboratory (GRAIL) mission was designed to map the structure of the lunar interior from crust to core and to advance the understanding of the Moon's thermal evolution by mapping the gravity field of the Moon globally at high resolution (Zuber et al., 2013). The mission consisted of two spacecraft with inter-satellite Ka-band range-rate (KBRR) tracking as the single science instrument (Asmar et al., 2013; Klipstein et al., 2013), using a modified version of the ranging system used on the Gravity Recovery and Climate Experiment (GRACE) mission (e.g., Tapley et al., 2004). The two spacecraft were also tracked from the Earth using the Deep Space Network (DSN). GRAIL had two separate mission phases. GRAIL's Primary Mission (PM) lasted from 1 March 2012 until 29 May 2012, where the spacecraft flew at
a mean altitude of 55 km above the lunar surface, with a lower average altitude of around 20 km in early March and late May. It was followed by an Extended Mission (XM), which lasted from 30 August until 14 December, where the spacecraft had a mean altitude of 23 km before 18 November, and between 11–20 km through 14 December.

Global gravity field models are conveniently expressed in spherical harmonics, which are global basis func-
tions. Initial analysis of the GRAIL PM data resulted in a global model of the lunar gravity field of spherical
harmonic degree and order 420 (spatial resolution of 13 by 13 km at the lunar equator; Zuber et al., 2013),
surpassing the original requirements of determining a degree and order 180 model (a resolution of 30 by
30 km). These initial mission results showed that gravity and topography at the Moon are highly correlated,
up to 98% over a wide range of wavelengths, reflecting the preservation of crater relief in a highly fractured
crust (Zuber et al., 2013). These results were also used to characterize the lunar crust, showing that its den-
sity is lower than what was generally assumed (Wieczorek et al., 2013). This low bulk density of the crust
allowed the determination of the thickness of the lunar crust in agreement with constraints from the Apollo
seismic data, finding an average crustal thickness between 34 and 43 km. In turn, the results for the thickness
of the lunar crust allowed a reassessment of the Moon's bulk chemical composition (Taylor & Wieczorek,
2014), showing that the Earth and the Moon are similar in their concentration of refractory elements, putting
strong constraints on the formation of the Moon. Gradients of the Bouguer gravity (free-air gravity corrected
for the attraction of the surface topography) showed the existence of linear anomalies that are interpreted
as giant magmatic dikes, indicating evidence of early expansion of the Moon (Andrews-Hanna et al., 2013).

The models of the gravity field of the Moon were GRAIL's level-2 science product, and in order to assure
high-quality products, two separate teams analyzed the GRAIL data with different software and analysis
strategies to determine gravity models, cross-validating their results for consistency; there was one group at
the Goddard Space Flight Center (GSFC) and one group at the Jet Propulsion Laboratory (JPL). Subsequent
analysis of the GRAIL data by these groups has resulted in gravity field models of increasing degree and order
(Lemoine et al., 2013, 2014; Konopliv et al., 2013, 2014). The most recent models have either a maximum
degree and order of 1200 (Goossens et al., 2016) or 1500 (Park et al., 2015), but the local resolution in cer-
tain areas can be much less than the maximum degree of the models. GRAIL data have also been analyzed
by different groups using different processing techniques, focusing on using the PM data and on validating
the GRAIL science team's results (Arnold et al., 2015; Klinger et al., 2014; Wirnsberger et al., 2019). In addi-
tion to these global models, local solutions have been derived as well (Goossens et al., 2014; Goossens et al.,
2018; Han, 2013). These high-resolution gravity models (globally and locally) have enabled many studies on
the properties of the lunar interior: revealing the structure of impact basins (Neumann et al., 2015), deter-
mining the structure of the Procellarum region on the nearside of the Moon (Andrews-Hanna et al., 2014),
determining density (Wieczorek et al., 2013; Besserer et al., 2014; Han et al., 2014; Jansen et al., 2017) and
porosity (Besserer et al., 2014; Ding et al., 2018; Soderblom et al., 2015; Wieczorek et al., 2013) variations in
the crust at various scales, determining the thickness of the mare basalts (Gong et al., 2016), and determin-
ing the structure of basins such as Mare Orientale (Andrews-Hanna et al., 2018; Jansen et al., 2019; Zuber
et al., 2016), among many other studies concerning the lunar crust and basins. Finally, GRAIL data have
also greatly increased knowledge of the deep interior of the Moon, constraining parameters pertaining to
the Moon's layered structure such as core and density and viscosity of the layers, using improved measure-
ments of the Moon's mass, moment of inertia, and tidal response (Matsuyama et al., 2016; Matsumoto et al.,
2015; Williams et al., 2014).

As the resolution of the GRAIL gravity models improved, high correlations between gravity and topography
were obtained at smaller and smaller scales. At high degrees (small spatial scales), gravity and topography
are expected to have correlations close to unity, as most surface topography is expected to be uncompensated
(e.g., Wieczorek, 2015). Laterally varying density of the uppermost crust together with loading effects results
in correlations close to but different from unity. These high correlations have been exploited to infer both
lateral and vertical variations in crustal density. Besserer et al. (2014) and Han et al. (2014) both investigated
models of density variations with depth, finding surface porosities of 20% or higher, and density variations
in the crust due to compaction. In addition, Besserer et al. (2014) found a decrease in density, instead of an
increase, beneath the mare areas. Han et al. (2014) limited their analysis of gravity and topography to cover
spherical harmonic degrees 110 to 360 due to their use of models based on PM data. Besserer et al. (2014)
used the degree range 250 to 550, using models that included XM data with a maximum degree and order of
900. They were able to globally map lateral and vertical density variations. Yet at smaller scales (for degrees
<table><caption>**Table 1**
Summary of Data Used and Data Weights Applied in Our Analysis of GRAIL DATA</caption><thead><tr><th rowspan="2">Data type</th><th colspan="2">Mission phase</th></tr><tr><th>Primary</th><th>Extended</th></tr></thead><tbody><tr><td>DSN</td><td>Number of data points</td><td></td></tr><tr><td>KBRR</td><td>437,194</td><td>530,932</td></tr><tr><td></td><td>1,499,380</td><td>4,388,039</td></tr><tr><td></td><td>Number of arcs</td><td></td></tr><tr><td></td><td>40</td><td>58</td></tr><tr><td></td><td>Data weights</td><td></td></tr><tr><td>DSN</td><td>0.12 mm/s</td><td>0.12 mm/s</td></tr><tr><td>KBRR</td><td>0.03 µm/s</td><td>0.05 µm/s</td></tr></tbody></table>

larger than the range 600–700), the correlations between gravity and topography for these models decrease rapidly, and this limits the use of these models for geophysical analysis to only those degrees where the correlations are still deemed high enough.

Gravity models are affected by factors such as the geographically varying sensitivity of the data, resulting in noise in the models, especially at the higher degrees. This introduces spurious signals in the model, and as a result, the equation system that governs the estimation of gravity field parameters from satellite tracking data is typically ill-conditioned. Additional information is needed to stabilize the system and obtain a solution, which can be realized with constraints. In gravity field determination, the Kaula rule (Kaula, 1963, 1966) is often used as a regularization constraint. It assumes that coefficients have a nominal value of 0, and prescribes a variance of $\hat{\beta}/n^2$ around that zero value, where n is the spherical harmonic degree and $\hat{\beta}$ is a constant, which for the Moon is often either $25 \cdot 10^{-5}$ or $36 \cdot 10^{-5}$ (e.g., Konopliv et al., 2013; Lemoine et al., 2013). Such a spectral constraint is isotropic and often applied uniformly, which can lead to an underestimation of peak amplitudes (e.g., Floberghagen, 2002; Konopliv et al., 1999).

In Goossens et al. (2017), we introduced a different kind of constraint. We aimed to construct a constraint that would leave part of the solution determined by the data alone, yet that would use information from topography to improve correlations, with the goal to improve the recovery of fine-scale features in the gravity model. We showed that from a vector of M gravity model parameters $\mathbf{x}_a$ with non-zero elements, we can construct an *a priori* (pseudo-)inverse-covariance matrix such that for large values of the weight factor $\lambda$ applied to the constraint (indicating how strongly, relatively to the data, the constraint is to be applied in a standard constrained least-squares solution), the solution $\mathbf{x}_{k+1}$ for the k-th iteration of a standard Gauss-Newton iteration (e.g., Seber & Wild, 1989) would approach $\alpha\mathbf{x}_a$, and the correlation between gravity and topography would approach unity. We note that for such large values of the weight factor, the solution with the Kaula constraint would be **0**. We will show results for various values of the weight factor $\lambda$ in section 3.2. Following our earlier results, we will explore values between $\lambda$=0.1 and 100, while always keeping the relative weight on the data matrix at 1 (noting that the data matrix has data weights built in, and these are discussed in section 2.1 and shown in Table 1). The constructed constraint assigns infinite variance in the direction of $\mathbf{x}_a$, whereas in all directions orthogonal to $\mathbf{x}_a$, the preferred state is **0** (the same as for the Kaula constraint). The factor $\alpha$ is completely determined by the data. This leaves one degree of freedom, namely $\alpha$, and so we named our constraint *rank minus one* (RM1). We also showed that if we chose $\mathbf{x}_a$ to be the coefficients of an expansion in spherical harmonics of topography-induced gravity, following the procedure outlined in Wieczorek and Phillips (1998), then $\alpha$ can be interpreted as a scale factor for the (constant) bulk density of the crust.

To obtain the value of the density of the crust, it is assumed that the estimated gravity $\mathbf{g}_{\text{obs}}$ is related to the gravity induced by topographic relief $\mathbf{g}_{\text{relief}}$ of a crust of constant density $\rho_{\text{estim}}$ are as follows:

$$ \mathbf{g}_{\text{obs,nm}} = \rho_{\text{estim}} \mathbf{g}_{\text{relief,nm}} + I_{nm}, \qquad (1) $$

where n and m are spherical harmonic degree and order respectively. The contribution I is that part of the estimated gravity not modeled by topography, and it is assumed to be a random variable that is uncorrelated with $\mathbf{g}_{\text{relief}}$ and which has a zero mean. Following Wieczorek et al. (2013), an unbiased estimate of the density of the crust per degree n, $\rho_{\text{eff}}(n)$, also called *effective density*, can then be obtained from:

$$ \rho_{\text{eff}}(n) = S_{\mathbf{g}_{\text{obs}} \mathbf{g}_{\text{relief}}}(n) / S_{\mathbf{g}_{\text{relief}} \mathbf{g}_{\text{relief}}}(n), \qquad (2) $$

where crosspower $S_{ab}(n)$ (or autopower, if both quantities a and b are the same) is defined as:

$$ S_{ab}(n) = \sum_{m=-n}^{m=n} a_{nm} b_{nm}. \qquad (3) $$

Spherical harmonic degree is related to spatial scales, and thus, the effective density represents crustal density at different spatial scales. In Wieczorek et al. (2013), effective density was used to determine the bulk and laterally varying density of the crust. This effective density was also employed in (Besserer et al., 2014) and
Han et al. (2014) to infer vertical variations in addition to the lateral variations. In Goossens et al. (2017), we showed that applying this constraint to a pre-GRAIL data system results in a bulk density estimate for the Moon close to that obtained from GRAIL data, 2550 kg m⁻³ (Wieczorek et al., 2013). The effective density spectrum for the RM1-constrained solution presented in that work stays stable for its entire degree range. Extensive details and examples, as well as the application of the constraint to Mars, can be found in Goossens et al. (2017) and accompanying supporting information, so we will not repeat this here.

In this work, we apply our RM1 constraint to our latest analysis of the entire GRAIL data set. We determine a new degree and order 1200 model, and we generate solutions using either the standard Kaula constraint or our RM1 constraint. We then compare correlation and effective density spectra and use the RM1 model to determine lateral and vertical density variations of the crust of the Moon, with an emphasis on extending the results from Besserer et al. (2014) to smaller spatial scales. This work is structured as follows: we briefly discuss the GRAIL data analysis and solution methods in section 2. We present the gravity model results in section 3 and the density results in section 4. Finally, we present a discussion of the results in section 5 and the conclusions in section 6.

## 2. Data Processing and Method of Solution

Our analysis of the GRAIL data is the same as that described in publications presenting our primary and extended mission data models (Lemoine et al., 2013, 2014). We refer to these publications for all the details, and we give only a summary here, including any changes or updates we applied to our processing. First, we describe the processing of the GRAIL radio tracking data, and then we describe our solution strategy to estimate a gravity field of degree and order 1200 in spherical harmonics.

### 2.1. GRAIL Data Processing Methods

We use the NASA GSFC GEODYN II Orbit Determination and Geodetic Parameter Estimation software (Pavlis & Nicholas, 2017), which has been used extensively to analyze tracking and altimetric data from various missions orbiting different planets such as Mercury (e.g., Genova et al., 2018; Mazarico et al., 2014), Mars (e.g., Genova et al., 2016; Lemoine et al., 2001), and of course the Moon (e.g., Goossens et al., 2011; Lemoine et al., 1997; Mazarico et al., 2010) and Earth (e.g., Luthcke et al., 2006). The force and measurement models we use in our analysis are in general, the same as we used in our previous GRAIL work. For the model presented here, we reprocessed the GRAIL data using our latest publicly available degree and order 1200 model, GRGM1200A model (Goossens et al., 2016). This model already incorporated some of the changes listed below, when compared to the processing for the earlier GRGM900C model (Lemoine et al., 2014). For this work, we update our relativity modeling following Genova et al. (2018) for the GRAIL satellites, but we find no significant improvements in, for example, KBRR data fit. We use the DE430 ephemerides (Folkner et al., 2014) to obtain the positions of all planets and to define the lunar reference system in which we estimate the gravity field. In contrast to our processing that resulted in the GRGM900C model, our arc length for XM data is now also on average 2.5 days, the same as it was for PM data. We also lower the elevation cut-off angle for DSN data used in our analysis to 5°, with the elevation being measured as the angle between the local horizon at the DSN tracking site and the line-of-sight to the satellite orbiting the Moon. While low elevation data might be more affected by atmospheric signals, we do not find a degraded fit for the DSN data. We applied the same media corrections as in our earlier GRAIL analysis (Lemoine et al., 2013). Lowering the elevation cut-off increases the amount of DSN data in our model.

We list the number of arcs, number of data points, and data weights in Table 1. PM KBRR data have a sampling interval of 5 s, while XM KBRR data have a sampling interval of 2 s. The count interval of the Doppler data is 10 s in both phases. When compared to the number of data used for GRGM900C (see Lemoine et al., 2014's Table S1), we note that there are indeed more DSN data used, but also fewer KBRR data in this solution. This is because of a much more aggressive KBRR data editing that was undertaken for this analysis. This is the result of noticing occurrences of KBRR residuals that show a sudden, high-frequency behavior without having any apparent correlation with a lunar surface feature. We illustrate such an occurrence in Figure 1, where we show KBRR residuals (with respect to an earlier degree 900 model) for a span of 5 min in an XM arc. While the KBRR residuals typically oscillate, there are two instances in this example where their amplitude suddenly increases. After an in-depth analysis of the KBRR data, it was found that these occurrences very often coincide with certain reaction wheels on both GRAIL satellites reaching a rotation
<figure><img src="image_3.png" /><figcaption>**Figure 1.** KBRR residuals during an arc in the Extended Mission, showing an example of the sudden, high-frequency noise that was found to be related to reaction wheel speeds onboard the satellites.</figcaption></figure>

frequency of about 100 Hz (N. Harvey, JPL, priv. comm., 2014). The high frequency of these KBRR residual signals make it perhaps unlikely that they will be mapped into gravity features, yet on the other hand, as can be seen from Figure 1, their duration is easily several tens of seconds, during which the satellites have traveled a distance larger than the model resolution. Hence, we deleted KBRR data that show such oscillations.

The parameters we estimate can be conveniently divided into local and global parameters, where local parameters are those that influence only the data in an arc, and global parameters are those that influence all data. The local parameters we estimate are the position and velocity of each satellite at the start of the arc, one solar radiation pressure parameter per arc per satellite, and a large set of empirical accelerations. These consist of one constant acceleration and sine and cosine amplitudes of an acceleration with a period of one orbital revolution (which is about 110 min). We estimate these for both the along-track and cross-track directions (also called transverse and normal, respectively). We did not include the radial direction as we found that those accelerations adversely affect the recovery of low-degree terms and the potential degree 2 Love number. We estimate such sets of parameters (6 per satellite) every quarter orbit, and we apply time-correlation to ensure smoothness of the parameters (see Lemoine et al., 2013 for how these correlations are applied). The *a priori* sigma on these acceleration parameters is $10^{-9}$ m/s² for those in the PM, and $5 \times 10^{-9}$ m/s² for those in the XM. We chose larger uncertainties for the XM because of the lower altitude: when first processing the XM data, we found their inclusion improved correlations with topography (indicating that the larger uncertainty does not mean that the accelerations absorb gravity signal). We now also include the estimation of biases on the Doppler data during the XM, to account for remaining unmodeled signal in the Doppler data. The estimated values are in general $10^{-4}$ Hz or smaller, and their inclusion improves the Doppler data fit. We found slight improvements in KBRR fit as well when biases are included, but overall, the improvements are below the reported KBRR noise level. Whereas we previously estimated time biases on the KBRR data per arc, we are now using the final version of the KBRR data (release 4) for which this is no longer necessary. The global parameters we estimate are the spherical harmonic coefficients of the gravity field model up to degree and order 1200. These spherical harmonic coefficients are the standard geodesy $4\pi$ normalized coefficients (e.g., Kaula, 1966), and the reference radius we use for the gravity field model is 1738 km. We also included the estimation of the Moon's gravitational parameter *GM*, the product of the gravitational constant *G* and the Moon's mass *M*, and its Love number $k_2$, but the results do not change significantly from those for GRGM900C, so they are not included here.

## 2.2. Gravity Field Estimation Method

To generate the partial derivatives of the measurements with respect to the estimation parameters, we first iterate each arc by estimating the local parameters only. We use a batch least-squares approach (e.g., Montenbruck & Gill, 2000; Tapley et al., 2004) based on normal matrices for this, with constraints applied to those parameters that have listed an *a priori* sigma in the previous section. Each arc is deemed to be converged when the relative difference of two subsequent iterations' residuals root-mean-square (RMS) is less than or equal to 0.5 %. Then, we generate the partials with respect to all parameters. We note that a degree and order 1200 gravity field model has 1,442,397 coefficients, and for our processing, we have thus turned to the NASA Center for Climate Simulation (NCCS) supercomputer Discover, located at NASA GSFC.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** RMS power of various gravity field models. Thin lines denote the associated error curves for each model. The error curves have their calibration factors (see text) applied. The RMS power of the coefficient differences (see also Figure 3) between GRGM1200B and its RM1 variant with $\lambda = 1$ is also shown. Figure 2a shows the full degree range whereas Figure 2b shows a close-up of the higher degrees.</figcaption></figure>

We applied a method based on QR factorization (Golub & van Loan, 1989) to solve the resulting system. We make use of the Square Root Information Filter (SRIF; Bierman, 1977) in our implementation, and for details, we refer to Lemoine et al. (2013). We generate solutions with the standard Kaula constraints, where $\hat{\beta}$ is $36 \cdot 10^{-5}$, and separately, we generate solutions with our new RM1 constraint. For the latter, we choose $\mathbf{x}_a$ to be the expansion in spherical harmonics of the gravitational potential induced by surface topography, using a reference “unit” density of 1000 kg m$^{-3}$, following Wieczorek and Phillips (1998), where we used a series expansion of topography up to the 20th power (Wieczorek, 2015). We use the spherical harmonic expansion of topography generated from the Lunar Reconnaissance Orbiter's (LRO) Lunar Orbiter Laser Altimeter (LOLA) instrument (Smith et al., 2016), expressed in the Moon's principal axis system. For our RM1 constraint, we apply various weight factors $\lambda$, where higher values indicate a stronger constraint. Both the Kaula and RM1 constraints are applied for degrees larger than 600. We name the Kaula-constrained solution GRGM1200B. RM1 solutions will be denoted as GRGM1200B + RM1, together with the weight factor $\lambda$ that was used for a particular solution, since we investigate several cases: $\lambda = 0.1, 1, 10$, and 100. In Goossens et al. (2017), we expressed our RM1 constraint in terms of normal equation systems. We have now updated this to be expressed in the form of SRIFs.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Absolute values of the differences in coefficients between GRGM1200B and the RM1 model with $\lambda = 1$. Negative orders refer to the $\bar{S}_{nm}$ coefficients, positive ones to the $\bar{C}_{nm}$ coefficients.</figcaption></figure>

### 3. Gravity Model Results

We now discuss the results of the determination of a gravity field model of degree and order 1200 using either the Kaula constraint or our new RM1 constraint. We discuss the power spectrum of the model, the resulting (global) correlations with topography and effective density spectrum, Bouguer anomalies, data fit, and degree strength.

#### 3.1. Power Spectrum and Coefficient Differences

We compute the RMS power per degree spectrum $\sigma(n)$ of a model expressed in spherical harmonics as follows:

$$ \sigma(n) = \sqrt{\frac{\sum_{m=0}^{m=n} (\bar{C}_{nm}^2 + \bar{S}_{nm}^2)}{2n + 1}}, \quad (4) $$

where $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are the $4\pi$-normalized spherical harmonic coefficients. We show the power spectra of several models in Figure 2, including the associated error spectra, which are computed following equation (4) with the coefficient values replaced by the formal errors from the inversion for each coefficient. We also show the power of the coefficient differences between GRGM1200B and its RM1 version with $\lambda = 1$. The power spectra for GRGM1200A and GRGM1200B are nearly the same. Both models use the same data weights and have the same Kaula rule applied. The power for the RM1 models follows the power for topography-induced gravity closely (not shown in Figure 2), as expected. The spectra for GRGM1200B and those for the RM1 models start to diverge around degree $n=650$, and more sharply so around $n=750$.

When comparing the power spectrum of constrained models with that of the unconstrained model, differences are already apparent at degree $n=600$. For comparison, we also included the power spectrum of the unconstrained version of GRGM900C in Figure 2. By extending the resolution of the model from 900 to 1200, the unconstrained model's power increases sharply at a lower degree for the model with the larger maximum degree, which is likely due to the uneven spatial data coverage: as we will show in section 3.5, there are fewer areas with high-degree strength values.

In order to obtain a more realistic estimate of the error in our models, we calibrate the errors in our solutions by applying a calibration factor obtained from the final (constrained) SRIF (see Lemoine et al., 2014 for
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Correlations between the topography-induced gravity and various gravity solutions (a,c), and effective density spectra (b,d). The bottom row (c,d) shows close-ups between degrees 500 and 900 of the top row results to clearly show the various spectra around *n*=600, the degree from which both Kaula and the RM1 constraint have been applied.</figcaption></figure>

details; low degrees might have larger calibration factors, e.g., Williams et al., 2014), but we only apply one factor to all degrees. The calibration factors that we obtain are 1.635 for GRGM1200A, 1.555 for GRGM1200B, 1.529 for GRGM1200B+RM1 $\lambda$ = 0.1, 1.561 for GRGM1200B+RM1 $\lambda$ = 1, and 1.620 for GRGM1200B+RM1 $\lambda$ = 10. The calibration factor for the unconstrained system is 1.606. Calibration factors closer to 1 indicate that the *a priori* formal errors for that system are closer in agreement to those observed. The calibration factor for the RM1 models is close to that of GRGM1200B itself, and it increases as $\lambda$ increases, indicating that as the RM1 constraint is applied stronger, the formal errors will become increasingly optimistic. The error spectra for the RM1 models show a discontinuity at degree *n*=600, the degree from which the constraint is applied, and the discontinuity is stronger as $\lambda$ increases (and the formal errors thus become more and more optimistic). It is hardly present for $\lambda$ = 0.1. The power spectra themselves show no such discontinuity.

The power in the coefficient differences between GRGM1200B and the RM1 $\lambda$ = 1 model crosses the RM1 model's (calibrated) error curve around degree *n* = 660. This is also the degree around which the differences in power spectrum become apparent. For the lower degrees, the differences are well below the various error curves, and they always stay below GRGM1200B's (calibrated) error curve. Around degree *n* = 360, the difference curve shows a break, which it does again at *n* = 600 (beyond which the constraint is applied). This indicates that the RM1 constraint has an influence on coefficients well below *n* = 600. We show the coefficient differences between these two models in more detail in Figure 3, in a triangle plot that shows the absolute value of the difference per coefficient. The differences are largest for the high degree, high-order terms, as the RM1 constraint is much more able to influence individual coefficient values by driving them to those from topography-induced gravity (with a scale factor remaining). Streaks in both $\bar{C}_{nm}$ and $\bar{S}_{nm}$ coefficients are visible as well, starting around *n* = 360, corresponding to the break in the difference power spectrum. This is the degree around which resonances can occur (see Floberghagen, 2002, and discussion of our earlier models in Lemoine et al., 2013), and similar streaks are (weakly) visible around degrees 540 and 720.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Correlations between various gravity solutions and topography-induced gravity (a) and the corresponding effective density spectra (b). Here we focus on the results for the unconstrained GRGM1200B model, comparing them with the results for the Kaula-constrained GRGM1200B model and one model with the RM1 constraint.</figcaption></figure>

## 3.2. Correlations with Topography, Effective Density Spectrum, and Average Bulk Crustal Density

With the crosspower per degree $S_{ab}(n)$ given in equation (3), the correlation per degree $\gamma(n)$ between two sets of spherical harmonic coefficients can be expressed as

$$ \gamma(n) = \frac{S_{ab}(n)}{\sqrt{S_{aa}(n)S_{bb}(n)}}. \qquad (5) $$

The effective density spectrum for various models is computed using equation (2). We present the correlations with gravity-from-topography and the effective density spectra in Figure 4.

Correlations for GRGM1200B are only slightly better than those for GRGM1200A, and both show a sharp drop in correlations around the same degree, $n \approx 600$, which is also the degree from which the Kaula rule is applied. Correlations for an unconstrained 1200 model are shown in Figure 5, and they drop sharply at an even lower degree, $n \approx 500$, so the Kaula constraint helps to extend the correlations. Correlations for the RM1 models are much higher for the high degrees, which is as expected because the RM1 constraint was designed that way. For comparison, we also include the correlations for the degree and order 1500 model GL1500E (Park et al., 2015). This model was processed at JPL with different software (Konopliv et al., 2013, 2014), and a different constraint: for this model, a constant constraint of $8 \cdot 10^{-10}$ was applied starting at $n = 701$. The correlations for this model are better than those for GRGM1200A and B, and almost the same as those for the weakly constrained RM1 model with $\lambda = 0.1$.

With increasing $\lambda$, the correlations for the RM1 models approach unity rapidly. Figure 4 also shows a close-up of correlations between degrees 500 and 900, and it can be seen that at degree $n=600$, a jump in correlations occurs for values $\lambda \ge 1$. Such a jump does not occur for the model with $\lambda = 0.1$. For this model, the correlations are smooth, and they drop off at a higher degree than the Kaula-constrained model. The power spectra in Figure 2 do not show these jumps and neither do the effective density spectra in Figure 4.

As $\lambda$ becomes larger, the effective density spectrum will become flatter, showing fewer variations because the solution will approach $\alpha x_g$, which means only one effective density for the entire degree range, prescribed by $\alpha$ (Goossens et al., 2017). In Figure 4, this starts to be visible for the solution using $\lambda = 100$ but not for the solutions using smaller values for $\lambda$.

The effective density spectrum for the unconstrained model (Figure 5) shows large variations that are physically unrealistic, but we note that they do not drop immediately as they do for the Kaula constrained model, but rather stay somewhat stable until degree $n \approx 800$. For GL1500E (Figure 4), the effective density spectrum makes a sharp drop at $n = 701$, which is caused by the constant constraint applied. For degrees smaller than 700, the effective density spectrum of GL1500E follows that of the RM1 models (Figure 4d). Yet even though the RM1 $\lambda = 0.1$ correlations are similar to those from GL1500E around degree $n = 700$, the RM1 model's effective density spectrum is stable to higher degrees.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Average bulk crustal density as obtained by averaging the effective density spectrum versus the maximum degree used in the averaging. We show two results: one that starts the averaging at *n*=150 and one that starts at *n*=250.</figcaption></figure>

The application of a constraint may result in bias in the solution (e.g., Floberghagen, 2002; Xu, 1992). When introducing the RM1 constraint in Goossens et al. (2017), we used the constraint on a pre-GRAIL data set and compared it to the GRGM900C model, showing that we obtain similar results with minimal bias. While we cannot do a similar test on our degree and order 1200 models, we do not have direct indications that our solutions are heavily biased. The RM1 constraint can be seen in the framework of a transformation from spherical harmonics into a different basis, where one direction is spanned by gravity-from-topography (our choice for x_a). As stated, the RM1 constraint assigns infinite variance to the x_a direction, and for all other directions, the preferred state is zero. For GRAIL, this undamped mode is dominant, as seen by the high correlations between gravity and gravity-from-topography. While for large constraint factors λ it is inevitable that the solution projected in the x_a direction is affected, the effects of bias on it are mitigated: the combination of parameters that the RM1 constraint does not damp corresponds to a strong physical mode, and because the remaining signal is likely small, this protects the dominant mode from being too biased. For instance, in Figure 4d, the effective density spectra for solutions using different values of λ are clustered together and do not show different trends (see also section 5.1).

The average bulk crustal density can be estimated from the effective density spectrum by averaging this
spectrum over a range of degrees (Wieczorek et al., 2013). In order to account for uncertainties in the gravity
field, we generated 1,001 clone models, which are an ensemble of solutions of the same statistical family
(Sabaka et al., 2014). Following Goossens et al. (2017), we average the effective density spectra of all clones
of the RM1 model with $\lambda = 1$ to obtain the bulk crustal density and its uncertainty. Using the degree range
$n=150-310$, which was used in Wieczorek et al. (2013), we obtain a crustal density of $2515 \pm 25 \text{ kg m}^{-3}$.
If we increase the degree range over which we average, the bulk density will become smaller because the
effective density is smaller at higher degrees. We show this in Figure 6. By including higher degrees, more
weight is given to shallower depths, thus skewing the average to lower values. If we fit an analytical model
to the global effective density spectrum to obtain a global density-depth relationship (see section 4.1 for the
expressions), we can then use the results to obtain an average crustal density by integrating the density-depth
relationship. If we limit the integration to the upper 10 km of the crust, which corresponds to the variations
of the relief that the Bouguer correction removes, we obtain values of $2460 \text{ kg m}^{-3}$ for the linear model and
$2497 \text{ kg m}^{-3}$ for the exponential model, close to the estimate from the effective density spectrum. We discuss
density-depth relations in detail in Section 4. For results here, when computing Bouguer anomalies and
power spectra in the next sections, we will use a density of $2500 \text{ kg m}^{-3}$, consistent with our average from
the effective density spectrum and consistent with values used by Wieczorek et al. (2013) and Neumann
et al. (2015).

**3.3. Bouguer Anomalies**

The Bouguer field is computed by taking the difference between the free-air gravity as determined from the tracking data and the contribution by topography. For the latter, we again use topography-induced gravity as was used in our RM1 constraint. We use a constant density of 2500 kg m⁻³ (see above) for the crust. We then subtract the two sets of spherical harmonics and expand them to compute the radial derivative of gravity at the reference surface, which is called the Bouguer anomaly. We note that as $\lambda \to \infty$, the Bouguer field would become zero, if the correct density was used to compute topography-induced gravity, which would be $\alpha$ times the reference density used. We show the Bouguer anomalies for the RM1 model with $\lambda = 1$ in Figure 7, where the anomalies are expanded between degrees 7 and 900, in order to limit the influence from noise from the high degrees and signals from the deep interior from the low degrees, in particular, the degree-1 term of the Bouguer anomaly that compensates the degree-1 surface topography.

The Bouguer anomalies in Figure 7 show variations for the mascon basins and several other areas such as the Procellarum Kreep Terrane and around South Pole-Aitken, indicating densities different from the assumed value of 2500 kg m⁻³. Compared to the Bouguer map from Neumann et al. (2015), the map in Figure 7 is
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Bouguer anomalies for the GRGM1200B RM1 model with $\lambda = 1$. We used a constant crustal density of 2500 kg m$^{-3}$ to compute the Bouguer anomalies, and the anomalies were expanded between degrees $n = 7-900$. The map is in Mollweide projection centered on the farside.</figcaption></figure>

expanded to higher degree, yet conclusions for the basins from Neumann et al. (2015) remain unchanged. Our new Bouguer map can especially give more information on the density of smaller scale features, by virtue of the higher resolution.

In addition to the full Bouguer anomaly map shown in Figure 7, we also present band-filtered Bouguer maps. In Figure 8, we show the Bouguer anomalies filtered between degrees 150 and 1080, for the GRGM1200B standard Kaula-constrained model, the GL1500E model (which used a constant constraint), and the RM1 model with $\lambda = 1$. We have been using these filtered Bouguer anomaly maps to indicate processing issues which can show up as very localized features. This becomes especially clear for these band-filtered anomalies since most of the dominant signals will be filtered out. We use a cut-off here of $n=1080$ instead of 900 to emphasize the high-frequency part of the gravity field. As can be seen in Figure 8, the band-filtered Bouguer anomalies for the RM1 model are much smoother than those for the standard models GRGM1200B and GL1500E (by design). The minimum and maximum values of the Bouguer anomalies for the RM1 model are an order of magnitude smaller than those for the other models. And although GRGM1200B's extremes are larger, the statistics in terms of mean and RMS of the Bouguer anomalies for GRGM1200B and GL1500E are the same for this degree range.

The anomalies for the standard models show several locations with severe striping. GRGM1200B and GL1500E have striping in different areas, highlighting differences in processing between the two groups. The striping is mostly removed in the RM1 model: enforcing improved correlations with topography removes the striping by weighting the influence of KBRR residuals against the topography information. In that sense, the effect is the same as the neighbor smoothing constraints (Rowlands et al., 2010; Sabaka et al., 2010) that we applied in a separate, local analysis of GRAIL data (Goossens et al., 2014), showing that the neighbor constraints helped to improve the correlations with topography.

Apart from the striping, GL1500E (to be clear, here also expanded up to $n = 1080$ despite its larger size) shows more variations with larger anomalies in many places, likely due to its constant constraint. The standard models also show strong signals over parts of the farside. These are not present in the anomalies for the RM1 model. This indicates that for the larger part, correlations between gravity and topography are driven by the signal over the farside: the improved correlations for the RM1 model result in smoother Bouguer anomalies over the area with the strongest signal. Still, the RM1 model shows small scale differences in the farside area, indicating density differences with respect to the chosen (constant) density for the Bouguer correction, which we will investigate in section 4.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Band-filtered Bouguer anomalies, filtered between *n*=150–1080, for the GRGM1200B model (a), for GL1500E (b), and for the RM1 model with $\lambda$ = 1 (c). The assumed crustal density was again 2500 kg m⁻³. The minimum, maximum, and RMS values of the Bouguer anomalies are indicated beneath each map. The projection is the same as in Figure 7.</figcaption></figure>

In Figure 9, we show the Bouguer RMS power spectra plots, computed in the same way as those in Figure 2, for GRGM1200B and the RM1 models with $\lambda$ = 0.1, 1 and 10. The Bouguer spectrum for GRGM1200B is continuous and starts to increase again slightly before degree *n* = 600, as the power in the model (as can be seen in Figure 2) starts to decrease because of the Kaula constraint applied. The power for the RM1 model with $\lambda$ = 0.1 is close to that of GRGM1200B, but slightly lower at the higher degrees, due to the topography-based constraint. For the RM1 models with $\lambda$ = 1 and 10, there is a discontinuity at *n* = 600, the degree from which the constraint is applied. This discontinuity is stronger for the model with the stronger constraint (and it disappears for the model with $\lambda$ = 0.1). After degree *n* = 600, the RM1 free-air spectra follow power from topography closely, so the Bouguer power should be less than that for a model like GRGM1200B that quickly starts to drop in terms of its free-air power at that degree. The discontinuity in the power spectrum however could indicate that we could apply the constraint from an even earlier degree in order to have a smoother transition. We did not try this, and note that the free-air power spectrum is smooth, as is the effective density spectrum.

### 3.4. Data Fit

We also evaluate the new models for their performance in fitting the tracking data. We compute the *a* posteriori fit for the start model GRGM1200A, for the Kaula constrained model GRGM1200B, and for the RM1
<figure><img src="image_11.png" /><figcaption>**Figure 9.** RMS power spectra for the Bouguer expansion for GRGM1200B and its RM1 variants using $\lambda = 0.1, 1$ and $10$, and the GL1500E model. We used a constant density of 2500 kg m$^{-3}$ when computing the Bouguer models. The inset shows a close-up of the higher degrees.</figcaption></figure>

model using $\lambda = 1$. We present the RMS per arc of the fit to the Doppler and KBRR data for the extended mission part in Figure 10. We present only the XM data fits because of the lower satellite altitudes, meaning that these data have the best sensitivity to small scale gravity features. For most arcs, the improvements for GRGM1200B when compared to GRGM1200A are minor. The biggest improvements in Doppler fit occur for the late XM arcs, where the satellites reached their lowest altitudes. KBRR fit behaves similarly, with GRGM1200B showing a smoother data fit overall, as some small peaks in data fit during the earlier XM phase are reduced. One arc in particular in early December 2012 showed a larger than average fit for the GRGM1200A model, and GRGM1200B now performs much better for that arc.

The fit for the RM1 model is almost the same as that for the GRGM1200B model. For most arcs, the KBRR fit is slightly higher, but the differences are below noise level. We also inspected the values of the estimated empirical parameters (not shown), as they can also be used to indicate the performance of a model: smaller estimated values for these accelerations are interpreted as indications of an improved model. We found that the differences between the parameters for the two models are small. On average, most of the accelerations for the RM1 model are somewhat smaller than those for the GRGM1200B model, but the difference is much smaller than the actual level of the accelerations themselves. This indicates that the use of the RM1 constraint has only a small effect on the data fit. We note here that the constrained estimation problem of gravity field determination is described by a cost function that is to be minimized with two components: one describing the fit to the data and one describing the constraint, and their relative contribution is given by the relative weights (as described by $\lambda$). Indeed, fits improve somewhat (but again only a small amount) for RM1 models with smaller values of $\lambda$. Yet we want to stress that fitting the data is only a part of the goal to estimate a model for geophysical interpretation. We deem the stable effective density spectrum of Figure 4 more substantial than small differences in levels of data fit.

### 3.5. Anomaly Errors and Degree Strength

The associated errors of the models can be represented spatially by propagating the formal error covariance into a map of anomaly errors. However, for models of high degree and order, the determination and propagation of the error covariance matrix presents a considerable computational burden. Hence, an alternative method was developed for GRAIL, the details of which can be found in Sabaka et al. (2014). This method
<figure><img src="image_12.png" /><figcaption>**Figure 10.** RMS of Doppler (a, top) and KBRR (b, bottom) residuals per arc for GRGM1200A, GRGM1200B, and the RM1 model using $\lambda = 1$, for the extended mission phase.</figcaption></figure>

is based on the generation of *clone* models, see also section 3.2. If enough clones are used, they capture the same information as a full propagation of the error covariance. We have generated 1,001 clones for each model and used those to compute the anomaly error maps. Earlier tests have shown that this number of clones is enough to capture the details of the desired anomaly maps (Lemoine et al., 2014). In addition to the anomaly error maps, we also computed the degree strength of the models using these clones (generated from the constrained covariance), which was first derived for Venus gravity models by Konopliv et al. (1999). The degree strength at each location is computed by determining at which degree the anomaly error equals the expected signal from a given Kaula rule. In other words, it presents, spatially, the degrees at which the signal-to-noise ratio becomes 1. For the sake of convenience, we assume the same Kaula rule, $36 \cdot 10^{-5}/n^2$, when computing the degree strength for the RM1 models, even though these models are not derived with a Kaula rule. This likely does not affect the degree strength value significantly, as this is mostly determined by the level of the associated errors. If anything, the use of a Kaula rule for the RM1 degree strength computations will result in a conservative estimate, as Figure 2 clearly shows lower power in the Kaula-constrained model (which follows the prescribed Kaula rule at higher degrees), which means that the crossing with the RM1 error spectrum will occur at a lower degree value than if some estimate of the signal covariance for the RM1 model were used. We also account for the calibration factors as listed in section 3.1 when computing the anomaly errors and the degree strength.

We show the anomaly errors and degree strength in Figure 11. It is immediately clear that the RM1 model has smaller anomaly errors and higher degree strength values when compared to GRGM1200B, following the RM1 model's lower errors as shown in the spectrum in Figure 2. Despite calibration factors applied, it is difficult to determine if the models are properly calibrated because of the presence of systematic errors, although efforts to calibrate gravity models have been undertaken in the past (e.g., Lerch, 1991). We therefore focus more on the spatial differences between Figures 11a and 11b than on their amplitude differences. For the RM1 model, the (relatively) larger errors are concentrated around the equator, possibly because of larger spacing between satellite ground tracks. The equatorial band shows longitudinal variations in the maximum errors, the locations of which compare well with the locations for the largest errors for GRGM1200B. These larger errors are mostly driven by variations in satellite altitude, as is shown for example when Figure 11 is
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Anomaly errors from covariance up to degree *n* = 1200 (top, a and b), and degree strength (bottom, c and d). Results are shown for GRGM1200B and the GRGM1200B model with the RM1 constraint with *λ* = 1 (models indicated in the panels). The projection is the same as in Figure 7.</figcaption></figure>

compared with Figure 1 from Konopliv et al. (2014), which shows the minimum spacecraft altitude above lunar topography during the entire GRAIL mission. We do note here that the degree strength over Mare Orientale is only slightly higher (around degree 800) than that for its surroundings, despite the low-altitude campaign conducted towards the end of the GRAIL mission (e.g., Zuber et al., 2016). Those data are of course included here. Because this was not a repeated fly-over of the area but a campaign towards the end of the mission, the cross-track spacing is about 1° (a consequence of the Moon's rotation and the spacecraft's orbital period), and this may be responsible for a smaller increase of degree strength than may be expected, despite the much lower altitudes over the area.

Overall, the RM1 model shows a more uniform error map because the constraint utilizes the full topography spherical harmonic model with different coefficients for each degree and order, instead of a constraint that is applied degree-wise. This is consistent with the smaller band-filtered Bouguer anomalies from Figure 8 for the RM1 model. The constraint also removes striping from the solutions, and this also shows in the anomaly error map. The anomaly map for the Kaula-constrained GRGM1200B model also shows several isolated areas with locally larger errors. These areas coincide with locations where spurious KBRR residuals were deleted because of the high-frequency signal shown in Figure 1. The error map for the RM1 model shows these areas as well, but again at smaller amplitudes and due to the general smoothness of that map they are less visible.

The degree strength plots in Figures 11c and 11d are more alike, with higher degree strength values for the RM1 model because it has smaller errors. The degree strength results show that spatial variations are present in both models, despite the consistently high global correlations with topography for the RM1 model over the entire degree range. This will be further explored in the next section by means of localization.

**4. Lateral and Vertical Density Variations Using the RM1 Model**

With the improved correlations and stable effective density spectrum over the entire degree range, we now revisit the determination of lateral and vertical density variations in the lunar crust. We follow the original approach of Besserer et al. (2014), who used the earlier GRGM900C model to produce localized effective density spectra and then used these spectra to fit parameters of models of how density varies with depth. We will use our RM1 models to explore the impact of increased resolution on these estimated density structures.
In this section, we will first give a short summary of the theoretical models and localization used. We will then show results using our RM1 $\lambda$ = 10 model. We choose this model as it shows an effective density spectrum with smaller variations (see Figure 4) without being overconstrained. We will discuss the impact of choosing this model over others with smaller values of $\lambda$ later in section 5.

## 4.1. Theoretical Effective Density Spectra

Besserer et al. (2014) derive expressions for the effective density spectrum for two different models of vertical density variations: a model where density increases linearly with depth $z$, $\rho_{\text{lin}}(z)$ and a model where it increases exponentially, $\rho_{\text{exp}}(z)$. These one-dimensional density structures are written as:

$$ \rho_{\text{lin}}(z) = \rho_{\text{s,lin}} + az, \qquad (6) $$

$$ \rho_{\text{exp}}(z) = \rho_{\text{depth}} - (\rho_{\text{depth}} - \rho_{\text{s,exp}})e^{-z/d}, \qquad (7) $$

where $\rho_{\text{s,lin}}$ and $\rho_{\text{s,exp}}$ are surface densities for the linear and exponential models, respectively, $a$ is the density gradient for the linear model, $\rho_{\text{depth}}$ is the density at depth for the exponential model, and $d$ is a length scale (also called e-fold depth, $d > 0$) for the exponential model that indicates how quickly the value of $\rho_{\text{depth}}$ is reached (quicker for smaller values of $d$). Besserer et al. (2014) then derive the effective density spectra for each model as:

$$ \rho_{\text{eff,lin}}(k) = \rho_{\text{s,lin}} + \frac{a}{k}, \qquad (8) $$

$$ \rho_{\text{eff,exp}}(k) = \rho_{\text{s,exp}} + \frac{\rho_{\text{depth}} - \rho_{\text{s,exp}}}{1 + kd}. \qquad (9) $$

The variable $k$ is related to the spherical harmonic degree $n$ through

$$ k = \frac{\sqrt{n(n+1)}}{R}, \qquad (10) $$

where $R$ is the average radius of the planet. These expressions are derived assuming a Cartesian geometry. Han et al. (2014) used a spherical geometry and derived the effective density spectrum for any density-with-depth model $\rho(z)$ as:

$$ \rho_{\text{eff,spherical}}(n) = \sum_{p=1}^{P_{\text{max}}} \Delta\rho_p \left(\frac{r_p}{r_1}\right)^{n+2}, \qquad (11) $$

where $p$ is a counter for the number of interfaces, or layers (up to $P_{\text{max}}$), that are used for the crust, $\Delta\rho_p$ is defined as $\rho_p - \rho_{p-1}$ with $\rho_p$ the density at depth $z_p$, and $r_p$ is the corresponding radius (where $r_p < r_{p-1}$). The surface density is given by $\rho_1$ and $r_1$ is the reference radius. The density contrast $\Delta\rho_p$ can then be computed from equations (6) or (7), and by summing over the interfaces, the final effective density spectrum can be constructed. We used 1,000 layers in our implementation. Note that in the spherical case, the thickness of the crust is used when it is divided into $P_{\text{max}}$ layers to compute $z_p$ and $r_p$, whereas no crustal thickness is used as input in the Cartesian case (Besserer et al., 2014 derive the effective density expressions from an integral to infinity). We tested both expressions and for the value of crustal thickness adopted in this work of 40 km (Wieczorek et al., 2013), both spectra are very close, especially at higher degrees $n$, as expected. For the exponential model, we take $\rho_{\text{depth}}$ to be the (spatially varying) grain density derived in Huang and Wieczorek (2012) and used in Wieczorek et al. (2013) and taken from the GRAIL Crustal Thickness Archive (Wieczorek, 2012). This assumes explicitly that the porosity at depth is zero. The grain density of Huang and Wieczorek (2012) was determined by using a mineralogical norm approach, where an empirical relation between rock density and iron and titanium oxides was derived. The grain density map was then determined using Lunar Prospector gamma ray spectrometer iron and titanium abundance maps. After discussion of results and associated errors for both theoretical models, we will briefly comment on the differences between the two in section 5.4.

## 4.2. Estimation Parameters and Localization

Following Besserer et al. (2014) and Han et al. (2014), we can fit theoretical models to measure effective density with equations (8) and 9 or equation (11) and determine the parameters of interest: $\rho_{\text{s,lin}}$, $\rho_{\text{s,exp}}$, $a$,
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Maximum degree values used to fit the localized spectra, plotted geographically, for localization using a cap radius of 15° (a, left) and using a cap radius of 7.5° (b, right). The maximum degree is the degree at which the localized correlations drop below 0.9. The projection is the same as in Figure 7.</figcaption></figure>

and d. Lateral density variations in the surface densities and gradient or e-fold depth are obtained by means of localization. Following the work of Wieczorek & Simons (2005, 2007), local effective density spectra can be estimated by multiplying the data by one or more specially constructed localization windows and then expanding the result in spherical harmonics. We make use of their localization windows that concentrate optimally their power within a spherical cap of a given radius and with a maximum spherical harmonic bandwidth $L_{win}$. The windows are then rotated to the longitude and latitude of interest on the sphere. If more than one localization window is used, the individual effective density spectra are averaged. Besserer et al. (2014) estimated local spectra using a cap radius of 15°, with $L_{win} = 58$ to obtain 30 tapers (the windowing functions) with a concentration factor of each taper ≥ 0.99 (see Wieczorek & Simons, 2005). They then used the degree range $n = 250–550$ of the resulting effective density spectrum to estimate the density-depth parameters, by computing theoretical effective density spectra for a set of parameters and by comparing the computed spectra with the measured ones. By repeating this procedure for different locations on the surface, they thus obtain a spatial map of values for $\rho_{s,lin}$ and a, or $\rho_{s,exp}$ and d. These lateral variations come solely from the moving localization: the density-depth relationships themselves have a uniform density laterally since they only describe vertical variations. The size of the localization cap is then important, as for geologically complex areas, the assumption of laterally uniform properties will not hold if the area is too large, in which case forward-modeling efforts may be applied (e.g., Deutsch et al., 2019).

Our goal is to use the RM1 models to extend the resolution of this analysis. This can be done in different ways: by considering larger degrees $n$ when fitting the theoretical models and by considering smaller radii of the spherical caps used in the localization. For both, models with improved resolution, such as the RM1 models are necessary: effective density spectra that are stable up to higher degrees will allow a larger degree range for the fit; thus, increasing the number of observations (effective density at each degree $n$) for the fits, and smaller radii for the localization will result in larger values for the windowing degree $L_{win}$ if the same number of tapers will be used. Because windowing results in the coefficients being only reliable up to degree $n_{max} - L_{win}$ (Wieczorek & Simons, 2005), having a larger $n_{max}$ will help extend the resolution, especially if $L_{win}$ is larger. For the RM1 models, we will show that we can extend the analysis to an $n_{max}$, which increases more than the increase in $L_{win}$ needed for a smaller cap radius (see Figure 12).

Multitaper analysis is preferred in general because averaging the results from multiple orthogonal tapers produces a better estimate, as the energy coverage of the data over the region of interest will be more uniform (Wieczorek & Simons, 2005). By exploiting the increase in resolution of our models in this analysis, we expect improved estimates of especially the surface densities: both theoretical models of effective density taper off towards the surface density values $\rho_{s,lin}$ or $\rho_{s,exp}$ for large degrees $n$ (and as a result, for relatively larger values of $k$), as can be seen from equations (8) and (9). Thus, having more data at higher degrees is expected to improve the surface density estimates.

In the following parts, we present results for lateral and vertical density variations, where we extend the resolution of the analysis in various ways. We start with the Besserer et al. (2014) localization parameters: a cap radius of 15° with a spectral bandwidth of $L_{win}=58$, which allows us to use 30 tapers to estimate the localized effective density spectrum at a given location (in latitude and longitude). We compute these spectra
at intervals of 5° spatially. We used the freely available SHTOOLS software package for the localized analysis (Wieczorek & Meschede, 2018). We fit theoretical spectra to determine values ρ<sub>s,lin</sub>, ρ<sub>s,exp</sub>, a, and d, by initially using the degree range *n* = 250–650, extending the range from the value of 550 used by Besserer et al. (2014). We note that because L<sub>win</sub> = 58, the spectra cover the degree range *n* = 59–1142. Next, we decrease the cap radius to 7.5°. In order to use the same number of tapers with a concentration factor ≥ 0.99, we need to increase L<sub>win</sub> to 116, which does not pose a problem for the RM1 models (in terms of narrowing the usable degree range too much). For this set of localized spectra (where, because of the chosen value of L<sub>win</sub> = 116, the spectra now cover the range *n* = 117–1084), we again first fit theoretical spectra over the fixed degree range *n* = 250–650. In addition, we also vary the maximum degree for the fit according to each measured localized correlation spectrum (but we keep the minimum always at *n* = 250): we determine *n*<sub>max</sub> to be the degree value where the localized correlations drop below 0.9. This corresponds to a signal-to-noise ratio of about 4.26, if it is assumed that deviations from unity correlations are due to gravity signals uncorrelated with topography (Wieczorek, 2008). This results in a set of spatially varying maximum degree values, which we show in Figure 12. In this Figure, we interpolated the results from the spectra spaced by 5° into a map of resolution 1° by 1°, centered in this case on the farside. We show two maps, since the different radii for the spherical caps used result in different localized correlation spectra.

Both maps are similar, with the map for the cap radius of 7.5° having in general slightly lower maximum degree values. As expected, the largest maximum degree values can be found over the farside highlands, the area where correlations between gravity and topography are highest and which are closest to where the measurements were made by the orbiting spacecraft. There, we greatly exceed previously used maximum degree values, and the increase in L<sub>win</sub> does not have an effect. Effectively, we can use the entire degree range for the farside area. Values lower than n<sub>max</sub> = 650 are found in areas on the nearside (mostly because of the presence of the mare regions) and in South Pole-Aitken (SPA). For the latter, correlations are lower because of decreased data sensitivity, as the distance from the spacecraft to the surface is larger over SPA. In the interior of SPA, the maximum degree values are close to 560. These variations in maximum degree are also likely related to the assumption of uniform properties in each spatially windowed area as mentioned above: this assumption is more likely to be true for the battered highlands than for the geologically complex mare areas. Finally, the results in Figure 12 depend to some extent on the value of λ for the RM1 model, as models with large values for λ would reach maximum resolution everywhere. The variations in Figure 12 thus show that the chosen model using λ = 10 does not appear to be overly constrained.

The maximum degree maps presented here can also be compared to the degree strength maps as presented in Figure 11. The degree strength maps are much more related to data coverage and spacing, with the largest values towards the poles, and in fact lower values over the central farside. While signal in Figure 12 correlates with specific regions on the Moon (mare, farside highlands, SPA), this is not the case for the degree strength, which is dominated by long-wavelength signals.

Finally, we obtain values for ρ<sub>s,lin</sub>, ρ<sub>s,exp</sub>, a, and d by means of least-squares. The exponential model is non-linear for variable *d* (see equation (9)), and so, we used iterated least-squares; we select starting values for the estimated parameters, then compute the difference between the measured and computed spectra and adjust the parameters by setting up a linearized normal equation system. This procedure is iterated until the relative change in the RMS of the misfit between subsequent iterations is less than or equal to 0.01%. In the case of the exponential model, we use a step-size of 0.1 in our least-squares solutions to assure stability. This means that for each iteration, we update the parameters with 0.1 times the adjustment (instead of using the full adjustment, or a step-size of 1). We found that while this increases the number of iterations, the solutions are more stable as with a step-size of 1 they often would not converge. We tested our implementation using the same parameter settings as Besserer et al. (2014) and obtained the same results.

**4.3. Results for the Linear Model**

We present the results for the surface density ρ<sub>s,lin</sub> and the gradient value a in Figure 13. Results for the extended degree range *n* = 250–650 using the same localization as Besserer et al. (2014) are very close to those original results, thus establishing the baseline. Density inversions, where density decreases with depth, can be found in the mare regions from the negative values for the gradient term a. As the resolution is increased by using smaller caps but still the same degree range (Figures 13c and 13d), the maps show more
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Results for the surface density (left column) and gradient (right column) for the linear model. Different localization and maximum degree values for the fit were applied: a cap of radius 15° with the degree range *n* = 250–650 (a,b), a cap of 7.5° with the degree range *n* = 250–650 (c,d), and a cap of 7.5° with the varying degree range from Figure 12 (e,f). Boundaries of the mare regions are included in the maps. The projection is the same as in Figure 7.</figcaption></figure>

variations. Several locations now show lower surface densities. These remain consistent as the resolution is further enhanced by varying the maximum degree values (Figures 13e and 13f). SPA, in general, has a higher density than its surroundings, as a result of a more iron-rich composition (see also Wieczorek et al., 2013 and Besserer et al., 2014), yet there is a portion in the northern part with lower densities. As such, the SPA area stands out from its surroundings in the map. These results for SPA hold for the varying maximum degree map, despite the fact that the values in SPA are lower than 650, as stated in the previous section. Mare Orientale also stands out in Figure 13c, as having a larger surface density, yet this area is somewhat smaller in the varying maximum degree map, possibly owing to slightly lower maximum degree values at roughly the center of the basin, as seen in Figure 12.

The gradient values in Figures 13b, 13d, and 13f show a similar grouping as the surface density maps. Areas of density inversion are more fragmented in the mare regions but overall consistent with earlier results. As resolution is increased, a few new areas with density inversions appear, most notably in the area of the Apollo basin in SPA, and on the edge of the farside highlands on the northern hemisphere. Overall, the gradients in the highlands are higher than those for the earlier results.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Results for the surface density (left column) and e-fold depth (right column) for the exponential model. The same localization as for the panels in Figure 13 applies. The map projections are also the same, and the mare boundaries are again included. Gray areas in the e-fold maps indicate locations where in the linear model $a < 0$, or where a reasonable fit for $d$ was not possible (see text for further explanations).</figcaption></figure>

### 4.4. Results for the Exponential Model

Results for the exponential model's surface density $\rho_{s,exp}$ and the e-fold depth value $d$ are presented in Figure 14. The results using the same localization as Besserer et al. (2014) (a cap of 15° radius) but extending the degree range to $n = 250–650$ again match the original results closely, with only a few departures.

We note that Besserer et al. (2014) did not show the surface densities for the exponential model. Surface densities for the exponential model as shown in Figures 14a, 14c, and 14e are similar to those for the linear model, with the exception that the surface densities for the exponential model are in general lower than those for the linear one. This can be traced to the equations for the effective density as given by equations 8 and 9. It can be shown that if both models describe the same measured spectrum with the same slope, then the surface densities for the models can be related to each other by:

$$ \rho_{s,lin} = \rho_{s,exp} + \frac{a}{dk^2}. \qquad (12) $$

Thus, for spectra where $a > 0$ and $d > 0$ (the latter should always be the case), the linear surface density will be larger.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Effective density spectrum for the RM1 model using $\lambda = 1$ (black). To show variations in effective density caused by variations in the gravity field coefficients, we include the spectra of 1,001 clones for the model using $\lambda = 1$ (dark gray) and spectra of 1,001 clones for the models using $\lambda = 0.1$ (light gray). We also show the result of fitting linear functions with two different slopes to the spectra, with the best fit for the lower degree range in red and that for the higher range in green (see also text, and Figure 16).</figcaption></figure>

The surface densities for the exponential model show larger variations than the linear model does. For both cases, SPA stands out, with larger densities in parts of the interior of the basin, but the exponential model shows a larger extent for the surrounding lower densities. Both the exponential and linear surface densities for Mare Orientale are larger than the densities of the surroundings. The exponential model also shows a band of lower densities on the northern farside. To some extent, this is present in the linear surface densities as well, but much less pronounced and, again, at higher values.

In Figure 14, gray areas in the e-fold maps indicate areas for which the linear models obtain negative gradients, $a < 0$. This is close to the masking used in Besserer et al. (2014), where values with $a < 5$ were excluded, which mainly affected the mare regions, as the exponential model is assumed to be a better fit for non-mare areas. The exponential model was designed for density increasing with depth. Density inversions can in principle be described by simply having a surface density larger than the grain density at depth, $\rho_{\text{depth}}$, yet we find issues with the least-squares fit, in that we cannot find reasonable density inversions. Even when choosing a start value in the least-squares fit of $\rho_{\text{s,exp}} > \rho_{\text{depth}}$, this results in either negative values or abnormally large values for $d$. This instability that sometimes occurs when fitting the exponential model is also the reason why we use a step-size of 0.1 in the iterated least-squares fits, as stated in section 4.2. We note that Gong et al. (2016) used a simple two-layer model to invert for mare thickness.

In our results with the extended degree range and smaller localization cap radius, the Apollo basin again stands out with its values for the surface density and the e-fold depth. SPA itself stands out as an area with low values for $d$, as originally indicated by Besserer et al. (2014), but other parts of the farside now reach similar low $d$ values. Overall, the signatures in the linear gradient and exponential $d$ maps are very similar, and we note that higher linear gradients are in principle similar to low $d$ values, as both reach higher densities at relatively shallow depths.

## 5. Discussion

We now take a closer look at the results presented for the RM1 models. First, we discuss in more detail the global effective density spectrum, by looking at variations arising from the errors in the gravity field coefficients. We also discuss the apparent change of slope in the global effective density spectrum. We discuss the influence of the choice of the value for $\lambda$ for the RM1 model we used in our analyses: do the results for surface density and gradients change when we choose a model that is constrained less strongly? We then present an analysis of the expected errors in surface density and trends. Following these analyses, we briefly comment on the differences between the linear and exponential models. We show the resulting surface porosities for the linear model, and finally, we discuss differences in overall surface density for three distinct regions: the farside highlands, SPA, and the mare regions.

### 5.1. Global Effective Density Spectrum: Clones and Slope

We introduced the clone models in section 3.5, where they were used to assess anomaly errors, circumventing the need to propagate the entire covariance matrix. Here we use the same clone models to assess the variations in the global effective density spectrum as presented in Figures 4b and 4d. We use 1,001 clones for two RM1 models: with $\lambda = 0.1$ and 1. We use the less constrained models to obtain larger variations than the more tightly constrained model using $\lambda = 10$ would show. The results are shown in Figure 15.

The variations in effective density are relatively small up to degree $n = 750$, and they grow larger for higher degrees, as expected since the associated errors for those coefficients are larger. The variations for the model using $\lambda = 0.1$ are higher than those for the model using $\lambda = 1$, due to the looser constraint. As stated before, despite decreasing correlations with topography for higher degrees for the $\lambda = 0.1$ model, its effective density spectrum remains stable (Figure 4d). The spectra for the clones also remain stable. A spike in the effective density spectrum can be seen around degree $n = 1080$, especially for the $\lambda = 1$ model. This degree value is
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Histograms of the values of the estimated slopes from 1,001 clones (a), separated into a slope for a lower degree range ($n = 250 - n_{break}$, red) and a slope for a higher degree range ($n = n_{break} - 1000$, green), and histograms of the degree $n_{break}$ where the break occurs (b).</figcaption></figure>

a multiple of 360, the degree around which resonances can occur, as stated in section 3.1. Despite streaks in the coefficient differences around degrees 360, 540, and 720, as shown in Figure 3, we do not see similar spikes in the variations of the effective density spectrum, although some level of variation in the effective density spectrum is associated with these degree values.

The effective density spectrum shows signs of having a break in the slope, where the higher degrees appear to have a shallower slope than lower degrees (this was of course already noticeable in Figure 4). To investigate whether this slope is indeed present in the spectrum, we fit a linear function with two separate slopes to the effective density spectrum, where we also estimate the value of the degree of the breakpoint, $n_{break}$, using iterative least-squares. We note that the nature of the constraint can easily bias the results obtained here, as stronger constraints result in flat effective density spectra (Goossens et al., 2017), thus guaranteeing that there will be a difference in slope, even for less strongly constrained solutions. We will use the RM1 model with $\lambda = 1$ here, to ensure that the constraint influence is not too strong, and we will discuss results for models with $\lambda = 0.1$ and 10.

We fit the functions using the effective density values between degrees $n = 250$ and 1000. The estimated parameters are the effective density at degree $n = 250$, the two slopes, and $n_{break}$. We choose the initial value of $n_{break}$ at 450, well below the degree $n > 600$ for which the RM1 constraint is applied, so as not to bias the results for this parameter. We performed this estimation for all 1,001 clones, and we show the histograms of the values for the slopes and the breakpoint in Figure 16.

The slopes for the lower degree range and the higher degree are clearly separated, for all clones. This indicates that a separate slope is indeed present. We also further confirmed this by applying an F-test that showed that the two-slope model indeed has a significantly better fit than a one-slope model. The histogram for the slope for the higher degree range has a wider spread because the effective density spectrum shows larger variations for that range due to the increased uncertainty of the gravity field coefficients. The value for $n_{break}$
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Results for surface density when using the RM1 model with $\lambda = 1$ instead of 10 (left column, a and c), and differences with the surface densities from the RM1 $\lambda = 10$ model (right column, b and d). The top row shows results for the linear model, the bottom row for the exponential model. The localization used a cap size of 7.5° and the fit used the degree range $n = 250–650$. The projection is the same as in Figure 7.</figcaption></figure>

shows a wide range, with values between $n_{break} = 610$ and 789. We note that care needs to be taken when interpreting these results, considering that the constraint is applied for $n > 600$. As discussed above, the effective density spectrum becomes flat when the constraint factor $\lambda$ is large, since the lateral density variations disappear and only the scale factor $\alpha$ on a constant density remains. In such cases, one would expect a change in slope, and as the constraint is applied stronger, it will be exactly at the degree from which the constraint is applied. The spread for $n_{break}$ however is wide, and starts always after $n=600$, showing it is not necessarily confined to be close to the degree from which the constraint is applied. We also performed the same analysis on clones for the $\lambda = 0.1$ and 10 models. For both, the two slopes are still clearly separated. For the $\lambda=0.1$ model, the range of slopes for the higher degree part becomes wider, as expected since it shows more variations (as also indicated by the variations from the clones as shown in Figure 15). The values for $n_{break}$ also show a wider spread. For the $\lambda=10$ model the histograms are narrower. Yet even for this model, the slope for the higher degrees is not significantly closer to 0 (which we would otherwise expect if the constraint dominated).

We performed this analysis to investigate the behavior of the effective density and to investigate whether or not the trend observed in earlier GRAIL models would continue or if the spectrum will gradually flatten. One should expect the effective density spectrum to level off because higher degrees mostly describe the shallower parts of the crust. Hence, as discussed in section 3.2, we can also estimate a global density-depth relation, where the higher degree part of the effective density spectrum constrains the value for the surface density. For high degrees, the effective density should approach the density of the uppermost layer, and this is not likely to be extremely low, which would be the case if the slope in effective density from the lower degrees is continued. We thus interpret the separate slopes that we find as being indicative of a global average of the surface density.

## 5.2. Influence of Constraint Factor $\lambda$ on the Density Variations

In sections 4.3 and 4.4, we determined the surface density distribution and vertical gradients from the RM1 $\lambda = 10$ model. Here we investigate to what extent the results change if we choose a model with a looser constraint, by redoing the density analysis with the RM1 model using $\lambda = 1$. We compare results using the fixed degree range $n = 250–650$ with localization using a cap radius of 7.5°, assuming that we will see similar
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Assessment of the errors in the estimated parameters: surface density, gradient, and e-fold depth. The left column corresponds to results for the linear density-depth model, and the right column corresponds to the exponential model. The first row (a, b) shows the values of the misfit between the localized effective density spectra and the best-fit theoretical models. The middle row (c, d) shows the errors in the recovered surface density from the clones. The bottom row (e, f) shows the errors in the gradients: the linear gradient (e) and the e-fold depth (f). The projection is the same as in Figure 7.</figcaption></figure>

differences when using different localization or degree ranges. In Figure 17, we show the surface densities for the linear and exponential models for the RM1 model with $\lambda = 1$ and the differences with the $\lambda = 10$ model.

The differences for the surface densities for the linear model, as shown in Figure 17b, are small, indicating that the results for the linear model from both RM1 models are close together. Those for the exponential model are higher and have more (spatial) variation. The median of the differences is only $1 \text{ kg m}^{-3}$ for the linear model, and $4 \text{ kg m}^{-3}$ for the exponential model (although the spread in values is larger, especially for the exponential model). For areas with relatively large differences, such as occurring in SPA, the maximum differences are about $100 \text{ kg m}^{-3}$ for the linear model, which is still about 5% or less of the nominal value.
We found that the exponential model is more difficult to fit, despite using iterative least-squares, resulting in unrealistically large surface density differences in some areas (such as SPA), despite the small median value of differences. Small changes in start values, or in the localized effective density spectrum itself, can sometimes lead to relatively large changes in the results, and hence the wider range of density differences. While a grid-based search for the minimum misfit, as performed in the work of Besserer et al. (2014), would always produce some sort of minimum for the parameter space included, we find that a least-squares approach here is insightful in highlighting the differences between the depth-variation models, and in highlighting such shortcomings as unstable fits for the exponential model. Locations where the differences between the $\lambda = 1$ or 10 models are largest roughly coincide with locations where the e-fold depth could not be determined (see Figure 14).

Overall however, as stated before, the density patterns from the linear and exponential model are very similar, as can be seen when comparing Figures 17a and 17c, with the exponential model showing on average lower surface densities, as explained in section 4.4. The maps from Figure 17 ($\lambda = 1$) also compare well with those shown in Figures 13 (linear model, $\lambda = 10$) and 14 (exponential model, $\lambda = 10$). Differences for the gradients and e-fold depth are not shown in Figure 17, but they show the same pattern; they are small for the linear model except for in a few locations, and they are larger for the exponential model, but overall, the patterns are very similar. Concluding, the results for both constraints, $\lambda = 1$ or 10, are very similar, justifying our use of the $\lambda = 10$ model.

### 5.3. Error Analysis of Recovered Surface Density and Gradient Parameters Using Clone Models

We can assess the results of the least-squares fit in different ways: in terms of the misfit between measured and theoretical best-fit spectrum, and in terms of the expected errors of the recovered parameters. The latter can be computed from the resulting covariance of the least-squares fit. However, this would not take into account the effect of errors in the gravity field itself. Instead, we opt to use clone models in the same way as we did when computing the anomaly errors and degree strength in section 3.5. We perform the localized analysis to fit effective density spectra on a set of clone models, determining maps of surface density and gradient parameters for each clone. From the variations between the results for the clone models, we can compute the variations in the recovered parameters. In this case, we used only 46 randomly chosen clones instead of the full set of 1,001. Generating results for each model is time-consuming because of the localization involved. We compared the average surface density and gradients (or e-fold) maps from the 46 clones and found that they were very close to the map from the original model. This gives us confidence in the derived parameter variations, while noting that strictly speaking we only find some lower bound for the parameter errors, as the use of more clones is likely to increase the variations. For the results shown here, we localized the models using a cap size of 7.5°, and we fitted the effective density spectra between degrees $n = 250$ and 650.

In Figure 18, we show the results of this analysis with the clones. Several areas in the maps for the exponential model (the right column of the Figure) are gray because a fit could not be obtained for those spectra, as discussed in the previous section. The misfit maps (Figures 18a and 18b) are very similar for both density-depth models. The largest errors occur mainly in the mare areas, where density inversions occur. A noticeable difference between the linear and exponential misfit values can be seen for Mare Orientale, where the exponential model shows much larger misfit values. The maps for the errors for the surface densities (Figures 18c and 18d) also are similar, when compared between the theoretical models as well as when compared with the misfit maps. The linear model shows overall slightly lower misfit values when compared to the exponential model. The lowest values for both occur over the farside, where correlations with topography are the highest. Both show increased errors in SPA, likely owing to the relatively higher spacecraft altitude. The maps for the gradient (Figure 18e) and e-fold depth (Figure 18f) are similar as well.

The misfit and the surface density both use the unit for density. The misfit however is the misfit for the *effective* density and thus, indicates an average kind of density (throughout the crust). The error values for the surface density are relatively small, when compared to the level of the densities themselves, as shown in Figures 13 and 14: the mean is $11 \text{ kg m}^{-3}$ for the linear model and $17 \text{ kg m}^{-3}$ for the exponential model, respectively, whereas the surface densities are two orders of magnitude larger. The mean values for the gradients/e-fold depths are similar, at $2.3 \text{ kg m}^{-3} \text{ km}^{-1}$ for the linear model, and $2.8 \text{ km}$ for the exponential model, respectively, but the spread is larger for the exponential model. Finally, the misfit has a mean of $22 \text{ kg m}^{-3}$ for the linear model and a mean of $25 \text{ kg m}^{-3}$ for the exponential model, respectively. These results
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Surface porosity for the model of linear density-depth variations, where we used a cap size of 7.5° for the localization, and the degree range 250–650 for the least-squares fit. Areas in gray are those where the density is larger than the grain density (porosity < 0). The projection is the same as in Figure 7.</figcaption></figure>

can be combined with those shown in Figures 13 and 14 to describe variations in the density structure of the lunar crust, together with (lower) error bounds.

## 5.4. Comments on the Linear Model Versus the Exponential Model for Depth Variations

Since we have compared results for the linear and exponential theoretical models for density-depth variations in the lunar crust, it is natural to ask which, if either, performs better. However, based on the performance of the theoretical models as shown in the previous section, we cannot readily choose one over the other. While the exponential model may have a slightly higher mean misfit value, the distribution of misfit values for both models is nearly identical. The slightly higher mean for the exponential model is merely a result of that model showing a longer tail towards higher misfit values, owing to our finding that it was more difficult to fit despite our use of iterative least-squares. The variations in the recovered parameters from the clones do look slightly better for the linear model, but the difference between the two appears to be small. It thus becomes more a matter of preference. The exponential model has a more physical meaning (e.g., Besserer et al., 2014; Han et al., 2014) as it is related to an increase in density with depth due to compaction, but it also shows lower surface density values when compared to the linear model. The surface densities of the linear model appear more reasonable. The linear model can also be used more readily to handle density inversions and to discriminate between mare and non-mare regions (Besserer et al., 2014).

## 5.5. Surface Porosity for the Linear Model

Knowledge of the (surface) porosity of the lunar crust is important, since it controls the gravity signature of impacts (e.g., Milbury et al., 2015). We can turn the values for surface density, as presented in either Figure 13 or Figure 14, into porosities in a straightforward way, since:

$$ \rho_s = (1 - \phi)\rho_{\text{grain}}, \qquad (13) $$

where $\rho_s$ is the surface density, $\phi$ the porosity, and $\rho_{\text{grain}}$ the grain density (which is equal to $\rho_{\text{depth}}$ from equation (7)). We again use the grain density from Huang and Wieczorek (2012). In this case, we only present the results for the linear model as shown in Figure 13e, where we used a cap size of 7.5° in the localization and where we used the degree range 250–650 in the least-squares fit. We choose to show this model because the exponential model results in lower surface densities, which would imply even higher porosities. We show the resulting surface porosities in Figure 19. We did not filter or mask any particular area, and areas with surface densities higher than the derived grain density (such as occurs in the mare regions) are shown here with a porosity of 0%.

Our resulting porosities are in general consistent with those presented in earlier works (Besserer et al., 2014; Han et al., 2014; Wieczorek et al., 2013), in terms of the spatial variations. The increased resolution of the models allows us to use a higher degree value in the fit, which in turn means a better sensitivity to shallower
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Histograms for the surface densities inside the farside highlands, SPA, and nearside mare regions on the Moon. We used the results of the linear model, with a spherical cap of 7.5° and the degree range 250–650. (see Figures 13c and 13d).</figcaption></figure>

features, and we thus obtain slightly lower surface densities than earlier studies. This results in slightly higher surface porosities, especially over the farside. We also note that Apollo in SPA and Mare Orientale show 0% porosity because of their higher surface densities.

**5.6. Density Parameters for Distinct Regions: The Farside Highlands, SPA, and Mare Areas**

We now investigate differences in density, gradient, and porosity, for three distinct regions on the Moon: the farside highlands, SPA, and the nearside mare regions. We use results from a spherical cap of 7.5° for the localization and using the degree range 250–650 to determine the model parameters. For each region, we select points well within that region for which we have localized effective density spectra. For the highlands, we used a circle with a radius of 20° centered on 180° E, 45° N, as this is an area in the middle of the farside highlands, and we thus take it to be representative for the highlands crust. For SPA, we used a circle with a radius of 15° centered on 191.1° E, 53.2° S (the center of the best fitting ellipse for SPA following Garrick-Bethell and Zuber (2009); the chosen circle fits well within this ellipse), and for the mare regions, we selected points that were at least 10° away from the mare boundaries in order to minimize the effects of having spectra that contain a mix of geologically distinct areas. We have on average 50 effective density spectra per area, from our original 5° spacing (see section 4.2). For each effective density spectrum at such a location, we take the estimated value for surface density and gradient and their associated errors, and then we generate 1,000 points per location from a Gaussian distribution, in order to populate a histogram. We perform this analysis for both the linear and exponential model. We show the resulting histogram for the densities of the linear model in Figure 20, and we list the mean and standard deviation for each parameter (obtained from the histogram values) in Table 2 for both models.

The results show that, as expected, the highlands densities are well-constrained. The results for the exponential model are in agreement with those from Besserer et al. (2014). The linear model finds a higher surface density value for SPA than the exponential model; however, we note that the spread for the surface densities for the exponential model is exceptionally large. In addition, the small e-fold value for SPA for the exponential model means that density quickly increases with depth. The signal for the mare regions is not as well-constrained as that for the other regions in the linear model, and we do not include results for the mare areas for the exponential model. Most of the effective density spectra in the mare regions are not well-determined, and correlations between gravity and topography tend to be lower than 0.9 or even 0.8. In Figure 20, we only include results for spectra with negative slopes. If we expand the spectra used for Figure 20, we find another peak at lower densities. The current peak at 2300 kg m⁻³ stems from one spectrum with a low surface density. The mare areas are also flatter than the highlands, resulting in a weaker gravity signal (SPA suffers from this as well because of its low elevation). The mare regions are also more
<table><caption>**Table 2**
Mean Values and Standard Deviations of Surface Densities, Gradients, and Porosity for Different Regions</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">Linear model</th></tr><tr><th>Surface density [kg m<sup>-3</sup>]</th><th>Gradient [kg m<sup>-3</sup> km<sup>-1</sup>]</th><th>Porosity [%]</th></tr></thead><tbody><tr><td>Farside highlands</td><td>2324 ± 59</td><td>27 ± 11</td><td>19.7 ± 2</td></tr><tr><td>SPA</td><td>2411 ± 117</td><td>49 ± 23</td><td>19.7 ± 4</td></tr><tr><td>Nearside maria</td><td>2836 ± 261</td><td>-94 ± 78</td><td>3.7 ± 9</td></tr><tr><th rowspan="2">Region</th><th colspan="3">Exponential model</th></tr><tr><th>Surface Density [kg m<sup>-3</sup>]</th><th>E-fold depth [km]</th><th>Porosity [%]</th></tr><tr><td>Farside highlands</td><td>2264 ± 84</td><td>17 ± 15</td><td>21.8 ± 3</td></tr><tr><td>SPA</td><td>2035 ± 637</td><td>1.7 ± 0.9</td><td>32 ± 21</td></tr><tr><td>Nearside maria</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

geologically heterogeneous when compared to the highlands. The resulting porosity for the mare regions is about half that of the other areas, at 10% compared to 19% for the highlands and SPA.

## 6. Conclusions

We have determined models of the lunar gravity field of degree and order 1200 in spherical harmonics from GRAIL data. Apart from determining a model with the standard spectral constraint applied in gravity field determination (called the Kaula constraint), we also determined models with a constraint based on topography information. This constraint was introduced in Goossens et al. (2017), and it aims to improve correlations between gravity and topography. It has been designed such that a scale factor between gravity and topography is free to adjust. This scale factor is fully determined by the data instead of by the constraint, and the constraint is named RM1 (rank-minus-one). The scale factor is interpreted to be related to the average bulk density of the crust. Both the Kaula constraint and the RM1 constraint are applied for degrees larger than 600. We determined several RM1 models by applying a different weight factor $\lambda$ to the constraint system in the normal equation system used to solve the least-squares problem of determining a gravity field model from satellite tracking data. A larger value for $\lambda$ indicates a stronger constraint. We showed that applying this constraint results in models with high correlations between gravity and topography, as the constraint was designed to do. We have evaluated the RM1 models in terms of their power spectra, anomaly errors, degree strength, and fit to the data, comparing them with the Kaula-constrained model.

In addition, the RM1 models were shown to have stable effective density spectra over the entire degree range, whereas the Kaula-constrained model can only be used globally up to a degree of about $n = 700$, after which the correlations with topography, and the effective density spectrum, decrease rapidly. By averaging the effective density spectrum over the degree range $n=250-1200$, we find a bulk crustal density of $2357 \pm 82$ kg m$^{-3}$, slightly lower than reported before. Variations in the effective density spectrum can be used to determine lateral and vertical density profiles in the crust. We used a localized analysis to determine the effective density spectrum for areas on the Moon, and we then fitted theoretical models of density variations with depth to this localized spectrum. We used theoretical models of linear and exponential variations with depth and determined the surface density and linear gradient or e-fold depth from the localized spectra. Because of the increased resolution of our gravity field model and because the effective density spectrum is now stable over the entire degree range, we were able to improve the resolution of this analysis when compared to previous work in two ways: by using smaller areas in the localization (spherical caps with a radius of 7.5° instead of 15°), and by extending the degree range over which the fit was performed. From the localized correlation spectra, we determined the degree where correlations drop below a value of 0.9 (corresponding to a signal-to-noise ratio of about 4.26), and the results show that we can use the entire model for most of the farside. For parts of the nearside, as well as for SPA, the maximum resolution is around degree 550. The used density-depth relationships assume laterally uniform properties, and thus being able to reduce the size of the localization area is important. This uniform property assumption will work better in an area such as the highlands than in geologically complex areas such as the nearside maria.

We presented maps of surface density and density gradients in the lunar crust. Our analysis is broadly consistent with earlier findings, with more small-scale variations due to the increased resolution. Results for
the linear model are very similar in pattern to those from the exponential model (for the surface densities as well as for the vertical profiles). We show that the exponential model mostly results in lower surface density values than the linear model, if both models fit the same slope of the effective density spectrum. In terms of the resulting surface densities, the interior of SPA is more pronounced, with higher surface densities than its surroundings, with lower densities in the northern part. The SPA area is distinct from its surroundings. Mare Orientale also stands out as having a higher surface density. In the mare regions, areas of density inversions are more scattered than found before. As we increase the resolution, a few new areas with inversions appear, most notably in the vicinity of the Apollo basin in SPA, and on the edge of the farside highlands in the northern hemisphere. The Apollo basin has mare materials in its interior, and a negative gradient here is consistent with Besserer et al. (2014) and Gong et al. (2016). The gradients in the highlands are larger than those for the earlier results. An error analysis using clone models (an ensemble of solutions of the same statistical family, which allows us to estimate errors without having to propagate the entire covariance matrix) shows that the surface densities and vertical profiles are determined best over the farside highlands, where uncertainties as low as a few kg m⁻³ are obtained. The misfit between measured and modeled localized spectra results in a mean of about 25 kg m⁻³. These results are to be interpreted as lower bounds on the errors because we only used 46 clones due to computational constraints. Nonetheless, the resulting average surface density maps from these clones is very close to that of the model, noting that the expected value of the clones is the model itself by definition. Our results also indicate differences in surface density, density gradient, and porosity between three different regions on the Moon: the farside highlands, SPA, and the mare regions. We find that the highlands surface density is close to the global surface density and that SPA has a slightly higher density. Both regions have a similar porosity of 19%. The mare regions have a density that in general is higher, and their porosity is around 10% (but with large uncertainty).

The theoretical spectra are not linear in their dependence on spherical harmonic degree, and they approach the surface density for high degrees. Inspection of the global effective density spectrum shows an apparent change in slope for the higher degrees, where the slope is smaller. To characterize a possible change in slope, we fitted a function with two slopes and a breakpoint to the global effective density spectrum. We repeated this for a set of 1,001 clone models, and the results show that indeed the two slopes are significantly different, with the slope for the higher degrees smaller than that for the lower degrees, and with the breakpoint occurring after the degree from which the constraint is applied (n = 600). This result suggests that the effective density spectrum reaches a global average of the surface density. It should be noted that the RM1 constraint, when applied strongly, will result in an effective density spectrum with zero slope for degrees larger than the degree from which the constraint is applied, and there is thus the danger of circular reasoning. We tested models using λ = 0.1 and 10 and both still show a separation of slopes, with non-zero slopes for the higher degrees. In addition, we showed the effective density spectra for clone models to assess the uncertainty in effective density. While the spectra for the clones from the λ = 0.1 models show larger variations than those for the λ = 1 model, they are stable over the entire degree range and do not collapse like the spectra for Kaula-constrained models do.

Locally higher resolution models may be extracted from the GRAIL data, owing to the relatively higher fit to the inter-satellite tracking data for the final part of the mission with very low altitudes above the surface, which indicates that there is still gravity information present. However, it might be difficult to accommodate such signal in a global model. We have shown that our RM1 models can overcome this, resulting in models with stable effective density spectra up to high degrees. That makes these models especially valuable in (localized) spectral studies of the structure of the lunar crust.

References

Andrews-Hanna, J. C., Asmar, S. W., Head, J. W. III, Kiefer, W. S., Konopliv, A. S., Lemoine, F. G., et al. (2013). Ancient igneous intrusions and early expansion of the Moon revealed by GRAIL gravity gradiometry. *Science*, **339**, 675–678. https://doi.org/10.1126/science.1231753

Andrews-Hanna, J. C., Besserer, J., Head, J. W. III, Howett, C. J. A., Kiefer, W. S., Lucey, P. J., et al. (2014). Structure and evolution of the lunar Procellarum region as revealed by GRAIL gravity data. *Nature*, **514**, 68–71. https://doi.org/10.1038/nature13697

Arnold, D., Bertone, S., Jäggi, A., Beutler, G., & Mervart, L. (2015). GRAIL gravity field determination using the celestial mechanics approach. *Icarus*, **261**, 182–192. https://doi.org/10.1016/j.icarus.2015.08.015

Asmar, S. W., Konopliv, A. S., Watkins, M. M., Williams, J. G., Park, R. S., Kruizinga, G., et al. (2013). The scientific measurement system of the Gravity Recovery and Interior Laboratory (GRAIL) Mission. *Space Science Reviews*, **178**(1), 25–55. https://doi.org/10.1007/s11214-013-9962-0

Acknowledgments

The data used in this study are all
available in the Geosciences Node of
the NASA Planetary Data System. The
lunar gravity field models presented in
this study will be made available at the
PDS, as well as on our own data portal
at https://pgda.gsfc.nasa.gov, along
with auxiliary data such as the results
presented in several figures. This work
was supported by the GRAIL Project
and the NASA Discovery Program, and
we acknowledge the entire GRAIL
Science Team for many constructive
discussions about the gravity field
models. We sincerely thank the NASA
Center for Climate Simulation (NCCS)
at NASA Goddard Space Flight Center
for the use of their computational
resources and their support. We thank
Francis Nimmo (UC Santa Cruz) for
discussions about the density-depth
models, and for his formal review that
further improved this paper. We thank
the editor, David Baratoux, and an
anonymous reviewer: their comments
helped improve this paper. We
extensively used the freely available
SHTOOLS software (Wieczorek &
Meschede, 2018) for the localization.
All figures were made with the freely
available GMT software (Wessel
et al., 2013).
Besserer, J., Nimmo, F., Wieczorek, M. A., Weber, R. C., Kiefer, W. S., McGovern, P. J., et al. (2014). GRAIL gravity constraints on the vertical and lateral density structure of the lunar crust. *Geophysical Research Letters*, *41*, 5771–5777. https://doi.org/10.1002/2014GL060240

Bierman, G. J. (1977). *Factorization methods for discrete sequential estimation*. New York: Academic Press, Inc.

Deutsch, A. N., Neumann, G. A., Head, J. W., & Wilson, L. (2019). GRAIL-identified gravity anomalies in Oceanus Procellarum: Insight into subsurface impact and magmatic structures on the Moon. *Icarus*, *331*, 192–208. https://doi.org/10.1016/j.icarus.2019.05.027

Ding, M., Soderblom, J. M., Bierson, C. J., Nimmo, F., Milbury, C., & Zuber, M. T. (2018). Constraints on lunar crustal porosity from the gravitational signature of impact craters. *Journal of Geophysical Research: Planets*, *123*, 2281–2294. https://doi.org/10.1029/2018JE005654

Floberghagen, R. (2002). *Lunar gravimetry* (Vol. 273). Dordrecht, The Netherlands: Kluwer Academic Publishers.

Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. (2014). The planetary and lunar ephemerides DE430 and DE431 (IPN Progress Report 42-196): Jet Propulsion Laboratory.

Garrick-Bethell, I., & Zuber, M. (2009). Elliptical structure of the lunar South Pole-Aitken basin. *Icarus*, *204*, 399–408. https://doi.org/10.1016/j.icarus.2009.05.032

Genova, A., Goossens, S., Lemoine, F., Mazarico, E., Neumann, G., Smith, D., & Zuber, M. (2016). Seasonal and static gravity field of Mars from MGS, Mars Odyssey and MRO radio science. *Icarus*, *272*, 228–245. https://doi.org/10.1016/j.icarus.2016.02.050

Genova, A., Mazarico, E., Goossens, S., Lemoine, F. G., Neumann, G. A., Smith, D. E., & Zuber, M. T. (2018). Solar system expansion and strong equivalence principle as seen by the NASA MESSENGER mission. *Nature Communications*, *9*(1), 289. https://doi.org/10.1038/s41467-017-02558-1

Golub, G. H., & van Loan, C. F. (1989). *Matrix computations* (2nd ed.). Baltimore, MD, USA: The Johns Hopkins University Press.

Gong, S., Wieczorek, M. A., Nimmo, F., Kiefer, W. S., Head, J. W., Huang, C., et al. (2016). Thicknesses of mare basalts on the Moon from gravity and topography. *Journal of Geophysical Research: Planets*, *121*, 854–870. https://doi.org/10.1002/2016JE005008

Goossens, S., Lemoine, F. G., Sabaka, T. J., Nicholas, J. B., Mazarico, E., Rowlands, D. D., et al. (2016). A Global Degree and Order 1200 Model of the Lunar Gravity Field Using GRAIL Mission Data. In *47th Lunar and Planetary Science Conference, abstract 1484*, The Woodlands, Texas, USA.

Goossens, S., Matsumoto, K., Liu, Q., Kikuchi, F., Sato, K., Hanada, H., et al. (2011). Lunar gravity field determination using SELENE same-beam differential VLBI tracking data. *Journal of Geodesy*, *85*, 205–228. https://doi.org/10.1007/s00190-010-0430-2

Goossens, S., Sabaka, T., Fernández Mora, A., & Heijkoop, E. (2018). A High-Resolution Global Map of Lunar Gravity from Patched Local Solutions Using GRAIL Data. In *49th Lunar and Planetary Science Conference, abstract 2477*.

Goossens, S., Sabaka, T. J., Genova, A., Mazarico, E., Nicholas, J. B., & Neumann, G. A. (2017). Evidence for a low bulk crustal density for Mars from gravity and topography. *Geophysical Research Letters*, *44*, 7686–7694. https://doi.org/10.1002/2017GL074172

Goossens, S., Sabaka, T. J., Nicholas, J. B., Lemoine, F. G., Rowlands, D. D., Mazarico, E., et al. (2014). High-resolution local gravity model of the south pole of the Moon from GRAIL extended mission data. *Geophysical Research Letters*, *41*, 3367–3374. https://doi.org/10.1002/2014GL060178

Han, S. C. (2013). Determination and localized analysis of intersatellite line of sight gravity difference: Results from the GRAIL primary mission. *Journal of Geophysical Research: Planets*, *118*, 2323–2337. https://doi.org/10.1002/2013JE004402

Han, S. C., Schmerr, N., Neumann, G., & Holmes, S. (2014). Global characteristics of porosity and density stratification within the lunar crust from GRAIL gravity and LOLA topography data. *Geophysical Research Letters*, *41*, 1882–1889. https://doi.org/10.1002/2014GL059378

Head, J. W., Johnson, B. C., Keane, J. T., Kiefer, W. S., McGovern, P. J., Neumann, G. A., et al. (2018). Ring faults and ring dikes around the Orientale basin on the Moon. *Icarus*, *310*, 1–20. https://doi.org/10.1016/j.icarus.2017.12.012

Huang, Q., & Wieczorek, M. A. (2012). Density and porosity of the lunar crust from gravity and topography. *Journal of Geophysical Research*, *117*, E05003. https://doi.org/10.1029/2012JE004062

Jansen, J. C., Andrews-Hanna, J. C., Li, Y., Lucey, P. G., Taylor, G. J., Goossens, S., et al. (2017). Small-scale density variations in the lunar crust revealed by GRAIL. *Icarus*, *291*, 107–123. https://doi.org/10.1016/j.icarus.2017.03.017

Jansen, J. C., Andrews-Hanna, J. C., Milbury, C., Head, J. W. III, Li, Y., Melosh, H. J., & Zuber, M. T. (2019). Radial gravity anomalies associated with the ejecta of the Orientale basin. *Icarus*, *319*, 444–458. https://doi.org/10.1016/j.icarus.2018.09.034

Kaula, W. M. (1963). Determination of the Earth's gravitational field. *Reviews of Geophysics and Space Physics*, *1*, 507–551. https://doi.org/10.1029/RG001i004p00507

Kaula, W. M. (1966). Theory of satellite geodesy, applications of satellites to geodesy. Waltham, MA, USA: Blaisdell Publishing Company.

Klinger, B., Baur, O., & Mayer-Gürr, T. (2014). GRAIL gravity field recovery based on the short-arc integral equation technique: Simulation studies and first real data results. *Planetary and Space Science*, *91*, 83–90. https://doi.org/10.1016/j.pss.2013.12.001

Klipstein, W., Arnold, B., Enzer, D., Ruiz, A., Tien, J., Wang, R., & Dunn, C. (2013). The lunar Gravity Ranging System for the Gravity Recovery and Interior Laboratory (GRAIL) Mission. *Space Science Reviews*, *178*(1), 57–76. https://doi.org/10.1007/s11214-013-9973-x

Konopliv, A. S., Banerdt, W. B., & Sjogren, W. L. (1999). Venus gravity: 180th degree and order model. *Icarus*, *139*, 3–18. https://doi.org/10.1006/icar.1999.6086

Konopliv, A. S., Park, R. S., Yuan, D.-N., Asmar, S. W., Watkins, M. M., Williams, J. G., et al. (2013). The JPL lunar gravity field to spherical harmonic degree 660 from the GRAIL Primary Mission. *Journal of Geophysical Research: Planets*, *118*, 1415–1434. https://doi.org/10.1002/jgre.20097

Konopliv, A. S., Park, R. S., Yuan, D.-N., Asmar, S. W., Watkins, M. M., Williams, J. G., et al. (2014). High-resolution lunar gravity fields from the GRAIL Primary and Extended Missions. *Geophysical Research Letters*, *41*, 1452–1458. https://doi.org/10.1002/2013GL059066

Lemoine, F. G., Goossens, S., Sabaka, T. J., Nicholas, J. B., Mazarico, E., Rowlands, D. D., et al. (2013). High-degree gravity models from GRAIL primary mission data. *Journal of Geophysical Research: Planets*, *118*, 1676–1698. https://doi.org/10.1002/jgre.20118

Lemoine, F. G., Goossens, S., Sabaka, T. J., Nicholas, J. B., Mazarico, E., Rowlands, D. D., et al. (2014). GRGM900C: A degree-900 lunar gravity model from GRAIL primary and extended mission data. *Geophysical Research Letters*, *41*, 3382–3389. https://doi.org/10.1002/2014GL060027

Lemoine, F. G., Smith, D. E., Rowlands, D. D., Zuber, M. T., Neumann, G. A., Chinn, D. S., & Pavlis, D. E. (2001). An improved solution of the gravity field of Mars (GMM-2B) from Mars Global Surveyor. *Journal of Geophysical Research*, *106*(E10), 23,359–23,376. https://doi.org/10.1029/2000JE001426

Lemoine, F. G. R., Smith, D. E., Zuber, M. T., Neumann, G. A., & Rowlands, D. D. (1997). A 70th degree lunar gravity model (GLGM-2) from Clementine and other tracking data. *Journal of Geophysical Research*, *102*(E7), 16,339–16,359. https://doi.org/10.1029/97JE01418

Lerch, F. (1991). Optimum data weighting and error calibration for estimation of gravitational parameters. *Bulletin Géodesique*, *65*, 44–52. https://doi.org/10.1007/BF00806341
Luthcke, S. B., Rowlands, D. D., Lemoine, F. G., Klosko, S. M., Chinn, D., & McCarthy, J. J. (2006). Monthly spherical harmonic gravity field solutions determined from GRACE inter-satellite range-rate data. *Geophysical Research Letters*, **33**, L02402. https://doi.org/10.1029/2005GL024846

Matsumoto, K., Yamada, R., Kikuchi, F., Kamata, S., Ishihara, Y., Iwata, T., et al. (2015). Internal structure of the Moon inferred from Apollo seismic data and selenodetic data from GRAIL and LLR. *Geophysical Research Letters*, **42**, 7351–7358. https://doi.org/10.1002/2015GL065335

Matsuyama, I., Nimmo, F., Keane, J. T., Chan, N. H., Taylor, G., Wieczorek, M. A., et al. (2016). GRAIL, LLR, and LOLA constraints on the interior structure of the Moon. *Geophysical Research Letters*, **43**, 8365–8375. https://doi.org/10.1002/2016GL069952

Mazarico, E., Genova, A., Goossens, S., Lemoine, F. G., Neumann, G. A., Zuber, M. T., et al. (2014). The gravity field, orientation, and ephemeris of Mercury from MESSENGER observations after three years in orbit. *Journal of Geophysical Research: Planets*, **119**, 2417–2436. https://doi.org/10.1002/2014JE004675

Mazarico, E., Lemoine, F., Han, S. C., & Smith, D. (2010). GLM3-1, a degree 150 lunar gravity model from the historical tracking data of NASA Moon orbiters. *Journal of Geophysical Research*, **115**, E05001. https://doi.org/10.1029/2009JE003472

Milbury, C., Johnson, B. C., Melosh, H. J., Collins, G. S., Blair, D. M., Soderblom, J. M., et al. (2015). Preimpact porosity controls the gravity signature of lunar craters. *Geophysical Research Letters*, **42**, 9711–9716. https://doi.org/10.1002/2015GL066198

Montenbruck, O., & Gill, E. (2000). *Satellite Orbits*. Berlin: Springer-Verlag.

Neumann, G. A., Zuber, M. T., Wieczorek, M. A., Head, J. W., Baker, D. M. H., Solomon, S. C., et al. (2015). Lunar impact basins revealed by Gravity Recovery and Interior Laboratory measurements. *Science Advances*, **1**, 9. https://doi.org/10.1126/sciadv.1500852

Park, R. S., Konopliv, A. S., Yuan, D.-N., Asmar, S. W., Watkins, M. M., Williams, J. G., et al. (2015). A high-resolution spherical harmonic degree 1500 lunar gravity field from the GRAIL mission. In *American Geophysical Union Fall Meeting*, San Francisco, California, USA. Paper 67307, abstract G41B-01.

Park, R. S., Konopliv, A. S., Yuan, D.-N., Asmar, S. W., Watkins, M. M., Williams, J. G., et al. (2010). Global mass flux solutions from GRACE: A comparison of parameter estimation strategies - Mass concentrations versus Stokes coefficients. *Journal of Geophysical Research*, **115**, B01403. https://doi.org/10.1029/2009JB006546

Pavlis, D. E., & Nicholas, J. B. (2017). GEODYN II system description (Vols. 1-5) (*contractor report*). Greenbelt, MD: SGT Inc.

Sabaka, T. J., Nicholas, J. B., Goossens, S., Lemoine, F. G., & Mazarico, E. (2014). Error propagation for high-degree gravity models developed from the GRAIL mission. Retrieved from http://pds-geosciences.wustl.edu/grail/grail1-1-lgrs-5-rdr-v1/grail1001/extras/clones/clonedescription.pdf

Sabaka, T. J., Rowlands, D. D., Luthcke, S. B., & Boy, J. P. (2010). Improving global mass flux solutions from Gravity Recovery and Climate Experiment (GRACE) through forward modeling and continuous time correlation. *Journal of Geophysical Research*, **115**, B11403. https://doi.org/10.1029/2010JB007533

Seber, G. A. F., & Wild, C. J. (1989). *Nonlinear regression*. New Jersey, USA: Wiley-Interscience.

Smith, D. E., Zuber, M. T., Neumann, G. A., Mazarico, E., Lemoine, F. G., Head, J. W. III, et al. (2016). Summary of the results from the Lunar Orbiter Laser Altimeter after seven years in orbit. *Icarus*, **283**, 70–91. https://doi.org/10.1016/j.icarus.2016.06.006

Soderblom, J. M., Evans, A. J., Johnson, B. C., Melosh, H. J., Miljkovic, K., Phillips, R. J., et al. (2015). The fractured Moon: Production and saturation of porosity in the lunar highlands from impact cratering. *Geophysical Research Letters*, **42**, 6939–6944. https://doi.org/10.1002/2015GL065022

Tapley, B., Bettadpur, S., Ries, J., Thompson, P., & Watkins, M. (2004). GRACE measurements of mass variability in the Earth system. *Science*, **305**, 503–505. https://doi.org/10.1126/science.1099192

Tapley, B., Schutz, B., & Born, G. (2004). *Statistical orbit determination*. Burlington, MA, USA: Elsevier Academic Press.

Taylor, G. J., & Wieczorek, M. A. (2014). Lunar bulk chemical composition: A post-Gravity Recovery and Interior Laboratory reassessment. *Philosophical Transactions of the Royal Society of London Series A*, **372**, 20130242. https://doi.org/10.1098/rsta.2013.0242

Wessel, P., Smith, W. H. F., Scharroo, R., Luis, J., & Wobbe, F. (2013). Generic mapping tools: Improved version released. *EOS, Transactions American Geophysical Union*, **94**(45), 409–410. https://doi.org/10.1002/2013EO450001

Wieczorek, M. A. (2008). Constraints on the composition of the martian south polar cap from gravity and topography. *Icarus*, **196**, 506–517. https://doi.org/10.1016/j.icarus.2007.10.026

Wieczorek, M. (2012). GRAIL Crustal Thickness Archive [Data set]. Zenodo https://doi.org/10.5281/zenodo.997347

Wieczorek, M. A. (2015). Gravity and topography of the terrestrial planets, (2nd ed.), *Treatise on geophysics* (pp. 153–193). Oxford: Elsevier. https://doi.org/10.1016/B978-0-444-53802-4.00169-X

Wieczorek, M. A., & Meschede, M. (2018). SHTools: Tools for working with spherical harmonics. *Geochemistry, Geophysics, Geosystems*, **19**, 2574–2592. https://doi.org/10.1029/2018GC007529

Wieczorek, M. A., Neumann, G. A., Nimmo, F., Kiefer, W. S., Taylor, G. J., Melosh, H. J., et al. (2013). The crust of the Moon as seen by GRAIL. *Science*, **339**(6120), 671–675. https://doi.org/10.1126/science.1231530

Wieczorek, M. A., & Phillips, R. J. (1998). Potential anomalies on a sphere: Applications to the thickness of the lunar crust. *Journal of Geophysical Research*, **103**(E1), 1715–1724. https://doi.org/10.1029/97JE03136

Wieczorek, M. A., & Simons, F. J. (2005). Localized spectral analysis on the sphere. *Geophysical Journal International*, **162**, 655–675. https://doi.org/10.1111/j.1365-246X.2005.02687.x

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. https://doi.org/10.1007/s00041-006-6904-1

Williams, J. G., Konopliv, A. S., Boggs, D. H., Park, R. S., Yuan, D.-N., Lemoine, F. G., et al. (2014). Lunar interior properties from the GRAIL Mission. *Journal of Geophysical Research: Planets*, **119**, 1546–1578. https://doi.org/10.1002/2013JE004559

Wirnsberger, H., Krauss, S., & Mayer-Gürr, T. (2019). First independent Graz lunar gravity model derived from GRAIL. *Icarus*, **317**, 324–336. https://doi.org/10.1016/j.icarus.2018.08.011

Xu, P. (1992). The value of minimum norm estimation of geopotential fields. *Geophysical Journal International*, **111**, 170–178.

Zuber, M. T., Smith, D. E., Lehman, D. H., Hoffman, T. L., Asmar, S. W., & Watkins, M. M. (2013). Gravity Recovery and Interior Laboratory (GRAIL): Mapping the lunar interior from crust to core. *Space Science Reviews*, **178**(1), 3–24. https://doi.org/10.1007/s11214-012-9952-7

Zuber, M. T., Smith, D. E., Neumann, G. A., Goossens, S., Andrews-Hanna, J. C., et al. (2016) Gravity field of the Orientale basin from the Gravity Recovery and Interior Laboratory Mission. *Science*, **354**, 438–441. https://doi.org/10.1126/science.aag0519

Zuber, M. T., Smith, D. E., Watkins, M. M., Asmar, S. W., Konopliv, A. S., Lemoine, F. G., et al. (2013). Gravity Field of the Moon from the Gravity Recovery and Interior Laboratory (GRAIL) Mission. *Science*, **339**, 668–671. https://doi.org/10.1126/science.1231507