High-resolution local gravity model of the south pole of the Moon from GRAIL extended mission data.

We estimated a high-resolution local gravity field model over the south pole of the Moon using data from the Gravity Recovery and Interior Laboratory's extended mission. Our solution consists of adjustments with respect to a global model expressed in spherical harmonics. The adjustments are expressed as gridded gravity anomalies with a resolution of 1/6° by 1/6° (equivalent to that of a degree and order 1080 model in spherical harmonics), covering a cap over the south pole with a radius of 40°. The gravity anomalies have been estimated from a short-arc analysis using only Ka-band range-rate (KBRR) data over the area of interest. We apply a neighbor-smoothing constraint to our solution. Our local model removes striping present in the global model; it reduces the misfit to the KBRR data and improves correlations with topography to higher degrees than current global models. KEY POINTS: We present a high-resolution gravity model of the south pole of the Moon Improved correlations with topography to higher degrees than global models Improved fits to the data and reduced striping that is present in global models.

1. Introduction

The Gravity Recovery and Interior Laboratory (GRAIL) mission was designed to map the structure of the lunar interior through high-precision global gravity mapping [Zuber et al., 2013a]. The mission was composed of two spacecraft with Ka-band intersatellite tracking as the single science instrument, complemented by tracking from Earth using the Deep Space Network (DSN) [Asmar et al., 2013]. The mission consisted of two phases: a primary mission which lasted from 1 March 2012 to 29 May 2012, during which the spacecraft had a mean altitude of 55 km above lunar surface, and an extended mission which lasted from 30 August 2012 to 14 December 2012. In the extended mission, the spacecraft had a mean altitude of 23 km above lunar surface until 18 November. After that, they were lowered to 20 km and finally to 11 km on 6 December. Initial analysis of the primary mission data resulted in a gravity field model in spherical harmonics of degree and order 420 [Zuber et al., 2013b], showing high correlations with topography up to high degrees and improving on the pre-GRAIL lunar gravity models such as those by Konopliv et al. [2001], Matsumoto et al. [2010], and Goossens et al. [2011], which were based on Lunar Prospector, Selenological and Engineering Explorer (SELENE), and other historical data. Subsequent analysis of the primary mission data resulted in models of degree and order 660 [Konopliv et al., 2013; Lemoine et al., 2013]. Owing to the much lower average altitudes above lunar surface in the extended mission, recent analysis has resulted in further improvements, in the form of models of degree and order 900 as presented by Konopliv et al. [2014] and Lemoine et al. [2014]. Both models have variable spatial resolution between degrees 600–680 over the central farside of the Moon to degree 900 at the poles. This variable spatial resolution of the gravity field is due to spacecraft altitude and ground track spacing [Konopliv et al., 2014]. Aliasing and relatively higher fits for the low-altitude part of the extended mission indicate that there is remaining signal in the Ka-band range-rate (KBRR) data even beyond degree 900. However, global spherical harmonics are not optimal when the spatial data coverage is heterogeneous, which is the case here at short wavelengths.

In instances where global models require a large number of parameters to describe high-resolution features, the estimation of these parameters may become unstable, and local representations can be of great utility: they have the potential to extract high-resolution information from the data, while not demanding the computational burden of determining a full global model of similar resolution. In satellite geodesy there is a long history of local representations that have been applied successfully to different data types from satellites orbiting various planetary bodies. Examples include analysis for Venus [Barriot and Balmino, 1992], Mars [Beuthe et al., 2006], and Jupiter's moon Ganymede [Anderson et al., 2004; Palguta et al., 2005], without even mentioning the many applications to Earth-orbiting satellites. Prior to the SELENE and GRAIL missions the tracking data distribution for the Moon was extremely asymmetrical with no farside data available, and thus, the Moon was especially suitable for the application of local methods, mostly using Lunar Prospector data: Sugano and Heki [,] and Goossens et al. [,] estimated anomalies while Han [2008] and Han et al. [2011] used localized spherical harmonic functions to model the nearside gravity field. These same functions were used by Goossens et al. [2012] to estimate the gravity field over the South Pole-Aitken basin on the farside of the Moon using SELENE data. In addition, local representations allow one to use constraints other than the spectral constraints often used in global gravity modeling [e.g., Rowlands et al., 2010; Sabaka et al., 2010].

