Efficient Forward Modeling of Gravitational Fields in Spherical Harmonic Domain with Application to Lunar Topography

1. Introduction

Forward modeling of gravitational fields plays an important role in data processing and geological interpretation in gravity exploration. It has been wildly used to calculate the gravity responses of topography, sedimentary layer, crystalline layer and Moho interface [1,2,3,4,5]. Then these effects can be removed from the free-air gravity anomaly to get the residual gravity anomaly that only produced by the inhomogeneity of the upper mantle of the lithosphere [6,7,8,9]. In addition, gravity forward modeling, as the core of iteration of 3D density inversion, is typically employed in mineral resources exploration [10,11,12] and lithospheric density structures of planetary interiors [13,14,15]. However, traditional gravity forward methods are mostly based on Cartesian coordinate system [16,17]. In recent years, with the development of satellite gravity, more and more high-precision and resolution gravity field models have been applied to the study of 3D density structures and evolution history of the Earth [18,19], the Moon [20,21,22] and Mars [23,24]. Therefore, gravity forward methods and relative methods of data processing and inversion should be changed from Cartesian coordinate system to spherical coordinates.

Gravity forward modeling in spherical coordinates can be carried out in spatial domain and spherical harmonic domain (i.e., frequency domain). In spatial domain, we usually divide the study region into some spherical prisms, i.e., tesseroids [25]. The gravity responses of these tesseroids can be obtained by the accumulation of gravity signal generated by each tesseroid. However, the 3D Newton’s integral of the tesseroid is usually hard to be solved analytically except for the case that the source regions are sphere and sphere shell with constant density [26]. Therefore, numerical integration methods are alternatively more common to tackle this issue, e.g., the Taylor series expansion [27,28] and the 2D/3D Gauss-Legendre quadrature (GLQ) method [29,30,31]. However, these numerical methods always lead to absolutely wrong results when the observation points are close to tesseroids. Uieda et al. [32] proposed an adaptive discretization method that has greatly increased the computational precision ensuring the relative error less than 0.1%. In addition, the computational efficiency will be another severe issue for large-scale or global gravity forward and inversion especially when the discretization is fine. Recently, some strategies, e.g., the kernel matrix equivalent storage [4,33] and the fast kernel-vector multiplication [5,34] have been proposed which decreased the computing time by three orders of magnitude compared with traditional method. Though the issues on computational efficiency and accuracy of these methods in spatial domain have been well settled, global and regional gravitational fields are now more commonly released in the form of spherical harmonic coefficients which makes it difficult for these spatial domain methods to simulate the gravity responses with specified degrees and orders. However, this will be very simple for the spherical harmonic domain methods.

Spherical harmonics are special functions defined on the surface of a sphere that form an orthogonal basis set for square-integrable functions, by which non-singular spectra expressions for gravity anomalies and gradients are derived based on the sums of spherical harmonic coefficients of gravitational field models [35,36]. Similar to the Parker’s formula [37] and the interface iterative inversion method [38] in Cartesian coordinate, Wieczorek and Phillips [39] derived the formula for calculating gravity anomaly of single undulating interface with constant density and the iterative inversion method of density interface in spherical harmonic domain, which were used to lunar topography correction and estimation of lunar crustal thickness [40,41]. This method is high-efficient and suitable for gravity correction using high-order and high-resolution terrain models, and gravity inversions of a single density interface. However, the binomial series expansion in this method has great effect to computational accuracy and convergence in inversion. Root et al. [42] proposed a correction to mitigate this erroneous behavior, by which it is applied to gravity forward modeling of Earth’s density structures. Instead of implementing the binomial expansion of the product of topographic relief and 2D density distribution in Wieczorek and Phillips [39], Šprlák et al. [43,44] present an innovative approach to construct the gravitational potential spectra in spherical coordinates using the discrete Fourier transform. In addition to topography corrections, these spherical harmonic domain methods have been also widely used in the calculation of gravity disturbances of lithosphere from crustal models with variable density distributions [2,8,45]. Although forward modeling methods of the gravity anomaly have been extensively studied and applied, studies on spherical harmonic spectra of gravitational gradient tensors remain comparatively rare. Besides, there are few studies on gravity forward and inversion modeling of 3D density anomalies by using the spherical harmonic methods [46].

To adapt to the requirement of gravity forward and inversion of 3D density anomalies with arbitrary degrees and orders, this study proposes an efficient forward modeling method of gravitational field in spherical harmonic domain, in which the gravity anomalies and gradient components are expressed as a form of spherical harmonic synthesis of spherical harmonic coefficients of arbitrary 3D density distributions. Two synthetic models demonstrate that the proposed method has improved the computational efficiency by four orders of magnitude and with a similar level of computational accuracy compared with the optimized 3D GLQ method [4,5]. Finally, we apply this method to calculated the topography correction and Bouguer gravity anomaly of the Moon.

