Terrain Corrections in Gravity and Gradiometry

Abstract

Before the computation of short-wavelength and long-wavelength components of the geoid undulations from terrestrial data and the two latest satellite missions, i.e. gravity (GRACE mission) and gradiometry (GOCE mission) measurements, the terrain corrections must be determined. Since the corrections enter the first of the three steps of the Remove-Compute-Restore (RCR) procedure for applying Stokes’s integral, this study focuses on determining these corrections. Formulation of the effects introduced and the effects are computed over high elevated topography in Ireland using Helmert’s second condensation method. Finally, the effects of topography on geoid height determinations are presented.

1. Introduction

The geoid is defined as an equipotential surface along which the Earth’s gravity potential (W) is constant and equal to a reference value $W_0.$ This datum is chosen such that the geoid coincides with a mean level of the oceans and can be mathematically extended over the continents. As a result of the unequal distribution of masses in the Earth’s interior, the geoid is irregularly shaped. It describes the figure of the Earth by a physical quantity, the gravity potential, in contrast to the idealized geometrical figure of a reference ellipsoid. The separations between the two surfaces are called the geoid undulation N, or geoidal heights.

The local gravity potential, $W_{\mathit{\text{local}}},$ value is derived from gravimetric $\left(N_{\mathit{\text{grav}}}\right),$ and geometric $\left(N_{\mathit{geo}}\right)$ geoid undulations. The geometric geoid undulations are obtained from GNSS and leveling data. In contrast, the gravimetric geoid undulations are computed based upon their long-wavelength components from a Global Geopotential Model (GGM), and their short-wavelength components using the terrestrial data through the so-called Remove-Compute-Restore (RCR) approach from gravity measurements.¹

1.1 The Remove-Compute-Restore procedure (RCR)

The RCR procedure is a method that fulfils Stokes’ requirements² for computing the geoid from:

• Terrestrial gravity measurements or,

• Satellites gradiometry³ measurements (e.g. GRACE, GOCE).

A gravimeter measures magnitude of gravity,

$$g=\mid gradV\mid ,$$E1

whereas the gradiometer measures the components of the gradient of gravity. We will focus on the rr component of the gradiometric tensor,

$$V_{\mathit{rr}}={\left(gradgradV\right)}_{\mathit{rr}}.$$E2

1.1.1 Removing the gravitational effect of the residual topographical masses

Subtracting the gravitational effect of the residual topographical masses, $\mathit{\delta V},$ from the actual anomalous gravitational potential T creates potential $T^h$ that is harmonic outside the geoid,

$$T^h=T-\mathit{\delta V}.$$E3

The gravity attraction of the residual topographical masses is then

$$\mathit{\delta A}≔\frac{\partial \mathit{\delta V}}{\partial r}$$E4

at the point of the gravity measurements, and

$$\mathit{\delta E}≔\frac{{\partial }^2\mathit{\delta V}}{\partial r^2}$$E5

at point of gradiometric measurements.

To make a potential harmonic in a space above the geoid, these effects have to be calculated and removed from the observations:

$$\mathrm{\Delta }{g}^h=\mathrm{\Delta }{g}^{\mathit{obs}}-\mathit{\delta A}$$E6

$$\mathrm{\Delta }{V}_{\mathit{rr}}^h=\mathrm{\Delta }{V}_{\mathit{rr}}^{\mathit{obs}}-\mathit{\delta E}.$$E7

1.1.2 Computation of the residual geoid (co-geoid)

The first step of the computation is continuing $\mathrm{\Delta }{g}^h$ and $\mathrm{\Delta }{V}_{\mathit{rr}}^h$ from the surface or satellite elevation to the geoid. This is usually performed by a harmonic downward continuation (DWC) method [1].

Assuming that the data have been continued down to the geoid, the two basic requirements are met and the computation of the geoidal heights can now be carried out. However, the computation gives the height $N^h$ of the equipotential surface of $T^h$ that is called co-geoid.

$$N^h=\frac{R}{4{\mathit{\pi \gamma }}_Q}{\int }_{{\mathrm{\Omega }}_0}\mathrm{\Delta }{g}^{h}S\left(\psi \right)d\mathrm{\Omega },$$E8

where $\psi$ is the spherical solid angle between the computation point and an integration point and $S\left(\psi \right)$ is the Stokes’s function ([2], p. 104) integrated over the full solid angle ${\mathrm{\Omega }}_0$

$$S\left(\psi \right)=\frac{1}{sin\frac{\psi }{2}}-6sin\frac{\psi }{2}+1-5cos\psi -3cos\psi ln\left(sin\frac{\psi }{2}+{sin}^2\frac{\psi }{2}\right).$$E9

Likewise, the heights of the co-geoid can be determined from gradiometric data as

$$N^h=\frac{R^2}{4{\mathit{\pi \gamma }}_Q}{\int }_{{\mathrm{\Omega }}_0}\mathrm{\Delta }{V}_{\mathit{rr}}^h{K}_{\mathit{rr}}\left(\psi \right)d\mathrm{\Omega },$$E10

where the kernel $K_{\mathit{rr}}\left(\psi \right)$ is the Green’s function as given by Martinec [3]

$$K_{\mathit{rr}}\left(\psi \right)=-3+6sin\frac{\psi }{2}+\left(1-3cos\psi \right)ln\left(\frac{sin\frac{\psi }{2}}{1+sin\frac{\psi }{2}}\right).$$E11

1.1.3 Adding the contribution of the topography to the solution

Finally, the heights of the geoid will be obtained by adding the difference $\mathit{\delta N}$ shown in Figure 1 which is called the primary indirect effect on geoid.

Figure 1.

co-geoid and geoidal heights.

To find the expression for $\mathit{\delta N},$ we start with geoidal height N is derived from Bruns’ formula

$$N=\frac{T}{{\gamma }_Q},$$E12

where ${\gamma }_Q$ is the normal gravity on a reference ellipsoid. Substituting for T from Eq. (60) we obtain Eq. (3)

$$N=\frac{T^h+\mathit{\delta V}}{{\gamma }_Q}.$$E13

Applying Bruns’ formula to the undulation of the co-geoid, $N^h=T^h/{\gamma }_Q,$ the geoidal height N is

$$N=N^h+\mathit{\delta N},$$E14

where

$$\mathit{\delta N}=\frac{\mathit{\delta V}}{{\gamma }_Q},$$E15

and $\mathit{\delta V}$ must be taken at a point on the geoid $P_g$.

GRACE and GOCE are the two latest satellite missions for precise, long-wavelength geoid determination from gravity and gradiometry measurements, respectively. Prior processing GRACE and GOCE data, the terrain corrections have to be determined.

The RCR procedure for gravity and gravimetry is summarized in Figure 2.

Figure 2.

Summary of the RCR procedure for gravity and gradiometry.

It shows the motivation for correcting the data by the effects of the topographical masses. These are the direct topographical effects $\mathit{\delta A}$ computed at surface or satellite altitude, the direct topographical effect on gradiometry $\mathit{\delta E}$ computed at satellite altitude and the indirect topographical effect on geoid $\mathit{\delta N}$ computed on the geoid.

In the intervening period of the theory of physical geodesy, a wide range of methods have been developed, e.g. Airy-isostatic reduction method [4], Residual Terrain Model (RTM) scheme [5], Helmert’s first and second method of condensation [4, 6, 7, 8]. The choice of each method is area-dependent, and the correlation between the spatial resolution of DEMs, and the elevation of computation points or data surrounding the computation points show the most suitable spatial resolution of DEMs that provide intended geodetic accuracies. This is investigated numerically over high elevated topography in Ireland using Helmert’s second condensation method.

The next section is devoted to express them mathematically and prepare them for numerical realization.