The goal of the research presented here is to extract more information from the tracking data than can be modeled reasonably using global spherical harmonics. In the case of GRAIL KBRR data, our interest is especially in areas where the altitude of the spacecraft was low and/or where the coverage was dense. Both Konopliv et al. [2014] and Lemoine et al. [2014] show that the degree strength of their models is at the maximum degree at the poles because of the confluence of tracks. Therefore, we estimated a local gravity field model for the south pole of the Moon covering a cap of 40° radius, as adjustments to a global background model, at a resolution of 1/6° by 1/6° in both longitude and latitude, which is equivalent to the resolution of a degree and order 1080 model in spherical harmonics (spatial block size = 5 km).

2. Data and Methods

GRAIL data consist of both DSN tracking (both in two-way and one-way modes) and intersatellite KBRR. While DSN data mostly determine the absolute positioning of the satellites, it is the contribution of KBRR data acquired at low altitudes that has resulted in the unprecedented GRAIL model sizes. KBRR data in the extended mission have a sampling of 2 s. Because we aim to estimate local parameters in one area of the Moon, we choose to use only KBRR data over the area of interest, excluding DSN tracking altogether. It is well known that each KBRR observation from the Gravity Recovery and Climate Experiment (GRACE) and GRAIL missions is approximately directly proportional to the difference in gravity potential at the two satellite locations. This makes KBRR data ideal for estimating local parameters from a subset of the entire mission data set [Rowlands et al., 2005]. KBRR data outside the region do not contribute, and data inside the region do not contain significant information from outside the region.

Our analysis of the GRAIL extended mission KBRR data is based on the technique of Rowlands et al. [2002]. They showed that a short-arc analysis using KBRR data only is capable of recovering the underlying gravity field model. They used a baseline representation of the 12 parameters describing the difference vector between the two satellites and the state vector of the midpoint of the two satellites, instead of using the Cartesian inertial state vectors of the two satellites, all at the initial epoch. They showed that KBRR data are only sensitive to a subset of these baseline parameters, facilitating the coestimation of gravity parameters and orbit parameters from short arcs. This approach has been used successfully in the determination of monthly gravity solutions using GRACE data, for example Rowlands et al. [2005]; Luthcke et al. [2006], and later publications.

Since we adjust gravity parameters on a spherical cap at the south pole, and since the GRAIL satellites were in polar orbits, each orbit passes over the area and each arc covering the area is approximately 25 min long. Using the orbits from a global gravity field model (GRGM900A, a predecessor to the GRGM900C model discussed in Lemoine et al. [2014]), we create short arcs over the south pole area with initial conditions and other parameters from these orbits and then adjust only the baseline parameters to the KBRR data. We emphasize that in this baseline analysis, we do not use the DSN data. The baseline parameters we adjust are those given in Rowlands et al. [2002] as being most sensitive to the KBRR data: the pitch of the GRAIL A-B position baseline, the magnitude of the GRAIL A-B velocity baseline, and the pitch of the GRAIL A-B velocity baseline. We use the NASA Goddard Space Flight Center's (GSFC) GEODYN II Orbit Determination and Geodetic Parameter Estimation Program [Pavlis et al., 2013] for the short-arc analysis. We adjust only the three baseline parameters for the short arcs, and once these arcs are converged we use GEODYN again to generate the partial derivatives of the KBRR data points with respect to the adjusted baseline parameters and the chosen gravity parameters.

We represent local gravity as gravity anomalies on a regular grid of 1/6° by 1/6°. This resolution is equivalent to that of a degree and order 1080 expansion in spherical harmonics, and it was chosen in order to assess whether or not we could increase the resolution of our current global models (which are degree and order 900). For a cap of 40° radius (with a 1/6° of a degree cap on the south pole itself), this leads to the estimation of a total of 518,401 anomalies over the selected area, for which we use 992,080 KBRR data points from the extended mission. This means that the ratio of observations and parameters is less than two and that constraints will be important for the solution. We discuss the constraints below, and we discuss their influence on the solutions in section 3. Despite the fact that the longitudinal size of the grid cells physically decreases toward the pole, we did not compensate for this by decreasing the number of grid cells. Results in section 3 will show we do not see artifacts in this dense south pole grid.