2. Theory and Methods

2.1. Spatial Domain Forward Modeling Method

In spherical coordinates, the global spherical shell is typically discretized into tesseroids [25], and observation points are on a spherical surface aligning with the tesseroids mesh, as shown in Figure 1.

Here we define a local north-oriented frame (LNOF) where z-axis points to upward in geocentric radial direction, x-axis towards the north and y-axis directs to the east. We suppose that P (r, φ, λ) is an arbitrary observation point and Q (rʹ, φʹ, λʹ) is the center position of an arbitrary tesseroid bounded by (r₁, r₂), (φ₁, φ₂), and (λ₁, λ₂) in three directions, where r, φ and λ (also for rʹ, φʹ and λʹ) are geocentric radius, latitude, and longitude, respectively. Therefore, the gravitational potential on point P by tesseroid Q with homogeneous density ρ can be written as Newton’s integral [28,31,32]:

$$V\left(r,\phi ,\lambda \right)=G\rho {\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}{\displaystyle \underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}\frac{{r^{\prime }}^2\cos {\phi }^{\prime }}{\ell }d{r}^{\prime }d{\phi }^{\prime }d{\lambda }^{\prime }}}}$$

(1)

where G is the gravitational constant; l is the distance between P and Q:

$$\begin{array}{l}\ell =\sqrt{{r^{\prime }}^2+r^2-2r^{\prime }r\cos \psi }\\ \cos \psi =\sin \phi \sin {\phi }^{\prime }+\cos \phi \cos {\phi }^{\prime }\cos \left({\lambda }^{\prime }-\lambda \right)\end{array}$$

(2)

According to the coordinate transformation relationship between the geocentric spherical coordinates and the local rectangular coordinates [27,31], it is easy to obtain their partial derivatives of V to r, φ and λ, then the gravitational accelerations g_α and gradient tensors g_αβ can be obtained by their combinations [47,48], where α, β∈(x, y, z):

$$\left[g_{\alpha }\right]={\left[\begin{array}{ccc}{g}_x& g_y& g_z\end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial V}{\partial \phi }}& {\displaystyle \frac{1}{r\cos \phi }}{\displaystyle \frac{\partial V}{\partial \lambda }}& {\displaystyle \frac{\partial V}{\partial r}}\end{array}\right]}^{\mathrm{T}}$$

(3)

$$\left[g_{\alpha \beta }\right]=\left[\begin{array}{ccc}{\displaystyle \frac{1}{r^2}}\left({\displaystyle \frac{{\partial }^{2}V}{\partial {\phi }^2}}+r{\displaystyle \frac{\partial V}{\partial r}}\right)& {\displaystyle \frac{1}{r^2\cos \phi }}\left({\displaystyle \frac{{\partial }^{2}V}{\partial \phi \partial \lambda }}+\tan \phi {\displaystyle \frac{\partial V}{\partial \lambda }}\right)& {\displaystyle \frac{1}{r}}\left({\displaystyle \frac{{\partial }^{2}V}{\partial \phi \partial r}}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial V}{\partial \phi }}\right)\\ {\displaystyle \frac{1}{r^2\cos \phi }}\left({\displaystyle \frac{{\partial }^{2}V}{\partial \phi \partial \lambda }}+\tan \phi {\displaystyle \frac{\partial V}{\partial \lambda }}\right)& {\displaystyle \frac{1}{r^2{\cos }^2\phi }}\left({\displaystyle \frac{{\partial }^{2}V}{{\partial }^2\lambda }}+r{\cos }^2\phi {\displaystyle \frac{\partial V}{\partial r}}-\cos \phi \sin \phi {\displaystyle \frac{\partial V}{\partial \phi }}\right)& {\displaystyle \frac{1}{r\cos \phi }}\left({\displaystyle \frac{{\partial }^{2}V}{\partial r\partial \lambda }}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial V}{\partial \lambda }}\right)\\ {\displaystyle \frac{1}{r}}\left({\displaystyle \frac{{\partial }^{2}V}{\partial \phi \partial r}}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial V}{\partial \phi }}\right)& {\displaystyle \frac{1}{r\cos \phi }}\left({\displaystyle \frac{{\partial }^{2}V}{\partial r\partial \lambda }}-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial V}{\partial \lambda }}\right)& {\displaystyle \frac{{\partial }^{2}V}{{\partial }^{2}r}}\end{array}\right]$$

(4)

The gravitational potential, three components of anomalies, and gradient tensors can be further rearranged as [32]:

$$V\left(r,\phi ,\lambda \right)=G\rho {\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}{\displaystyle \underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}\frac{1}{\ell }\kappa \cdot d{r}^{\prime }d{\phi }^{\prime }d{\lambda }^{\prime }}}},$$

(5)

$$g_{\alpha }\left(r,\phi ,\lambda \right)=G\rho {\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}{\displaystyle \underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}\frac{{\Delta }_{\alpha }}{{\ell }^3}\kappa \cdot d{r}^{\prime }d{\phi }^{\prime }d{\lambda }^{\prime }}}},$$

(6)

$$g_{\alpha \beta }\left(r,\phi ,\lambda \right)=G\rho {\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{{\phi }_1}{\overset{{\phi }_2}{\int }}{\displaystyle \underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}\left(\frac{3{\Delta }_{\alpha }{\Delta }_{\beta }}{{\ell }^5}-\frac{{\delta }_{\alpha \beta }}{{\ell }^3}\right)\kappa \cdot d{r}^{\prime }d{\phi }^{\prime }d{\lambda }^{\prime }}}},$$

(7)

where δ_αβ is Kronecker’s delta (δ_αβ = 1, if α = β and δ_αβ = 0, if α ≠ β), and

$$\begin{array}{l}{\Delta }_x=r^{\prime }\left[\cos \phi \sin {\phi }^{\prime }-\sin \phi \cos {\phi }^{\prime }\cos \left({\lambda }^{\prime }-\lambda \right)\right] \\ {\Delta }_y=r^{\prime }\cos \phi \sin \left({\lambda }^{\prime }-\lambda \right)\\ {\Delta }_z=r^{\prime }\cos \psi -r\\ \kappa ={r^{\prime }}^2\cos {\phi }^{\prime }\end{array}$$

(8)

There is no analytical solution of Equations (5) – (7) since they contain elliptical integrals about φ' and λ'. Therefore, we should adopt numerical integration methods to solve it. Since the 3D GLQ method is numerically better than the Taylor series expansion method [31], we choose the 3D GLQ to represent the spatial domain method. For simple, taking one-dimensional GLQ as an example, a line integral of continuous function f(x) in [a, b] interval can be expressed as:

$${\displaystyle \underset{a}{\overset{b}{\int }}f\left(x\right)}dx\approx {\displaystyle \frac{b-a}{2}}{\displaystyle \sum _{k=0}^N{W}_k}f\left({\displaystyle \frac{b-a}{2}}{x}_k+{\displaystyle \frac{a+b}{2}}\right)$$

(9)

where N is the number of Gaussian nodes; x_k and W_k are the k-th Gaussian node and Gaussian coefficient in [1,1] interval, respectively, where the commonly used Gaussian nodes and coefficients are given by Wild-Pfeiffer [31]. In the following comparison, we use two Gaussian nodes to fit the trade-off between computational accuracy and efficiency.

We note that if we simply use this 3D GLQ method to gravity forward modeling of tesseroids, its computational efficiency and accuracy will be very low, especially when observation points are close to the top of tesseroids and when the source region is discretized finely. Therefore, we will employ the optimal methods of adaptive discretization by Uieda et al. [32] and kernel matrix equivalent storage by Zhao et al. [5] for comparison with the new proposed method in the next section.

2.2. Forward Method in Spherical Harmonic Domain

The gravitational field outside the Earth satisfies Laplace equation, so it can also be expressed as spherical harmonic expansion. Based on this theory, here we firstly derive the equations of the gravitational field of a single-layer spherical shell model in the form of spherical harmonic expressions, and then the expressions of 3D source region are obtained by accumulation of each layer. The gravitational potential in Equation (1) can be reshaped as:

$$V(r,\phi ,\lambda )=G{\displaystyle \underset{S}{\iint }{\displaystyle \underset{r^{\prime }}{\int }\frac{\rho (r^{\prime },{\phi }^{\prime },{\lambda }^{\prime })}{\ell (r,\phi ,\lambda ;r^{\prime },{\phi }^{\prime },{\lambda }^{\prime })}}}dSd{r}^{\prime }$$

(10)

where S denotes the range of spherical integral.

The gravitational potential V and 3D density distribution ρin Equation (10) can be written in the form of normalized spherical harmonic functions according to the spherical harmonic expansion theory [49]:

$$V(r,\phi ,\lambda )={\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l\left[{}^1\overline{V}_l^m(r){\overline{R}}_l^m(\phi ,\lambda )+{}^2\overline{V}_l^m(r){\overline{T}}_l^m(\phi ,\lambda )\right]}}$$

(11)

