1. Introduction
===============

The twin spacecraft of the NASA Discovery Gravity Recovery and Interior Laboratory (GRAIL) mission mapped the gravity field of the Moon in two phases: during the Primary Mission (1 March to 29 May 2012) from a mean altitude of 55 km and during the Extended Mission (30 August to 14 December 2012) from a mean altitude of 23 km \[*Zuber et al.*, [@b28], [@b29]\]. The main data produced by GRAIL include the Ka-Band Range-Rate (KBRR) data between the two satellites, acquired by the GRAIL Lunar Gravity Ranging System and the Doppler data from the antennae of the Deep Space Network (DSN) that tracked the spacecraft in lunar orbit from the Earth \[*Asmar et al.*, [@b1]; *Klipstein et al.*, [@b9]\].

Initial analysis of the GRAIL Primary Mission data led to spherical harmonic models complete to degree and order 420 \[*Zuber et al.*, [@b30]\] that improved on the pre-GRAIL gravity models derived from Lunar Prospector, Kaguya, and other historical data \[*Konopliv et al.*, [@b11]; *Matsumoto et al.*, [@b18]\] by a factor of 3 to 4 in spatial resolution, and 3 to more than 5 orders of magnitude in quality. Subsequently, more detailed analyses of the Primary Mission data by the GRAIL gravity teams at the Jet Propulsion Laboratory (JPL), and the NASA Goddard Space Flight Center (NASA GSFC), produced models to degree 660 in spherical harmonics, GL0660B \[*Konopliv et al.*, [@b12]\] and GRGM660PRIM \[*Lemoine et al.*, [@b17]\]. This paper describes GRGM900C, a lunar gravity solution to degree and order 900 in spherical harmonics (resolving spatial blocks of 6 km) obtained from the analysis of the Primary and Extended Mission data by the NASA GSFC GRAIL gravity team. The JPL GRAIL gravity team solution to degree 900 is described by *Konopliv et al.* \[[@b13]\]. The degree 900 models developed by the two gravity teams differ in terms of the orbit determination software used, a priori models, data editing, and parameter estimation strategies.

2. Extended Mission Summary
===========================

During GRAIL\'s Extended Mission, the GRAIL spacecraft orbits were maintained within a tight band of apoapse and periapse altitudes. From 30 August to 17 November 2012 the minimum periapse altitude ranged from 10 to 20 km and the apoapse altitude from 25 to 35 km. From 17 November 2012 through the end of the mission, the periapse altitude generally ranged from 5 to 10 km while the apoapse altitude ranged from 15 to 25 km. The low altitude of the GRAIL spacecraft made the mission operations more complex, requiring sometimes as many as three maneuvers per week to maintain the desired orbit \[*Sweetser et al.*, [@b24]; *Ryne et al.*, [@b22]\]. The end-of-mission operations from 4 December 2012 onward included orbits where the GRAIL spacecraft were targeted to fly within 5--10 km over part of Orientale \[*Zuber et al.*, [@b29]\]. The orbit maneuvers also maintained the spacecraft separation between 40 and 60 km in order to reduce KBRR signal multipath, whereas during the Primary Mission the spacecraft separation ranged from 70 to 200 km. The Extended Mission included effectively 3 1/2 mapping cycles, compared to three for the Primary Mission. A mapping cycle is a complete rotation of the Moon with respect to the ascending node of the GRAIL spacecraft orbits. The GRAIL Extended Mission ground tracks were interleaved to improve the longitudinal sampling of the lunar gravity field. However, the lunar gravity field was not sampled uniformly in terms of the periapse altitude. Whereas a swath of longitudes just east of Orientale was sampled at 5--10 km periapse height, both the South Pole Aitken Basin (due partially to its depth) and a swath of longitudes in the central nearside (∼335°E to ∼5°E) are sampled at about 20 km of periapse altitude. (See Figure S1 in the supporting information.)

3. Data, Modeling, and Method of Solution
=========================================

3.1. Data Overview
------------------