Our gravity anomalies Δg are related to the disturbing potential T (from W = U + T, where W is the gravity potential and U the normal potential of a reference field) in the following way: in a spherical approximation, the anomaly can be expressed as a function of T through Δg =− ∂T/∂r − 2T/r, with r the spherical coordinate for radial distance. The potential T at a point P can be expressed in terms of anomalies Δg(Q) at points Q on an exterior boundary σas where R is the reference radius of the surface (1738 km in our case) and the kernel S(P,Q) is the Stokes-Pizetti kernel [Heiskanen and Moritz, 1984]. If the disturbing potential is expressed in spherical harmonics, and the anomalies are global, then both anomalies and standard spherical harmonics are completely equivalent, and anomalies can be expressed in spherical harmonics, with the anomalies carrying a factor (l − 1)in their expansion, where l is the degree.

The partial derivatives of KBRR data with respect to the anomalies can be constructed from the acceleration on the satellite, by differentiation of equation 1. When computing the acceleration contribution per anomaly (from which the partial derivative is also computed), we approximate the integral in equation 1 by a summation. For the integration over each separate anomaly in equation 1, each block (of size 1/6° by 1/6°) is divided into four subblocks (1/12° by 1/12°; except for 36 subblocks in the shape of wedges for the polar cap), and the integration is then approximated by the sum of the kernel function S(P,Q)evaluated at each center point of the subblock multiplied by the area of the subblock. This assumes that there are no anomaly contributions from outside of the chosen grid, which, in general, leads to boundary effects on the local solutions. This is further discussed in section 3.

Once the partial derivatives of the data with respect to all parameters are generated, the normal equation system is formed. Due to the size of the problem, we use the supercomputers of the NASA Center for Climate Simulation (NCCS) at NASA GSFC. We apply a neighbor-smoothing constraint to our solutions following Rowlands et al. [2010] and Sabaka et al. [2010]. These spatial constraints are straightforward: each anomaly pair is forced to be equal, but each constraint equation is only given a weight W that decreases with increasing distance between the ith and jth anomalies, according to W= exp(1 − d/D) with d being the distance between the two anomalies (computed from their center points), and D being the correlation distance, which was here set equal to the latitudinal size of each grid cell (which is thus the same for all anomalies because the grid is regular). We discuss these constraints further in section 3. We apply the smoothing constraint to the full anomaly field, i.e., on Δgglobal+Δgadj where Δgglobal is the vector with anomalies corresponding to the global spherical harmonics starting model and Δgadj is the vector with gravity anomalies for the adjustment. We apply the constraint to the full field (not just to the estimated gravity anomaly parameters) because the anomalies from GRGM900A show stripes (cf. Figure S1 in the supporting information). Once the constraint system is formed, it is added to the data system with a scaling factor μ, resulting in the following equation: where Ais the matrix with the partial derivatives of the data with respect to the parameters, W is the data weight matrix, Dis the matrix with partial derivatives of the constraints, W the constraint weight matrix with elements W as defined above, r is the vector with data residuals, and r is the vector with constraint residuals which is expressed as r=−DΔgglobal. The partial derivatives A and the data residuals r are evaluated using orbits determined with the global field GRGM900A. The anomalies for the global field are computed from its spherical harmonics expansion. The baseline parameters are not included in Δgadj but are accounted for by a modification of the data weight matrix W (by taking into account their Schur complement), following Sabaka et al. [2010]. The baseline parameters are unconstrained and free to adjust.

Figure 1

Maps in stereographic projection centered at the south pole for (a and b) free-air anomalies (up to 50°) and (c and d) Bouguer disturbances (up to 70°S). The local model is shown in Figure 1a, the difference between the local model and GRGM900A is shown in Figure 1b, GRGM900A (l = 7–900) is shown in Figure 1c, and a spherical harmonic transform of the local model (l = 7–900) is shown in Figure 1d.

3. Results

We processed the GRAIL extended mission data as outlined in the previous section. We had 1388 short arcs covering the selected south pole area. Our global starting gravity field model was GRGM900A. We weighted KBRR data at 0.1 μm/s throughout the entire extended mission because that is the level of fit with our current global models, but as in the processing of those models [Lemoine et al., 2014], we noted improvements when the low-altitude data were relatively downweighted. A factor of 0.5 on the normal system for the December arcs was found to give the best results in terms of data fit, anomaly patterns, and correlations with topography (all of which will be shown and explained below). We also experimented with the weight of the constraint system, denoted as μ in equation 2, and we found that a value of μ = 1·10−4 leads to the best overall results. We discuss this further at the end of this section.

