Comments and corrections to: “Ellipsoidal spectral properties of the Earth’s gravitational potential and its first and second derivatives” by Bölling and Grafarend (2005) in J. Geod. 79(6-7):300-330

Abstract

In this contribution, we comment on and thoroughly identify errors which occur in the research article titled “Ellipsoidal spectral properties of the Earth’s gravitational potential and its first and second derivatives” by Bölling and Grafarend (2005) published in J.Geod. 79(6-7):300-330. The major issues affecting the main theoretical results by Bölling and Grafarend (J Geod 79(6-7):300-330, 2005) are as follows: (1) The tangential-tangential tensor spherical harmonics are not orthogonal for different harmonic degrees and the same harmonic orders. Consequently, the corresponding solution of the spherical gradiometric boundary value problem and the spherical integral formulas relating the tangential-tangential components of the incremental gravity gradient tensor are invalid. (2) The tangential-tangential, tangential-shear, tangential-dilation, tangential-normal, and normal-normal tensor ellipsoidal harmonics are not orthogonal for different harmonic degrees and the same harmonic orders. As a result, the solutions of the ellipsoidal gradiometric boundary value problems and the ellipsoidal integral formulas transforming the components of the incremental gravity gradient tensor among themselves are incorrect.

1 Introduction

Gravity field is an important property of any planetary body and represents one of the three main pillars of the modern geodesy (Plag et al. 2009). Countless articles have discussed theoretical, numerical, or practical matters of this phenomenon in the geodetic, geophysical, and planetary science literature. One of the studies on the gravity field is titled “Ellipsoidal spectral properties of the Earth’s gravitational potential and its first and second derivatives” by Bölling and Grafarend (2005) (herein referred to as BG05). This theoretical work introduces the incremental gravity potential (Sect. 2*), sets up the incremental gravity gradient and the incremental gravity gradient tensor (Sect. 3*), solves the gradiometric boundary value problem (BVP, Sect. 4*), and derives integral transforms of the incremental gravity potential onto the incremental gravity gradients and those transforming the incremental gravity gradients among themselves (Sect. 5*).

All these topics are presented in both the spherical approximation and its ellipsoidal (oblate spheroidal) counterpart. The spherical treatment employs the spherical coordinates, quantities are related to the spherical local reference frame, and the preferred mathematical parametrizations are the external spherical harmonic expansions. The ellipsoidal approximation makes use of the Jacobi ellipsoidal coordinates, quantities refer to the ellipsoidal local reference frame, and the external ellipsoidal harmonic series are the matching mathematical representation.

We have thoroughly studied the topics presented in BG05 and highly value the authors’ systematic approach. Nevertheless, we found fundamental mathematical issues that significantly affect the main theoretical outputs. In our article, we identify errors in BG05. In particular, we discuss two decompositions of the incremental gravity gradient tensor (Sect. 2), investigate orthogonality properties of the tensor ellipsoidal harmonics (Sect. 3), and summarize consequences for BG05 (Sect. 4). Appendices and the Electronic Supplementary Material (ESM) contain related (mathematical and formal) aspects.

All references from BG05, e.g., to sections, tables, or equations, are appended by an asterisk, while those without any sign can be found within our study. The notation and nomenclature are consistent with BG05 to the highest extent. Mathematical formulas in BG05 use the notation with metric coefficients. Here, we prefer explicit forms that are completely equivalent with those in BG05 after inserting the metric coefficients. A theoretical minimum for understanding our explanations and arguments can be found in Appendix A. Nevertheless, we strongly recommend the reader to scrutinize the paper by BG05 before reading our article.

2 Comments on the two decompositions of the incremental gravity gradient tensor

The incremental gravity gradient tensor is a central quantity in BG05. Its representation in the local spherical reference frame is introduced in Sect. 3.1*, see Eq. (A1). The tensor components are the boundary conditions when solving the spherical gradiometric BVP in Sect. 4.1*. Spherical integral formulas relating the incremental gravity potential to the tensor components including those relating the tensor components among themselves are derived in Sect. 5.1*.

An important step for both, the solution of the spherical gradiometric BVP in Sect. 4.1* and the derivation of the integral formulas in Sect. 5.1*, is an appropriate decomposition of the incremental gravity gradient tensor. A four-part decomposition is presented in Table 6*, see Eq. (A3), which divides the tensor into the Tangential-Shear (TS), Tangential-Dilation (TD), Tangential-Normal (TN), and Normal-Normal (NN) parts. Table 10*, on the other hand, contains a three-part decomposition of the incremental gravity gradient tensor into the Tangential-Tangential (TT), TN, and NN parts, see Eq. (A2). While the four-part decomposition is not further employed upon its introduction in BG05, the three-part version is the starting point for the solution of the gradiometric BVP in Sect. 4.1*.