The Deep Space Network (DSN) tracking data (two-way tracking at S band and one-way tracking at X band) and the intersatellite Ka-Band Range-Rate (KBRR) data obtained during the Extended Mission are similar to those obtained in the Primary Mission \[*Asmar et al.*, [@b1]; *Lemoine et al.*, [@b17]\]. We used the Level 1B tracking data, which are a derived product based on the dual one-way Ka-Band ranging between the two satellites \[*Klipstein et al.*, [@b9]\]. The GRAIL Level 1B data are time-tagged in Barycentric Dynamic Time using a procedure outlined by *Kruizinga et al.* \[[@b14]\]. The time-tagging of the data is dependent on the quality of the reconstructed orbits, and so the Level 1B team produced several versions of the intersatellite tracking data, as improved lunar gravity models were developed \[*Kruizinga et al.*, [@b15]\]. The GRGM900C gravity solution discussed in this paper was based on the processing of the version 3 Level 1B KBRR data. The Primary Mission Level 1B KBRR data had a sampling of 5 s, whereas the Extended Mission data had a sampling of 2 s. As a consequence, the noise floor of the Primary Mission KBRR data is thought to be 0.03 μm/s, but for the Extended Mission the noise floor is thought to be 0.05--0.07 μm/s (N. Harvey, Jet Propulsion Laboratory, unpublished data, 2013).

3.2. Force Modeling
-------------------

The gravity field potential due to the Moon is modeled as a spherical harmonics series \[*Kaula*, [@b8]; *Heiskanen and Moritz*, [@b7]\],

where *G* is the Universal constant of gravitation, *M* is the mass of the Moon, are the normalized associated Legendre polynomials, *R*~*e*~ is the reference radius (1738 km), *ϕ*, *θ* and *r* are the longitude, latitude, and radius at the evaluation point. and are the normalized Stokes coefficients. Our solution, GRGM900C, estimates the lunar *GM* and the normalized Stokes coefficients to degree and order 900. In addition, we also estimated the degree 2 Love number *k*~2~. The model, GRGM900C, is based on the DE421 lunar and planetary ephemerides \[*Williams et al.*, [@b27]\].

Other aspects of the force modeling follow the strategy used for the Primary Mission \[*Lemoine et al.*, [@b17]\]. The modeling of the solar radiation pressure for GRAIL requires the detailed modeling of the Sun\'s shadow function, since the shape of the topography can alter this function over that predicted from a spherical or ellipsoidal Moon, especially for near-full-Sun orbits where grazing transitions are long \[*Lemoine et al.*, [@b17]\]. We derived the shadow function for the Extended Mission, taking into account the GRAIL spacecraft trajectories and the lunar topography determined by the altimeter data from the Lunar Reconnaissance Orbiter, Lunar Orbiter Laser Altimeter (LOLA) \[*Smith et al.*, [@b23]\], archived in the Planetary Data System (LOLA/Lunar Reconnaissance Orbiter Radio Science Planetary Data System Data Node, <http://imbrium.mit.edu>). The planetary radiation pressure, including the reflected solar radiation due to the lunar albedo, and the radiation from the lunar thermal emission are modeled as described in *Lemoine et al.* \[[@b17]\]. The relative magnitudes of these radiation pressure accelerations on GRAIL are discussed by *Park et al.* \[[@b19]\].

3.3. Solution Development
-------------------------

We used the NASA GSFC GEODYN II Orbit Determination and Geodetic Parameter Estimation Program \[*Pavlis et al.*, [@b20]\] to analyze the GRAIL tracking data and to compute the partial derivatives of the adjusted parameters. Because GEODYN uses a fixed-step Cowell integrator, the only change introduced with respect to *Lemoine et al.* \[[@b17]\] was to reduce the integration step size from 1.0 to 0.5 s, in order to provide an adequate sampling of the gravitational potential field of the Moon for the degree and order 900 fields, given the mean spacecraft orbital velocity of ∼1.7 km/s.