2. Topographical effects

The term terrain effect is used to express the gravitational effects of topographical masses on gravity anomalies, deflection of the vertical and other observed quantity. It can be classified according to a location of anomalous masses. Topographical effects are the direct influence of the visible topography in mountainous areas; isostatical effects account for a hypothesized isostatic compensation, whilst the residual terrain model (RTM) effects account for short-wavelength topographic irregularities referring topographic elevations to a smooth mean elevation surface, which may be defined, for instance, by spherical harmonic expansion of topographic heights.

2.1 Topographical masses and the Bouguer plate

2.1.1 Topographical masses

The topographical masses are the masses outside the geoid and below the topographical surface. The gravitational potential $V^t$ generated by the topographical masses is

$$V^t\left(r,\mathrm{\Omega }\right)=G{\int }_{{\mathrm{\Omega }}_0}{\int }_{r^{\prime }=r_g\left({\mathrm{\Omega }}^{\prime }\right)}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}\frac{\varrho \left(r^{\prime },{\mathrm{\Omega }}^{\prime }\right)}{L\left(r,\psi ,r^{\prime }\right)}{r^{\prime }}^{2}d{r}^{\prime }d{\mathrm{\Omega }}^{\prime }.$$E16

where G is the Newton’s gravitational constant, $G=6.67\cdot {10}^{-11}{\mathrm{m}}^3.{\mathrm{kg}}^{-1}.{\mathrm{s}}^{-2},$ $\varrho \left(r^{\prime },{\mathrm{\Omega }}^{\prime }\right)$ is the mass density inside the Earth’s interior located at $P^{\prime }\left(r^{\prime },{\mathrm{\Omega }}^{\prime }\right),L\left(r,\psi ,r^{\prime }\right)$ is the distance between P and $P^{\prime }$ and $\psi$ is the angular distance between the geocentric directions $\mathrm{\Omega }=\left(\vartheta ,\phi \right)$ and ${\mathrm{\Omega }}^{\prime }=\left({\vartheta }^{\prime },{\phi }^{\prime }\right),$ see Figure 3, i.e. $\psi \left(\mathrm{\Omega },{\mathrm{\Omega }}^{\prime }\right)$ is the spherical solid angle between P and $P^{\prime },$

Figure 3.

Spherical coordinates of the computation point $P\left(r,\mathrm{\Omega }\right),$ an integration point $P^{\prime }\left(r^{\prime },{\mathrm{\Omega }}^{\prime }\right),$ the distance L and solid angle $\psi$ between P and $P^{\prime }$.

$$cos\psi =cos\vartheta cos{\vartheta }^{\prime }+sin\vartheta sin{\vartheta }^{\prime }cos\left(\phi -{\phi }^{\prime }\right).$$E17

The argument notation in $L\left(r,\psi ,r^{\prime }\right)$ is used to emphasize the fact that L depends on radial distances r and $r^{\prime },$ and the angular distance $\psi .$

$$L\left(r,\psi ,r^{\prime }\right)≔\sqrt{{r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2}.$$E18

To abbreviate notations, we introduced the symbol $\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)$ for an indefinite radial integral of the Newton kernel,

$$\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)≔{\int }_{r^{\prime }}\frac{{r^{\prime }}^2}{L\left(r,\psi ,r^{\prime }\right)}d{r}^{\prime }.$$E19

Assuming that the density of the topographical masses does not vary in radial direction, that is $\varrho \left(r^{\prime },{\mathrm{\Omega }}^{\prime }\right)=\varrho \left({\mathrm{\Omega }}^{\prime }\right),$ and substituting Eq. (19) in Eq. (16), the Newton’s volume integral for the gravitational potential $V^t$ becomes

$$V^t\left(r,\mathrm{\Omega }\right)=G{\int }_{{\mathrm{\Omega }}_0}\varrho \left({\mathrm{\Omega }}^{\prime }\right){\left.\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=r_g\left({\mathrm{\Omega }}^{\prime }\right)}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}d{\mathrm{\Omega }}^{\prime }.$$E20

2.1.2 Bouguer plate

The Bouguer plate, used as an approximate model in gravity and gravity anomaly computations accounts for the bulk of topographical effects.

In Cartesian geometry Figure 4a, the topography around the gravity station P is approximated by an infinite plate of thickness $H_P$ and the masses between the geoid and the Earth’s surface have a constant density $\varrho$ equal to mean topographical density ${\varrho }_0=2670\mathrm{kg}.{\mathrm{m}}^{-3}$.

Figure 4.

The Bouguer plate in Cartesian and spherical geometry. (a) Infinite Bouguer plate. (b) Spherical Bouguer layer.

In spherical geometry Figure 4b, the Bouguer plate is regarded as a spherical layer of thickness $H_P$ and density ${\varrho }_0.$ The gravitational potential of the spherical Bouguer layer is

$$V^B\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}{\left.\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=r_g\left(\mathrm{\Omega }\right)}^{r_t\left(\mathrm{\Omega }\right)}d{\mathrm{\Omega }}^{\prime }.$$E21

For evaluating this integral, the geoid is approximated by a sphere of radius R

$$r_g\left(\mathrm{\Omega }\right)=R\mathrm{so}{r}_t\left(\mathrm{\Omega }\right)=R+H\left(\mathrm{\Omega }\right).$$E22

For R we consider the mean radius of the Earth, R = 6371 km. The integral Eq. (21) can be evaluated analytically, e.g. Wichiencharoen [9],

