1 Motivation

Most gravity surveys are now positioned using Global Navigation Satellite Systems (GNSSs), which deliver the 3D geodetic coordinates of each gravity observation after corrections for offsets between the GNSS antenna reference point and the gravity sensor. These can then be used to compute the gravity disturbance

$$\begin{aligned} \delta g_S =g_S -\gamma _S \end{aligned}$$
(1)

where \(\gamma _S\) is normal gravity at the same 3D position as the gravity observation \(g_S\) that has been reduced to gravity datum and corrected for instrumental drift and tides. To a second-order approximation near the Earth’s surface, \(\gamma _S\) can be computed analytically from the GNSS-derived ellipsoidal height h and geodetic latitude \(\varphi \) of the observation point S using (e.g., Heiskanen and Moritz 1967, p 79)

$$\begin{aligned} \gamma _S = \gamma \left[ {1-2(1+f+m-2f\sin ^{2}\phi ) \frac{h}{a}+3\frac{h^{2}}{a^{2}}}\right] \end{aligned}$$
(2)

with

$$\begin{aligned} \gamma =\gamma _e \frac{1+k\sin ^{2}\varphi }{\sqrt{1-e^{2} \sin ^{2}\varphi }} \end{aligned}$$
(3)

where f is the geometrical flattening of the normal ellipsoid, m is the ratio of gravitational and centrifugal accelerations at the equator of the normal ellipsoid, a is the equatorial radius of the normal ellipsoid, k is the normal gravity constant, \(\gamma _{e}\) is normal gravity acceleration at the equator, and e is the first numerical eccentricity of the normal ellipsoid. Numerical values of these parameters for the GRS80 normal ellipsoid are given in Moritz (1980) and reprinted in the Geodesist’s Handbook.

While the downward continuation of gravity disturbances is beyond the scope of this article, it nevertheless requires that: (1) \(\delta g_S\) are downward-continued to the geoid to give \(\delta g\) before being convolved with Hotine’s kernel or some modification thereof; and (2) the gravity disturbances must correspond to a harmonic disturbing potential in the solution domain. The availability of these downward-continued gravity disturbances \(\delta g\) then allows for computation of the geoid N by Hotine’s integral (Hotine 1969, Chap 29), which is a solution to a fixed or the second boundary-value problem of potential theory in spherical approximation (e.g., Heiskanen and Moritz 1967, p 36).

One advantage of using gravity disturbances over gravity anomalies to compute the geoid is that they are not adversely affected by, e.g., uncertain or ambiguous realisations of vertical datums and their associated height systems. In Australia, for instance, the vertical datum contains a confirmed tilt with respect to the geoid (Featherstone and Filmer 2012), regional distortions [e.g., Featherstone et al. (2011) and the citations therein], and uses a height system that is not compatible with the geoid (Filmer et al. 2010). These will cause errors in the computed terrestrial gravity anomalies (cf. Heck 1990) that can propagate into the combined geoid solution unless modelled and/or filtered out (cf. Vaníček and Featherstone 1998; Wang et al. 2011; Featherstone et al. 2011).

Modern regional geoid models are often computed from an Earth gravitational model (EGM) augmented by one, more or all of land, airborne, ship-borne and altimeter-derived marine gravimetry, depending on data availability and the region of interest. GNSS is used to coordinate and navigate ship-borne gravity surveys, which results in better Eötvös corrections and horizontal locations so that crossover adjustments are more effective (cf. Wessel and Watts 1988). As the mean sea surface departs from the geoid by up to \({\sim }2\) m, account also has to be made for mean dynamic topography (MDT) in the derivation of gravity disturbances from satellite altimetry (Zhang 1998). Land and airborne gravity surveys are also coordinated with GNSS, thus facilitating the direct computation of gravity disturbances. A prior geoid model would be needed to determine gravity anomalies from GNSS-coordinated gravity surveys, resulting in a circular argument (cf. Vaníček et al. 1992), and reinforcing the benefit of directly using gravity disturbances in Hotine’s integral.

Hotine’s integral, or its inverse, has been used for: (1) geoid determination from GNNS-positioned airborne gravimetry (e.g., Schwarz and Li 1996; Kearsley et al. 1998; Forsberg et al. 2000; Novák and Heck 2002; Novák 2003; Novák et al. 2003; Alberts and Klees 2004; Serpas and Jekeli 2005; Sjöberg and Eshagh 2009) or land gravimetry (e.g., Kirby 2003), (2) marine gravity field and MDT estimation from satellite radar altimetry (e.g., Rapp 1980; Zhang and Blais 1993; Rapp and Wang 1994; Zhang and Sideris 1996; Zhang 1998), and (3) global Earth and planetary gravity field modelling (e.g., Sjöberg 1989; Barriot and Balmino 1992).

However, the spatial coverage of GNSS-coordinated gravity data is currently limited, so there is a need to modify Hotine’s integral to reduce the truncation error that results from the omission of gravity disturbances in the far zones beyond the area of interest. Admittedly, this truncation error can be reduced by the inclusion of a high-degree EGM (Sect. 2.2). However, the kernel modification and cap radius can be used additionally as a filter to tune the relative data contributions (cf. Vaníček and Featherstone 1998; Kern et al. 2003; Featherstone 2003a). It is this property that will be emphasised more than only reduction of the truncation error.

Early modifications to Stokes’s kernel (e.g., Molodensky et al. 1962) were formulated to reduce the truncation error alone because of limited spatial coverage of terrestrial gravity anomalies at that time. Many of the subsequent modifications also consider EGMs (Appendix A). The limited spatial coverage of gravity disturbances coordinated by GNSS is probably the same now as it was in the years soon after the portable gravimeter was developed, so the motivation for modifications to Hotine’s kernel is similar now to as it was then for Stokes’s kernel. However, the advent of new EGMs derived from dedicated satellite gravimetry (e.g., Pail et al. 2011) and EGM2008 (Pavlis et al. 2012) has changed the approaches to regional geoid computation for a couple of reasons: (1) the truncation and approximation errors are lessened considerably when a very high-degree EGM is used (Sect. 2.2), and (2) a satellite-derived EGM that is more reliable in the low and medium degrees allows for higher degrees of modification so as to place more reliance on the geoid derived from satellite gravimetry (e.g., Sects. 3.1 and 5).