The three-part decomposition, however, leads to an incorrect result for this reason. When solving the spherical gradiometric BVP, a crucial property of the tensor spherical harmonics is their orthogonality, i.e., that the integral from a multiplication of two tensor basis functions over the sphere vanishes for different harmonic degrees and orders. For the TT tensor spherical harmonics of the three-part decomposition, the orthogonality is taken for granted when solving the spherical TT gradiometric BVP in Eqs. (67*)-(74*) in Tab. 11*, but not asserted anywhere in BG05. Importantly, we prove in Appendix B that the TT tensor spherical harmonics are not orthogonal for all different harmonic degrees and orders.

Overall, the spherical gradiometric BVP should have been solved for the TS, TN, and NN parts of the four-part decomposition, as the corresponding tensor spherical harmonics are orthogonal by Eq. (35*) in Tab. 6*. The solution of the gradiometric BVP for the TD part is omitted, as it is identical to the one for the NN part by the Laplace equation, see (e.g., Martinec 2003).

The incremental gravity gradient tensor in the ellipsoidal local reference frame is also decomposed into three (TT, TN, NN) and four (TS, TD, TN, NN) parts, see Tabs. 12* and 9* and Eqs. (A5) and (A6). The three-part decomposition is the initial point for the solution of the ellipsoidal gradiometric BVP in Sect. 4.2*, but the four-part decomposition is not employed again. We could eventually expect similar disputes as already discussed for the spherical approximation. However, the ellipsoidal approximation suffers from even more serious errors, see below.

3 Orthogonality properties for the tensor ellipsoidal harmonics

In this section, we investigate orthogonality properties for the tensor ellipsoidal harmonics. To the best of our knowledge, a proof of such properties has not been presented in the literature. Therefore, a detailed treatment of this topic is provided here supplementing the general presentation in BG05.

3.1 Basic relations

The tensor ellipsoidal harmonics of degree k, order l, and type i are basis functions for the mathematical parametrization of the incremental gravity gradient tensor. These are denoted ${\textbf{T}}_{k l}^{i}$, $i \in \{\textrm{TT}, \textrm{TS}, \textrm{TD}, \textrm{TN}, \textrm{NN}\}$, and result from the application of the corresponding differential operators in Eqs. (A5) and (A6) to $e_{kl}(\varphi ,\lambda )$, i.e., the $4\pi$-normalized ellipsoidal harmonics of degree k and order l. BG05 provide explicit formulas for TS, TD, TN, and NN tensor ellipsoidal harmonics in Eqs. (C.22*)-(C.25*).

The orthogonality properties of the five tensor ellipsoidal harmonics can be concluded from the values of the integrals:

$$\begin{aligned}&\big \langle {\textbf{T}}_{k_1l_1}^{i}\big (\lambda , \varphi , u\big ), {\textbf{T}}_{k_2l_2}^{i}\big (\lambda , \varphi , u\big )\big \rangle \nonumber \\&\quad = \frac{1}{S} \int _{S} w(\varphi ,u)\ {\textbf{T}}_{k_1l_1}^{i}\big (\lambda , \varphi , u\big ) : {\textbf{T}}_{k_2l_2}^{i}\big (\lambda , \varphi , u\big )\ {\textrm{d}} S, \nonumber \\&\qquad i \in \{\textrm{TT}, \textrm{TS}, \textrm{TD}, \textrm{TN}, \textrm{NN}\}. \end{aligned}$$

(1)

The symbol  :  indicates the tensor product, and S is the area of the ellipsoid:

$$\begin{aligned} S= & {} S(u) \nonumber \\= & {} 4 \pi \sqrt{u^2 + \varepsilon ^2} \bigg [\frac{u^2}{4\varepsilon } \ln \frac{\sqrt{u^2 + \varepsilon ^2} + \varepsilon }{\sqrt{u^2 + \varepsilon ^2} - \varepsilon } + \frac{\sqrt{u^2 + \varepsilon ^2}}{2}\bigg ], \end{aligned}$$

(2)

and the infinitesimal surface area is as follows:

$$\begin{aligned} \textrm{d}S = \textrm{d}S(u,\varphi ) = \sqrt{u^2 + \varepsilon ^2}\, \sqrt{u^2 + \varepsilon ^2 \sin ^2 \varphi } \cos \varphi \ \textrm{d}\varphi \ \textrm{d}\lambda .\nonumber \\ \end{aligned}$$

(3)

The weight function w is defined by Eq. (39*) (for $u = b$) and is generalized for $u \ne b$ as:

$$\begin{aligned} w(\varphi ,u)= & {} \frac{1}{\sqrt{u^2 + \varepsilon ^2 \sin ^2 \varphi }}\nonumber \\{} & {} \times \bigg [\frac{u^2}{4\varepsilon } \ln \frac{\sqrt{u^2 + \varepsilon ^2} + \varepsilon }{\sqrt{u^2 + \varepsilon ^2} - \varepsilon } + \frac{\sqrt{u^2 + \varepsilon ^2}}{2}\bigg ]. \end{aligned}$$

(4)