$$V^B\left(r,\mathrm{\Omega }\right)=\left\{\begin{array}{ll}4\mathit{\pi G}{\varrho }_0\frac{1}{r}\left[R^{2}H\left(\mathrm{\Omega }\right)+{\mathit{RH}}^2\left(\mathrm{\Omega }\right)+\frac{1}{3}{H}^3\left(\mathrm{\Omega }\right)\right]& ,r\ge r_t\left(\mathrm{\Omega }\right),\\ 2\mathit{\pi G}{\varrho }_0\left[{\left(R+H\left(\mathrm{\Omega }\right)\right)}^2-\frac{2}{3}\frac{R^3}{r}-\frac{1}{3}{r}^2\right]& ,R\le r\le r_t\left(\mathrm{\Omega }\right),\\ 4\mathit{\pi G}{\varrho }_0\left[\mathit{RH}\left(\mathrm{\Omega }\right)+\frac{1}{2}{H}^2\left(\mathrm{\Omega }\right)\right]& ,r\le R.\end{array}\right.$$E23

2.1.3 Terrain roughness

Since the actual Earth’s surface deviates from the Bouguer sphere, there are deficiencies and abundances of topographical masses with respect to the mass of the Bouguer plate Figure 5. These contribute to the topographical potential $V^t$ through the term $V^R$ as

$$V^t\left(r,\mathrm{\Omega }\right)=V^B\left(r,\mathrm{\Omega }\right)+V^R\left(r,\mathrm{\Omega }\right).$$E24

Figure 5.

Roughness of the terrain.

The terrain roughness term $V^R$ is expressed by the Newton integral

$$V^R\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}\left[{\left.\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}-{\left.\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{r_t\left(\mathrm{\Omega }\right)}\right]d{\mathrm{\Omega }}^{\prime }.$$E25

2.2 Compensated masses and the Helmert condensation layer

2.2.1 Compensation of the gravitational effects of topographical masses

The equipotential surfaces of $V^t$ undulate by several hundreds of meters with respect to a level ellipsoid. The fact that the known undulations of the geoid are much smaller than those induced by potential $V^t$ indicates that there must exist a compensation mechanism which reduces the gravitational effect of topographical masses [10].

As the masses are compensated in some way [2], we can introduce the gravitational potential of compensated masses $V^c$ as an approximation of the topographical potential $V^t.$ The difference between $V^t$ and $V^c$ is called the residual topographical potential $\mathit{\delta V}:$

$$\mathit{\delta V}≔V^t-V^c.$$E26

Two extremely idealized isostatic compensation models (see Figure 6) were proposed to approximate the effect of topographical abundances from surface gravity observations.

Figure 6.

topography and compensation layers of Pratt-Hayford and Airy-Heiskanen models ([2], Fig 3.16).

The Pratt-Hayford model was outlined by J.H Pratt in 1854 and put into a mathematical form by J.F Hayford. According to Pratt, the mountains have risen from the underground somewhat like a fermenting dough [2, sect.3]. The topographical masses are compensated by varying density distribution within the layer of a constant thickness D = 100 km under the geoid and a density ${\varrho }_c\left(\mathrm{\Omega }\right)$.

The Airy-heiskanen model, proposed by G.B Airy in 1855 and formulated for geodetic purposes by W.A Heiskanen, assumes that the mountains are floating on a fluid lava of higher known density ${\varrho }_1$ (somewhat like an iceberg floating on water), so that the higher the mountain, the deeper it sinks [2, sect.3]. The topographical masses are compensated by varying thickness $t\left(\mathrm{\Omega }\right)$ of a compensation layer which surface is situated at $T=30$ km deep. The density of the compensation layer is considered constant and equal to the density difference ${\varrho }_c={\varrho }_1-{\varrho }_0\left(>0\right).$

2.2.2 Helmert condensation layer

In the limiting case, the topographical masses may be compensated by a thin mass layer located on the geoid (somewhat like a glass sphere made over very thin but very robust glass [2]). As shown Figure 7, the topographical masses are condensed as a surface mass layer on the geoid. This kind of compensation is called Helmert 2nd condensation [11] approximating the actual potential of the topographical masses $V^t$ by the potential of a single layer $V^c$ described by Newton’s surface integral as

Figure 7.

Helmert condensation layer of density $\sigma$.

$$V^c\left(r,\mathrm{\Omega }\right)={\mathit{GR}}^2{\int }_{{\mathrm{\Omega }}_0}\sigma \left({\mathrm{\Omega }}^{\prime }\right)L^{-1}\left(r,\psi ,R\right)d{\mathrm{\Omega }}^{\prime },$$E27

where $\sigma \left(\mathrm{\Omega }\right)$ is a surface density of Helmert’s condensation layer, and $L^{-1}$ is the reciprocal distance 1/L. Analogously to Eq. (24), we can rewrite Eq. (27) as

$$V^c\left(r,\mathrm{\Omega }\right)=V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)+V^{\sigma ,R}\left(r,\mathrm{\Omega }\right),$$E28

where

$$V^{\sigma ,R}\left(r,\mathrm{\Omega }\right)={\mathit{GR}}^2{\int }_{{\mathrm{\Omega }}_0}\left[\sigma \left({\mathrm{\Omega }}^{\prime }\right)-\sigma \left(\mathrm{\Omega }\right)\right]L^{-1}\left(r,\psi ,R\right)d{\mathrm{\Omega }}^{\prime }.$$E29

The symbol $V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)$ denotes the gravitational potential of a spherical layer with density $\sigma \left(\mathrm{\Omega }\right)$ and radius R,

$$V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)≔{\mathit{GR}}^2\sigma \left(\mathrm{\Omega }\right){\int }_{{\mathrm{\Omega }}_0}{L}^{-1}\left(r,\psi ,R\right)d{\mathrm{\Omega }}^{\prime },$$E30

that may be evaluated analytically by

$$V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)=\left\{\begin{array}{ll}4\mathit{\pi G\sigma }\left(\mathrm{\Omega }\right)\frac{R^2}{r}& ,r>R,\\ 4\mathit{\pi G\sigma }\left(\mathrm{\Omega }\right)R& ,r\le R.\end{array}\right.$$E31

The condensation density $\sigma \left(\mathrm{\Omega }\right)$ can be chosen in a variety of ways. For this study, it will be chosen as in Martinec and Vaníček and Martinec [10] such that the principle of conservation of topographical masses [9] is respected. Assuming that

$$V^B\left(r,\mathrm{\Omega }\right)=V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)\text{for}r=r_t\left(\mathrm{\Omega }\right),$$E32

and substituting Eqs.(23) and (31) into (32), we find

$$\sigma \left(\mathrm{\Omega }\right)={\varrho }_0\tau \left(\mathrm{\Omega }\right),$$E33

with

$$\tau \left(\mathrm{\Omega }\right)=H\left(\mathrm{\Omega }\right)\left(1+\frac{H\left(\mathrm{\Omega }\right)}{R}+\frac{H^2\left(\mathrm{\Omega }\right)}{3R^2}\right).$$E34

Thus, Eq. (29) will be written

$$V^{\sigma ,R}\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{R}^2{\int }_{{\mathrm{\Omega }}_0}\left[\tau \left({\mathrm{\Omega }}^{\prime }\right)-\tau \left(\mathrm{\Omega }\right)\right]L^{-1}\left(r,\psi ,R\right)d{\mathrm{\Omega }}^{\prime }.$$E35

2.3 Indirect topographical effect

2.3.1 Indirect topographical effect on potential

The primary indirect topographical effect is the residual potential $\mathit{\delta V}=V^t-V^c$ evaluated at a point $P_g\left(R,\mathrm{\Omega }\right)$ on the geoid. Considering Eq. (24) for the topographical potential $V^t$, Eq. (28) for the condensation potential $V^c$ and replacing r by R, we may split $\mathit{\delta V}$ into two terms as

$$\mathit{\delta V}\left(R,\mathrm{\Omega }\right)=\delta V^B\left(R,\mathrm{\Omega }\right)+\delta V^R\left(R,\mathrm{\Omega }\right),$$E36

where the Bouguer term and the terrain roughness term are respectively given by

$$\delta V^B\left(R,\mathrm{\Omega }\right)=V^B\left(R,\mathrm{\Omega }\right)-V^{\sigma ,B}\left(R,\mathrm{\Omega }\right),$$E37

$$\delta V^R\left(R,\mathrm{\Omega }\right)=V^R\left(R,\mathrm{\Omega }\right)-V^{\sigma ,R}\left(R,\mathrm{\Omega }\right).$$E38

The subtraction of Eqs. (23) and (31) at $r\le R$ leads to

$$\delta V^B\left(R,\mathrm{\Omega }\right)=-2\mathit{\pi G}{\varrho }_{0}H\left(\mathrm{\Omega }\right)\left(1+\frac{2}{3}\frac{H\left(\mathrm{\Omega }\right)}{R}\right),$$E39

and the subtraction of Eqs. (25) and (35) gives