Figure 1 shows free-air anomalies and Bouguer disturbances (the latter of which are related to the disturbing potential through δg =− ∂T/∂r, with a factor of (l + 1) in their spherical harmonic expansion, as opposed to the (l − 1)factor for anomalies) for our preferred local solution (Figures 1a and 1d for free air and Bouguer, respectively), Bouguer disturbances for GRGM900A (Figure 1c), and anomaly differences between the local solution and GRGM900A (Figure 1b). The free-air anomalies are shown up to 50°S, whereas the Bouguer disturbances are shown up to 70°S in order to focus on features close to the south pole. We also refer to Figures S2–S5 which show enlargements of these anomaly maps. The local model shows far fewer stripes than the global model, as the stripes in the difference map (Figure 1b) negate stripes in the anomaly map of GRGM900A (cf. Figure S1). We attribute this to our short-arc analysis where we adjust the orbit for each track, and, more importantly, to the neighbor smoothing that was applied to the full anomaly field, thus including the anomalies from the starting model. We also tested solutions with either a damping constraint (a diagonal constraint system, thus without neighbor smoothing) or with the neighbor smoothing applied to the adjustments only, and both solutions show more stripes in the anomalies. We also tested a neighbor-smoothing constraint with a varying correlation distance D (see section 2), where D was multiplied by the cosine of the latitude of the grid cell center point to account for varying size toward the pole. While this may slightly improve anomalies further away from the pole, the resulting small correlation distance at the pole leads to more striping in the solution there. An optimum between a constant and varying correlation distance may be found,but this was outside the scope of our current analysis. We also note that despite the decreasing physical size in longitude of the grid cells toward the pole, there are no apparent artifacts, likely thanks to the constant correlation distance applied.