Mostly in analogy to those already proposed and used for Stokes’s integral and gravity anomalies (Appendix A), this article will present: (1) deterministic modifications to Hotine’s integral, where the user can control the errors in some prescribed ways; (2) stochastic modifications to Hotine’s integral, where error spectra are embedded in an attempt to control the balance amongst data and/or truncation errors; and (3) their band-limited hybrid combinations, where they are combined according to the perceived relative benefits of each. For instance, the ‘user’ may have a good understanding of the data errors in parts of the geopotential spectrum so can use a stochastic modifier in those bands, but say use a deterministic modifier in other bands (cf. Sect. 5). The options are many, but the hybrid combinations do provide much more flexibility for the ‘user’.

Previous authors who have investigated modifications to Hotine’s kernel comprise Jekeli (1979, 1980b); Guan and Li (1990); Sjöberg and Nord (1992); Vaníček et al. (1992); Zhang and Blais (1993); Zhang (1998); Novák (2003); Novák et al. (2003); Alberts and Klees (2004); and Sjöberg and Eshagh (2009), so their results will only be summarised as part of the review component of this paper. Importantly, most of these authors conclude that Hotine’s integral with gravity disturbances can be superior to Stokes’s integral with gravity anomalies. The other modifications to Hotine’s kernel presented herein will be based on adaptations of the principles previously applied to Stokes’s kernel, as well as some new formulations that can be applied back to Stokes’s or other kernels.

2 Basics

2.1 Spherical Hotine integral

In terms of spherical polar coordinates of spherical distance \(\psi \) and azimuth \(\alpha \) centred on each computation point, Hotine’s integral in spherical approximation reads

$$\begin{aligned} N=\frac{r}{4\pi \gamma }\int _0^{2\pi } {\int _0^\pi {H \delta g \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(4)

and the spherical Hotine kernel is

$$\begin{aligned} H&= \csc {(\psi /2)}-\ln \left( {1+\csc {(\psi }/{2)}} \right) \nonumber \\&= \sum _{n=0}^\infty {\frac{2n+1}{n+1}} \;P_n (\cos \psi ) \end{aligned}$$
(5)

where \(P_n (\cos \psi )\) is the Legendre polynomial of degree n, and r is the radius to the surface of the normal ellipsoid, which can reduce the ellipsoidal correction for the spherical approximation to a manageably small value (cf. Claessens 2006, Chap 6), and \(\gamma \) is normal gravity on the surface of the normal ellipsoid (Eq. 3) as is demanded by Bruns’s formula (e.g., Heiskanen and Moritz 1967, p 85).

When the integration domain in Eq. (4) is truncated to a spherical cap of radius \(\psi _0\) centred on each computation point, this results in an approximation of the geoid height

$$\begin{aligned} \widehat{N}_1 =\frac{r}{4\pi \gamma }\int _0^{2\pi } {\int _0^{\psi _0} {H \delta g \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(6)

where the corresponding truncation error (\(N=\widehat{N_1}+\Delta N\)) is

$$\begin{aligned} \Delta N=\frac{r}{2\pi }\sum _{n=0}^\infty {\left[ {\int _{\psi _0}^\pi {H P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }} \right] } \;\delta g_n \end{aligned}$$
(7)

and \(\delta g_n\) is the \(n\)-th degree spherical harmonic of the gravity disturbance, and the integral in square parentheses yields the truncation coefficients. Recursion formulas for Eq. (7) are given in, e.g., the Appendices of Jekeli (1979). In the remainder of this article, only the integral forms will be presented instead of cluttering the presentation with [too many] recursions.

2.2 Inclusion of an EGM

One very simple way to reduce \(\Delta N\) in Eq. (7) is by subtracting the gravity disturbances computed from an EGM to spherical harmonic degree L (\(\delta g_L\)) from the observed and downward-continued \(\delta g\), and then add back the geoid contribution of the same EGM to the same degree. In physical geodesy, this is commonly referred to as the remove– compute–restore technique.

Albeit well known for Stokes’s integral (e.g., Vincent and Marsh 1974; Rapp and Rummel 1975), this strategy for Hotine’s integral seems first-attributable to Rapp (1980). This gives the residual gravity disturbance

$$\begin{aligned} \delta g^{L}=\delta g-\sum _{n=0}^L {\delta g_n} \end{aligned}$$
(8)

Equation (6) then becomes

$$\begin{aligned} \widehat{N_2}=N_L +\frac{r}{4\pi \gamma }\int _0^{2\pi } {\int _0^{\psi _0} {H \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(9)

where \(N_L\) is the geoid undulation given by an EGM to degree L, and the truncation error becomes

$$\begin{aligned} \Delta N^{L}=\frac{r}{2\pi }\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }} \right] } \;\delta g_n \end{aligned}$$
(10)

If \(\left\| {\delta g^{L}} \right\| <\left\| {\delta g} \right\| \) over the far zones beyond the spherical cap (\(\psi _0 <\psi \le \pi \)), then \(\delta N^{L}\) should be \({<}\delta N\), indicating that this strategy can reduce the truncation error (cf. Vaníček and Sjöberg 1991). However, the practical implementation is not as simple as the theory may imply. Any EGM is imperfect, and if the terrestrial gravity data contain low-frequency errors, there will be leakage of errors during the evaluation of the last term in Eq. (9) (cf. Vaníček and Featherstone 1998). Therefore, it is important to consider the kernel modification not only as means to reduce the truncation error, but also as a filter to reduce leakage of any low-frequency errors from the terrestrial data (cf. Omang and Forsberg 2002; Featherstone et al. 2011; Wang et al. 2011).

Another benefit of including an EGM in is that the residual geoid height computed from the last term in Eq. (9)—as well as its variants presented in the remainder of this article—is smaller in magnitude (a metre or so vs. up to 100 m), so are less subject to approximation errors. Also, as EGMs become more homogeneously reliable—notably those derived from satellite gravimetry—the motivation for using the remove–compute–restore approach will become stronger.

In order to remain as general as possible, and in line with current widespread practice of utilising an EGM in regional geoid computations under the remove–compute–restore scheme, the following modifications will be applied only to Eqs. (5) and (9). The constants outside the integral terms will be abbreviated to \(\kappa =r{/}{(4\pi \gamma })\) and \(c=r{/}{(2\pi }).\)

3 Deterministic and hybrid modifications

3.1 Remove Legendre polynomials (modification D1)

A simple deterministic kernel modification is to subtract polynomial terms from Eq. (3). For Stokes’s kernel, this approach is generally attributed to Wong and Gore (1969), though de Witte (1967) alluded to it.

Vaníček et al. (1992) and Sjöberg and Nord (1992) have applied this approach to Hotine’s kernel, thus Eq. (9) becomes what herein is termed the spheroidal Hotine kernel

$$\begin{aligned} \widehat{N_{D1}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D1} \left( M \right) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \text{ for} M\le L\nonumber \\ \end{aligned}$$
(11)

with the subscript D1 denoting this as the first deterministic modification and so on

$$\begin{aligned} H_\mathrm{D1} (M)&= H-\sum _{n=0}^M {\frac{2n+1}{n+1}} \;P_n (\cos \psi ) \nonumber \\&= \sum _{n=M+1}^\infty {\frac{2n+1}{n+1}} \;P_n (\cos \psi ) \end{aligned}$$
(12)

and

$$\begin{aligned} \Delta N_\mathrm{D1} = c\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{D1} \left( M \right) P_n \left( {\cos \psi } \right) \sin \psi \,\mathrm{d}\psi }} \right] } \;\delta g_n\nonumber \\ \end{aligned}$$
(13)