$$\begin{array}{l}\delta V^R\left(R,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}\left[{\left.\widetilde{L^{-1}}\left(R,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}-{\left.\widetilde{L^{-1}}\left(R,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{r_t\left(\mathrm{\Omega }\right)}\right.\\ \left.-R^2\left[\tau \left({\mathrm{\Omega }}^{\prime }\right)-\tau \left(\mathrm{\Omega }\right)\right]L^{-1}(R,\psi ,R)\right]d{\mathrm{\Omega }}^{\prime }.\end{array}$$E40

Therefore, substituting Eqs. (39) and (40) in Eq. (36), we obtain the expression of the primary topographical indirect effect on potential:

$$\begin{array}{l}\mathit{\delta V}\left(R,\mathrm{\Omega }\right)=-2\mathit{\pi G}{\varrho }_0{H}^2\left(\mathrm{\Omega }\right)\left(1+\frac{2}{3}\frac{H\left(\mathrm{\Omega }\right)}{R}\right)\\ +G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}\left[{\left.\widetilde{L^{-1}}\left(R,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}-{\left.\widetilde{L^{-1}}\left(R,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{r_t\left(\mathrm{\Omega }\right)}\right.\\ \left.-R^2\left[\tau \left({\mathrm{\Omega }}^{\prime }\right)-\tau \left(\mathrm{\Omega }\right)\right]L^{-1}(R,\psi ,R)\right]d{\mathrm{\Omega }}^{\prime }.\end{array}$$E41

The unit of this effect is ${\mathrm{m}}^2.{\mathrm{s}}^{-2}.$

2.3.2 Indirect topographical effect on the geoid

To correct the geoidal heights N by this effect, $\mathit{\delta V}$ is divided by the normal gravity ${\gamma }_Q$ (see Eq. (15)),

$$\mathit{\delta N}\left(R,\mathrm{\Omega }\right)=\frac{\mathit{\delta V}\left(R,\mathrm{\Omega }\right)}{{\gamma }_Q},$$E42

which is the primary indirect effect on the geoid. Where the normal gravity can be taken as the mean value of the gravity of the Earth ${\gamma }_Q=9.81\mathrm{m}.{\mathrm{s}}^{-2}.$ The unit of $\mathit{\delta N}$ is meters.

2.3.3 The secondary indirect topographical effect (SITE) on gravity

This effect is expressed by means of the Primary Indirect Topographical Effect (PITE) on gravity, $\delta V_{P_g},$ at a point on the geoid multiplied by 2/R,

$$\mathit{\delta S}\left(\mathrm{\Omega }\right)=\frac{2}{R}\delta V_{P_g}\left(\mathrm{\Omega }\right).$$E43

The unit of this effect is mGal.

Notice that the radial derivative of the Newton surface and volume integrals, which is required for computing $\mathit{\delta A}\left(\mathrm{\Omega }\right),$ is

$$\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}={\int }_{r^{\prime }}{r^{\prime }}^2\frac{\partial L^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r}d{r}^{\prime },$$E44

and using the derivation of ${\left(u^{\alpha }\right)}^{\prime }=\alpha u^{\prime }{u}^{\alpha -1}$ in Eq. (18),

$$\frac{\partial L^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r}=-\frac{r-r^{\prime }cos\psi }{{\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{3/2}},$$E45

where substituting the expression (45) in (44), the analytical expression is given by,

$$\begin{array}{l}\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}=\left(3r^{2}cos\psi +r\left(r^{\prime }-6cos{\psi }^2{r}^{\prime }\right)+cos\psi {r^{\prime }}^2\right)L^{-1}\left(r,\psi ,r^{\prime }\right)\\ +r\left(3cos{\psi }^2-1\right)ln\mid L\left(r,\psi ,r^{\prime }\right)-rcos\psi +r^{\prime }\mid +C.\end{array}$$E46

2.4 Direct topographical effect (DTE)

2.4.1 Direct topographical effect on surface gravity

Differentiating the residual topographical potential $\mathit{\delta V}$ with respect to r and evaluating the result at the point on the topography $P_t\left(r_t\left(\mathrm{\Omega }\right),\mathrm{\Omega }\right)$, we obtain the gravitational attraction caused by the direct topographical effect on surface gravity:

$$\mathit{\delta A}\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial \mathit{\delta V}\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)}.$$E47

Substituting for $\mathit{\delta V}$ from Eq. (26), we can write

$$\mathit{\delta A}\left(r,\mathrm{\Omega }\right)=A^t\left(r,\mathrm{\Omega }\right)-A^c\left(r,\mathrm{\Omega }\right),$$E48

where

$$A^t\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial V^t\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)}\text{and}{A}^c\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial V^c\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)},$$E49

are the radial components of the gravitational attraction induced by the topographical and compensated masses a the point on the Earth’s surface, respectively. Considering Eq. (24) for $A^t$ and Eq. (28) for $A^c$, the attraction change $\mathit{\delta A}$ may be split into two terms

$$\mathit{\delta A}\left(r,\mathrm{\Omega }\right)=\delta A^B\left(r,\mathrm{\Omega }\right)+\delta A^R\left(r,\mathrm{\Omega }\right),$$E50

where $\delta A^B$ represents the Bouguer term and $\delta A^R$ the terrain roughness term.

Let us start with the determination of the Bouguer term

$$\delta A^B\left(r,\mathrm{\Omega }\right)=A^B\left(r,\mathrm{\Omega }\right)-A^{\sigma ,B}\left(r,\mathrm{\Omega }\right),$$E51

where

$$A^B\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial V^B\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)}\text{and}{A}^{\sigma ,B}\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)}.$$E52

Subsituting the radial derivative of Eq. (23) at $r\ge r_t\left(\mathrm{\Omega }\right)$ for $A^B$

$$\frac{\partial V^B\left(r,\mathrm{\Omega }\right)}{\partial r}=-4\mathit{\pi G}{\varrho }_0\frac{1}{r^2}\left[R^{2}H\left(\mathrm{\Omega }\right)+{\mathit{RH}}^2\left(\mathrm{\Omega }\right)+\frac{1}{3}{H}^3\left(\mathrm{\Omega }\right)\right],$$E53

similarly, the radial derivative of Eq. (31) at $r>R$ for $A^{\sigma ,B}$

$$\frac{\partial V^{\sigma ,B}\left(r,\mathrm{\Omega }\right)}{\partial r}=-4\mathit{\pi G\sigma }\left(\mathrm{\Omega }\right)\frac{R^2}{r^2}.$$E54

Using Eqs. (33)–(34) for $\sigma \left(\mathrm{\Omega }\right),$ we obtain $\delta A^B\left(r,\mathrm{\Omega }\right)=0.$ Consequently, the direct topographical effect $\mathit{\delta A}$ only consists of the terrain roughness contribution $\delta A^R$ that is

$$\delta A^R\left(r,\mathrm{\Omega }\right)=A^R\left(r,\mathrm{\Omega }\right)-A^{\sigma ,R}\left(r,\mathrm{\Omega }\right),$$E55

where

$$A^R\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial V^R\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)}\text{and}{A}^{\sigma ,R}\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial V^{\sigma ,R}\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_t\left(\mathrm{\Omega }\right)}.$$E56

Substituting the radial derivations of Eq. (25) for $A^R$ and Eq. (35) for $A^{\sigma ,R}$ into Eq. (55), we obtain the expression of the direct topographical effect on surface gravity as

$$\begin{array}{l}\mathit{\delta A}\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}\left[-{\left.\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}\right|}_{r^{\prime }=R}^{r_t\left(\mathrm{\Omega }\right)}\right.\\ -{\left.R^2\left[\tau \left({\mathrm{\Omega }}^{\prime }\right)-\tau \left(\mathrm{\Omega }\right)\right]\frac{\partial L^{-1}\left(r,\psi ,R\right)}{\partial r}\right]}_{r=r_t\left(\mathrm{\Omega }\right)}d{\mathrm{\Omega }}^{\prime }.\end{array}$$E57

2.4.2 Direct topographical effect on satellite gravity

Differentiating $\mathit{\delta V}$ with respect to r and evaluating the result at the point of measurement on satellite $\left(r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right),\mathrm{\Omega }\right)$, we obtain the change of the gravitational attraction caused by the direct topographical effect on satellite gravity

$$\mathit{\delta A}\left(r,\mathrm{\Omega }\right)={\left.\frac{\partial \mathit{\delta V}\left(r,\mathrm{\Omega }\right)}{\partial r}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}.$$E58