Table 1 Numerical values of the integrals (1.6)–(1.10) in the ESM ($a = 1.0$ m, $b = 0.5$ m, $u = 0.6$ m)

Table 2 Numerical values of the integrals (1.6)–(1.10) in the ESM ($a = 6378137.0$ m, $b = 6356752.3141$ m, $u = 6606752.3141$ m)

Note that BG05 prefer employing the general bra-ket notation, i.e., the left-hand side of Eq. (1). The explicit integral representation after the first equality sign can be deduced from Eq. (86*) in Tab. 13*.

BG05 list orthogonality relationships for TS, TD, TN, and NN tensor ellipsoidal harmonics in Eqs. (61*) and (C.26*), which state that the integral (1) equals $\delta _{k_1k_2}\, \delta _{l_1l_2}$ for $i \in \{\textrm{TS}, \textrm{TD}, \textrm{TN}, \textrm{NN}\}$. An equivalent orthogonality relationship is assumed for the TT tensor ellipsoidal harmonics, when solving the TT ellipsoidal gradiometric BVP in Eqs. (83*)-(90*) in Tab. 13*, but is not declared anywhere in BG05.

For convenience, explicit formulas of Eq. (1) for the TT, TS, TD, TN, and NN tensor ellipsoidal harmonics are listed in the ESM.

3.2 Numerical proof of orthogonality properties

The integrals (1.6)–(1.10) in the ESM have to vanish $\forall k_1 \ne k_2$ and be nonzero if $k_1 = k_2$ for the tensor ellipsoidal harmonics being orthogonal, as assumed in BG05. An analytical proof of non-orthogonality of the tensor ellipsoidal harmonics for $k_1 = 2$, $k_2 = 0$, and $l = 0$ is presented in the ESM. However, for arbitrary non-negative integer harmonic degrees $k_1$ and $k_2$, and order l, use of analytical tools may be extremely challenging and cumbersome. Therefore, and for efficiency and simplicity of implementation, we preferred a numerical proof to its analytical counterpart.

A computer program for the numerical calculation of the integrals (1.6)–(1.10) in the ESM was coded in the programming language C (Kernighan and Ritchie 1988). We implemented a Gauss-Legendre quadrature algorithm for the numerical integration (Press et al. 2002). Numerical tests showed that the quadrature order 10,000 was the most optimal and thus adopted in our experiments. The associated Legendre functions of the first kind ${\bar{P}}_{kl}$ were calculated by the forward column recursion (e.g., Holmes and Featherstone 2002). For computation of the associated Legendre functions of the second kind $Q_{kl}$ and their derivatives with respect to u, we employed the representation in terms of the Gauss hypergeometric function ${}_{2} F_{1}$ (Hobson 1965, p. 108) and the routine gsl_sf_hyperg_2F1 implemented in the GNU Scientific Library (Galassi et al. 2013).

To reveal and eliminate any numerical issues or systematic errors in the C program, particular attention was devoted to independent checks by:

• (1) The closed formulas (2.1)–(2.5) in the ESM for the case $k_1 = 2$, $k_2 = 0$, and $l = 0$.

• (2) The closed formulas of integrals (1.7)–(1.10) in the ESM being $(n-1)n(n+1)(n+2) \delta _{k_1 k_2} \delta _{l_1 l_2}$, $(n+1)^2(n+2)^2 \delta _{k_1 k_2} \delta _{l_1 l_2}$, $2n(n+1)(n+2)^2 \delta _{k_1 k_2} \delta _{l_1 l_2}$, and $(n+1)^2(n+2)^2 \delta _{k_1 k_2} \delta _{l_1 l_2}$, respectively, when setting $a = b = u = 1.0$, i.e., $\varepsilon ^2 = 0$ (e.g., Rummel 1997; Martinec 2003).

• (3) A routine evaluating the identical integrals written in the mathematical software Mathcad 15 (Larsen 2010).

All three validation means confirmed correctness of the numerical implementation in the C program for different input parameters, see below.

We evaluated the integrals (1.6)–(1.10) in the ESM for the major semi-axis $a = 1.0$ m, minor semi-axis $b = 0.5$ m, and the size of the minor semi-axis of the computational point $u = 0.6$ m (i.e., above the ellipsoid). This oblate ellipsoid was chosen for convenience. Table 1 lists the numerical values for harmonic degrees $k_1$; $k_2 = k_1, k_1 - 2, k_1 - 4, \ldots , 0$, and orders l up to 5. These reach approximately three orders of magnitude from 10$^{0}$ up to 10$^{3}$ and are nonzero. These integrals were also calculated for the harmonic degrees $k_2 = k_1 - 1, k_1 - 3, k_1 - 5, \ldots , 1$ (not listed in Table 1) and their magnitudes ranged from 10$^{-16}$ up to 10$^{-13}$. As all computations were performed in the double precision, such values are numerically negligible and indicate that the integrals vanish in this case.

