1 Introduction

Geodetic measurements of a planetary body’s reaction to tidal forces, expressed by the tidal Love numbers \(k_2\) and \(h_2\), provide a key constraint on its interior. The Love numbers \(k_2\) and \(h_2\) characterize the response of the gravity field and the shape to the tidal potential, respectively. Measurements of \(k_2\) have contributed to the determination of the Moon’s core size (Williams et al. 2014), the detection of a subsurface ocean on Titan (Iess et al. 2012), and the rheological properties of Mercury’s mantle (Padovan et al. 2014; Genova et al. 2019). Measurements of \(h_2\) would reveal the size of Mercury’s solid inner core (van Hoolst et al. 2007; Steinbrügge et al. 2018a) and the thicknesses of Europa’s and Ganymede’s ice shells (Wahr et al. 2006; Kamata et al. 2016).

Fig. 1
figure 1

Love number \(h_2\) results and their accuracies from this study and literature values. Williams et al. (2014) do not report an accuracy, but we assume a linear propagation of the 0.9% error of their observed \(k_2\) value. For planets such as the Moon, with a small core and fairly uniform density in the mantle, \(h_2\) and \(k_2\) are rather tightly related

Full size image

Our understanding of the interior of the Moon has recently greatly improved through a re-analysis of seismic data from the Apollo era (Weber et al. 2011; Garcia et al. 2011) and the Gravity Recovery and Interior Laboratory (GRAIL) mission which determined the gravity field, the moment of inertia, and \(k_2\) with unprecedented accuracy (Zuber et al. 2013). Furthermore, the lunar \(h_2\) has been determined by Lunar Laser Ranging (LLR, Williams et al. 2013; Pavlov et al. 2016; Viswanathan et al. 2018 and by laser altimetry (Mazarico et al. 2014), making the Moon the only planetary body other than Earth, for which measurements of both \(h_2\) and \(k_2\) exist (Fig. 1).

Williams et al. (2014) created a family of lunar interior models which satisfy all GRAIL-based geodetic parameters, including in particular the value of the \(k_2\) Love number with a quoted uncertainty of 0.9%. At the same time, they retained the density profile derived from seismic constraints by Weber et al. (2011), but adjusted inner and outer core radii and the radius of the low-velocity zone in the lower mantle. All members of this model family give \(h_2 = 0.0424\). Because \(h_2\) and \(k_2\) are closely related, in particular for an object with a small core and rather uniform mantle density, this value represents arguably the best available estimate for the lunar \(h_2\). Three groups which determined \(h_2\) from LLR (Williams et al. 2013; Pavlov et al. 2016; Viswanathan et al. 2018) each used slightly different data sets and their own ephemerides. Williams et al. (2013) first incorporated the new gravity field and \(k_2\) results from the GRAIL mission into their DE430 ephemerides to find \(h_2 = 0.0476\pm 0.0064\). Pavlov et al. (2016) notably compared two different Earth tidal models and found \(h_2 = 0.0430\pm 0.0010\) and \(h_2 = 0.0410\pm 0.0010\) for the DE430 model and a model recommended by the International Earth Rotation and Reference Systems Service (IERS), respectively. Viswanathan et al. (2018) first presented results including new infrared LLR data and found \(h_2 = 0.0439\pm 0.0002\). Mazarico et al. (2014) used altimetric data and a method where the radial displacement is determined through interpolation of intersecting ground tracks at crossover points which yields \(h_2 = 0.0371\pm 0.0033\). There are thus significant differences between the \(h_2\) value derived from \(k_2\) measurements (Williams et al. 2014; Williams and Boggs 2015) and the results of LLR (Williams et al. 2013; Pavlov et al. 2016; Viswanathan et al. 2018) on the one hand, and the \(h_2\) value derived from laser altimetry using the crossover method (Mazarico et al. 2014) on the other hand, as also previously noted by Matsumoto et al. (2015) and Harada et al. (2016).

This study uses laser altimetry to determine the lunar \(h_2\), but instead of relying on crossover points it applies a different method which simultaneously solves for \(h_2\) and coefficients which represent the global topography of the Moon on an equirectangular grid, as originally proposed by Koch et al. (2008, 2010). The parameterization of the topography into 2D cubic B-splines provides an advantage over the method of Koch et al. (2010) who only used cubic splines in longitude direction.

For the retrieval of \(h_2\), we use data from the Lunar Orbiter Laser Altimeter (LOLA). LOLA is a multi-beam laser altimeter with a ground pattern of five spots with a distance of approximately \(25\,\mathrm {m}\) from each other. The primary objective of LOLA is the generation of topographic maps of the lunar surface with appropriate resolution and accuracy for future robotic and human exploration (Smith et al. 2010a). Laser pulses are fired at a frequency of \(28\,\mathrm {Hz}\). While the measurement precision on a flat surface is \(10\,\mathrm {cm}\) (Smith et al. 2017), knowledge of the spacecraft position and orientation limit the radial accuracy to about \(1\,\mathrm {m}\) (Mazarico et al. 2018). For a review of LOLA’s achievements see Smith et al. (2017). LOLA is a payload of the Lunar Reconnaissance Orbiter (LRO), which has been orbiting the Moon since June 2009. From September 2009 to December 2011, it was in a near-circular \(50\,\mathrm {km}\) mapping orbit and entered an elliptical orbit afterward (Mazarico et al. 2018).

On other solar system bodies, such as Mercury (van Hoolst et al. 2007; Steinbrügge et al. 2018a), Ganymede (Moore and Schubert 2003; Jara-Orué and Vermeersen 2016; Kamata et al. 2016; Kimura et al. 2019), and Europa (Moore and Schubert 2000; Wahr et al. 2006; Steinbrügge et al. 2018b), \(h_2\) might be retrievable with a higher accuracy due to stronger tides, with amplitudes on the order of one to several meters compared to \(15\,\mathrm {cm}\) at the Moon. With increasing accuracy, tighter constraints can be emplaced on the respective interior structures. Due to the lack of seismic data and highly accurate \(k_2\) measurements from dedicated gravity missions like GRAIL on those bodies, a determination of \(h_2\) is crucial. Future missions strive to detect the body tide of Mercury (BepiColombo; Koch et al. 2008, 2010; Thor et al. 2020) and Ganymede (Jupiter Icy Moons Explorer Steinbrügge et al. 2015, 2019) by laser altimetry.

The goal of this study is thus twofold: Further constrain the lunar \(h_2\); and demonstrate the utility of the applied method for future application to other solar system bodies.

Fig. 2
figure 2

The tidal potential at selected locations on the Moon as a function of time, computed from Eq. 1 using the DE421 ephemerides (Williams et al. 2008)

Full size image
Fig. 3
figure 3

The range of the tidal potential exerted on the Moon by the Earth at each location on the surface over the period of the LRO circular orbit phase (September 2009 to December 2011), computed using Eq. 1 and the DE421 ephemerides (Williams et al. 2008). The determination of the range uses a time discretization of \(10^5\,\mathrm {s}\). The equivalent radial tidal displacements are computed with Eq. 9 assuming \(h_2 = 0.0371\) (Mazarico et al. 2014)

Full size image

2 Theory

The tidal potential at any point on the surface of the Moon due to the Earth is given by (Murray and Dermott 1999)

$$\begin{aligned} V_{2, \mathrm {tot}}(\theta , \lambda , t) = -\dfrac{\mu _{\mathrm {Earth}} R^2}{2r^3(t)}(3\cos ^2\psi (\theta , \lambda , t) - 1), \end{aligned}$$
(1)

where \(\mu _{\mathrm {Earth}} = 398{,}600.435436 \,\mathrm {km^3s^{-2}}\) (Folkner et al. 2014) is the gravitational parameter of the Earth, \(R = 1737.4 \,\mathrm {km}\) (Archinal et al. 2018) is the radius of the Moon, r is the distance between the centers of mass of the Earth and the Moon, \(\psi \) is the Moon-centric angle between the point on the surface and the center of mass of the Earth, \(\theta \) is colatitude, \(\lambda \) is longitude, and t is time. Eq. 1 only considers the second degree of the potential because each higher degree is smaller by a factor \(R/r \approx 220\). Using precise ephemerides, r(t), \(\psi (\theta , \lambda , t)\), and subsequently \(V_{2, \mathrm {tot}}(\theta , \lambda , t)\) can be determined to high accuracy for any given location on the Moon’s surface. Figures 2 and 3 show the tidal potential as a function of time and its range over the period of the LRO circular orbit phase, from which measurements are used in this study, considering also the potential exerted by the Sun from the analog description to Eq. 1 using \(\mu _{\mathrm {Sun}} = 132{,}712{,}440{,}041.9394\, \mathrm {km^3\,s^{-2}}\). The potential exerted by the Sun is about 70 times weaker than the potential exerted by the Earth. The asymmetry in Fig. 3 is due to the integration time of 2.24 years, which is short compared to the axial precession period of 18.6 years.

Kaula (1961, 1964) transformed Eq. 1 into Keplerian elements. Then, using the fact that the Moon is in a 1:1 spin–orbit resonance with Earth, the tidal potential can be written as

$$\begin{aligned} V_{2, \mathrm {tot}}(\theta , \lambda , t)= & {} -\dfrac{\mu _{\mathrm {Earth}} R^2}{a^3}\left( f_1 P_{20}(\cos \theta ) \right. \nonumber \\&\left. +\,f_2 P_{22}(\cos \theta )\cos (2\lambda ) \right. \nonumber \\&\left. + \,f_3 P_{22}(\cos \theta ) \sin (2\lambda )\right) , \end{aligned}$$
(2)

where

$$\begin{aligned} f_1= & {} -\dfrac{1}{2}\sum _{q = -\infty }^\infty G_{21q}(e)\cos (qM(t)) \end{aligned}$$
(3)
$$\begin{aligned} f_2= & {} \dfrac{1}{4}\sum _{q = -\infty }^\infty G_{20q}(e)\cos (qM(t)) \end{aligned}$$
(4)
$$\begin{aligned} f_3= & {} -\dfrac{1}{4}\sum _{q = -\infty }^\infty G_{20q}(e)\sin (qM(t)), \end{aligned}$$
(5)

a is the semi-major axis, e is the eccentricity, \(P_{lm}\) are the unnormalized associated Legendre polynomials of degree l and order m, \(G_{lpq}(e)\) are eccentricity polynomials given by Kaula (1961), Eqs. 23–26, and M(t) is the mean anomaly of the Earth in its orbit around the Moon. The derivation of Eq. 2 neglects the \(6.7^\circ \) obliquity of the Moon with respect to its orbital plane (Ward 1975).

Fig. 4
figure 4

The static component of the tidal potential computed using the analytical approximation of Eq. 6

Full size image

Because of the locked rotation of the Moon, the tidal potential can be decomposed into a static and a dynamic part. The static part can be obtained either from a spectral analysis of a numerical evaluation of the tidal potential over long periods of time or from regarding the \(q=0\) component of Eqs. 2–5. The latter gives a time-independent component

$$\begin{aligned} V_{2, \mathrm {stat}}(\theta , \lambda , t)= & {} -\dfrac{\mu _{\mathrm {Earth}} R^2}{a^3}\left( -\dfrac{1}{2}G_{210}(e) P_{20}(\cos \theta ) \right. \nonumber \\&\left. + \,\dfrac{1}{4}G_{200}(e) P_{22}(\cos \theta )\cos (2\lambda )\right) . \end{aligned}$$
(6)

This analytic approximation to the static component of the lunar tidal potential is shown in Fig. 4 for \(a = 384{,}400 \,\mathrm {km}\) and \(e = 0.0554\) (Standish 2001). Slightly different solutions are obtained from numerical evaluations, depending on analysis method and evaluation period.

The dynamic part of the lunar tidal potential

$$\begin{aligned} V_{2, \mathrm {dyn}} = V_{2, \mathrm {tot}} - V_{2, \mathrm {stat}} \end{aligned}$$
(7)

derives from the eccentricity of the Moon’s orbit, its obliquity, its non-uniform rotation, and solar tides. Its frequency spectrum lies in the range of days to several years, with the highest amplitude at the \(27.2\,\mathrm {d}\) monthly period (Williams and Boggs 2015).

The Love numbers characterize the response of the planet to the tidal potential exerted on it by other celestial bodies. The \(k_2\) Love number describes the resulting secondary gravitational potential

$$\begin{aligned} V' (\theta , \lambda , t) = k_2 V_2 (\theta , \lambda , t) \end{aligned}$$
(8)

caused by internal mass re-distribution due to tidal forcing. The \(h_2\) Love number describes the radial displacement of the surface

$$\begin{aligned} u_\mathrm {r} (\theta , \lambda , t) = h_2\dfrac{V_2 (\theta , \lambda , t)}{g}, \end{aligned}$$
(9)

where \(g = \mu _{\mathrm {Moon}}/R^2 = 1.62 \,\mathrm {ms^{-2}}\) is the gravitational acceleration at the surface and \(\mu _{\mathrm {Moon}} = 4902.800066 \,\mathrm {km^3s^{-2}}\) is the gravitational parameter of the Moon (Folkner et al. 2014). The Love numbers are functions of the spatial distribution of material properties in the planet’s interior and can be computed from radial profiles of density, shear modulus, and viscosity (Segatz et al. 1988; Moore and Schubert 2000). Equations 8 and 9 are the simplest description of tidal deformation and neglect viscoelasticity, lateral heterogeneities, and the frequency dependence of the Love numbers. A more detailed description would use different Love numbers for each of the frequency components in the tidal potential. However, most of the power of the tidal potential is concentrated within a tight band around the \(27.2\,\mathrm {d}\) monthly period, as is evident from Eq. 2, when considering that terms with \(|q|>1\) are second-order effects in eccentricity. Therefore, the descriptions in Eqs. 8 and 9 with a single Love number, which have been frequently applied in previous studies (e.g., Mazarico et al. 2014), are good approximations to the real secondary gravitational potential and radial surface displacement.

A special case of the frequency dependence of the Love numbers is the reaction at zero frequency. When inserting the dynamic potential \(V_{2, \mathrm {dyn}}\) as \(V_2\) in Eq. 9, the resulting radial displacement is the tidal deformation with a maximum peak-to-peak amplitude \(< 30\,\mathrm {cm}\) (Fig. 3). When inserting the total potential as \(V_2\) in Eq. 9, while using the same \(h_2\) valid for the tidal frequency band, the resulting radial displacement would contain an additional static bulge with a maximum deformation of \(\sim 50\) cm, proportional to the potential shown in Fig. 4. The actual deformation of the Moon in response to the static component of the tidal potential is of course much larger, because this forcing acts at a much larger time scale, where the Moon reacts as a fluid. It cannot be described by the same \(h_2\) Love number used to express deformation in the tidal frequency band. Furthermore, the actual static tidal bulge of the Moon was frozen in early in its history when the Moon was on a closer orbit and had a weaker lithosphere (Keane and Matsuyama 2014; Qin et al. 2018).

3 Methods

We simultaneously extract the lunar solid body tide and the static global topography. A single LOLA observation of the topographic elevation consists of the static, time-invariant topography \(T_\mathrm {stat}\) at that location, the radial displacement \(u_\mathrm {r}\) of the surface due to tides at time \(t_k\) (Eq. 9), and measurement and model errors, contained in the term \(e_k\):

$$\begin{aligned} T_k (\theta _k, \lambda _k, t_k) = T_\mathrm {stat}(\theta _k, \lambda _k) + u_\mathrm {r}(\theta _k, \lambda _k, t_k) + e_k. \end{aligned}$$
(10)

Here, \(k = 1, ..., K\) is the index of the individual topographic observations, of which there can be up to 5 per laser shot. Modeling the topography only serves the purpose of removing it from the measured signal and limiting the size of the error term \(e_k\). Accurate models for the topography of the Moon that are suitable for further analysis are available from LOLA data (Smith et al. 2010b) and from LOLA and photogrammetric data (Barker et al. 2016).

The static part of the topography is usually simply referred to as the shape of the Moon. Here, we parameterize the static topography as an expansion in 2D cubic B-spline basis functions. The expansion can be written as

$$\begin{aligned} T_\mathrm {stat}(\theta _k, \lambda _k) = \sum _{i=1}^I\sum _{j=1}^J c_{ij} S_{ij} (\theta _k, \lambda _k) \end{aligned}$$

where \(S_{ij}\) are the spline basis functions depending on position, \(c_{ij}\) are their coefficients, and I and J are the number of splines that are used in latitude and longitude direction, respectively, with \(N = I\cdot J\) being their total number. For the definitions of the spline functions \(S_{ij}\), see Koch et al. (2010), Eq. 9, 11, 14–17. 2D cubic B-spline basis functions have the property that at any point \((\theta _k, \lambda _k)\) on the surface, only 16 functions are nonzero. The splines are defined on the equirectangular projection of the spherical surface onto a \((\theta , \lambda )\)-plane. Each spline function is centered at one grid point (i, j) on the map and is nonzero in the two adjacent grid intervals in longitude and latitude in both directions. The projections of the grid cells back onto the sphere are squares at the equator (\(J = 2I\)) and evolve over trapezoids in higher latitudes to triangles right at the poles. Approaching the poles, the splines become shorter in longitudinal direction, while staying the same length in latitudinal direction. However, the usage of an equirectangular projection is necessary to accommodate the 2D cubic B-spline basis functions on a sphere and to provide an efficient computation scheme (Koch et al. 2010). Any topographic signal with a wavelength smaller than the grid cell size \(360^\circ /J\) cannot be modeled by the splines and contributes to the model error term \(e_k\). Our approach is to assume that the topography at smaller scales is independent and identically distributed (i.i.d.). The parameterization using 2D cubic B-splines is advantageous over, e.g., a spherical harmonic expansion because locality is retained while still allowing for a high smoothness of the solution. Additionally, an expansion in spherical harmonics would be computationally too expensive (Koch et al. 2010). Our method is an improvement over that developed by Koch et al. (2010), who used cubic splines only in longitudinal direction and step functions in latitudinal direction. Their justification for using step functions in latitudinal direction was the denser data coverage in that direction and a reduction in computational cost. However, seeing that using cubic splines instead of step functions in longitudinal direction provided a crucial improvement for their retrieval accuracy of \(h_2\), they recommended doing the same in latitudinal direction in a further study. As the involved computational challenge has become manageable by now, we apply cubic splines, which are smooth enough to model planetary topography well, in both coordinate directions.

We gather the spline coefficients \(c_{ij}\) and \(h_2\) in a parameter vector \({\varvec{x}}\) and use a least-squares adjustment to solve the observation equation (Eq. 10) for \({\varvec{x}}\). The resulting normal equation system is highly sparse, band-structured, and has \(N+1\) equations and \(51N + 1\) nonzero elements, which can be reduced to \(26N + 1\) elements due to symmetry. Because of the inhomogeneous coverage of the lunar surface with measurements, it can occur that some grid cells are poorly sampled or not sampled at all. Such data gaps cause instabilities in the solution. The linear inverse ill-posed problem must therefore be regularized, minimizing

$$\begin{aligned} ({\varvec{Ax}} - {\varvec{T}})^\top ({\varvec{Ax}} - {\varvec{T}}) + \alpha {\varvec{x}}^\top {\varvec{R}}{\varvec{x}} \end{aligned}$$

in a least-squares sense, where the vector \({\varvec{T}}\) contains the K observations \(T_k\), \({\varvec{A}}\) is the design matrix derived from Eq. 10, \(\alpha \) is a regularization parameter, and \({\varvec{R}}\) is the regularization matrix which constrains the adjustment by minimizing the second derivative of the topography at the grid points \((\theta _i, \lambda _j)\)

$$\begin{aligned}&{\varvec{\nabla }}\cdot {\varvec{\nabla }}T_\mathrm {stat}(\theta _i, \lambda _j) \nonumber \\&\quad = \sum \limits _{r = i - 1}^{i + 1}\sum \limits _{s = j - 1}^{j + 1}c_{rs}{\varvec{\nabla }}\cdot {\varvec{\nabla }}S_{rs}(\theta _i, \lambda _j)\nonumber \\&\quad = \sum \limits _{r = i - 1}^{i + 1}\sum \limits _{s = j - 1}^{j + 1}c_{rs}\left( \dfrac{\partial ^2}{\partial \theta ^2}S_{rs}(\theta _i, \lambda _j) \right. \nonumber \\&\qquad \left. + \,\dfrac{1}{\sin ^2\theta }\dfrac{\partial ^2}{\partial \lambda ^2}S_{rs}(\theta _i, \lambda _j)\right) . \end{aligned}$$
(11)

In Eq. 11, we have applied a Cartesian Laplace operator with a correction term for the change of cell width as a function of latitude for simplicity. The estimated parameters are given by

$$\begin{aligned} {\varvec{{\hat{x}}}} = ({\varvec{A}}^\top {\varvec{A}} + \alpha {\varvec{R}})^{-1}{\varvec{A}}^\top {\varvec{T}} \end{aligned}$$
(12)

and the adjustment residuals are

$$\begin{aligned} {\varvec{{\hat{e}}}} = {\varvec{T}}-{\varvec{A{\hat{x}}}}. \end{aligned}$$
(13)

The regularization enables a trade-off between smoothing the topography solution and minimizing the inevitable biases in the results. The bias can be approximated as (Xu 1992)

$$\begin{aligned} \mathrm {bias}({\varvec{{\hat{x}}}}) \approx ({\varvec{A}}^\top {\varvec{A}} + \alpha {\varvec{R}})^{-1}\alpha {\varvec{R{\hat{x}}}}, \end{aligned}$$
(14)

where \({\varvec{{\hat{x}}}}\) is the biased parameter solution. We apply randomized generalized cross-validation (Kusche and Klees 2002) to determine the regularization parameter causing the lowest bias, and find that it is always the smallest possible regularization parameter that still stabilizes the normal equation matrix enough for it to be solvable. The value of the smallest possible \(\alpha \) depends on both, the number of observations K, and the number of parameters N. We empirically derived the relation \(\alpha \approx 10^{-8} K / N\), which simplifies the cumbersome procedure of determining an optimal regularization for each pair of N and K. Computing the estimated bias from Eq. 14 for each solution assures us that the chosen regularization was not too strong. To compromise between obtaining an unbiased \(h_2\) result and a smooth topography solution, we set \(\alpha = 10^{-3} K / N\). While this represents a stronger regularization than required for a stable solution, it does not cause any significant bias on the final \(h_2\) result (see Sect. 5.4). A parallel direct sparse solver solves the normal equation system.

To determine the formal accuracy of the resulting \(h_2\) value, we estimate its variance from the adjustment residuals. The variance of unit weight (e.g., Xu et al. 2006) can be estimated from the adjustment residuals as

$$\begin{aligned} {\hat{\sigma }}^2=\dfrac{1}{K-(N+1)}{\varvec{{\hat{e}}}}^\top {\varvec{{\hat{e}}}}. \end{aligned}$$

The estimated covariance of estimated parameters is then given by

$$\begin{aligned} \hat{{\varvec{Q}}}_{\hat{{\varvec{x}}}}={\hat{\sigma }}^2({\varvec{A^\top A}})^{-1}. \end{aligned}$$
(15)

The normal equation matrix actually used in the computation of the covariance augments \({\varvec{A^\top A}}\) by the regularization term \(\alpha {\varvec{R}}\). Since this regularization is small, its influence on the variance is negligible.

In the context of this study, it is irrelevant for the estimation of \(h_2\), whether the dynamic or total potential (or, in fact, the dynamic potential plus any constant term) is used in Eq. 9. This is because \(h_2\) is the only parameter in our model that is sensitive to temporal variations. Whether the definition of the tidal deformation \(u_\mathrm {r}\) includes some of the static topography (namely the aforementioned \(\sim 50\) cm bulge) or not, only affects the parameters that model the static topography, but the result for \(h_2\) will be the same in both cases. The static effect contained in the time-dependent topography is an extremely smooth global signal (see Fig. 4) and therefore can be described precisely by the high-resolution parameterization of the static topography \(T_\mathrm {stat}\). In this study, we have used both the total and the dynamic potential. Like Mazarico et al. (2014), we have computed the total potential caused by the Earth and the Sun from the DE421 ephemerides and masses as given by Spacecraft, Planet, Instrument, Camera-matrix, Events (SPICE) kernels (Williams et al. 2008; Acton et al. 2018). To compute the dynamic potential, we have subtracted the analytic approximation of Eq. 6 shown in Fig. 4. Both approaches have yielded identical results for \(h_2\), but topography results which differ by the \(\sim 50\) cm time-invariant part. In order to compare the static topography derived from the total potential with other topographic models, one would have to correct for this missing part by adding the displacement caused by the static part of the tidal potential first.

Only data which are located in a region containing measurements from different tidal phases contribute to the determination of \(h_2\). These regions usually, but not necessarily, contain one or multiple crossovers. With finer grid resolution, the size of these regions, and thereby the amount of contributing measurements, decreases. Measurements outside of such regions do not contribute to the determination of \(h_2\), but they do not bias it either, because they have sufficient freedom to fit any value of \(h_2\). Only in the theoretical case where no sufficiently small regions with measurements of different tidal phases exist, will the topography coefficients attempt to model the tidal signal and cause a bias in the \(h_2\) result. This would be the case when the grid cell size is smaller than the typical shot-to-shot distance.

4 Data

We work with a subset of the over \(7\times 10 ^ 9\) range measurements that have been recorded and published on the Planetary Data System (PDS, Neumann 2009a). As a first filtering step, we only include LOLA shots within a pre-defined topographic range of \(\pm 13\,\mathrm {km}\) and observed at solar phase angles of \(<90^\circ \). This ensures that extreme outliers in height and LOLA shots that were recorded during the night where the so-called LOLA thermal blanket anomaly occurs do not enter the evaluation (Smith et al. 2010b, 2017). As a second step, we divide the surface of the Moon in 16 tiles and for each derive a coarse \(400\,\mathrm {m}/\mathrm {pixel}\) LOLA Digital Terrain Model (DTM) compiled from all available tracks from all mission phases. Each track is then co-registered to the 16 tiles (Gläser et al. 2013, 2018). At each tile, we sort out the entire track segment if the mean height difference or the standard deviation to the DTM tile is \(>100\,\mathrm {m}\) or if less than 500 points were co-registered to the tile. As a third step, to ensure a homogeneous spatial distribution, we select only measurements from the near-circular orbit phase of LRO extending over the 27 months between September 2009 and December 2011 (orbit numbers 1005–11,403), when LOLA achieved global data coverage. As a fourth step, we run a first adjustment as described in Sect. 3. Visual inspection of the spatial distribution of adjustment residuals \({\varvec{{\hat{e}}}}\) (Eq. 13) reveals five further outlier orbits (orbit numbers 1803, 2351, 7756, 10,302, 11,226) which we removed for a remaining total of \(K = 3{,}686{,}466{,}983\) measurements from 10,016 orbits.

5 Results

5.1 Complete data set

First, we use the complete data set to assess the sensitivity of the \(h_2\) estimate on the resolution for the static topography by varying the number N of base spline functions. Results for resolutions between 5 and 31 grid points per degree are shown in Fig. 5. For example, a resolution of 16 grid points per degree needs \(N = 2 I^2 = 2\cdot (16 \cdot 180) ^ 2 = 16,588,800\) splines and is equivalent to a cell size of \(1.9\,\mathrm {km}\cdot 1.9\,\mathrm {km}\) at the equator. Results for \(h_2\) from low-resolution models scatter widely, but this scatter decreases drastically at high resolutions. The maximum resolution is limited by the computational expense. Any topographic signal at smaller scales than the distance between two grid points cannot be modeled and contributes to the model error \(e_k\) of a specific measurement \(T_k\). A resolution of five grid points per degree corresponds to \(6.1\,\mathrm {km}\cdot 6.1\,\mathrm {km}\) at the equator. The increase in the model error with decreasing resolution is also reflected by the increase in the formal error of \(h_2\) (Eq. 15), from \(3\times 10^{-4}\) at 31 grid points per degree to 0.0014 at 5 grid points per degree. These formal errors are clearly smaller than the variation of the results with changing resolution. They are therefore to be understood as an indication of precision rather than of accuracy. From the decrease in both scatter and formal uncertainty with increasing resolution, it becomes clear that modeling the topography at high resolution is necessary to achieve reliable results.

Fig. 5
figure 5

Love number \(h_2\) from complete data set (blue) and reduced data set, only using every sixth data point (red), as a function of resolution of the topographic grid. For better visibility, data points are plotted with a slight offset from the integer value for resolution. The error bars indicate the formal error of the adjustment. The blue horizontal line and shaded area show the weighted mean and uncertainty of the results from the complete data set

Full size image
Fig. 6
figure 6

Top: Topography overlain by a hillshade. Middle: RMS of the adjustment residuals at a resolution of 16 grid cells per degree, computed over one grid cell. The global RMS residual is 24.8 m at this resolution. Bottom: Difference between the topography of the LOLA Gridded Data Record (GDR) V2 data set (Neumann 2014) and the topography model derived in this study at each \(1/16^\circ \cdot 1/16^\circ \) grid cell. The color gray indicates grid cells without measurements. For clarity, only a region on the lunar surface is shown

Full size image
Fig. 7
figure 7

Like Fig. 6, but for different region on the Moon centered at \(35^\circ \,\mathrm {N}, 20^\circ \, \mathrm {E}\)

Full size image
Fig. 8
figure 8

Histogram of the \(K = 3{,}686{,}466{,}983\) adjustment residuals \({\varvec{{\hat{e}}}}\) for different resolutions of the topographic model

Full size image

To obtain an additional check on the results, we compute the residuals of the adjustment \({\varvec{{\hat{e}}}}\) (Eq. 13). We then compute the root mean square (RMS) of the residuals over one grid cell and display the result in Figs. 6 and 7 for two regions of the lunar surface. Regions with a smooth topography, such as maria and crater plains, show small adjustment residuals, whereas high values are found in rougher regions, such as the highlands in general, crater rims, and small craters within the maria. The residuals depend on the resolution of the topography. The RMS residual decreases from 111 m at a resolution of 5 grid points per degree to 9.1 m at a resolution of 31 grid points per degree (see also Fig. 8). Dedicated DTMs obtain lower residuals with respect to the LOLA measurements. For example, SLDEM (Barker et al. 2016) achieves RMS vertical residuals of \(\sim 3\) m. This is possible because of their much higher resolution of 512 pixels per degree. With a \(\sim 16\) times coarser resolution, we achieve a residual that is only \(\sim 3\) times larger. This indicates that our spline model for the topography fits the LOLA measurements very well.

Fig. 9
figure 9

Histogram of the difference between the topography of the LOLA Gridded Data Record (GDR) V2 data set (Neumann 2014) and the topography model derived in this study at each of the \(N = 16{,}588{,}800\) \(1/16^\circ \cdot 1/16^\circ \) grid cells

Full size image

We also compare our topographic solution with a DTM derived from laser altimetry (Figs. 6, 7). Again, we notice a correlation with the smoothness of the terrain. The two terrain models agree well in smooth terrain and have larger disagreements in rougher terrain, where the impact of the choice of interpolation methods is large. The global RMS difference between the two models is \(61\,\mathrm {m}\) and the mean difference is \(46\,\mathrm {cm}\). We emphasize that the topographic model produced in this study is merely a by-product. It serves to show a good agreement with the actual topography of the Moon. For further geological or geophysical analysis, we recommend the usage of dedicated terrain models such as SLDEM (Becker et al. 2016) or the LOLA GDR (Neumann 2014).

5.2 Reduced data set

To carry out tests on synthetic data for further characterization of the error sources, we first define a reduced set of real data containing only every sixth measurement. It consists of 614,411,186 measurements while maintaining a representative geometry and distribution of used receiver channels. Again, we solve for \(h_2\) and vary the resolution of the topographic grid (Fig. 5), revealing a difference with respect to the complete data set of \(\le 1.5\times 10^{-4}\) for \(h_2\) and a formal error that is larger by a factor of approximately 2.5. The detection of tidal displacement with an amplitude \(\lesssim 10\) cm in the presence of much larger adjustment residuals and RMS differences to topography models is only possible due to the large amount of observations and the time-dependence of the tidal signal.

Fig. 10
figure 10

Power spectra of the topography of the Moon from the LOLA SHADR model (red), from the SLDEM complemented with the LDEM (blue), and of one realization of a randomly generated topography from \(l = 300\) to \(l = 7999\) (yellow). The simulated spectrum is constrained by a power law which fits the LOLA SHADR power spectrum from \(l = 300\) to \(l = 2050\), the maximum degree of the spherical harmonic expansion of the LOLA SHADR. The difference between the LOLA SHADR and SLDEM spectra is barely discernible below \(l\approx 1000\)

Full size image

5.3 Synthetic data set

We synthetically generate 100 datasets that use the actual footprint positions and epochs of the reduced dataset and simulated topography measurements. This allows the evaluation of the actual measurement geometry in the most direct way. Since the difference between the results of the complete and the reduced datasets is small, conclusions drawn from a reduced synthetic dataset can be applied to the complete dataset. Each simulated measurement consists of four parts:

  1. 1.

    A spherical harmonic model of the observed lunar topography up to degree 299 (equivalent to a resolution of 18.2 km at the equator). We use the LOLA SHADR model exclusively using LOLA data (Neumann 2009b; Smith et al. 2017). We also calculated a power spectrum from the SLDEM by complementing it with the LDEM (Neumann 2014) at latitudes polewards of \(60^\circ \). Their power spectra are very similar up to \(l\approx 1000\) (Fig. 10).

  2. 2.

    A synthetic spherical harmonic model of the lunar topography from degree 300 to 7999 (equivalent to a resolution of 680 m at the equator) according to a power law \(al^b\). The coefficients of the power law are chosen as \(a = 3\times 10^9\,\mathrm {m^2}\) and \(b = -2.8\) to approximate the observed topography at higher degrees and to allow at the same time to vary its detailed structure for test purposes (Fig. 10). Each spherical harmonic coefficient of degree l is randomly generated according to a normal distribution with zero mean and variance \(\sigma ^2 = al^b (2l + 1) ^ {-1}\). We do not extend the synthetic model to higher degrees than 7999 due to the computational expense. The spherical harmonic model is transformed into the spatial domain (Schaeffer 2013) and sampled at the laser spot coordinates using second-order Lagrangian interpolation.

  3. 3.

    The lunar topography at the smallest scales generated using Gaussian process regression. Each measurement’s expectancy and variance are computed based on the topography values of previous nearby measurements and reflect the measurement error and the power contained in the topography at smaller scales than in the spherical harmonic model. The expectancy and variance are given as \({\varvec{S}}_{UK}{\varvec{S}}_{KK}^{-1}{\varvec{T}}_\mathrm {ss}\) and \(\sigma _\mathrm {ss}^2 + \sigma _\mathrm {re}^2 - {\varvec{S}}_{UK}{\varvec{S}}_{KK}^{-1}{\varvec{S}}_{UK}^\top \), respectively, where \({\varvec{T}}_\mathrm {ss}\) are the small-scale topography values of previous observations, \({\varvec{S}}_{KK}\) is the covariance matrix of previous observations, \({\varvec{S}}_{UK}\) is a vector containing the covariances between the new and all previous observations, \(\sigma _\mathrm {ss}^2\) is the variance of the small-scale topography, and \(\sigma _\mathrm {re}\) is the range error. The covariance is assumed isotropic and following a Matérn model. The Matérn class of covariance functions provides a flexible description (Matérn 1986; Handcock and Stein 1993; Guttorp and Gneiting 2006; Guinness and Fuentes 2016). It is given by

    $$\begin{aligned} \mathrm {cov}(\psi )={\left\{ \begin{array}{ll} \sigma _\mathrm {ss}^2 + \sigma _\mathrm {re}^2,&{} \text {if } \,\,\psi = 0\\ \sigma _\mathrm {ss}^2\dfrac{2^{1-\nu }}{\varGamma (\nu )} \left( \dfrac{2\sqrt{\nu }}{\rho }\psi \right) ^\nu {\mathcal {K}}_\nu \left( \dfrac{2\sqrt{\nu }}{\rho }\psi \right) ,&{} \text {if }\,\, \psi >0 , \end{array}\right. } \end{aligned}$$

    where \(\psi \) is the spherical distance between the two observations, \(\nu \) is the smoothness parameter, \(\rho \) is the decorrelation distance, \(\varGamma \) is the gamma function, and \({\mathcal {K}}_\nu \) is the modified Bessel function of the second kind of order \(\nu \). Its spectral form

    $$\begin{aligned} M_l = \sigma _\mathrm {ss}^2 \dfrac{\varGamma (\nu + \frac{1}{2})}{\sqrt{\pi } \varGamma (\nu )(2l + 1)} \left( \dfrac{4 \nu }{\rho ^2}\right) ^\nu \left( \dfrac{4 \nu }{\rho ^2} + l^2\right) ^{-\nu -\frac{1}{2}} \end{aligned}$$

    is asymptotic to a power law with exponent \(b = -2 (\nu + 1)\) for large l (Guinness and Fuentes 2016). We assume \(\sigma _\mathrm {re} = 1\,\mathrm {m}\) for the LOLA range error (Mazarico et al. 2018). The topographic power at small scales is obtained by an integration of the power law from \(l = 8000\) to infinity and amounts to \(\sigma _\mathrm {ss} = 12.6\,\mathrm {m}\). We derive the smoothness parameter \(\nu = 0.4\) from the power law exponent \(b = -2.8\) to be consistent with the assumption made for the large-scale topography. For the decorrelation distance, the resolution of the spherical harmonic model in the latitudinal direction \(0.01125^\circ \) (\(341\,\mathrm {m}\)) is assumed. We only consider previous measurements of the same ground track to significantly reduce the computational load by ignoring the correlation of the small-scale topography between different ground tracks. Due to the linear nature of the locations of known footprints that are considered, the expectancy and variance of the new measurement mostly depends on a handful of previous measurements. We consider all previous measurements with a covariance \(\mathrm {cov}(\psi ) \ge 10^{-3} \sigma _\mathrm {ss}^2\), which is reached for at most 47 observations.

  4. 4.

    The radial displacement due to tides (Eq. 9) using an a priori \(h_2\) of 0.04.

Fig. 11
figure 11

RMSE and bias of \(h_2\) obtained from 100 synthetically generated data sets using the footprint positions of the reduced data set, as a function of the resolution of the topographic grid

Full size image

We determine \(h_2\) from each of the 100 synthetically generated sets of 614,411,186 measurements, which are based on the reduced data set. From the 100 results, we compute the root-mean-square error (RMSE) which describes the spread around the a priori value for \(h_2\). The RMSE generally decreases with increasing resolution (Fig. 11), from 0.0238 at 5 grid points per degree and reaching a minimum of 0.0014 at 30 grid points per degree. We also compute the bias, defined as the difference between the mean of the results from the 100 realizations and the a priori value for \(h_2\), and find that it ranges from \(-\,0.0017\) to 0.0023. However, the fact that the bias spreads evenly around zero for different resolutions, with a mean of \(-\,0.0003\) and a standard deviation of 0.0009, indicates that it is caused rather by the specific geometry associated with a certain grid resolution, than by the data themselves. The observed bias likely originates from ignoring the autocorrelation of topography between measurements of different tracks, as this causes mismatches of several meters within the crossover region. We tested this statement by simulating both cases, considering and ignoring the mutual correlation of different tracks, for a single crossover. The recovery of a constant offset between the two tracks, representing a time-dependent signal, yields a bias which is \(\sim 30\) times larger in the case without mutual correlation after averaging the results of 100 random realizations. This supports that the observed retrieval bias from synthetic data is only an artifact of the pessimistically simplified simulation and should not be present in the results obtained from real data.

The formal errors of the complete dataset, which indicate the precision of the measurement, always remain smaller than the RMSE obtained from the synthetic data set. Due to the simplified simulation of small-scale topography and neglect of other potential error sources, the RMSE still provides an overly optimistic measure of accuracy of the results obtained from real data. However, it gives an indication of the dependence of the \(h_2\) accuracy on the resolution of the topographic grid.

5.4 Final \(h_2\) result

We use the inverse square of the RMSE as weights to compute a weighted mean and weighted standard deviation of the \(h_2\) values obtained from the complete data set for each resolution from 5 to 31 grid points per degree, resulting in \(h_2 = 0.0387\pm 0.0016\), where the error bar indicates one standard deviation. This value captures the accuracy indicated by the scattering \(h_2\) results for different resolutions. We add to this weighted standard deviation the standard deviation of the bias obtained from the synthetic data set, 0.0009, in order to incorporate the additional uncertainty this implies. This results in a final \(h_2 = 0.0387\pm 0.0025\). This result is robust with respect to the chosen topography resolutions. For example, considering resolutions from 11 to 31 grid points per degree leads to a final \(h_2 = 0.0387\pm 0.0021\).

6 Discussion

Figure 1 compares the result of this study with results from previous studies. The result obtained in this work is in agreement with a previous study using LOLA data and the crossover method (Mazarico et al. 2014). This demonstrates that the method presented in this study is capable of reproducing results achieved using the same data set.

However, neither of the LOLA-based values agree with the result from lunar interior structure models constrained by the Love number \(k_2\) (\(h_2 = 0.0424\), Williams et al. 2014) within their respective error bars. On the one hand, this could indicate deficiencies in the LOLA data or their processing. Unaccounted perturbations due to thermal or instrumental effects acting at the tidal frequency of \(27.2\,\mathrm {d}\) may cause biases in the results. While the method presented in this study excels at filtering out random noise due to the large amount of measurements, such systematic errors would bias the result. On the other hand, the discrepancy could also indicate problems with the GRAIL measurements of \(k_2\) or the modeling.

LLR as a different range measurement technique cannot alleviate this problem either. Three different studies have used different sets of LLR data, ephemerides, and tidal models to retrieve \(h_2\) (Williams et al. 2013; Pavlov et al. 2016; Viswanathan et al. 2018), but only the result of Williams et al. (2013) with its large error bar and the result of Pavlov et al. (2016) using the Earth tidal model of Williams et al. (2013) agree with the modeling result based on \(k_2\) (Williams et al. 2014). In total, the differences between the various \(h_2\) results from LLR are larger than the uncertainty of \(h_2\) retrieved from LOLA data. Viswanathan et al. (2018) needed to treat low-degree coefficients of the lunar gravity field as free parameters in order to be able to fit a lunar interior model to their LLR measurements. This may also indicate a low reliability of the \(k_2\) result by GRAIL, which could potentially explain the discrepancy between the \(h_2\) result of this study and the modeled result of Williams et al. (2014).

Laser altimetry is more sensitive than LLR to the radial tidal displacement itself (Mazarico et al. 2014). One possible issue with the retrieval of \(h_2\) from LLR is that all measurements are taken on the near side of the Moon. If the hemispheric dichotomy of the Moon’s surface (e.g., Wieczorek et al. 2006) extends into the elastic properties in its interior, a Love number determined using only near-side measurements could be biased (Zhong et al. 2012). LLR determination of solid body tides may also be biased by thermal expansion of retroreflectors on the Moon and the regolith which acts on monthly frequencies and has an amplitude in the millimeter range (Williams and Boggs 2015).

The crossover method employed by Mazarico et al. (2014) has previously been the only method to successfully determine \(h_2\) of any planetary body from laser altimetry and has also found wide application in simulations (Mazarico et al. 2010, 2015; Steinbrügge et al. 2015, 2018c; Hussmann et al. 2016). However, a disadvantage of that method is that in a near-polar orbit, crossovers occur mainly at high latitudes where the tidal signal is weaker than at the equator. Both on the Moon and on Mercury, the tidal displacements at the poles are up to \(\sim 70\%\) weaker than at low latitudes. Furthermore, if a single-beam laser altimeter were to be used, the grazing angles at which the few crossovers at lower latitudes occur would be unfavorable for an accurate retrieval (Koch et al. 2008; Mazarico et al. 2014). Fortunately, this limitation does not affect multi-beam laser altimeters like LOLA. However, the upcoming BepiColombo Laser Altimeter (BELA) at Mercury (Thomas et al. 2007; Hussmann et al. 2018; Thomas et al. 2019) and Ganymede Laser Altimeter (GALA, Kimura et al. 2019; Hussmann et al. 2019) will be single-beam laser altimeters. Furthermore, the amount of crossovers of these missions will be limited by larger orbital periods and short mission durations due to hostile radiation environments. In the case of BELA, another limitation is the lower pulse repetition rate of only \(10\,\mathrm {Hz}\) (Kallenbach et al. 2013). Using the reduced data set, we have achieved a result that differs by less than 1% when omitting 83.3% of the data, which is promising with respect to future application of the presented method to experiments which provide smaller data sets than LOLA.

The resolution of the topography achieved by an expansion in spherical harmonics is significantly lower than the resolution achieved here. While Koch et al. (2008) and Steinbrügge et al. (2019) solved to spherical harmonic degree \(L = 64\) and \(L = 60\), equivalent to a lunar topographic resolution of 85 or 91 km, respectively, the splines used in this study can resolve topography down to scales of 0.98 km. However, the crossover method interpolates topography over the even smaller scale of the shot-to-shot distance, which is approximately 25 m (Mazarico et al. 2010).

For the solution, we had to make the assumptions that the error term \(e_k\), which is dominated by the unmodeled small-scale topography, is uncorrelated and has equal variance. In practice, these assumptions are violated because topography at the typical spot-to-spot distance of \(\sim 25\) m is correlated and the variance is expected to be significantly lower in regions with smoother topography because the spline model can better fit smooth topography. However, when generating synthetic measurements, we take the correlation of small-scale topography into account, and still obtain a high retrieval accuracy for the \(h_2\) results. Small-scale topography down to a wavelength of \(\sim 680\) m is synthetically generated as correlated topography by the spherical harmonic model with degree \(L = 7999\). Gaussian process regression models the autocorrelation of topography at even smaller scales. The bias obtained when evaluating synthetic measurements is likely due to the incomplete consideration of the correlation of small-scale topography, which is a pessimistic assumption. Since neglecting the autocorrelation of topography in the solution produces satisfactory results using synthetic data, the same should be the case using real data.

7 Conclusions

Koch et al. (2008, 2010) developed a method for the retrieval of \(h_2\) from laser altimetry data which does not use crossovers directly, but instead solves simultaneously for the global topography and \(h_2\) in a least-squares adjustment. This study advances this method further by implementing 2D splines as basis functions and making the solution strategy more robust. Koch et al. (2008, 2010) tested their method only with synthetic data. Here, we apply it for the first time to actual data, choosing the enormous dataset produced by LOLA. The result of \(h_2 = 0.0387\pm 0.0025\) agrees with a previous result from the same data set using a crossover method within its standard deviation (Mazarico et al. 2014), but is approximately 10% smaller than the probably most reliable value based on a lunar structure model that satisfies the observed value of the \(k_2\) Love number, which is known with an estimated accuracy of 1% (Williams et al. 2014). Nonetheless, the results in this work suggest that the method is capable of retrieving the \(h_2\) Love number of other Solar System objects with much larger tidal amplitudes than the Moon, such as Mercury and Ganymede, by laser altimeters to sufficient accuracy (Steinbrügge et al. 2019). Further studies must examine the influence of systematic orbit and pointing errors on the retrieved \(h_2\) and should target the retrieval of additional geophysical parameters such as the tidal lag, forced libration amplitudes, or regionally varying elastic properties.