Analogically to the direct topographical effect on surface gravity Eq. (57), where the radius of the computation point is $r=r_t\left(\mathrm{\Omega }\right)$, the direct topographical effect on satellite gravity will be:

$$\begin{array}{l}\mathit{\delta A}\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}\left[{\left.\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}-{\left.\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}\right|}_{r^{\prime }=R}^{r_t\left(\mathrm{\Omega }\right)}\right.\\ -{\left.R^2\left[\tau \left({\mathrm{\Omega }}^{\prime }\right)-\tau \left(\mathrm{\Omega }\right)\right]\frac{\partial L^{-1}\left(r,\psi ,R\right)}{\partial r}\right]}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}d{\mathrm{\Omega }}^{\prime },\end{array}$$E59

where $r_{\mathit{sat}}\left(\mathrm{\Omega }\right)$, is:

$$r_{\mathit{sat}}\left(\mathrm{\Omega }\right)=R+H_{\mathit{sat}},$$E60

and $H_{\mathit{sat}}$ is the flying altitude of the satellite⁴ which performs the gravity measurements. For instance, the value of $H_{\mathit{sat}}\left(\mathrm{\Omega }\right)$ for the satellite GRACE is $H_{\mathit{sat}}\left(\mathrm{\Omega }\right)=400$ km.

2.4.3 Direct topographical effect on gradiometry

Differentiating $\mathit{\delta V}$ twice with respect to r and evaluating the result at the point of measurement on satellite $\left(r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right),\mathrm{\Omega }\right)$, we obtain the rr component of the gradient of gravity caused by the direct topographical effect on gradiometry and Eq. (5) becomes

$$\mathit{\delta E}\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2\mathit{\delta V}\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}.$$E61

Substituting for the residual topographical potential $\mathit{\delta V}$ from Eq. (26), we can write

$$\mathit{\delta E}\left(r,\mathrm{\Omega }\right)=V_{\mathit{rr}}^t\left(r,\mathrm{\Omega }\right)-V_{\mathit{rr}}^c\left(r,\mathrm{\Omega }\right),$$E62

where

$$V_{\mathit{rr}}^t\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2{V}^t\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}\text{and}{V}_{\mathit{rr}}^c\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2{V}^c\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)},$$E63

are the rr components of the gradiometric tensor induced by the topographical and compensated masses a the point at satellite altitude, respectively. Considering Eq. (24) for $V_{\mathit{rr}}^t$ and Eq. (28) for $V_{\mathit{rr}}^c,\mathit{\delta E}$ may be split into two terms

$$\mathit{\delta E}\left(r,\mathrm{\Omega }\right)=\delta E^B\left(r,\mathrm{\Omega }\right)+\delta E^R\left(r,\mathrm{\Omega }\right),$$E64

where $\delta E^B$ represents the Bouguer term and $\delta E^R$ the terrain roughness term.

Let us start with the determination of the Bouguer term

$$\delta E^B\left(r,\mathrm{\Omega }\right)=V_{\mathit{rr}}^B\left(r,\mathrm{\Omega }\right)-V_{\mathit{rr}}^{\sigma ,B}\left(r,\mathrm{\Omega }\right),$$E65

where

$$V_{\mathit{rr}}^B\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2{V}^B\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}\text{and}{V}_{\mathit{rr}}^{\sigma ,B}\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2{V}^{\sigma ,B}\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}.$$E66

Taking the radial second derivative of Eq. (23) at $r\ge r_t\left(\mathrm{\Omega }\right)$ for $V_{\mathit{rr}}^B$

$$\frac{{\partial }^2{V}^B\left(r,\mathrm{\Omega }\right)}{\partial r^2}=8\mathit{\pi G}{\varrho }_0\frac{1}{r^3}\left[R^{2}H\left(\mathrm{\Omega }\right)+{\mathit{RH}}^2\left(\mathrm{\Omega }\right)+\frac{1}{3}{H}^3\left(\mathrm{\Omega }\right)\right],$$E67

the radial second derivative of Eq. (31) at $r>R$ for $V_{\mathit{rr}}^{\sigma ,B}$

$$\frac{{\partial }^2{V}^{\sigma ,B}\left(r,\mathrm{\Omega }\right)}{\partial r^2}=8\mathit{\pi G\sigma }\left(\mathrm{\Omega }\right)\frac{R^2}{r^3},$$E68

and using Eqs. (33)–(34) for $\sigma \left(\mathrm{\Omega }\right)$, we obtain $\delta E^B\left(r,\mathrm{\Omega }\right)=0$. Consequently, as for the direct topographical effect on gravity $\mathit{\delta A}$, the direct topographical effect on gradiometry $\mathit{\delta E}$ only consists of the terrain roughness contribution $\delta E^R$ that is

$$\delta E^R\left(r,\mathrm{\Omega }\right)=V_{\mathit{rr}}^R\left(r,\mathrm{\Omega }\right)-V_{\mathit{rr}}^{\sigma ,R}\left(r,\mathrm{\Omega }\right),$$E69

where

$$V_{\mathit{rr}}^R\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2{V}^R\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}\text{and}{V}_{\mathit{rr}}^{\sigma ,R}\left(r,\mathrm{\Omega }\right)={\left.\frac{{\partial }^2{V}^{\sigma ,R}\left(r,\mathrm{\Omega }\right)}{\partial r^2}\right|}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}.$$E70

Therefore, differentiating $V^R$ (Eq. (25)) and $V^{\sigma ,R}$ (Eq. (35)) twice with respect to r for $A^R$ and $A^{\sigma ,R}$, respectively, we obtain the expression for the direct topographical effect on gradiometry

$$\begin{array}{l}\mathit{\delta E}\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}\left[{\left.\frac{{\partial }^2\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)}-{\left.\frac{{\partial }^2\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}\right|}_{r^{\prime }=R}^{r_t\left(\mathrm{\Omega }\right)}\right.\\ -{\left.R^2\left[\tau \left({\mathrm{\Omega }}^{\prime }\right)-\tau \left(\mathrm{\Omega }\right)\right]\frac{{\partial }^2{L}^{-1}\left(r,\psi ,R\right)}{\partial r^2}\right]}_{r=r_{\mathit{sat}}\left(\mathrm{\Omega }\right)}d{\mathrm{\Omega }}^{\prime }.\end{array}$$E71

Here $H_{\mathit{sat}}$ is the flying altitude of the satellite which takes the gradiometry measurements. For instance, the satellite GOCE with $H_{\mathit{sat}}\left(\mathrm{\Omega }\right)=250$ km.

2.5 Computations of the topographical effects

2.5.1 The integral Newton kernels for numerical computation

The distance $L\left(r,\psi ,r^{\prime }\right)$, defined Eq. (18), can be written with the form $\sqrt{X}$, where X is the rational function

$$X={\mathit{ax}}^2+\mathit{bx}+c,\text{where}x=1,a=1,b=-2rcos\psi ,c=r^2.$$E72

Thus, we can write the Newton kernel $\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)$ Eq. (19) with the form

$$\int \frac{x^2}{\sqrt{X}}\mathit{dx}.$$E73

To solve this integral, we use the following equations Gradshteyn and Ryzhik [12]

$$\int \frac{x^2\mathit{dx}}{\sqrt{X}}=\left(\frac{x}{2a}-\frac{3b}{4a^2}\right)\sqrt{X}+\frac{3b^2-4\mathit{ac}}{8a^2}\int \frac{\mathit{dx}}{\sqrt{X}}$$E74

$$\int \frac{\mathit{dx}}{\sqrt{X}}=\frac{1}{\sqrt{a}}ln\mid 2\sqrt{\mathit{aX}}+2\mathit{ax}+b\mid +C,$$E75

for $a>0$ and where C is a constant value.