$$\rho (r^{\prime },{\phi }^{\prime },{\lambda }^{\prime })={\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l\left[{}^1\overline{\rho }_l^m(r^{\prime }){\overline{R}}_l^m({\phi }^{\prime },{\lambda }^{\prime })+{}^2\overline{\rho }_l^m(r^{\prime }){\overline{T}}_l^m({\phi }^{\prime },{\lambda }^{\prime })\right]}}$$

(12)

where ${}^1\overline{V}_l^m(r)$ and ${}^2\overline{V}_l^m(r)$ are fully normalized spherical harmonic coefficients of the gravitational potential with degree l and order m on a sphere surface with radius r; ${}^1\overline{\rho }_l^m(r^{\prime })$ and ${}^2\overline{\rho }_l^m(r^{\prime })$ are fully normalized spherical harmonic coefficients of the 2D density distribution on a layer with radius rʹ; L is the maximum spherical harmonic degree; ${\overline{R}}_l^m(\phi ,\lambda )$ and ${\overline{T}}_l^m(\phi ,\lambda )$ are fully normalized spherical harmonic functions which can be written as [49]:

$$\begin{array}{l}{\overline{R}}_l^m(\phi ,\lambda )={\overline{P}}_l^m(\sin \phi )\cos m\lambda \\ {\overline{T}}_l^m(\phi ,\lambda )={\overline{P}}_l^m(\sin \phi )sinm\lambda \end{array}$$

(13)

where ${\overline{P}}_l^m(\sin \phi )$ is the 4π fully-normalized associated Legendre function.

The Green function (1/l) in Equation (10) is the key to spherical harmonic domain methods which can be expressed as [49]:

$${\displaystyle \frac{1}{\ell }}={\displaystyle \frac{1}{r}}{{\displaystyle \sum _l^{\infty }\left(\frac{r^{\prime }}{r}\right)}}^l{\displaystyle \frac{1}{2l+1}}{\displaystyle \sum _{m=0}^l\left[{\overline{R}}_l^m(\phi ,\lambda ){\overline{R}}_l^m({\phi }^{\prime },{\lambda }^{\prime })+{\overline{T}}_l^m(\phi ,\lambda ){\overline{T}}_l^m({\phi }^{\prime },{\lambda }^{\prime })\right]}, r>r^{\prime }$$

(14)

We should note that the infinite series in Equation (14) only converge uniformly for any observation points that are r > rʹ, i.e., the observation points need to be outside of the source region. Substituting Equations (11), (12) and (14) into (10), and employing the orthogonal characteristic of spherical harmonic functions, we get:

$$\begin{array}{l}{}^1\overline{V}_l^m(r)=G{\displaystyle \underset{r^{\prime }}{\int }\frac{4\pi r^{\prime }}{2l+1}}{\left({\displaystyle \frac{r^{\prime }}{r}}\right)}^{l+1}{}^1\overline{\rho }_l^m(r^{\prime })d{r}^{\prime }\\ {}^2\overline{V}_l^m(r)=G{\displaystyle \underset{r^{\prime }}{\int }\frac{4\pi r^{\prime }}{2l+1}}{\left({\displaystyle \frac{r^{\prime }}{r}}\right)}^{l+1}{}^2\overline{\rho }_l^m(r^{\prime })d{r}^{\prime }\end{array}$$

(15)

Equation (15) shows that spherical harmonic coefficients of gravitational potential are calculated according to the radial integration of spherical harmonic coefficients of 3D density distribution, which is the core of gravity forward modeling in spherical harmonic domain. By substituting Equation (15) into (10), we can get spherical harmonic expressions of gravitational potential on a sphere with a given radius r outside the source region:

$$V(r,\phi ,\lambda )=4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{r^{\prime }}{2l+1}{\left(\frac{r^{\prime }}{r}\right)}^{l+1}}\left[{}^1\overline{\rho }_l^m(r^{\prime }){\overline{R}}_l^m(\phi ,\lambda )+{}^2\overline{\rho }_l^m(r^{\prime }){\overline{T}}_l^m(\phi ,\lambda )\right]}}d{r}^{\prime }$$

(16)

According to Equations. (3) and (4), we can obtain the spherical harmonic expressions of gravitational acceleration and gradient tensors:

$$g_x(r,\phi ,\lambda )=4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{1}{2l+1}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}\left[{}^1\overline{\rho }_l^m(r^{\prime })\cos m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })sinm\lambda \right]}}{\displaystyle \frac{\partial {\overline{P}}_l^m(\sin \phi )}{\partial \phi }}d{r}^{\prime }$$

(17)

$$g_y(r,\phi ,\lambda ) =4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{1}{\left(2l+1\right)\cos \phi }{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}m\left[-{}^1\overline{\rho }_l^m(r^{\prime })\sin m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })\cos m\lambda \right]}}{\overline{P}}_l^m(\sin \phi )d{r}^{\prime }$$