This modification makes Eq. (11) less sensitive to low-frequency errors in the terrestrial gravity disturbances by partially filtering them out (cf. Omang and Forsberg 2002; Wang et al. 2011). However, the filtering is only ever partial because the truncated integration domain allows for spectral leakage (cf. Vaníček and Featherstone 1998).

In the context of band-limited kernel modifications, the summation in Eq. (12) does not necessarily have to start at degree \(n=0\) (cf. Featherstone 2003a). Such a band-limited modification to Stokes’s kernel was used by Li and Sideris (1994) and to Hotine’s kernel by Novák and Heck (2002) and Novák et al. (2003); also see Colombo (1977). Extending this yet further, and not violating the restriction \(M \le L\) (otherwise components of the combined geoid model will be omitted), the band-limited modification can be applied over multiple bands of the user’s choice. This modification is straightforward to implement in existing software as recursion routines for Legendre polynomials are widely available (e.g., Press et al. 2007, Chap 6.7).

Figure 1 shows an example of the D1-modified Hotine kernel (Eq. 12) in relation to the spherical Hotine kernel (Eq. 5). The degree of D1 modification has been chosen arbitrarily at \(M=50\). As the degree of modification increases, the D1-modified kernel oscillates more rapidly. Thus, high degrees of modification can slow the numerical integration because more nodes are needed to determine the integral mean value of the modified kernel over each element (cf. Hirt et al. 2011).

Fig. 1
figure 1

Solid line The spherical Hotine kernel (Eq. 5); dotted line a spheroidal Hotine kernel (Eq. 12) for \(M=50\)

Full size image

3.2 Set kernel to zero at \(\psi _0\) (modification D2)

Another easy-to-implement deterministic modification is to set the kernel to zero at the truncation radius \(\psi _0\) by subtraction. There is some conjecture as to whether this type of modification, albeit for Stokes’s integral, should be first-attributed to Ostach (1970) or Meissl (1971), but the latter author has gained wider acceptance in the literature on modifications to Stokes’s kernel.

Rapp (1980) applied this strategy to Hotine’s integral, then attributing it to Jekeli (1980b), presumably because both papers were under consideration at around the same time; this is

$$\begin{aligned} \widehat{N_\mathrm{D2}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D2} (\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(14)

with

$$\begin{aligned} H_\mathrm{D2} (\psi _0)=H-H(\psi =\psi _0) \end{aligned}$$
(15)

and

$$\begin{aligned} \Delta N_\mathrm{D2} =c\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{D2} (\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }} \right] } \;\delta g_n\nonumber \\ \end{aligned}$$
(16)

Recursion formulas for the integral term in Eq. (16) are given in Guan and Li (1990) and Jekeli (1980b).

By forcing the kernel to zero at the truncation radius accelerates the convergence rate of the truncation error from \(\mathcal O \left( n^{{-}1}\right) \text{ to} \mathcal O \left( n^{{-}2}\right) \) (cf. Jekeli 1980a, 1981; Featherstone et al. 1998). However, faster convergence of a series does not necessarily guarantee smaller values of its coefficients. Since the truncation error is the sum of all terms (Eq. 16), the truncation error may even increase, as shown by Jekeli (1980a, 1981) for Stokes’s kernel. However, Guan and Li (1990) claim that this Ostach–Meissl-type modification to the Hotine kernel does decrease the truncation error.

3.3 Deterministic hybrid (modification D3)

The deterministic modifications in Eqs. (12) and (15) can be combined, as first proposed by Heck and Grüninger (1987) for Stokes’s integral, such that the Legendre polynomials are removed to some degree M such that (s.t.) the zero-crossing point of the D1-modified Hotine kernel in Eq. (12) coincides with the truncation radius \(\psi _0\); this is

$$\begin{aligned} \widehat{N_\mathrm{D3}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D3} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(17)

with

$$\begin{aligned} H_\mathrm{D3} (M,\psi _0)\!=\!H(M)\;\text{ s}.\text{ t}.\;H_\mathrm{D3} (M,\psi _0)\!=\!0\;\text{ at}\;\psi \!=\!\psi _0 \end{aligned}$$
(18)

and

$$\begin{aligned} \Delta N_\mathrm{D3} \!=\!c\sum _{n=L\!+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{D3} (M,\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }} \right] } \;\delta g_n \nonumber \\ \end{aligned}$$
(19)