Replacing Eq. (75) and “our” notations in Eq. (74), we obtain the analytical expression:

$$\begin{array}{l}\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)=\frac{1}{2}\left[\left(3rcos\psi +r^{\prime }\right)\sqrt{{r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2}+\right.\\ \left.+r^2\left(3{cos}^2\psi -1\right)ln|2\left(\sqrt{{r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2}-rcos\psi +r^{\prime }\right)|\right]+C\end{array}$$E76

Hence, the analytical expression $\widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)$ from the list of integrals of Gradshteyn and Ryzhik [12]

$$\begin{array}{l}\widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)=\frac{1}{2}\left[\left(3rcos\psi +r^{\prime }\right)L(r,\psi ,r^{\prime })+\right.\\ \left.+r^2\left(3{cos}^2\psi -1\right)ln|L(r,\psi ,r^{\prime })-rcos\psi +r^{\prime }|\right]+C_1\left(r,\psi \right),\end{array}$$E77

where $C_1\left(r,\psi \right)$ a “constant” which may depend only on the variables r and $\psi$. The computation of $\mathit{\delta N}$ uses the expression of $\widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)$.

2.5.2 First radial derivative of the Newton kernel

The radial derivative of the Newton kernel gives

$$\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}={\int }_{r^{\prime }}{r^{\prime }}^2\frac{\partial L^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r}d{r}^{\prime }.$$E78

By Eq. (18), since

$${\left(u^{\alpha }\right)}^{\prime }=\alpha u^{\prime }{u}^{\alpha -1},$$

we readily get

$$\frac{\partial L^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r}=-\frac{r-r^{\prime }cos\psi }{{\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{3/2}}.$$E79

By substituting the expression (79) in (78), and integrating the equation, we obtain the analytical expression

$$\begin{array}{l}\frac{\partial \widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}=\left(3r^{2}cos\psi +r\left(r^{\prime }-6cos{\psi }^2{r}^{\prime }\right)+cos\psi {r^{\prime }}^2\right)L^{-1}\left(r,\psi ,r^{\prime }\right)\\ +r\left(3cos{\psi }^2-1\right)ln\mid L\left(r,\psi ,r^{\prime }\right)-rcos\psi +r^{\prime }\mid +C_2\left(r,\psi \right).\end{array}$$E80

The computation of $\mathit{\delta A}$ uses the expression of $\partial \widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)/\partial r$.

2.5.3 Second radial derivative of the Newton kernel

Expression of ${\partial }^2\widetilde{L_a^{-1}/\partial r^2}$ emphasizes the second radial derivative of the Newton kernel, which has to be derived two times.

$$\frac{{\partial }^2\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}={\int }_{r^{\prime }}{r^{\prime }}^2\left(\frac{{\partial }^2{L}^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}\right)d{r}^{\prime }$$E81

where

$$\frac{{\partial }^2{L}^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}=\frac{\partial }{\partial r}\left(\frac{\partial L^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r}\right)$$E82

To obtain the second derivative of $L^{-1}\left(r,\psi ,r^{\prime }\right)$, we derive over r, its first derivative, Eq. (79), using the property

$$\left(\frac{u}{v}\right)=\frac{u^{\prime }v-u{v}^{\prime }}{v^2}.$$E83

We finally obtain:

$$\frac{\partial }{\partial r}\left(\frac{\partial L^{-1}\left(r,\psi ,r^{\prime }\right)}{\partial r}\right)=-\frac{1}{{\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{3/2}}+\frac{3{\left(r-r^{\prime }cos\psi \right)}^2}{{\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{5/2}}$$E84

Including Eq. (84) in Eq. (81), the second derivative of $\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)$ gives the followed integral:

$$\frac{{\partial }^2\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}={\int }_{r^{\prime }}-\frac{{r^{\prime }}^2}{{\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{3/2}}+\frac{3{r^{\prime }}^2{\left(r-r^{\prime }cos\psi \right)}^2}{{\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{5/2}}d{r}^{\prime }$$E85

After deriving a second time the reciprocal distance $L^{-1}$ and integrating this equation, we obtain the analytical expression

$$\begin{array}{l}\frac{{\partial }^2\widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}=\left[3r^{3}cos\psi \right.+r^2\left(1-12{cos}^2\psi \right)r^{\prime }\\ \left.+2rcos\psi \left(6{cos}^2\psi +1\right){r^{\prime }}^2+2\left(1-4{cos}^2\psi \right){r^{\prime }}^3\right]L^{-3}\left(r,\psi ,r^{\prime }\right)\\ +\left(3{cos}^2\psi -1\right)ln\mid L\left(r,\psi ,r^{\prime }\right)-rcos\psi +r^{\prime }\mid +C_3\left(r,\psi \right),\end{array}$$E86

where

$$L^{-3}\left(r,\psi ,r^{\prime }\right)={\left({r^{\prime }}^2-2r{r}^{\prime }cos\psi +r^2\right)}^{-3/2}.$$E87

The computation of $\mathit{\delta E}$ uses the expression of ${\partial }^2\widetilde{L_a^{-1}}\left(r,\psi ,r^{\prime }\right)/\partial r^2$.

3. Numerical studies

The determination of topographical effects from DEMs is a very time-consuming process, particularly when computations are required for large areas, such as a country or a continent. With a fine grid resolution, for instance, 50 m Quadratic Grid Resolutions (QGR) a unified spatial data structure, computations are beyond what a multi-processor computer can accomplish within a reasonable time-frame for an area such as Ireland. Numerical investigations have resulted that increasing the spatial resolution of DEM by a factor of two increases CPU computational time by a factor of fourteen.

Although the computational time is a factor to be taken into account, it is less important since it is the spatial resolution of DEM which is critical for improving accuracy required when a precise geoid is determined.

The Fast Fourier Transform (FFT), which relies on linearization and series expansions of the non-linear terrain effect integrals, provides a reduction in computational time by several orders of magnitude, compared to space domain integration methods [13]. In the FFT, the higher-order terms of the radially integrated Newton kernel expressed by the Taylor series expansion are neglected. In addition, the reciprocal distance 1/L is approximated by the planar distance $1/{\mathrm{\ell }}_0$ i.e. ${\mathrm{\ell }}_0=L\left(R,\psi ,R\right)$ [14], Eq. (15)], which is a sufficiently good only if ${\mathrm{\ell }}_0>H\left({\mathrm{\Omega }}^{\prime }\right)$. This condition is violated for computing topographical effects in rough terrain.

3.1 The study area

Topographical effects in high elevated topography in Ireland are computed from 50 m QGR, Ireland with a relatively flat terrain surrounded by the ocean. The study area is $1^{\circ }\times {0.7}^{\circ }$ limited between longitudes $-{5.8}^{\circ }$ and $-{6.8}^{\circ }$ and latitudes ${52.7}^{\circ }$ and ${53.3}^{\circ }.$ A map of topographic heights of the study areas is presented in Figure 8a, and the distribution of topography in meters, counted in the percentage of the number of points in 50 m QGR is demonstrated in histograms Figure 8b.

Figure 8.

Topographic heights over Co-Wicklow Ireland Figure 8a and distribution of elevations in meter, counted in percentage of the number of points from 50 $m^2$ DEM. (a) Topography - County Wicklow Ireland. (b) Distribution of topography.

3.2 The bound of integration area

The gravitational potential of Topographic Masses of finite thicknesses behaves like the potential of a thin layer when it is observed from a larger distance. This is explained by the behavior of integration kernels generating the potential of the gravitational potential $V^t\left(r,\mathrm{\Omega }\right)$ and $V^c\left(r,\mathrm{\Omega }\right).$ Figure 9 illustrates the behavior of the two kernels for the determination of PITE and DTE relative to maximum elevation in Ireland, which is the Carrauntoohil elevation of 1039 m [15]. In the determination of the primary indirect topographical effect, when the angular distance Ψ between the computation point and integration point increases, the integration kernel,

Figure 9.

The behavior of $K_t$ and $K_c$ in the dependence of varying the angular distance $\psi$ in the determination of PITE and DTE in Ireland. The elevation of computation and integration point are fixed to 1039 m and 1 m respectively.

$$K_{\mathit{\text{pite}}}^c\left(R,\psi ,R\right)=R^2\left(\frac{\tau \left({\mathrm{\Omega }}^{\prime }\right)}{L\left(R,\psi ,R\right)}\right)$$E88

generating the potential of Helmert’s condensation layer $V^c\left(r,\mathrm{\Omega }\right)$ approaches the integration kernel

$$K_{\mathit{\text{pite}}}^t\left(R,\psi ,H\left({\mathrm{\Omega }}^{\prime }\right)\right)={\left.\widetilde{L^{-1}}\left(R,\psi ,r^{\prime }\right)\right|}_{r^{\prime }=R}^{R+H\left({\mathrm{\Omega }}^{\prime }\right)}$$E89

generating the gravitational potential $V^t\left(R,\mathrm{\Omega }\right).$ Thus, the differences between two kernels $K_{\mathit{\text{pite}}}^t-K_{\mathit{\text{pite}}}^c=\delta K_{\mathit{\text{pite}}}$ decrease (see Figure 9-PITE kernel). In determination of DTE Eq. (57) a similar decrease occurs for the difference between two kernels $\left(K_{\mathit{dte}}^t-K_{\mathit{dte}}^c=\delta K_{\mathit{dte}}\right),$ when the integration kernel

$$K_{\mathit{dte}}^c\left(H\left(\mathrm{\Omega }\right),\psi ,R\right)=R^2\tau \left(\mathrm{\Omega }\right)\frac{\partial L^{-1}\left(r,\psi ,R\right)}{\partial r},$$E90

generating the potential of Helmert’s condensation layer $V^c\left(r,\mathrm{\Omega }\right)$ approaches the integration kernel

$$K_{\mathit{dte}}^t(H\left(\mathrm{\Omega }\right),\psi ,H\left({\mathrm{\Omega }}^{\prime }\right))=\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}{|}_{r^{\prime }=R}^{R+H\left({\mathrm{\Omega }}^{\prime }\right)},$$E91

generating the gravitational potential $V^t\left(r,\mathrm{\Omega }\right)$ of the TMs (see Figure 9-DTE). This also means that the magnitudes of $\delta K_{\mathit{\text{pite}}}$ or $\delta K_{\mathit{dte}}$ are largest in the immediate neighborhood of the computation point.

The choice of varying angular distance and fixing the elevation of integration point to 1039 m and computation point to 1 m or vice versa, enables us to determine the most attainable differences (maximum or minimum) between the kernels in question.

Furthermore, a numerical examination of $\delta K_{\mathit{\text{pite}}}$ and $\delta K_{\mathit{dte}},$ controlled by varying topographic height $H\left({\mathrm{\Omega }}^{\prime }\right)$ at a fixed angular distance $\psi ={0.0001}^{\circ },$ in the immediate neighborhood of the computation point shows that the larger the height of the integration point, the larger the difference between these kernels, and therefore the topographical effects are stronger (see Figure 10).

Figure 10.

Magnitudes of $\mid \delta K_{\mathit{\text{pite}}}\mid$ and $\mid \delta K_{\mathit{dte}}\mid$ in an immediate neighborhood $\left(\psi ={0.0001}^{\circ }\right)$ of the computation point with $H\left(\mathrm{\Omega }\right)=1$ and varying height of integration point.

Eqs. (36) and (57) assume the compensation is strictly local, which means $\mathit{\delta A}\left(\mathrm{\Omega }\right)$ consists of the terrain roughness contribution only. The longitudinal profiles (see Figure 11) show that the Bouguer components of PITE have a larger contribution to topographical effects than the terrain roughness components. Since DTE does not contain a Bouguer component, the correlation between DTE and DEM is in general smaller than that for PITE.

Figure 11.

Correlation between topographical effects and elevation of topography from $50m^2$ QGRs.

3.3 Direct topographical effect on gravity

3.3.1 On the Earth’s surface - residual effect $\delta A_{\mathit{\text{surf}}}$

The computation of the residual direct topographical effect on gravity on the Earth’s surface $\delta A_{\mathit{\text{surf}}}$ uses Eq. (57). Table 1 shows the minimum, mean, maximum and root mean square (rms) of $\delta A_{\mathit{\text{surf}}}$ in mGal.⁵

	Min	Mean	Max	rms
Topographic Heights (m)	0	58.600	920.000	$\pm 200.000$
$\delta A_{\mathit{\text{surf}}}$ (mGal)	−17.4	0.006	5.7	$\pm 0.7$
$\delta A_{\mathit{\text{surf}}}^t$ (mGal)	−27.5	−0.6	2.6	$\pm 1.5$

Table 1.

Direct topographical effect $\mathit{\delta A}$ on gravity.

Due to the low elevation of topography in the test area (maximum 920 m), $\delta A_{\mathit{\text{surf}}}$ is small. Numerical investigation with different grid resolutions shows that the flatter and lower the terrain, the smaller the reduction of the surface gravity observations in the remove step of the RCR procedure is. And a sparse griding reduces values of $\mathit{\delta A},$ so using as finer griding as possible is recommended.

Since the residual effect $\mathit{\delta A}$ is rather small, we compared it with the total effect to estimate the efficiency of reduction by the Helmert condensation using the following expression

$$A^t\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}{\left.\frac{\partial \widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r}\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)},$$E92

where $r=r_t\left(\mathrm{\Omega }\right)$ (see Table 1). The total direct topographical effect on gravity $A_{\mathit{\text{surf}}}^t$ which is the case when the Helmert condensation is not employed, is larger in amplitude than the residual effect $\delta A_{\mathit{\text{surf}}},$ and mostly negative. The values of $A_{\mathit{\text{surf}}}^t$ are prevailingly distributed around 0 mGal in flat areas and are highly correlated with the topography.

3.3.2 At satellite altitudes—residual effect $\delta A_{\mathit{sat}}$

The computation of the residual direct topographical effect on satellite gravity $\delta A_{\mathit{sat}}$ uses Eq. (59) for $r_{\mathit{sat}}=r_{\mathit{\text{GRACE}}}=R+400\mathrm{km}.$.

The residual effect on gravity $\delta A_{\mathit{sat}}$ computed at satellite altitude is of the order of ${10}^{-5}$ mGal, that is significantly smaller in comparison with $\delta A_{\mathit{\text{surf}}}$ which is of the order of 10’s to 100’s of mGals (see Table 2). The total direct topographical effect on satellite gravity $A_{\mathit{sat}}^t$ also uses Eq. (92) for $r=r_{\mathit{\text{GRACE}}}.$ Comparing the residual effect $\delta A_{\mathit{sat}}$ with the total effect $A_{\mathit{sat}}^t,$ we see that the total effect $A_{\mathit{sat}}^t$ has a long-wavelength feature. This shows that the gravitational signal of the topography is attenuated when going to satellite altitudes in such a way that short wavelengths of the gravitation are attenuated faster than long wavelengths. The fact that the total effect $A_{\mathit{sat}}^t$ is of the order of 10’s of mGals and the residual effect $\delta A_{\mathit{sat}}$ almost vanishes shows that the compensation of the masses by the Helmert condensation is an efficient way to process satellite gravity data. However, this is not the case for surface gravity data. Note that a coarser griding does not affect the results of $A_{\mathit{sat}}^t.$

	Min	Mean	Max	rms