Integrals (1.6)–(1.10) in the ESM were further evaluated for the realistic GRS80 reference ellipsoid (Moritz 2000). The lengths of the major and minor semi-axes are $a = 6378137.0$ m and $b = 6356752.3141$ m. Numerical values of the five integrals 250 km above the surface of the GRS80 reference ellipsoid, i.e., $u = b + 250$ km $= 6606752.3141$ m, are listed in Table 2. Magnitudes are of the order of $1 / a^4$ for harmonic degrees $k_1$; $k_2 = k_1, k_1 - 2, k_1 - 4, \ldots , 0$ and at least 13 orders of magnitude smaller for $k_2 = k_1 - 1, k_1 - 3, k_1 - 5, \ldots , 1$, when they vanish (not listed in Table 2).

Other numerical experiments were performed for both ellipsoids by changing the size of the minor semi-axes u and did not change our conclusions.

We have thus demonstrated that the tensor ellipsoidal harmonics are orthogonal for $k_2 = k_1 - 1, k_1 - 3, k_1 - 5, \ldots , 1$, but non-orthogonal for $k_2 = k_1, k_1 - 2, k_1 - 4, \ldots , 0$. This is in contradiction with BG05, who stated that the tensor ellipsoidal harmonics are orthogonal $\forall k_1 \ne k_2$.

4 Conclusions and consequences for (Bölling and Grafarend 2005)

In this contribution, we identified errors in the paper by Bölling and Grafarend (2005). Their main focus is to solve the gradiometric boundary value problem and to derive integral formulas of the incremental gravity potential onto the incremental gravity gradients and those mutually relating the incremental gravity gradients. These mathematical problems are treated in the spherical and ellipsoidal approximations.

The major issue for the spherical approximation stems from the fact that Bölling and Grafarend (2005) employed the three-part decomposition of the incremental gravity gradient tensor when solving the spherical gradiometric BVP. Nevertheless, two TT tensor spherical harmonics of the three-part decomposition are non-orthogonal for degrees $k_1$; $k_2 = k_1, k_1 - 2, k_1 - 4, \ldots , 0$, and the same orders greater than 0, as we proved in Appendix B. Consequently, the solutions of the TT spherical gradiometric BVP in the spectral and spatial domains, see Eqs. (70*) and (74*) in Tab. 11*, are incorrect. The solution of the TT spherical gradiometric BVP is further used when finding the integral formula between the TT part of the incremental gravity gradient tensor in Eqs. (101*)-(105*) in Table 15*. These expressions are also invalid.

For the ellipsoidal approximation, the fundamental error in (Bölling and Grafarend 2005) originates from the assertion that two TT, TS, TD, TN, and NN tensor ellipsoidal harmonics are orthogonal for different degrees and orders, see Eqs. (61*) and (C.26*). However, we demonstrated numerically in Sect. 3 that this is not true for degrees $k_1$; $k_2 = k_1, k_1 - 2, k_1 - 4, \ldots , 0$, and the same orders. As the orthogonality properties do not hold, the solutions of the ellipsoidal gradiometric BVPs in Tab. 13* are incorrect starting from Eq. (86*). In addition, the solutions of the ellipsoidal gradiometric BVPs are considered when deriving the integral formulas between the components of the incremental gravity gradient tensor. The corresponding formulas (131*)-(140*) in Tab. 19* are thus incorrect.

Correct solutions of the spherical gradiometric BVP were found in (e.g., van Gelderen and Rummel 2001; Martinec 2003; Tóth 2013). Similarly, spherical integral formulas mutually relating the incremental gravity gradients were derived in (e.g., Tóth et al. 2005, 2006; Šprlák et al. 2014). On the other hand, the ellipsoidal gradiometric BVP and ellipsoidal integral formulas transforming incremental gravity gradients among themselves have not been presented in the literature, except for the attempt by Bölling and Grafarend (2005). As they were proved invalid, these still remain open problems of the geodetic theory and the potential theory in general.

Appendices

Appendix A Two decompositions of the incremental gravity gradient tensor

As a theoretical minimum, we summarize basic equations for the incremental gravity gradient tensor. We represent this quantity in both the spherical and ellipsoidal local reference frames by three- and four-part decompositions.

1.1 The incremental gravity gradient tensor in the spherical local reference frame

The incremental gravity gradient tensor results from the application of the gradient operator to the incremental gravity potential $\delta w$ twice:

$$\begin{aligned}&\textrm{grad} \otimes \textrm{grad}\ \delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\quad = \Bigg ({\textbf{e}}_{\lambda }\, \frac{1}{\sqrt{g_{\lambda \lambda }}}\, \textrm{D}_{\lambda } + {\textbf{e}}_{\varphi }\, \frac{1}{\sqrt{g_{\varphi \varphi }}}\, \textrm{D}_{\varphi } + {\textbf{e}}_{r}\, \frac{1}{\sqrt{g_{rr}}}\, \textrm{D}_{r}\Bigg ) \nonumber \\&\qquad \otimes \Bigg ({\textbf{e}}_{\lambda }\, \frac{1}{\sqrt{g_{\lambda \lambda }}}\, \textrm{D}_{\lambda } + {\textbf{e}}_{\varphi }\, \frac{1}{\sqrt{g_{\varphi \varphi }}}\, \textrm{D}_{\varphi } + {\textbf{e}}_{r}\, \frac{1}{\sqrt{g_{rr}}}\, \textrm{D}_{r}\Bigg ) \delta w\big (\lambda , \varphi , r\big ). \end{aligned}$$

(A1)

The gradient operator is expressed in terms of the spherical coordinates $\lambda$ (spherical longitude), $\varphi$ (spherical latitude), and r (geocentric radius). The spherical local reference frame is defined by the vector basis $({\textbf{e}}_{\lambda }, {\textbf{e}}_{\varphi }, {\textbf{e}}_{r})$, $g_{\lambda \lambda } = r^2 \cos ^2 \varphi$, $g_{\varphi \varphi } = r^2$, and $g_{rr} = 1$ are the metric coefficients. The symbol $\otimes$ indicates the dyadic product. The dyadic products of two basis vectors of the spherical local reference frame ${\textbf{e}}_{o} \otimes \, {\textbf{e}}_{p}, o,p \in \{\lambda ,\varphi ,r\}$, are termed spherical dyadics. The incremental gravity gradient tensor in Eqs. (A1), (A2), and (A3) is a linear combination of the spherical dyadics.

By performing the differentiations according to the right-hand side of Eq. (A1) and using the continuity of $\delta w$, we can obtain:

$$\begin{aligned}&\textrm{grad} \otimes \textrm{grad}\ \delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\quad = \Bigg \{{\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\lambda }\, \bigg [- \frac{\tan \varphi }{r^2} \textrm{D}_{\varphi } + \frac{1}{r} \textrm{D}_{r} + \frac{1}{r^2 \cos ^2 \varphi } \textrm{D}_{\lambda }^2\Bigg ] \nonumber \\&\qquad + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\varphi }\, \frac{1}{r^2} \bigg [r\, \textrm{D}_{r} + \textrm{D}_{\varphi }^2\bigg ] \nonumber \\&\qquad + ({\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\varphi } + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\lambda })\, \frac{1}{r^2\, \cos \varphi } \left( \tan \varphi \, + \textrm{D}_{\varphi }\right) \textrm{D}_{\lambda }\Bigg \}\delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\qquad - \Bigg \{({\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{r} + {\textbf{e}}_{r} \otimes {\textbf{e}}_{\lambda })\, \frac{1}{r\, \cos \varphi }\left( \frac{1}{r} - \textrm{D}_{r}\right) \textrm{D}_{\lambda } \nonumber \\&\qquad + ({\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{r} + {\textbf{e}}_{r} \otimes {\textbf{e}}_{\varphi })\, \frac{1}{r}\left( \frac{1}{r} - \textrm{D}_{r}\right) \textrm{D}_{\varphi }\Bigg \}\ \delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\qquad + \Bigg \{{\textbf{e}}_{r} \otimes {\textbf{e}}_{r}\, \textrm{D}_{r}^2\Bigg \}\ \delta w\big (\lambda , \varphi , r\big ). \end{aligned}$$

(A2)

The last equation decomposes the tensor into the three parts:

• (1) Tangential-Tangential (TT), the first term on the right-hand side,

• (2) Tangential-Normal (TN), the second term on the right-hand side,

• (3) Normal-Normal (NN), the third term on the right-hand side.

The three-part decomposition of the incremental gravity gradient tensor summarized in Tab. 10* can be obtained by applying Eq. (A2) to the spherical harmonic expansion of Eq. (11*) for the incremental gravity potential.

Alternatively, the tensor can be split as follows:

$$\begin{aligned}&\textrm{grad} \otimes \textrm{grad}\ \delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\quad = \Bigg \{\frac{1}{2} \left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\lambda } - {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\varphi }\right) \nonumber \\&\qquad \times \bigg [- \frac{\tan \varphi }{r^2} \textrm{D}_{\varphi } + \frac{1}{r^2 \cos ^2 \varphi } \textrm{D}_{\lambda }^2 - \frac{1}{r^2} \textrm{D}_{\varphi }^2\bigg ] \nonumber \\&\qquad + \left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\varphi } + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\lambda } \right) \, \frac{1}{r^2\, \cos \varphi } \left( \tan \varphi + \textrm{D}_{\varphi }\right) \textrm{D}_{\lambda }\Bigg \}\delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\qquad + \Bigg \{\frac{1}{2} \left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\lambda } + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\varphi }\right) \nonumber \\&\qquad \times \bigg [- \frac{\tan \varphi }{r^2} \textrm{D}_{\varphi } + \frac{2}{r} \textrm{D}_{r} + \frac{1}{r^2 \cos ^2 \varphi } \textrm{D}_{\lambda }^2 + \frac{1}{r^2} \textrm{D}_{\varphi }^2\bigg ]\Bigg \}\delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\qquad - \Bigg \{\left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{r} + {\textbf{e}}_{r} \otimes {\textbf{e}}_{\lambda }\right) \, \frac{1}{r\, \cos \varphi }\left( \frac{1}{r} - \textrm{D}_{r}\right) \textrm{D}_{\lambda } \nonumber \\&\qquad + \left( {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{r} + {\textbf{e}}_{r} \otimes {\textbf{e}}_{\varphi }\right) \frac{1}{r}\left( \frac{1}{r} - \textrm{D}_{r}\right) \textrm{D}_{\varphi }\Bigg \}\ \delta w\big (\lambda , \varphi , r\big ) \nonumber \\&\qquad + \Bigg \{{\textbf{e}}_{r} \otimes {\textbf{e}}_{r}\, \textrm{D}_{r}^2\Bigg \}\ \delta w\big (\lambda , \varphi , r\big ). \end{aligned}$$