This hybrid modification simultaneously exploits the partial high-pass filtering properties of Eq. (12) and the accelerated rate of convergence of the truncation error from Eq. (16). However, this hybrid is less flexible to implement because M and \(\psi _0\) are inextricably linked in this case. For instance, a small \(\psi _0\) will dictate that M has to be large (cf. Fig 1) and possibly too much filtering will occur, and vice versa.

3.4 Deterministic hybrid (modification D4)

To counteract the above restriction, an alternative combination of Eqs. (12) and (15) can be implemented, where the Legendre polynomials are removed to a user-chosen degree M, and then the kernel is set to zero at the truncation radius \(\psi _0\) by subsequent subtraction (cf. Heck and Grüninger 1987; Featherstone et al. 1998).

This type of hybrid modified Hotine kernel was first introduced by Alberts and Klees (2004), but they did not elaborate upon its properties. In the notation adopted herein, this is

$$\begin{aligned} \widehat{N_\mathrm{D4}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D4} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(20)

with

$$\begin{aligned} H_\mathrm{D4} (M,\psi _0)=H(M)-H(M,\psi =\psi _0) \end{aligned}$$
(21)

and

$$\begin{aligned} \Delta N_\mathrm{D4}&= c\sum _{n=L+1}^\infty \nonumber \\&\times {\left[ {\int _{\psi _0}^\pi {H_\mathrm{D4} (M,\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\! \right] }\! \;\delta g_n\quad \end{aligned}$$
(22)

This hybrid offers the same properties as Eqs. (15) and (18) in terms of the accelerated convergence rate of the truncation error, but does not suffer the same restrictions, thus giving the ‘user’ much more control over the degree of filtering and integration domain chosen. Alberts and Klees (2004) chose \(M=20\) and \(\psi _0 =5^{\circ }\) from simulations and experiments with airborne gravity in their case study area. Importantly, the parameter values chosen depend on the study area, data sets used, their resolution and spatial extent; indeed, this applies to all the modifications.

3.5 Molodensky-type approach (modification D5)

(Molodensky et al. (1962), Chap 7) presented an approach to reduce the \(\text{ L}_{2}\) norm of the truncation error for the spherical Stokes’s kernel. This was later adapted for a higher-than-second-degree reference spheroid by Vaníček and Kleusberg (1987) and Vaníček and Sjöberg (1991); also see Martinec and Vanicek (1997). The Molodensky et al. (1962) approach for the spherical [not spheroidal] Hotine integral (Eq. 5) was presented in Sjöberg and Eshagh (2009), and Novák (2003) presented a band-limited version of the same kernel. However, neither of these modifications allow for the inclusion of an EGM, whereas the presentation herein does.

In analogy, the Molodensky-type modification to the spheroidal Hotine kernel (Eq. 12) was presented by Zhang (1998) as

$$\begin{aligned} \widehat{N_{\mathrm{D5}}}\!=\!N_L \!+\!\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D5} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(23)

with

$$\begin{aligned} H_\mathrm{D5} (M,\psi _0)\!=\!H(M)\!-\!\sum _{n=0}^M {\frac{2n\!+\!1}{2}} \;h_n (\psi _0) P_n (\cos \psi ) \end{aligned}$$
(24)

and

$$\begin{aligned}&\Delta N_\mathrm{D5} \nonumber \\&\!=\!c\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{D5} (M,\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }} \right] } \;\delta g_n \nonumber \\ \end{aligned}$$
(25)

Once the values of \(\psi _0\) and M have been chosen by the ‘user’, the \(h_n (\psi _0)\) modification coefficients are determined from the solution of a set of linear equations. However, the formulation of (Zhang (1998), Sect. 2.2) is not particularly intuitive for practical application, and also contains two typographical errors. Thus, a clearer formulation is given below.

By inference from Vaníček and Sjöberg (1991, Eq.18), the \(\text{ L}_{2}\) norm of the Hotine truncation error (Eq. 25) is minimised when

$$\begin{aligned} \int _{\psi _0}^\pi {H_\mathrm{D5} (M,\psi _0)P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi } =0, \quad 0\le n\le M\nonumber \\ \end{aligned}$$
(26)

Inserting Eqs. (24) and (12) in Eq. (26) gives

$$\begin{aligned} \int _{\psi _0}^\pi \left[ H(\psi )-\sum _{k=0}^M {\frac{2k+1}{k+1}} \;P_k (\cos \psi ) -\sum _{k=0}^M {\frac{2k+1}{2}} \;h_k (\psi _0) P_k (\cos \psi ) \!\right] \!P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi =0\nonumber \\ \end{aligned}$$
(27)

Using the abbreviation

$$\begin{aligned} e_{nk} (\psi _0)=\int _{\psi _0}^\pi {P_n (\cos \psi ) P_k (\cos \psi )\sin \psi \,\mathrm{d}\psi } \end{aligned}$$
(28)

leaves the desired system of (\(M+1\)) linear equations

$$\begin{aligned}&\sum _{k=0}^M {\frac{2k+1}{2}} \;h_k (\psi _0) e_{nk} (\psi _0)\nonumber \\&\quad =\int _{\psi _0}^\pi {H P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi } -\sum _{k=0}^M {\frac{2k+1}{k+1}} e_{nk} (\psi _0) \end{aligned}$$
(29)

A recursion formula for the second term in Eq. (29) (cf. Eq. 7) is derived in the Appendices of Jekeli (1979), and recursions for \(e_{nk} (\psi _0)\) (Eq. 28) are given in Paul (1973) or Hagiwara (1972, 1976). The right-hand-side of Eq. (29) is also the recursion used to compute the integral term in Eq. (13). For instance, the computer code from Featherstone (2003b) can be adapted to compute Eq. (29) and thence Eq. (24), as well as the other deterministic modifications presented herein.

3.6 Deterministic hybrid (modification D6)

The modified kernel in Eq. (24) can be forced to be zero at \(\psi _0\) by appropriate selections of \(h_n (\psi _0)\), e.g., by varying M and \(\psi _0\). This approach was first suggested for the Molodensky-modified spheroidal Stokes’s integral (cf. Sect. 3.5) by Featherstone et al. (1998), but it can also be applied to Eq. (24). In analogy to modification D3 (Sect. 3.3), \(h_n (\psi _0)\) are chosen such that the zero-crossing point of the modified Hotine kernel in Eq. (24) coincides with the truncation radius \(\psi _0\)