(18)

$$g_z(r,\phi ,\lambda )=-4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{l+1}{2l+1}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}\left[{}^1\overline{\rho }_l^m(r^{\prime }){\overline{R}}_l^m(\phi ,\lambda )+{}^2\overline{\rho }_l^m(r^{\prime }){\overline{T}}_l^m(\phi ,\lambda )\right]}}d{r}^{\prime }$$

(19)

$$\begin{array}{l}{g}_{xx}(r,\phi ,\lambda )=4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{1}{r\left(2l+1\right)}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}\left[{}^1\overline{\rho }_l^m(r^{\prime })\cos m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })sinm\lambda \right]}}\cdot \\ \left[{\displaystyle \frac{{\partial }^2{\overline{P}}_l^m(\sin \phi )}{\partial {\phi }^2}}-\left(l+1\right){\overline{P}}_l^m(\sin \phi )\right]d{r}^{\prime }\end{array}$$

(20)

$$\begin{array}{l}{g}_{xy}(r,\phi ,\lambda )=4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{1}{r\left(2l+1\right)}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}m\left[-{}^1\overline{\rho }_l^m(r^{\prime })\sin m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })\cos m\lambda \right]}}\cdot \\ \left[{\displaystyle \frac{\partial {\overline{P}}_l^m(\sin \phi )}{\cos \phi \partial \phi }}+{\displaystyle \frac{\tan \phi }{\cos \phi }}{\overline{P}}_l^m(\sin \phi )\right]d{r}^{\prime }\end{array}$$

(21)

$$g_{xz}(r,\phi ,\lambda )=-4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{l+2}{r\left(2l+1\right)}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}\left[{}^1\overline{\rho }_l^m(r^{\prime })\cos m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })sinm\lambda \right]}}{\displaystyle \frac{\partial {\overline{P}}_l^m(\sin \phi )}{\partial \phi }}d{r}^{\prime }$$

(22)

$$\begin{array}{l}{g}_{yy}(r,\phi ,\lambda )=-4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{1}{r\left(2l+1\right)}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}\left[{}^1\overline{\rho }_l^m(r^{\prime })\cos m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })sinm\lambda \right]}}\cdot \\ \left[m^2{\displaystyle \frac{{\overline{P}}_l^m(\sin \phi )}{{\cos }^2\phi }}+\left(l+1\right){\overline{P}}_l^m(\sin \phi )+{\displaystyle \frac{\tan \phi \partial {\overline{P}}_l^m(\sin \phi )}{\partial \phi }}\right]{\displaystyle \frac{\partial {\overline{P}}_l^m(\sin \phi )}{\partial \phi }}d{r}^{\prime }\end{array}$$

(23)

$$g_{yz}(r,\phi ,\lambda )=-4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{l+2}{r\left(2l+1\right)}{\left(\frac{r^{\prime }}{r}\right)}^{l+2}}m\left[-{}^1\overline{\rho }_l^m(r^{\prime })\sin m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })\cos m\lambda \right]}}{\displaystyle \frac{{\overline{P}}_l^m(\sin \phi )}{\cos \phi }}d{r}^{\prime }$$

(24)

$$g_{zz}(r,\phi ,\lambda )=4\pi G{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l{\displaystyle \underset{r^{\prime }}{\int }\frac{\left(l+1\right)\left(l+2\right)}{r^{\prime }\left(2l+1\right)}{\left(\frac{r^{\prime }}{r}\right)}^{l+3}}\left[{}^1\overline{\rho }_l^m(r^{\prime })\cos m\lambda +{}^2\overline{\rho }_l^m(r^{\prime })\sin m\lambda \right]}}{\overline{P}}_l^m(\sin \phi )d{r}^{\prime }$$

(25)

Equations (17) – (25) contain the singular terms of 1/cosφ and 1/cos²φ when φ is closing to the pole. In addition, some equations include the terms of derivatives of ${\overline{P}}_l^m(\sin \phi )$. These mathematical problems have been well solved by previous studies [35,36,50], so here we will not focus on it.

The above Equations (17) – (25) are only theoretical expressions in the form of continuous radial integral. It is hard to calculate spherical harmonic coefficients of continuous density distribution in actual forward modeling. Therefore, the 3D spherical shell model will be divided into several density layers in the depth direction, and each density layer will be approximated by numerical integration. Finally, the gravitational fields generated by all layers will be summed. Supposing that a spherical shell model with a radius interval of [r_a, r_b] is evenly divided into N_z layers with the thickness of each layer ∆r = (r_b − r_a)/N_z, the center radius of each layer of the spherical shell is:

$${r^{\prime }}_k=r_a+\left(k-0.5\right)\times \Delta r, k=1,2,\cdots ,N_z$$

(26)

Though the 3D source region has been divided into multi-layer spherical shells, ${}^1\overline{\rho }_l^m(r^{\prime })$ and ${}^2\overline{\rho }_l^m(r^{\prime })$ are still continuous functions of r'. Assuming that the radial subdivision interval ∆r is small enough, then the density distribution on the k-th layer can be seen as a constant equal to its center density. Therefore, the above equations can be numerically solved, taking the commonly used V, g_z and g_zz components as examples:

$$V(r,\phi ,\lambda )=4\pi G{\displaystyle \sum _{k=1}^{N_z}{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l\frac{\left[{}^1\overline{\rho }_l^m({r^{\prime }}_k){\overline{R}}_l^m(\phi ,\lambda )+{}^2\overline{\rho }_l^m({r^{\prime }}_k){\overline{T}}_l^m(\phi ,\lambda )\right]}{\left(2l+1\right)\left(l+3\right)r^{l+1}}\left[{\left({r^{\prime }}_k+\frac{1}{2}\Delta r\right)}^{l+3}-{\left({r^{\prime }}_k-\frac{1}{2}\Delta r\right)}^{l+3}\right]}}}$$

(27)

$$g_z(r,\phi ,\lambda )=-4\pi G{\displaystyle \sum _{k=1}^{N_z}{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l\begin{array}{l}\frac{\left(l+1\right)\left[{}^1\overline{\rho }_l^m({r^{\prime }}_k){\overline{R}}_l^m(\phi ,\lambda )+{}^2\overline{\rho }_l^m({r^{\prime }}_k){\overline{T}}_l^m(\phi ,\lambda )\right]}{\left(2l+1\right)\left(l+3\right)r^{l+2}}\times \\ \left[{\left({r^{\prime }}_k+\frac{1}{2}\Delta r\right)}^{l+3}-{\left({r^{\prime }}_k-\frac{1}{2}\Delta r\right)}^{l+3}\right]\end{array}}}}$$

(28)

$$g_{zz}(r,\phi ,\lambda ) =4\pi G{\displaystyle \sum _{k=1}^{N_z}{\displaystyle \sum _{l=0}^L{\displaystyle \sum _{m=0}^l\begin{array}{l}\frac{\left(l+1\right)\left(l+2\right)\left[{}^1\overline{\rho }_l^m({r^{\prime }}_k){\overline{R}}_l^m(\phi ,\lambda )+{}^2\overline{\rho }_l^m({r^{\prime }}_k){\overline{T}}_l^m(\phi ,\lambda )\right]}{\left(2l+1\right)\left(l+3\right)r^{l+3}}\times \\ \left[{\left({r^{\prime }}_k+\frac{1}{2}\Delta r\right)}^{l+3}-{\left({r^{\prime }}_k-\frac{1}{2}\Delta r\right)}^{l+3}\right]\end{array}}}}$$

(29)

According to Equations. (27) – (29), the steps of gravitational forward modeling of 3D spherical shell model in spherical harmonic domain are clear. We firstly expand density distributions of each layer into spherical harmonic coefficients which are multiplied by the coefficients of corresponding degree and order. Next, gravitational responses of each layer can be obtained by spherical harmonic synthesis of the result in the first step. Finally, these gravity effects of each layer are summed. For more specificity, Table 1 shows detailed procedures of our method with g_z component as example.