$\delta A_{\mathit{sat}}$ (${10}^{-4}$ mGal)	−3.1	−0.002	138.9	$\pm 6.7$
$\delta A_{\mathit{sat}}^t$ (mGal)	−0.87	−0.74	−0.51	$\pm 0.8$

Table 2.

Direct topographical effect $\mathit{\delta A}$ on satellite gravity.

3.4 Direct topographical effect on gradiometry—residual effect $\mathit{\delta E}$

The computation for the direct topographical effect on gradiometry $\mathit{\delta E}$ uses Eq. (71) for $r_{\mathit{sat}}=r_{\mathit{\text{GOCE}}}=R+250\mathrm{km}.$ Table 3 shows the minimum, mean, maximum and rms values of $\mathit{\delta E}$ in Eötvös.⁶ The values of $\mathit{\delta E}$ given in Table 3 are of the order of 0.1 $\mathrm{m}\stackrel{°}{\mathrm{E}}.$ The coarser topographical griding does not affect the results of the topographical effect $\mathit{\delta E}$ on gradiometry. Even though $\mathit{\delta E}$ has tiny values, these are still slightly correlated with the topography. As the measuring accuracy of the second derivatives by the GOCE gradiometer is 10 m$\stackrel{°}{\mathrm{E}},$ the effect of $\mathit{\delta E}$ should still be taken into account in the RCR procedure, especially in mountainous area.

	Min	Mean	Max	rms
$\mathit{\delta E}\left(m\stackrel{°}{E}\right)$	−2.13	−0.00001	0.7	$\pm 0.10$
$V_{\mathit{rr}}^t\left(\stackrel{°}{E}\right)$	0.01	0.7	0.11	$\pm 0.02$

Table 3.

Direct topographical effect $\mathit{\delta E}$ on gradiometry.

The total direct topographical effect $V_{\mathit{rr}}^t$ on gradiometry is computed by

$$V_{\mathit{rr}}^t\left(r,\mathrm{\Omega }\right)=G{\varrho }_0{\int }_{{\mathrm{\Omega }}_0}{\left.\frac{{\partial }^2\widetilde{L^{-1}}\left(r,\psi ,r^{\prime }\right)}{\partial r^2}\right|}_{r^{\prime }=R}^{r_t\left({\mathrm{\Omega }}^{\prime }\right)},$$E93

where $r=r_{\mathit{\text{GOCE}}}.$ Table 3 shows that the values of $V_{\mathit{rr}}^t$ are two orders in magnitude larger than $\mathit{\delta E}.$ The fact that the total effect $V_{\mathit{rr}}^t$ is much larger in amplitude than the residual effect $\mathit{\delta E}$ shows that the compensation of the masses by the Helmert condensation is again an efficient way to process satellite gravity gradient data as well as the satellite gravity data. The distribution of the total effect $V_{\mathit{rr}}^t$ between the maximum and minimum values is much more homogeneous compared to the distribution of the residual effect $\mathit{\delta E}$ with most values around 0$\stackrel{°}{\mathrm{E}}.$

3.5 Topographical effect on geoid heights

3.5.1 The effect of PITE on geoid heights

The primary indirect topographical effects on the geoid height determination are computed by Eq. (42),

$$\mathit{\delta N}\left(R,\mathrm{\Omega }\right)=\frac{\mathit{\delta V}\left(R,\mathrm{\Omega }\right)}{{\gamma }_Q}.$$

Figure 12b illustrates the effect of PITE in a high elevated topography at 50 m grid resolution.

Figure 12.

Illustration of DTE and PITE on geoid undulation for High elevated topography in Ireland at 50 meter grid resolutions. (a) Effects of DTE on geoid undulation. (b) Effects of PITE on geoid undulation.

3.5.2 The effect of DTE on geoid heights

The direct topographical effects on geoid heights, $\mathit{\delta D},$ are computed by applying Stokes’ integral to DTE as a known function $f\left(\mathrm{\Omega }\right)$ distributed in the near-zone spherical cap $C_{{\psi }_0}$

$$\delta D^{\mathrm{\ell },{\psi }_0}\left(\mathrm{\Omega }\right)=\frac{R}{4\pi }{\int }_{C_{{\psi }_0}}{f}^{\mathrm{\ell }}\left({\mathrm{\Omega }}^{\prime }\right)S^{\mathrm{\ell }}\left(\psi \right)d{\mathrm{\Omega }}^{\prime },$$E94

where $\delta D^{\mathrm{\ell },{\psi }_0}$ is a higher-degree (presented as superscript $\mathrm{\ell }$) correction to the geoid height, $C_{{\psi }_0}$ is a spherical cap of radius ${\psi }_0$ (for this study ${\psi }_0=1°$ is chosen) and $S^{\mathrm{\ell }}\left(\psi \right)$ is the spheroidal Stokes function, [16]. The effect of DTE on geoid heights is also presented in Figure 12a.

4. Conclusions

The RCR procedure for gravity and gravimetry is summarized in Figure 2. It shows the motivation for correcting the data by the effects of the topographical masses. These are the direct or indirect topographical effects computed at surface or satellite altitude. This study summarized the field corrections for the determinations of geoids by terrestrial data or the latest satellite missions. The integral of newton presented in Sec:2.1.1 in the form of the gravitational potential of topography Eq. (16), proves necessary to the expression of topographical effects. It made it possible to understand the observations made either $\partial V/\partial r$ or second derivatives ${\partial }^{2}V/\partial r^2.$ Furthermore, make it possible to fulfil the conditions necessary to determine a model of geoids whose centimetre precision is required.

Presentation of the RCR method allowed us to understand the interest in determining the effects of topographic masses in order to correct gravimetric and gradiometric measurements and to be able to apply the calculations to determine a geoid model.

In order to express the effects, we need to know the topographic masses defined by their height and density. We have seen that Newton’s integral for determining the geoid allows expressing the effects. The topography effects on the satellite measurements are calculated in a way similar to the effects on ground measurements.

The effects of topographic masses have much less impact on gravity measurements at the satellite level than on ground measurements.

Computing topographical effects for large areas is a very time-consuming process. Increasing the resolution of sampled DEM by a factor of 2 (e.g. from 100 m to 50 m quadrangle) increases the number of data by a factor of 4, and it increases the computational time by a factor of approximately 14. Thus, it is suggested to restrict the integration area to a small area of radius ${\psi }_0$ around the computation point.

A sparse grid size, particularly in rugged areas, is not sufficient to express the irregularities of the terrain and thus does not reveal properly the contribution to geoidal height due to terrain height variations.

With a tiny grid step size, the magnitude of the Bouguer component becomes comparable with that of the terrain roughness component, which reduces the correlation between DTE and PITE. Since DTE does not contain a Bouguer component, the correlation between DTE and DEM is generally smaller than that for PITE.

Numerical investigation shows that the Bouguer components of PITE have a larger contribution to topographical effects than the terrain roughness components.

Numerical examination of Kernel’s, controlled by varying topographic height at a fixed angular distance, in the immediate neighborhood of the computation point, shows that the larger the height of the integration point, the more significant the difference between these kernels, and therefore, the topographical effects are more substantial.

Comparing the results with different grid sizes shows (not shown here) an improvement in computation accuracy. Contrary to our expectations, it is not the case for calculations at satellite altitudes, so griding can be reduced, and a more refined grid does not change a long-wavelength feature of $V_{\mathit{rr}}^t.$

Acknowledgments

The authors would like to thank Dr. Pat Gill (TUS) and Dr. Nic Wilson (UCC) for their comments that have helped to improve the manuscript.