$$\begin{aligned} \widehat{N_\mathrm{D6}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D6} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(30)

with

$$\begin{aligned} H_\mathrm{D6} (M,\psi _0)=H_\mathrm{D5} (M,\psi _0)\; \text{ s.t} .\;H_\mathrm{D6} (M,\psi _0)\nonumber \\=0\;\text{ at}\;\psi =\psi _0 \end{aligned}$$
(31)

and

$$\begin{aligned} \Delta N_\mathrm{D6} =c\!\!\!\sum _{n=L+1}^\infty \!{\left[ \!{\int _{\psi _0}^\pi {H_\mathrm{D6} (M,\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\right] } \delta g_n\nonumber \\ \end{aligned}$$
(32)

However, the practical evaluation of Eq. (31) can be quite cumbersome because trial and error has to be used to meet the condition that this modified Hotine kernel is zero at \(\psi _0\). For some starting values of \(\psi _0\) and M, the kernel in Eq. (24) has to be evaluated (involving the inversion of Eq. 29), then plotted to see if it is zero at \(\psi _0\). If not, then the values of \(\psi _0\) and/or M have to be adjusted until it is. Evidently, this may involve a lot of work and is not so attractive given the following option.

3.7 Deterministic hybrid (modification D7)

A far simpler way to set the D5-modified kernel (Eq. 24) to zero at \(\psi _0\) is by subtraction (cf. Sects. 3.2 and 3.4). This strategy was suggested by Featherstone et al. (1998) for Stokes’s kernel. It allows for more user control over the values chosen for \(\psi _0\) and/or M in terms of the filtering properties of the kernel.

$$\begin{aligned} \widehat{N_\mathrm{D7}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{D7} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(33)

with

$$\begin{aligned} H_\mathrm{D7} (M,\psi _0)=H_\mathrm{D5} (M,\psi _0)-H_\mathrm{D5} (M,\psi =\psi _0) \end{aligned}$$
(34)

and

$$\begin{aligned} \Delta N_\mathrm{D7} =c\!\!\!\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{D7} (M,\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\!\right] } \;\delta g_n\nonumber \\ \end{aligned}$$
(35)

Practical implementation just involves the evaluation of Eqs. (29) and (24) for the user-chosen values of \(\psi _0\) and M and subtraction of the value of \(H_\mathrm{D5} (M,\psi _0)\) at \(\psi _0\) for \(0< \psi \le \psi _0 \)

4 Stochastic and hybrid modifications

Stochastic kernel modifications can be more subjective than deterministic modifications because reliable error spectra of the data involved are not always available or reliable (e.g., Sjöberg and Hunegnaw 2000). The error spectra of EGMs can be unrealistic if they are derived only from the diagonal of their variance covariance (VCV) matrices and are global estimates, so do not necessarily reflect the errors in a particular region (cf. Sjöberg 2005). Attempts are sometimes made to ‘calibrate’ these error spectra; nevertheless, they still remain global estimates.

The error spectra of terrestrial gravity data are even more problematic to estimate (e.g., Kern et al. 2003) and can vary quite considerably from region to region. Most often, simple covariance models are used (e.g., Ellmann 2005a), which render the stochastic modifiers more akin to least squares collocation and thus subject to the same simplifying assumptions such as stationarity and isotropy. As such, the physical acceptability of the stochastic modifiers is arguably less than for the deterministic modifiers.

4.1 Wenzel-type approach (modification S1)

Wenzel (1981, 1982, 1983) implemented a Wiener-type filter in Stokes’s integral, which can also be applied to Hotine’s integral to give

$$\begin{aligned} \widehat{N_\mathrm{S1}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{S1} \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(36)

with the subscript S1 denoting this as the first stochastic modification and so on

$$\begin{aligned} H_\mathrm{S1} =\sum _{n=0}^\infty {\frac{2n+1}{n+1}} \;w_n P_n (\cos \psi ) \end{aligned}$$
(37)

and

$$\begin{aligned} w_n =\frac{\sigma _n^2 \{\delta g_\mathrm{EGM} \}}{\sigma _n^2 \{\delta g_\mathrm{EGM}\}+\sigma _n^2 \{\delta g_\mathrm{T} \}}\;0\le n\le L \end{aligned}$$
(38)

where \(\sigma _n^2 \{\delta g_\mathrm{EGM}\}\) is the error degree variance of the gravity disturbances from the EGM and \(\sigma _n^2 \{\delta g_\mathrm{T}\}\) is the error degree variance of the terrestrial gravity disturbances. This modification is—by necessity—restricted to the degree L of EGM used in the combined solution for the geoid, such that

$$\begin{aligned} H_\mathrm{S1} =\sum _{n=0}^L {\frac{2n+1}{n+1}} \;w_n P_n (\cos \psi )\!+\!\!\!\sum _{n=L+1}^\infty {\frac{2n+1}{n+1}} \;P_n (\cos \psi )\nonumber \\ \end{aligned}$$
(39)

leaving

$$\begin{aligned} \Delta N_\mathrm{S1} {=}c\!\!\!\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\! \right] } \delta g_n =\Delta N\qquad \end{aligned}$$
(40)

to show that no specific attempt has been made to reduce the truncation error; it is just as large as for the truncated spherical Hotine integral (Eq. 7). Nevertheless, the truncation error is reduced already because of the use of the EGM to degree \(L\) (cf. Sect. 2.2).

4.2 Stochastic hybrid (modification S2)

Similar to the D1 modification (Sect. 3.1), the degree to which the Wiener-type filter is applied can be limited to any \(M \le L\) or any band(s) in that domain.

$$\begin{aligned} \widehat{N_\mathrm{S2}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{S2} (M) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(41)

with, e.g.,

$$\begin{aligned} H_\mathrm{S2} (M)&= \sum _{n=0}^M {\frac{2n+1}{n+1}} \;w_n P_n (\cos \psi )\nonumber \\&+\sum _{n=M+1}^\infty {\frac{2n+1}{n+1}} \;P_n (\cos \psi ) \end{aligned}$$
(42)

and \(\Delta N_\mathrm{S2} =\Delta N_\mathrm{S1} =\Delta N\), showing again that there is no reduction of the truncation error.

4.3 Stochastic hybrid (modification S3)

An Ostach–Meissl-type modification (cf. Sect. 3.2) can also be applied to Eq. (42), noting that if \(M=L\) it degenerates to Eq. (39) so can be implemented simply for both options.

$$\begin{aligned} \widehat{N_\mathrm{S3}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{S3} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(43)

with

$$\begin{aligned} H_\mathrm{S3} (M,\psi _0)=H_\mathrm{S2} (M)-H_\mathrm{S2} (M,\psi =\psi _0) \end{aligned}$$
(44)

and

$$\begin{aligned} \Delta N_\mathrm{S3} =c\!\!\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{S3} (M,\psi _0)P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\!\right] } \;\delta g_n\nonumber \\ \end{aligned}$$
(45)

This additional modification accelerates the rate of convergence of the truncation error.

4.4 Stochastic hybrid (modification S4)

This stochastic hybrid is an analogue of Heck and Grüninger (1987) and Featherstone et al. (1998) to achieve an accelerated rate of convergence of the truncation error without subtraction (Eq. 45), but by selecting the value of M for which the value of the kernel is zero at \(\psi _0\).

$$\begin{aligned} \widehat{N_\mathrm{S4}}=N_L +\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{S4} (M,\psi _0) \delta g^{L} \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} \end{aligned}$$
(46)

with

$$\begin{aligned} H_\mathrm{S4} (M,\psi _0)\!=\!H_\mathrm{S2} (M)\;\text{ s.t}\;H_\mathrm{S4} (M,\psi _0)\!=\!0\;\text{ at}\;\psi \!=\!\psi _0 \end{aligned}$$
(47)

and

$$\begin{aligned} \Delta N_\mathrm{S4} \!=\!c\!\!\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{S4} (M,\psi _0)P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\!\right] } \;\delta g_n\nonumber \\ \end{aligned}$$
(48)