Bouguer disturbances were computed by subtracting the contributions of topography to gravity, expressed in the same principal axis frame as gravity, using the formulation of Wieczorek and Phillips [1998]. We used Lunar Orbiter Laser Altimeter (LOLA) topography [Smith et al., 2010] archived in the Planetary Data System (LOLA/PDS Data Node, http://imbrium.mit.edu). We used a density of 2400 kg/m3 to account for the observed decrease in crustal density with higher degrees [Wieczorek et al., 2013; Han, 2013]. In order to compare two solutions of the same size so we can subtract the same topography-induced gravity potential, we generated global spherical harmonics for our local solution through a spherical harmonic transformation using Gauss-Legendre (GL) quadrature. We interpolated the adjustments Δgadj on a GL grid appropriate for a degree and order 1079 model, applied the spherical harmonic transformation, and constructed an updated global model by adding the resulting coefficients to the background global model. Our local model thus becomes a global model with a maximum size of degree and order 1079 which is only attained in the south pole area; elsewhere it is equivalent to the background model. For the Bouguer disturbances in Figure 1 we used both models up to degree 900. The Bouguer map for the local model shows fewer speckles and clearer features when compared to the map for the global model, which can help facilitate the geophysical interpretation of the structure of small-scale features [e.g., Zuber et al., 2014].

The free-air anomalies for the local solution in Figure 1a show a nonphysical speckle pattern at longitude 270°E, which is also where the adjustments are largest (Figure 1b). This is the area just below Mare Orientale where GRAIL attained its lowest altitudes above topography. We have indications that for these areas the current resolution of 1/6° by 1/6° of a degree is not sufficient to extract the full signal in the data [Zuber et al., 2014]. The speckle pattern is not as pronounced in the Bouguer disturbances in Figure 1, indicating that it appears beyond l = 900.

The fits to the KBRR data for the global and local models, using short arcs over the selected area for both, are shown in Figure 2. The fits for the local model were computed by adding the acceleration contributions of the adjustments Δgadj using equation 1 to the global model GRGM900A in the short-arc orbit determination. This shows that our local model fits the KBRR data better than the global starting model, sometimes by a considerable amount: on average, the fits improve by 12% for the whole extended mission and by 43% for the low December arcs only. There are still high fits for certain periods, most notably around 5 and 8 December, where the GRAIL satellites had their lowest periapsis heights. This indicates that information beyond degree and order 1080 is still present in the data, especially considering that the fit for those arcs is still well above the intrinsic noise levels of the data, which is thought to be between 0.05 and 0.07 μm/s (N. Harvey, Jet Propulsion Laboratory, Pasadena, CA, unpublished data, 2013). Figure 2 shows the fits for the last part of the GRAIL extended mission because the improvements are most clear there. The fits also improve for the earlier arcs in the extended mission, although there they are very close to the global model fits, at ∼0.12μm/s.

Figure 2

The root-mean-square (RMS) of fit to the KBRR data for GRGM900A and the local model, for the December arcs where the satellite altitude above topography was at its lowest.

Figure 3 shows correlations of various models (we again used the spherical harmonic transform for the local solution) with the gravity potential generated from topography. We show correlations localized over the south pole, following Wieczorek and Simons [2005], using a cap of 30°, with one taper and a concentration factor of 0.9999, resulting in a windowing function with Lwin=13. We chose the radius of 30° to avoid possible boundary effects at the edges of the solution, arising from the fact that we only use data over the area of interest, whereas data points close to the edge also contain information for unaccounted anomalies just beyond the grid. By testing small overlapping solutions, we found that the negative edge effects extend a few degrees at most, so our 10° cutoff is probably conservative.

Figure 3

Localized correlations with topography-induced gravity for various gravity field models. The localization used a cap with a radius of 30° centered on the south pole.

Figure 3 shows that the local model improves the correlations with topography considerably, and extends them to higher degrees. We also included the correlations for GRGM900C, a more recent model than GRGM900A, to show that the local solution also has higher correlations than that global model. Such high correlations are indicative of the pervasive fracturing of the lunar crust [Zuber et al., 2013b]. These results indicate that the smoothness of the maps in Figure 1 captures real gravity information and is not only the expression of the neighbor-smoothing constraint.

Our solutions depend strongly on the data weights W and the scaling factor μ for the constraint system; cf. equation 2. We generated solutions with various weights for the December data and with various scaling factors μ. We found that while anomaly maps for μ = 1·10−3, in general, look smoother (they do not show as many stripes below the Mare Orientale area), such solutions do not fit the data well, nor do they maintain correlations with topography as high as the μ = 1·10−4 solution. To assess the contribution of the constraints to the solution, we decomposed the gradient component of equation 2 into two parts, one describing the data (AWr) and one describing the constraints (μDWr). We then computed two separate solutions, each with only one of these components as a right-hand side to the equation system, which otherwise was the same for both solutions (thus including the constraint equation system μDWD). These solutions indicate the contribution of each component to the final solution. These solutions are shown in the Figure S6, showing that the data largely influence anomalies in areas where the spacecraft were low and that the constraints largely influence the anomalies where GRGM900A showed stripes. Finally, downweighting the December normal system reduces the striping in the area south of Mare Orientale, yet using a factor of 1.0 rather than 0.5 slightly increases the correlations with topography. This again indicates information present in the data beyond our current resolution. However, the KBRR fits for such a solution in arcs other than the December arcs, especially earlier ones at slightly higher altitudes, are not as good as that for a solution using a factor 0.5, and often even higher than the fits for the starting model. For these reasons, a solution using a factor of 0.5 on the December system and a scaling factor of μ = 1·10−4 is our preferred solution. For future solutions we might increase the grid resolution (in certain areas) to account for the information present in the low-altitude data.

4. Summary

We presented a local model of the gravity field of the Moon over the south pole estimated from GRAIL extended mission data. Our model is expressed as gridded gravity anomalies with a resolution of 1/6° by 1/6°, covering a cap over the south pole with a radius of 40°. We estimated the anomalies from a short-arc analysis using only KBRR data. We applied a neighbor-smoothing constraint to obtain our solution. Our local model removes striping present in the global models, and it improves fits to the KBRR data. We also transformed our local anomalies into spherical harmonics using Gauss-Legendre quadrature. This facilitates the comparison of Bouguer disturbances between the local and global models, and it allows us to compute correlations with topography. Our local model shows improved correlations with topography up to higher degrees than the global models. The spherical harmonic transformation from the original anomaly grid also has the benefit that it is easier to produce gravity disturbances and to compute gradients. Due to increased resolution and extended high correlations with topography, maps of the lunar gravity field such as from our local south pole model have the potential to improve the geophysical interpretation of small-scale features on the Moon.