Similar to the procedures of gravitational forward modeling in Cartesian coordinate system that based on FFT [51,52], the forward modeling in spherical harmonic domain need also the spherical harmonic forward and inverse transformations. Therefore, the computational efficiency of spherical harmonic transformations dominates the efficiency of spherical harmonic based methods. Wieczorek and Meschede [53] released a software package SHTOOLS (https://shtools.github.io/SHTOOLS/), in which the Fourier transform software FFTW [54] and the OpenMP parallel program under the “DUCC” library (https://github.com/litebird/ducc) are employed by default to greatly improved the computational efficiency of spherical harmonic expansion and synthesis. In this study, 3D gravity forward modeling in spherical harmonic domain is performed under the SHTOOLS software system as expressed in the steps 2 and 4 in Table 1.

3. Synthetic Forward Model Tests

3.1. Sphere Shell Model

Here we firstly design a simple spherical shell model to testify the efficiency of our method compared with 3D GLQ methods. A homogeneous spherical shell with density of 500 kg/m³ and thickness of 100 km ranging from 1638 to 1738 km, as shown in Figure 2, is divided into small tesseroids with varying intervals ranging from 0.25° to 20° and 1 layer in depth. The gravitational fields are calculated and compared at 10 km height above the top surface of the shell model with observation points aligning with the tesseroid mesh. The analytical solution of the gravitational field of the sphere shell is calculated by the difference of gravity responses of two spheres with different radii [27].

We calculate the g_z and g_zz at different mesh intervals (i.e., different mesh scales) using adaptive discretization method by Uieda et al. [32] and kernel matrix equivalent storage method by Zhao et al. [5], and this proposed spherical harmonic domain method, with their relative computational time shown in Figure 3. The absolute errors of g_z and g_zz by the spatial domain method and our proposed spherical harmonic domain method are shown in Figure 4. In addition, we also show the comparison of computational time and the maximum relative error of g_z and g_zz when the mesh interval is 0.25° (i.e., the discretization is 720 × 1440 × 1) in Table 2.

Figure 3 shows that the computational time increases rapidly with the decreases of mesh interval (i.e., with the increases of total number of tesseroids). Specifically, the computational efficiency by the proposed spherical harmonic domain method is increased by 3 and 6 orders of magnitude compared with the methods of Zhao et al. [5] and of Uieda et al. [32], respectively, for the g_z component when the mesh interval is less than 3° (Figure 3a). As for the g_zz component, the improvements of computational efficiency are more obvious, which are 4 orders and 7 orders of magnitude (Figure 3b), respectively. These improvements on the computational efficiency are quantified in Table 2, where the absolute computational times by the spatial domain methods are extremely long, while the computational time by the proposed spherical harmonic domain method is acceptable and suitable to global 3D gravity inversion.

As for computational accuracy, Figure 4 shows absolute errors of g_z and g_zz by both the spatial and spherical harmonic domain methods vary with latitude. Also, it is obviously that absolute errors by the proposed spherical harmonic domain method (Figure 4b and d) are lower than that by the spatial domain method (Figure 4a and c). Specifically, the maximum relative error of the g_z component by the proposed spherical harmonic domain method is decreased by 3 orders of magnitude compared with that by the spatial domain method, as shown in Table 2, while the improvement on the accuracy for the g_zz component is about two orders of magnitude.

3.2. Complex Synthetic Model

Here we design a more complex synthetic model to further test the effectiveness of the proposed method. The model is a spherical shell with a thickness of 100 km ranging from 1638 km to 1738 km, which is evenly discretized into 360 × 720 × 10 tesseroids (i.e., with a mesh interval of 0.5° in horizontal direction and ten layers in depths) with their density anomalies varying as:

$$\rho \left(r,\phi ,\lambda \right)=\left(1738-r\right)\cdot \left[\sin \left(2\phi \right)\cdot \sin \left(2\lambda \right)\right]$$

(30)

The density distribution of this model on a layer with its radius of r = 1682.5 km is shown in Figure 5. The gravitational fields are calculated on a spherical surface at 10 km above the top surface of this model. Figure 6 shows the g_z and g_zz of this model by the spatial domain method, and the differences between it and the proposed spherical harmonic domain method.

Although it is difficult for this complex model to have an analytical solution, the results by the spatial domain method with adaptive discretization strategy [32] can be regarded as standards for comparison. It can be seen in Figure 6 that the g_z and g_zz by the proposed spherical harmonic domain method are very close to (a) and (b) by the spatial domain method, with the biggest differences for g_z (c) and g_zz (d) between the two methods lower than 0.5 mGal and 5.0 E, respectively, which proves the effectiveness of our method.

4. Application to Lunar Topography Correction

Bouguer gravity anomaly is the comprehensive reflection of all density anomalies after removing the influence of topography from the free-air anomaly, which has been widely used in 3D gravity inversion to reveal the density structure and deformation mechanism of the lithosphere. In addition, Bouguer anomalies are also the basis for other gravity corrections, such as residual gravity anomalies used for Moho interface inversion by removing the gravity responses of sedimentary layers, crystalline crust layers and deep mantle effects from Bouguer anomalies. Here we employ the proposed method to calculate the gravity response of the lunar topography and then obtain the Bouguer gravity anomaly.

The latest lunar gravitational field model GL1500E [21] is used here to calculate the free-air gravity anomaly. This model is expressed in the form of spherical harmonic coefficients, where different degrees and orders represent gravity anomalies with different spatial resolutions [55]. In this study, the free-air gravity is obtained with truncated degrees of 450 equivalent to a horizontal resolution of 0.4° × 0.4° by using spherical harmonic synthesis methods [56]. Considering that the observation data are usually extended upwards to a certain height in inversions to reduce the influence of the noise in data [57], so we calculate the free-air gravity anomaly on a surface of r = 1,748 km which is 10 km height above the reference radius, as shown in Figure 7a.

Next, we evaluate gravity responses of the lunar topography. The lunar topography map (Figure 7b) is obtained by spherical harmonic synthesis of the model LRO_LTM05_2050 [58] using the same degrees of coefficients as the free-air gravity anomaly to keep consistency. Considering that our method is suitable for uniform layered shells which is inconsistent with the topography model, so it is necessary to discrete the terrain to plenty of tesseroids with constant thickness for each layer to reduce the error of fitting the terrain using tesseroids. Figure 7b shows the lunar topography ranges from −9.29 to 9.63 km, so we construct a spherical shell that slightly higher than the topography ranging from −10.00 to 10.00 km to cover the topography, and then this spherical shell is evenly divided into 400 thin layers with an interval of 0.05 km, as shown in Figure 8. Next, we need to judge the relationship between the depth of the geometric center of each tesseroid and the terrain. If the tesseroid is covered by the terrain, the correction density of this tesseroid is assigned as 2560 kg/m³ (for mountain) or −2560 kg/m³ (for basin), otherwise 0 kg/m³. The topography correction is also calculated on a surface of r = 1748 km consistent with the free-air gravity. Figure 9 shows the gravity effects of the lunar topography calculated by the spatial domain method, and the differences between it and the proposed method.

Figure 9 shows that the gravity effect of the lunar topography by the proposed method is similar to the result that by the spatial domain method which further proves the effectiveness and practicality of our method. The computational times of the topography corrections by our method and the spatial domain method are 88.69 s and 48789.98 s, respectively, indicating the high efficiency of the proposed method. However, the advantage of our method in terms of computational efficiency is slightly weakened in this application compared with that in the synthetic example in Section 3. This is because the premise of the proposed method is that all tesseroids of each layer have the same thickness. Therefore, when calculating terrain correction, the entire terrain needs to be divided into several thin layers (e.g., 400 layers in this example).

We remove gravity responses of topography from the free-air gravity to get the Bouguer gravity anomaly. Since the results of topography corrections by the spatial domain method and the proposed method are nearly close, here we just show Bouguer gravity anomaly that by the proposed method in Figure 10a. The Bouguer gravity anomaly in Figure 10a ranges from −500 mGal to 900 mGal, where the positive high anomalies are mainly distributed in most mascon basins and the South Pole-Aitken basin, while the negative anomalies are distributed on the Highland on the farside of the Moon. It can be found that the Bouguer gravity anomaly in Figure 10a is very close to the Bouguer result of Liang et al. [14], and the only difference between them is that the latter one used spherical harmonic coefficients with degrees and orders less than 60 to calculate the free-air gravity anomaly and the topography correction, with makes its resolution much lower than that of in this paper. However, the features of the Bouguer gravity anomaly in Figure 10a are very different from the results of Neumann et al. [56] and Zhao et al. [5]. It is obviously that the positive and negative Bouguer anomalies in Figure 10a display obvious asymmetries between nearside and farside of the Moon as well as between the northern and southern hemispheres. As for the results in Neumann et al. [56] and Zhao et al. [5] as shown in Figure 10b, positive high anomalies are significantly distributed beneath most mascon basins with less significant negative anomalies scattering around the positive anomalies (i.e., mascon basins). This is because they used 6−450 degrees and orders of spherical harmonic coefficients to calculate the free-air gravity anomaly and the topography correction, where the effect of the hemispheric asymmetry and the South Pole-Aitken impact has been naturally removed by discarding the signal with degrees and orders less than 6.

5. Conclusions

This study presented a high-efficient gravitational forward modeling method in spherical coordinate system based on spherical harmonic transform, in which the gravity anomalies and gradient tensors can be expressed as spherical harmonic synthesis form of spherical harmonic coefficients of 3D density distribution. Two synthetic forward examples demonstrated that the computational efficiency of the proposed method has improved by four orders of magnitude and with a similar level of computational accuracy compared with the optimized 3D GLQ method. The verified gravity forward method is then applied to calculate the topography correction and the Bouguer gravity anomaly of the Moon using the latest high-resolution gravitational field model and topography model. The result shows that the gravity effect of the lunar topography by the proposed method is similar to the result that by the spatial domain method. The computational efficiency of the topography correction is relatively lower than the synthetic forward examples because that the terrain should be evenly discretized into many thin layers to satisfy the requirement of the proposed forward method. In addition, the proposed method can calculate topography correction and Bouguer gravity anomaly with any desired degrees and orders, which is generally different for the spatial domain method. However, we should note that the proposed spherical harmonic domain method is only applicable to global gravity forward modeling, which may be analogized and applied to the spherical cap harmonic analysis for local gravity forward and inversion modeling.