(A3)

This four-part decomposition divides the tensor into:

• (1) Tangential-Shear (TS) part, the first term on the right-hand side,

• (2) Tangential-Dilation (TD) part, the second term on the right-hand side,

• (3) TN part, the third term on the right-hand side,

• (4) NN part, the fourth term on the right-hand side.

The four-part decomposition of the incremental gravity gradient tensor presented in Tab. 6* results from the application of Eq. (A3) to the spherical harmonic expansion of the incremental gravity potential in Eq. (11*).

1.2 The incremental gravity gradient tensor in the ellipsoidal local reference frame

The ellipsoidal coordinates are $\lambda$ (ellipsoidal longitude), $\varphi$ (reduced latitude), and u (minor semi-axis). The ellipsoidal local reference frame is given by a moving origin and the vector basis $({\textbf{e}}_{\lambda }, {\textbf{e}}_{\varphi }, {\textbf{e}}_{u})$. The dyadic products of two basis vectors of the ellipsoidal local reference frame ${\textbf{e}}_{o} \otimes \, {\textbf{e}}_{p}, o,p \in \{\lambda ,\varphi ,u\}$, represent ellipsoidal dyadics. The incremental gravity gradient tensor in Eqs. (A4), (A5), and (A6) is a linear combination of the ellipsoidal dyadics.

The incremental gravity gradient tensor referred to this frame is obtained from the twofold action of the gradient operator to $\delta w$:

$$\begin{aligned}&\textrm{grad} \otimes \textrm{grad}\ \delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\quad = \Bigg ({\textbf{e}}_{\lambda }\, \frac{1}{\sqrt{g_{\lambda \lambda }}}\, \textrm{D}_{\lambda } + {\textbf{e}}_{\varphi }\, \frac{1}{\sqrt{g_{\varphi \varphi }}}\, \textrm{D}_{\varphi } + {\textbf{e}}_{u}\, \frac{1}{\sqrt{g_{uu}}}\, \textrm{D}_{u}\Bigg ) \nonumber \\&\qquad \otimes \Bigg ({\textbf{e}}_{\lambda }\, \frac{1}{\sqrt{g_{\lambda \lambda }}}\, \textrm{D}_{\lambda } + {\textbf{e}}_{\varphi }\, \frac{1}{\sqrt{g_{\varphi \varphi }}}\, \textrm{D}_{\varphi } + {\textbf{e}}_{u}\, \frac{1}{\sqrt{g_{uu}}}\, \textrm{D}_{u}\Bigg )\delta w\big (\lambda , \varphi , u\big ), \end{aligned}$$

(A4)

with the metric coefficients $g_{\lambda \lambda } = (u^2 + \varepsilon ^2) \cos ^2 \varphi$, $g_{\varphi \varphi } = u^2 + \varepsilon ^2 \sin ^2 \varphi$, and $g_{uu} = (u^2 + \varepsilon ^2 \sin ^2 \varphi ) / (u^2 + \varepsilon ^2)$. $\varepsilon = \sqrt{a^2 - b^2}$ is the linear eccentricity given by the major semi-axis a and minor semi-axis b of the reference ellipsoid.

Differentiation by Eq. (A4) leads to the three-part decomposition in the form:

$$\begin{aligned}&\textrm{grad} \otimes \textrm{grad}\ \delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\quad = \Bigg \{{\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\lambda }\, \bigg [- \frac{\tan \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi } \nonumber \\&\qquad + \frac{u}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{u} + \frac{1}{(u^2 + \varepsilon ^2) \cos ^2 \varphi } \textrm{D}_{\lambda }^2\Bigg ] \nonumber \\&\qquad - {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\varphi }\, \frac{1}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \bigg [\frac{\varepsilon ^2 \sin \varphi \cos \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi } \nonumber \\&\qquad - \frac{u (u^2 + \varepsilon ^2)}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{u} - \textrm{D}_{\varphi }^2\bigg ] \nonumber \\&\qquad + ({\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\varphi } + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\lambda })\nonumber \\&\qquad \times \frac{1}{\sqrt{u^2 + \varepsilon ^2} \sqrt{u^2 + \varepsilon ^2 \sin ^2 \varphi }\, \cos \varphi } \left( \tan \varphi \, + \textrm{D}_{\varphi }\right) \textrm{D}_{\lambda }\Bigg \} \delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\qquad - \Bigg \{({\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{u} + {\textbf{e}}_{u} \otimes {\textbf{e}}_{\lambda })\, \frac{1}{\sqrt{u^2 + \varepsilon ^2 \sin ^2 \varphi }\, \cos \varphi }\nonumber \\&\qquad \times \left( \frac{u}{u^2 + \varepsilon ^2} - \textrm{D}_{u}\right) \textrm{D}_{\lambda } \nonumber \\&\qquad + ({\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{u} + {\textbf{e}}_{u} \otimes {\textbf{e}}_{\varphi })\, \frac{\sqrt{u^2 + \varepsilon ^2}}{u^2 + \varepsilon ^2 \sin ^2 \varphi }\nonumber \\&\qquad \times \bigg (\frac{u}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi } + \frac{\varepsilon ^2 \sin \varphi \cos \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{u} - \textrm{D}_{\varphi } \textrm{D}_{u}\bigg )\Bigg \}\delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\qquad + \Bigg \{{\textbf{e}}_{u} \otimes {\textbf{e}}_{u}\, \frac{1}{u^2 + \varepsilon ^2 \sin ^2 \varphi }\bigg [\frac{\varepsilon ^2 \sin \varphi \cos \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi } \nonumber \\&\qquad - \frac{u \varepsilon ^2 \cos ^2 \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{u} + (u^2 + \varepsilon ^2)\, \textrm{D}_{u}^2\bigg ]\Bigg \}\ \delta w\big (\lambda , \varphi , u\big ). \end{aligned}$$

(A5)

The first term on the right-hand side is the TT part, the second term is the TN part, and the third term is the NN part. The corresponding representation of the incremental gravity gradient tensor in Tab. 12* follows when the ellipsoidal harmonic expansion of the incremental gravity potential of Eq. (10*) is inserted into Eq. (A5).

Table 3 Numerical values of the integral (B1)

The four-part decomposition splits the tensor as follows:

$$\begin{aligned}&\textrm{grad} \otimes \textrm{grad}\ \delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\quad = \Bigg \{\frac{1}{2} \left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\lambda } - {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\varphi }\right) \nonumber \\&\qquad \times \bigg [- \frac{\tan \varphi \, (u^2 + \varepsilon ^2 - 2 \varepsilon ^2 \cos ^2 \varphi )}{(u^2 + \varepsilon ^2 \sin ^2 \varphi )^2} \textrm{D}_{\varphi } \nonumber \\&\qquad - \frac{u \varepsilon ^2 \cos ^2 \varphi }{(u^2 + \varepsilon ^2 \sin ^2 \varphi )^2} \textrm{D}_{u} \nonumber \\&\qquad + \frac{1}{(u^2 + \varepsilon ^2) \cos ^2 \varphi } \textrm{D}_{\lambda }^2 - \frac{1}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi }^2\bigg ] \nonumber \\&\qquad + \left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\varphi } + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\lambda } \right) \nonumber \\&\qquad \times \frac{1}{\sqrt{u^2 + \varepsilon ^2} \sqrt{u^2 + \varepsilon ^2 \sin ^2 \varphi }\, \cos \varphi } \left( \tan \varphi + \textrm{D}_{\varphi }\right) \textrm{D}_{\lambda }\Bigg \}\delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\qquad + \Bigg \{\frac{1}{2} \left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{\lambda } + {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{\varphi }\right) \bigg [- \frac{(u^2 + \varepsilon ^2) \tan \varphi }{(u^2 + \varepsilon ^2 \sin ^2 \varphi )^2} \textrm{D}_{\varphi } \nonumber \\&\qquad + \frac{u (2u^2 + 2\varepsilon ^2 - \varepsilon ^2 \cos ^2 \varphi )}{(u^2 + \varepsilon ^2 \sin ^2 \varphi )^2} \textrm{D}_{u} \nonumber \\&\qquad + \frac{1}{(u^2 + \varepsilon ^2) \cos ^2 \varphi } \textrm{D}_{\lambda }^2 + \frac{1}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi }^2\bigg ]\Bigg \}\delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\qquad - \Bigg \{\left( {\textbf{e}}_{\lambda } \otimes {\textbf{e}}_{u} + {\textbf{e}}_{u} \otimes {\textbf{e}}_{\lambda }\right) \nonumber \\&\qquad \times \frac{1}{\sqrt{u^2 + \varepsilon ^2 \sin ^2 \varphi }\, \cos \varphi }\left( \frac{u}{u^2 + \varepsilon ^2} - \textrm{D}_{u}\right) \textrm{D}_{\lambda } \nonumber \\&\qquad + \left( {\textbf{e}}_{\varphi } \otimes {\textbf{e}}_{u} + {\textbf{e}}_{u} \otimes {\textbf{e}}_{\varphi }\right) \nonumber \\&\qquad \times \frac{\sqrt{u^2 + \varepsilon ^2}}{u^2 + \varepsilon ^2 \sin ^2 \varphi }\bigg (\frac{u}{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi } + \frac{\varepsilon ^2 \sin \varphi \cos \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{u} \nonumber \\&\qquad - \textrm{D}_{\varphi } \textrm{D}_{u}\bigg )\Bigg \}\ \delta w\big (\lambda , \varphi , u\big ) \nonumber \\&\qquad + \Bigg \{{\textbf{e}}_{u} \otimes {\textbf{e}}_{u}\, \frac{1}{u^2 + \varepsilon ^2 \sin ^2 \varphi }\bigg [\frac{\varepsilon ^2 \sin \varphi \cos \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{\varphi } \nonumber \\&\qquad - \frac{u \varepsilon ^2 \cos ^2 \varphi }{u^2 + \varepsilon ^2 \sin ^2 \varphi } \textrm{D}_{u} + (u^2 + \varepsilon ^2)\, \textrm{D}_{u}^2\bigg ]\Bigg \}\ \delta w\big (\lambda , \varphi , u\big ). \end{aligned}$$

(A6)

Here, the first term on the right-hand side is the TS part, the second term is TD part, the third term is TN part, and the fourth term is the NN part. When the ellipsoidal harmonic expansion of the incremental gravity potential in Eq. (10*) is substituted into Eq. (A6) for $\delta w$, we get the four-part decomposition of the incremental gravity gradient tensor presented in Tab. 9*.

Appendix B Orthogonality properties of TT tensor spherical harmonics

Orthogonality of the TT tensor spherical harmonic functions can be deduced from the formula (1.1) in the ESM by setting $\varepsilon ^2 = 0$ and $u = 1.0$. Such integral vanishes $\forall l_1 \ne l_2$, because the integrals from the products $\cos l_1 \lambda \cos l_2 \lambda$ and $\sin l_1 \lambda \sin l_2 \lambda$ over the variable $\lambda$ are zero. Then, we have to consider only the case $l_1 = l_2 = l$ and the single variable integral over $\varphi$:

$$\begin{aligned}&\frac{1+\delta _{l,0}}{4} \int \limits _{-\pi /2}^{\pi /2} \bigg \{\bigg [\tan \varphi \, \textrm{D}_{\varphi } + (k_1+1) + \frac{l^2}{\cos ^2 \varphi }\bigg ]\, {\bar{P}}_{k_1 l}(\sin \varphi ) \nonumber \\&\quad \times \bigg [\tan \varphi \, \textrm{D}_{\varphi } + (k_2+1) + \frac{l^2}{\cos ^2 \varphi }\bigg ]\, {\bar{P}}_{k_2 l}(\sin \varphi ) \nonumber \\&\quad + \bigg [-(k_1+1) + \textrm{D}_{\varphi }^2\bigg ]\, {\bar{P}}_{k_1 l}(\sin \varphi )\ \nonumber \\ {}&\quad \times \bigg [-(k_2+1) + \textrm{D}_{\varphi }^2\bigg ] {\bar{P}}_{k_2 l}(\sin \varphi ) \nonumber \\&\quad + 4\ \frac{l}{\cos \varphi } \left( \tan \varphi + \textrm{D}_{\varphi }\right) {\bar{P}}_{k_1 l}(\sin \varphi )\, \nonumber \\ {}&\quad \times \frac{l}{\cos \varphi } \left( \tan \varphi + \textrm{D}_{\varphi }\right) {\bar{P}}_{k_2 l}(\sin \varphi ) \bigg \}\, \cos \varphi \, \textrm{d} \varphi . \end{aligned}$$

(B1)

Numerical values of the integral (B1) for harmonic degrees $k_1$; $k_2 = k_1, k_1 - 2, k_1 - 4, \ldots , 0$, and orders l up to 5 computed by the C program are listed in Table 3. Values smaller than $\pm 10^{-13}$ are numerically insignificant in comparison to other numbers and may safely be considered zero. When $l=0$, the integral (B1) vanishes for $k_1 \ne k_2$ and is nonzero otherwise. For $l>0$, on the other hand, the integral (B1) is nonzero not only for $k_1 = k_2$, but also for $k_2 = k_1 - 2, k_1 - 4, \ldots , 0$.

We calculated the integral (B1) for harmonic degrees $k_2 = k_1 - 1, k_1 - 3, k_1 - 5, \ldots , 1$ (not shown in Table 3). The numerical values are of the order 10$^{-13}$ and less, thus proving orthogonality only in this case.