The solutions were derived using the QR method, which allows for a numerically stable way of solving the normal equations in terms of a Cholesky square root matrix \[*Golub and van Loan*, [@b4]\]. The degree 900 solutions included estimates for 811,808 common parameters, and used up to 4080 processors and 6 terabytes of memory belonging to the supercomputers of the NASA Center for Climate Simulation (NCCS) at the NASA Goddard Space Flight Center. The arc lengths were 2.5 days for the data of the Primary Mission, where the arcs were bounded by spacecraft angular momentum desaturation events. For the Extended Mission data, we used 1 day arcs to limit the buildup of dynamical errors. Initially, these dynamical errors were mostly due to gravity mismodeling, especially in the first interim gravity solutions we developed with the Extended Mission data. We note that mismodeling of radiation pressure (both the direct solar radiation pressure, and the planetary radiation pressure) will also be accommodated by the estimation of empirical accelerations. The set of estimated arc parameters included the satellite initial state, a KBRR bias, and a KBRR time-tag bias per arc, and empirical accelerations. It is necessary to estimate a time-tag bias for the KBRR data to accommodate any residual time offsets between these data and the time tags of the DSN data \[*Kruizinga et al.*, [@b15]\]. The KBRR time-tag bias is stable over a data arc and typically of the order of a few milliseconds. We note that as updated versions of the KBRR Level1B data were made available, the values and the dispersion in these estimated KBRR time-tag biases were reduced \[e.g., *Fahnestock et al.*, [@b2]\]. We estimated empirical accelerations along track and cross track to the orbit, every quarter revolution, for each satellite. Acceleration parameters in time periods *t*~*i*~ and *t*~*j*~ are tied together with an exponentially decaying weight factor that has a time correlation of one-quarter revolution \[*Lemoine et al.*, [@b17], cf. equation (5)\]. An a priori sigma, *σ*~*A*~ of 1.0 ×10^−9^, was assigned to control the amplitude of the estimated accelerations. The number of observations and data weights used in the GSFC degree 900 solutions are summarized in Table S1 of the supporting information.