As for modifications D3 and D6 (Sects. 3.3 and 3.6), this requires cumbersome trial and error to determine the appropriate value of M, but is also complicated further by the choice of EGM used to provide \(\sigma _n^2 \{\delta g_\mathrm{EGM}\}\) and the model adopted for \(\sigma _n^2 \{\delta g_\mathrm{T}\}\) (i.e., via Eq. 38).

4.5 Sjöberg-type approach (modification S5)

Sjöberg (1980a, b, 1981, 1984a, b, 1986, 1991, 2003c) and Sjöberg and Hunegnaw (2000) have provided a series of incremental stochastic modifications to Stokes’s kernel, culminating in the variant in Sjöberg (2003b, Eq. 26). An attractive aspect of most of the Sjöberg-type modifications is that they attempt to simultaneously reduce the truncation error and errors originating from the EGM and terrestrial gravity data, or subsets thereof. The principal restriction is reliably estimating the error spectra of the terrestrial gravity disturbances \(\sigma _n^2 \{\delta g_T\}\), as well as the other caveats mentioned at the start of this Section.

Assuming that Sjöberg (2003b) gives the ‘final word’ on this class of modifications for Stokes’s integral, when applied to Hotine’s integral gives

$$\begin{aligned} \widehat{N_\mathrm{S5}}\!=\!\kappa \int _0^{2\pi } {\int _0^{\psi _0} {H_\mathrm{S5} (M,L) \delta g \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} +c\!\!\sum _{n=2}^L {b_n} \;\delta g_n \end{aligned}$$
(49)

with

$$\begin{aligned} H_\mathrm{S5} (M,L)=H-\sum _{n=0}^M {\frac{2n+1}{2}} \;s_n P_n (\cos \psi ) \end{aligned}$$
(50)

and

$$\begin{aligned} \Delta N_\mathrm{S5}=c\!\!\sum _{n=L+1}^\infty {\left[ {\int _{\psi _0}^\pi {H_\mathrm{S5} (M,L)P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\!\right] } \;\delta g_n\nonumber \\ \end{aligned}$$
(51)

Depending on the choices of M and L, the \(s_n\) modification coefficients are

$$\begin{aligned} 2\!\le \! n\!\le \! \min (M,L): s_n&\!=\!&\frac{2\sigma _n^2 \{\delta g_\mathrm{T} \}\left( {c_n \{\delta g\}+\sigma _n^2 \{\delta g_\mathrm{EGM} \}} \right) }{(n\!+\!1) d_n \{\delta g\}}\end{aligned}$$
(52)
$$\begin{aligned} (L\!+\!1)\le n\le M:\;s_n&= \frac{2\sigma _n^2 \{\delta g_\mathrm{T} \}}{(n\!+\!1)\left( {\sigma _n^2 \{\delta g_\mathrm{T} \}+c_n \{\delta g\}} \right) } \end{aligned}$$
(53)

where

$$\begin{aligned} d_n \{\delta g\}&= c_n \{\delta g\}\sigma _n^2 \{\delta g_{\mathrm{EGM}}\}\nonumber \\&\quad +\sigma _n^2 \{\delta g_\mathrm{T} \}\left( {c_n \{\delta g\}+\sigma _n^2 \{\delta g_\mathrm{{EGM}}\}}\right) \end{aligned}$$
(54)

and \(c_n \{\delta g\}\) is the degree variance of the gravity disturbances. For weighting the contribution of the EGM, the \(b_n\) coefficients are

$$\begin{aligned} 2\le n\le \min (M,L):\;b_n&= \frac{2\sigma _n^2 \{\delta g_\mathrm{T} \}c_n \{\delta g\}}{(n+1) d_n \{\delta g\}}\end{aligned}$$
(55)
$$\begin{aligned} n>\min (M,L):\;b_n&= 0 \end{aligned}$$
(56)

Ellmann (2005a, 2012) provides computer code that can be adapted to compute the Sjöberg-type modifiers to Hotine’s kernel, albeit only with an isotropic and stationary covariance model for \(\sigma _n^2 \{\delta g_\mathrm{T}\}\).

4.6 Stochastic hybrid (modification S6)

An Ostach–Meissl-type modification (cf. Sects. 3.2 and 4.3) can also be applied to Eq. (49) to give

$$\begin{aligned} \widehat{N_{S6}}{=}\kappa \!\!\int _0^{2\pi } {\int _0^{\psi _0} {H_{S6} (M,L,\psi _0) \delta g \sin \psi \,\mathrm{d}\psi \mathrm{d}\alpha }} {+}c\!\!\sum _{n=2}^L {b_n} \;\delta g_n\nonumber \\ \end{aligned}$$
(57)

with

$$\begin{aligned} H_{S6} (M,L,\psi _0)=H_\mathrm{S5} (M,L)-H_\mathrm{S5} (M,L,\psi =\psi _0) \end{aligned}$$
(58)

and

$$\begin{aligned} \Delta N_{S6} {=}c\!\!\sum _{n=L+1}^\infty \!{\left[ {\int _{\psi _0}^\pi \! {H_{S6} (M,L,\psi _0) P_n (\cos \psi )\sin \psi \,\mathrm{d}\psi }}\!\!\!\right] } \delta g_n\nonumber \\ \end{aligned}$$
(59)

The case of varying the parameters in stochastic modification S5 to achieve a zero crossing of the kernel at the truncation radius \(\psi _0\) is not considered on the grounds of practicality; there are simply too many parameters to trial to make it feasible versus the more pragmatic approach proposed here.

5 Hybrid and band-limited modifications: some suggestions

This section is restricted to a brief discussion of only a few of the options possible, though there are many others; the final choices are left ultimately to the ‘user’. This style of presentation is deliberate to encourage the ‘user’ to experiment with various combinations and permutations so as to tune the data combination for their data sources and area(s) of interest. As alluded to earlier, there is no specific requirement to use any single modification in isolation or to any particular degree or truncation radius, especially when treating them as band-pass filters to reduce errors in the combined solution for the geoid. This applies to Hotine’s, Stokes’s and many other geodetic integrals.

Modern EGMs, particularly those derived from GRACE and/or GOCE satellite gravimetry, are far superior at modelling the low-frequency geoid than terrestrial data alone. As such, it is logical to apply an as-strong-as-possible filter to the terrestrial data, e.g., to the degree that the EGM is considered reliable, so as to rely more upon the low-frequency geoid provided by that EGM. The D1 modification (Sect. 3.1) is the most powerful filter because it removes the low-degree polynomial terms altogether, but the amount of filtering also depends on the truncation radius used (cf. Vaníček and Featherstone 1998). The S1 and S2 modifications (Sects. 4.1 and 4.2) are less effective high-pass filters because they depend on the estimates of \(\sigma _n^2 \{\delta g_\mathrm{T}\},\) noting that if \(\sigma _n^2 \{\delta g_\mathrm{T}\}=0,\) they degenerate to the D1 modification.

One suggested strategy, but only in this author’s opinion (cf. Featherstone 2003a), is to use the D1 modification for the low degrees where the satellite-only EGM is superior to terrestrial data, then use other modifications in the bands where the satellite-only EGM starts to deteriorate, e.g., because of the attenuation of gravitation at satellite altitude. Assuming that GRACE and/or GOCE static gravity field models (e.g., Pail et al. 2011) are better than terrestrial gravity data below some degree L1, but start to deteriorate beyond this to degree \({L2}\,(\le L)\), e.g., the hybrid band-limited kernels from Eqs. (12), (39) and (41) can be combined to give

$$\begin{aligned} H_{\mathrm{D1}+\mathrm{S1}}&= H-\sum _{n=0}^{L1} {\frac{2n+1}{n+1}} \;P_n (\cos \psi )\nonumber \\&+\sum _{n=L1+1}^{L2} {\frac{2n+1}{n+1}} \;w_n P_n (\cos \psi ) \end{aligned}$$
(60)

Likewise, hybrid band-limited versions of the more sophisticated modifiers can be implemented together with the D1 modifier, which can reduce the truncation and other errors. Just as two other examples, combining Eqs. (12), (24) and (50) gives

$$\begin{aligned} H_{\mathrm{D1}+\mathrm{D5}}&\!=\!&H-\!\sum _{n=0}^{L1}\! {\frac{2n\!+\!1}{n\!+\!1}} \;P_n (\cos \psi )-\!\!\sum _{n=L1\!+\!1}^{L2}\! {\frac{2n\!+\!1}{2}} \;h_n P_n (\cos \psi )\quad \end{aligned}$$
(61)
$$\begin{aligned} H_{\mathrm{D1}+\mathrm{S5}}&\!=\!&H-\!\!\sum _{n=0}^{L1}\! {\frac{2n\!+\!1}{n\!+\!1}} P_n (\cos \psi )-\!\!\sum _{n=L1\!+\!1}^{L2}\! {\frac{2n\!+\!1}{2}} \;s_n P_n (\cos \psi )\quad \end{aligned}$$
(62)

The above three examples can be extended or simplified depending upon one’s confidence in the terrestrial gravity error spectra, say where the S5 modification is applied in bands where the error spectra are known and the D5 modification applied where they are not. Naturally, there are many more and alternative options than suggested here.

Figure 2 shows two examples of the hybrid band-limited modified Hotine kernels. Equation (60) is calculated using L1 = 50, where polynomial terms are removed completely in the band \(0\le n\le 50\), and the Wenzel-type modifier (S1) is computed to L2 = 2,160 (i.e., \(51\le n\le 2,160\)). EGM2008 (Pavlis et al. 2012) was used to provide the error degree variances of the EGM and \(\sigma _n^2 \{\delta g_\mathrm{T}\}\) are very crudely assumed to be 0.01 mGal\(^{2}\) for all degrees in a band-limited implementation of Eq. (38). The example presented for Eq. (61) also removes polynomials to L1 = 50 and applies a band-limited Vaníček–Kleusberg-type modifier (D5) to L2 = 100 (i.e., \(51\le n\le 100\)) computed for a spherical cap radius of \(\psi _0 =10^{\circ }\).

Fig. 2
figure 2

Solid line A hybrid band-limited D1+S1 modified Hotine kernel (Eq. 60) for \(L1=50\) and \(L2=2,160\) with \(\sigma _n^2 \{\delta g_{\mathrm{EGM}}\}\) computed from EGM2008 and an assumed \(\sigma _n^2 \{\delta g_\mathrm{T}\}\) of 0.01 mGal\(^{2}\); dotted line A hybrid band-limited D1 + S5-modified Hotine kernel (Eq. 61) for L1 = 50 and L2 = 100 and \(\psi _0=10^{\circ }\)

Full size image

As well as using band-limited modifiers, it is also possible to truncate the modified kernel to some degree L3 that is commensurate with the spatial resolution of the data (cf. Colombo 1977), and which can also avoid aliasing of high-frequency errors (cf. Kern et al. 2003; Novák et al. 2003). Thus, Eqs. (12), (60), (61) and (62) are adapted to omit the spherical Hotine kernel that contains infinite degrees (Eq. 5) and instead evaluate the kernel in bands that are driven by the perceived reliability of the EGM (L1), the types of modifications selected (L2) and the spatial resolution of the terrestrial gravity disturbances (L3). The overbar is used to distinguish these as the band-limited hybrid kernels.

$$\begin{aligned} \bar{{H}}_\mathrm{D1}&= \sum _{n=L1+1}^{L3} {\frac{2n+1}{n+1}} \;P_n (\cos \psi )\end{aligned}$$
(63)
$$\begin{aligned} \bar{{H}}_{\mathrm{D1}+\mathrm{S1}}&= \sum _{n=L1+1}^{L2} {\frac{2n+1}{n+1}} \;w_n P_n (\cos \psi )+\sum _{n=L2+1}^{L3} {\frac{2n+1}{n+1}} \;P_n (\cos \psi )\end{aligned}$$
(64)
$$\begin{aligned} \bar{{H}}_{\mathrm{D1}+\mathrm{D5}}&= \sum _{n=L1+1}^{L2} {\frac{2n+1}{2}} \;\;h_n P_n (\cos \psi )+\sum _{n=L2+1}^{L3} {\frac{2n+1}{n+1}} \;P_n (\cos \psi )\end{aligned}$$
(65)
$$\begin{aligned} \bar{{H}}_{\mathrm{D1}+\mathrm{S5}}&= \sum _{n=L1+1}^{L2} {\frac{2n+1}{2}} \;\;s_n P_n (\cos \psi )+\sum _{n=L2+1}^{L3} {\frac{2n+1}{n+1}} \;P_n (\cos \psi ) \end{aligned}$$
(66)

Naturally, multiple bands and modifiers are possible by combining the above; again, the choices are left to the ‘user’ depending on their data source(s) and area(s) of interest.

The rate of convergence of the truncation error can be accelerated by setting the Hotine kernels to zero at the truncation radius, which removes the discontinuity. This zero can be achieved by simple subtraction or appropriate choices of L1, L2, L3, \(w_n, h_n\) and/or \(s_n\), but noting that some require iteration and lessen the amount of control that the ‘user’ has over the filtering properties of the modifications. Another way to remove the discontinuity is to taper the kernel (cf. Forsberg et al. 2003), though this was not presented as a means to accelerate the convergence, but instead to avoid spectral discontinuities, which can lead to Gibbs fringing when transforming the kernel from the spatial to the spectral domain in FFT geoid computations.

In the most generic form that can be applied to any of the modified Hotine kernels, this is

$$\begin{aligned} H_*=\sum _{n=0}^{L3} {\alpha _n \frac{2n+1}{n+1}} \;m_n P_n (\cos \psi ) \end{aligned}$$
(67)

where \(H_*\) is the modified Hotine kernel, \(m_n\) are the modification coefficients (e.g., 1, \(w_n, h_n\) or \(s_n\)), and the tapering \(\alpha _n\) can be implemented by the following, but other methods of tapering can be used according to the user’s preference

$$\begin{aligned} \alpha _n = \left\{ {{\begin{array}{c@{\quad }cc} 1&{} \mathrm{for}&{} {0\le n<l_1} \\ {{(L3-n)}/{(L3-l_1)}}&{} \mathrm{for}&{} {l_1 \le n\le L3} \\ 0&{} \mathrm{for}&{} {n>L3} \\ \end{array}}} \right. \end{aligned}$$
(68)

where \(l_1\) is chosen by the ‘user’. Tapering could also be used if a combination of the modified kernels causes other discontinuities.

6 Closing remarks

As more and more terrestrial (land, marine and airborne) gravity observations are coordinated by GNSS, gravity disturbances are becoming available for regional geoid computation via Hotine’s integral. Until large surface areas are covered by these types of observations, the truncation error will be larger than desired. As such, this article has presented a suite of deterministic, stochastic, hybrid and band-limited modifications that can be trialled with gravity disturbances for regional geoid computation. Many are driven by modifications already well established for Stokes’s integral and gravity anomalies, some are new, and some can be applied back to Stokes’s and other geodetic integrals, particularly the band-limited and hybrid variants. However, the motivation for this work is not only the reduction of the truncation error, but also the optimal combination and filtering of the heterogeneous GNSS-based gravity data sources available for regional geoid computation.

These various modifications to Hotine’s kernel have been presented in a deliberately non-prescriptive manner so that the ‘user’ has total freedom to experiment with their combinations and permutations, provided that terms are not omitted, which will result in an incomplete spectral representation of the geoid model. It is suggested, but not necessarily recommended, that (1) D1 modifications (Sect. 3.1) are applied routinely because of the superior data now being provided by GRACE and/or GOCE EGMs, (2) the stochastic modifiers can be applied when the ‘user’ is confident with estimates of the error spectra of the data or wishes consider a stochastic interpretation of the gravity field, and (3) the deterministic modifiers can be applied when there is no reliable information of the error properties of the data or the ‘user’ does not wish to consider a stochastic interpretation of the gravity field.