Using the Primary and the Extended Mission data, we developed successively several interim models of varying size (*e.g.,*660 × 660, 720 × 720, and 900 × 900). Each time we obtained a new global solution, we reconverged the orbits and recomputed the partial derivatives before obtaining the subsequent model. The processing that resulted in the first degree and order 900 model, called GRGM900A, was based on a starting model of degree and order 720. Then GRGM900A was used to produce the set of partial derivatives that were used to develop the degree 900 models that are the subject of this paper, GRGM900B and GRGM900C. For all the solutions we applied a constraint on the total coefficient power per spherical harmonic degree, known as a Kaula power law constraint. For the degree 900 solutions we applied the Kaula constraint only above degree 600, after verifying the behavior of the unconstrained solution. We applied a power law constraint of 2.5 ×10^−4^/*ℓ*^2^ for GRGM900B, and a constraint of 3.6 ×10^−4^/*ℓ*^2^ for GRGM900C. Also, for GRGM900C we downweighted the data after 18 November 2012 (when the spacecraft had lower periapse altitudes), by applying a factor of 0.316 to the end-of-mission Cholesky square root matrix, which meant the data were effectively downweighted (both DSN and KBRR) by a factor of 10, from the nominal weights of 0.12 mm/s (for the DSN data) and 0.1 μm/s (for the KBRR data).

In addition to deriving the models, we also calibrate the a posteriori parameter error-covariance matrix by estimating a scalar multiplicative factor, *s*^2^, which corresponds to treating all data and constraints as a single statistical unit via the formula

where is the a posteriori weighted residual variance, *N* is the number of observations plus the Kaula constraint, and *M* is the number of parameters. The effect of calibration is to bring the formal errors into agreement with those actually observed \[cf. *Kusche*, [@b16]; *Lemoine et al.*, [@b17], sections 4.2.4 and section 5.4\].

4. Results
==========

4.1. Power Spectrum for GRAIL Gravity Solutions
-----------------------------------------------

We compute the root-mean-square (RMS) power of the coefficients, RMS~*ℓ*~(U), and RMS power of the coefficient standard deviations following *Kaula* \[[@b8]\]:

The RMS power for the GRAIL gravity solutions GRGM660PRIM and GRGM900C are shown in [Figure 1](#fig01){ref-type="fig"}. We show the coefficient power spectra as well as the calibrated spectrums of the solution standard deviations for GRGM900C and GRGM660PRIM. In addition, we show the RMS power of the coefficients for the GRGM900C unconstrained solution, where no Kaula constraint was applied. In Figure S2, we show the coefficient differences between the constrained and unconstrained solutions to degree 900, normalized by GRGM900C. Since the unconstrained solution begins to deviate from a smooth power law around degrees 630--650, we chose *ℓ* = 600 to be the degree above which we applied a Kaula constraint in GRGM900C. Note also that the power spectrum of GRGM660PRIM starts to drop more steeply after degree 400, indicating that there is not enough resolving power in the Primary Mission data for a global model of degree and order 660. The impact of the Extended Mission data can be seen in the improved error estimates of the midrange degrees (*ℓ*=200 to 500) by 2--3 orders over GRGM660PRIM.

![Power spectra of the Primary Mission GRAIL gravity solution, GRGM660PRIM \[*Lemoine et al.*, [@b17]\], GRGM900C, the new solution including Extended Mission data, and the unconstrained version of GRGM900C, where no Kaula constraint was applied. The power spectra of the coefficient standard deviations for GRGM660PRIM and GRGM900C are also shown.](grl0041-3382-f1){#fig01}

The calibrated errors of GRGM900C cross the coefficient power spectrum at *ℓ* = 790 ([Figure 1](#fig01){ref-type="fig"}), whereas the uncalibrated errors do not intersect the coefficient power spectrum. Referring to equation [2](#m2){ref-type="disp-formula"}, the calibration factors are 3.20 for GRGM900B and 1.93 for GRGM900C and directly scale the uncalibrated formal errors from the two solutions. They show that the downweighting of the end-of-mission data has reduced the systematic errors in the solution. This indicates the importance of data weighting and solution calibration on the estimation of the error spectrum.

The low-degree errors for GRGM900C are close to and even slightly above those of GRGM660PRIM. In addition, we report that the lunar gravitational parameter *GM* for GRGM900C is 4.90279997×10^12^±14147 m^3^/s^2^, and the degree 2 potential Love number *k*~2~ is 0.024116 ± 0.000108. For both these values the errors are actually higher than those reported for GRGM660PRIM in \[*Lemoine et al.*, [@b17]\]. This is mostly due to the 1 day arcs which we used in the processing of the Extended Mission data, which slightly degrade the sensitivity with respect to these parameters. Future models will extend the length of these arcs and increase the sensitivity with respect to the lower degrees.

4.2. Correlations With Topography
---------------------------------

The gravity solutions from the GRAIL Primary Mission showed high correlations of gravity and topography. We also evaluate this correlation for the new Extended Mission models. We compute the gravity induced by topography following *Wieczorek and Phillips* [@b25], assuming a uniform density of 2560 kg/m^3^\[*Wieczorek et al.*, [@b26]\], using the topography expansion to the ninth power. We see in [Figure 2](#fig02){ref-type="fig"} that the addition of the Extended Mission data has extended the range of degrees with high correlations. For GRGM660PRIM, the correlations begin to decrease below 0.98 after *ℓ*= 320, whereas for the two Extended Mission models, correlations ≥ 0.98 extend to *ℓ*≈ 638. The two Extended Mission models, GRGM900B and GRGM900C, differ in their treatment of the end-of-extended-mission data and the magnitude of the Kaula constraint. By relaxing the Kaula constraint and slightly downweighting the Extended Mission data after 18 November 2012, we modestly improve the global correlations, while maintaining the same levels of postfit residuals (see section 4.4). For this reason, our preferred solution is GRGM900C.

![Correlations between gravity solutions and gravity computed from LOLA topography, and Bouguer power spectrum. The top curves (referred to the right axis) show global correlations for the Primary Mission model, GRGM660PRIM, and the models including Extended Mission data, GRGM900B and GRGM900C. The bottom curves (referred to the left axis) show the spectrum of the Bouguer power for the three gravity models.](grl0041-3382-f2){#fig02}

4.3. Bouguer Gravity
--------------------

The Bouguer gravity, which highlights subsurface structures, is computed by subtracting the gravity predicted from topography, as just described, from the estimated geopotential coefficients. The RMS power is computed according to equation [3](#m3){ref-type="disp-formula"}.

The Bouguer power (see bottom curves of [Figure 2](#fig02){ref-type="fig"}) begins to increase around *ℓ*=300 for GRGM660PRIM, and around *ℓ*≈570 for GRGM900B and GRGM900C. These degrees represent, in global terms, where spurious power in the geopotential coefficients begins to affect the global geophysical interpretation based on the gravity solutions.

We also compute the Bouguer coefficient amplitude . This is illustrated in [Figure 3](#fig03){ref-type="fig"} for two gravity solutions, GRGM900C without application of a Kaula constraint and the final solution GRGM900C. The solution derived without Kaula illustrates which degrees are most reliably determined by the KBRR data. As pointed out by *Floberghagen* \[[@b3]\], since the Moon is a slow rotator, resonances occur only at high degrees (*ℓ*≈ 350--370, and *ℓ*≈700--740). In the unconstrained solution, it is these resonance and near-resonance orders that develop excess power in the Bouguer coefficients (see [Figure 3](#fig03){ref-type="fig"}a). In these GRAIL solutions, the lower-order coefficients (*m* less than 300) are better determined than the "wings" of the gravity model (high degrees and high orders). This phenomenon is also a characteristic of solutions for the Earth gravity derived from GRACE data \[*Rowlands et al.*, [@b21]\]. The Kaula constraint successfully reduces some of the spurious power in the coefficients ([Figure 3](#fig03){ref-type="fig"}b), however, a signature remains in the RMS Bouguer coefficient power spectrum (bottom curves of [Figure 2](#fig02){ref-type="fig"}).

![Amplitude of Bouguer gravity coefficients (a) based on GRGM900C-unconstrained (no Kaula constraint applied), and (b) based on GRGM900C, with a Kaula constraint of 3.6 × 10^−4^/*ℓ*^2^ applied for degree *ℓ* greater than 600.](grl0041-3382-f3){#fig03}

4.4. Data Fit
-------------

In addition to examining the correlations with topography and the power spectra of the solutions, we also computed the a posteriori orbital fits to the tracking data. The fits to GRGM900C are very similar to those of GRGM900B, even though we note an improvement in the topography correlations with GRGM900C (see section 4.2). The GRGM900C fits to the Primary Mission KBRR data are largely unchanged from those obtained with the model GRGM600PRIM and are in the range 0.04--0.06 μm/s. The GRGM900C a posteriori fits to the Extended Mission KBRR data through 18 November 2012 are in the range of 0.10--0.17 μm/s, whereas the fits to the data beyond 18 November 2012 range from 0.2 to 3.2 μm/s. We illustrate the Extended Mission KBRR data fits for GRGM900C in [Figure 4](#fig04){ref-type="fig"}, together with the periapse height for each data arc. The RMS fits to the DSN S band data (not shown) are in general close to the data weights of 0.12 mm/s, except for the latter part of the Extended Mission where they reach RMS fits of 0.5 mm/s. The low-periapse altitude of 5--10 km dramatically increases the gravity signal inherent in the satellite-tracking data. In certain regions and along some orbital tracks, the inherent resolution is more than can be accommodated in a global spherical harmonic model of degree and order 900.

![RMS residuals of the KBRR data per data arc in the Extended Mission with the gravity model GRGM900C, shown together with the periapse altitude for each orbital arc.](grl0041-3382-f4){#fig04}

4.5. Anomaly Errors and Degree Strength
---------------------------------------

The propagation of the gravity field errors into a map of anomaly errors is helpful for evaluating the model\'s performance geographically and aids the geophysical interpretation using that model. However, the propagation of a full covariance matrix of a degree and order 900 field is computationally prohibitive. Instead, we have applied a Monte Carlo approach which uses the Cholesky square root matrix from QR to compute an ensemble of statistically consistent members, or "clones," to replace the error-covariance matrix in these calculations. We have verified this procedure by comparing its gravity anomaly errors for the SGM150J model \[*Goossens et al.*, [@b5]\] with those obtained from using the full error-covariance matrix and have found the convergence rate to be the expected , where *K* is the ensemble size. An example of the use of Monte Carlo calculations to propagate error functions for large covariances can be found in *Gundlich et al.* \[[@b6]\]. In Figure S3 in the supporting information, we show the anomaly errors to *ℓ*=720 for GRGM900C. The total global RMS error to degree and order 720 is 8.49 mGal. This compares to an RMS global anomaly error from GRGM660PRIM to degree and order 660 of 30.49 mGal \[*Lemoine et al.*, [@b17]\].

We also seek to evaluate the spatial variability in the model resolution. We do this by computing the degree strength of the solution following *Konopliv et al.* \[[@b10]\]. First, we evaluate the anomaly error at each degree for each geographic location using the Monte Carlo approach just described, and then we determine at what degree the anomaly error matches the Kaula law signal strength. The result is a map that illustrates the resolution of GRGM900C (see [Figure 5](#fig05){ref-type="fig"}). We see that the lowest degree strength is 550--600 over a broad swath of the central farside. It is possible that the slight deweighting of the post 18 November 2012 data affects the degree strength over the farside. Over Orientale, the degree strength of the model is 700--725; however, we note that the KBRR data obtained from the lowest altitudes at the end of the mission over Orientale have been deweighted in this solution. The model has the highest resolution (degree strength of 900) over the polar regions due to the extra sampling that comes from the convergence of the orbital ground tracks.

![Degree strength of the GRGM900C gravity model. The map shows spatially at what degree the gravity field errors match the expected signal power from the applied Kaula power law. The map is a Mollweide projection centered on 270°E, with the lunar farside to the left of the figure center, and the lunar nearside to the right of the figure center.](grl0041-3382-f5){#fig05}

5. Summary
==========

Using the method of QR factorization, we have derived a gravity field of the Moon from the GRAIL primary and Extended Mission data complete to degree and order 900 in spherical harmonics. The tracking data, in particular the GRAIL KBRR data, completely resolve the coefficients below *ℓ* = 600 (9 km block size). A power law constraint of 3.6 ×10^−4^/*ℓ*^2^ was applied for *ℓ*≥600. The inclusion of the Extended Mission GRAIL data improves the lunar gravity field resolution by 2--3 orders of magnitude over the midrange degree bands, compared to the Primary Mission GRAIL model, GRGM660PRIM. The variable sampling of the Moon at low altitude (see Figure S1 in the supporting information) combined with the RMS of fit to the KBRR data that is higher than the intrinsic noise, especially for the data from the last month of the Extended Mission, suggest that a degree 900 model does not capture all the signal inherent in the tracking data. Future work should expand the solutions to higher degree to obtain improved resolution or should assess the utility of local solutions over bounded regions to fully extract the signal inherent in the KBRR tracking data.

The work described herein was supported by the GRAIL Project and the NASA Discovery Program. We acknowledge Gerhard Kruizinga and his team at the Jet Propulsion Laboratory (JPL) for their meticulous work in developing the GRAIL Level 1B data sets. We also acknowledge the NASA Center for Climate Simulation (NCCS) at the NASA Goddard Space Flight Center for their assistance in optimizing the performance of the supercomputers to obtain the GRAIL degree 900 solutions. The GRAIL Level 1B data used in this study are available from the Geosciences Node of the Planetary Data System <http://pds-geosciences.wustl.edu/missions/grail/default.htm>.

The Editor thanks Rune Floberghagen and an anonymous reviewer for their assistance evaluating this paper.

Key Points
==========

1.  New gravity model of the Moon obtained to degree and order 900

2.  The model is obtained from analysis of GRAIL intersatellite data

3.  The degree strength of the model ranges from 575--900

Supporting Information
======================

[^1]: **Citation**: Lemoine, F. G., et al. (2014), GRGM900C: A degree 900 lunar gravity model from GRAIL primary and extended mission data
