Topographic Effects in Geoid Determinations

Abstract

Traditionally, geoid determination is applied by Stokes’ formula with gravity anomalies after removal of the attraction of the topography by a simple or refined Bouguer correction, and restoration of topography by the primary indirect topographic effect (PITE) after integration. This technique leads to an error of the order of the quasigeoid-to-geoid separation, which is mainly due to an incomplete downward continuation of gravity from the surface to the geoid. Alternatively, one may start from the modern surface gravity anomaly and apply the direct topographic effect on the anomaly, yielding the no-topography gravity anomaly. After downward continuation of this anomaly to sea-level and Stokes integration, a theoretically correct geoid height is obtained after the restoration of the topography by the PITE. The difference between the Bouguer and no-topography gravity anomalies (on the geoid or in space) is the “secondary indirect topographic effect”, which is a necessary correction in removing all topographic signals. In modern applications of an Earth gravitational model (EGM) in geoid determination a topographic correction is also needed in continental regions. Without the correction the error can range to a few metres in the highest mountains. The remove-compute-restore and Royal Institute of Technology (KTH) techniques for geoid determinations usually employ a combination of Stokes’ formula and an EGM. Both techniques require direct and indirect topographic corrections, but in the latter method these corrections are merged as a combined topographic effect on the geoid height. Finally, we consider that any uncertainty in the topographic density distribution leads to the same error in gravimetric and geometric geoid estimates, deteriorating GNSS-levelling as a tool for validating the topographic mass distribution correction in a gravimetric geoid model.

1. Introduction

Stokes’ formula [1] is fundamental for gravimetric geoid determination. As it requires no masses outside the sphere of computation, traditionally ([2], Chapter 3) the topographic signal on gravity is removed or reduced by a compensation mass below or by a density layer at sea-level (direct topographic effect; DITE on gravity). Another topographic correction is the free-air correction, which provides a downward continuation (dwc) of gravity from the Earth’s surface to sea level. By subtracting normal gravity at the reference ellipsoid, a Bouguer type of gravity anomaly is obtained. Adding the so-called “secondary indirect topographic effect” (SITE), Stokes’ operator yields the regularized geoid or co-geoid, and finally, after adding also the primary indirect topographic effect (PITE) on the geoid for restoration of the topographic signal, the geoid height follows.

More or less the same topographic corrections are applied in one way or another also in the modern remove-compute-restore (RCR) technique (e.g., [3,4,5,6,7]). Some methods start from the classical gravity anomaly on the geoid (as above), while others (like [5]) start from Molodensky’s surface gravity anomaly [8]. The latter anomaly is also used in the KTH method, also called least squares modification of Stokes’ formula with additive corrections (LSMSA; [9,10]), where all topographic corrections (the DITE, SITE, PITE, as well as the dwc correction) are added as total effects to the modified Stokes’ formula.

Below, we discuss and compare some of the above principles for topographic corrections in geoid determination. To make the presentation more transparent, the methods are given without corrections for atmospheric masses and ellipsoidal effects (e.g., [11,12]).

2. The Traditional Solution to Stokes’ Formula

Stokes’ formula determines the geoid height N with respect to a selected reference ellipsoid from gravity anomalies $\Delta g$ by the global integral:

$$N=\frac{R}{4\pi {\gamma }_0}{\displaystyle \underset{\sigma }{\iint }S\left(\psi \right)}\Delta gd\sigma$$

(1)

where R is the radius of the sphere of computation (which approximates sea level), $\sigma$ is the unit sphere, $S\left(\psi \right)$ is Stokes’ function with argument $\psi$ being the geocentric angle between integration and computation points, and ${\gamma }_0$ is normal gravity at the reference ellipsoid.

As gravity is observed on the Earth’s surface, the application of Equation (1) requires that the gravitational effect of the topography must be removed from observed gravity. This is performed by applying the so-called direct topographic effect (DITE), which in the classical approach is solved by a simple or refined Bouguer gravity correction, the latter being:

$$-B=\frac{\partial V^T}{\partial r}$$

(2)

where $V^T$ is the potential of the topography and r is the radius of the observation point.

The gravity value must be downward continued to sea level (the dwc correction), which is accomplished by applying the free-air correction:

$$F=-H(\partial \gamma /\partial h)\approx 2{\gamma }_{0}H/R$$

(3)

where H is the orthometric height, to the observed anomaly (e.g., ([2], Chapter 3) and ([7], Chapter 8)).

An additional small correction, the secondary indirect topographic effect (SITE), (denoted $d{g}_I$) is also needed for the gravity anomaly to comply with the potential at the so-called co-geoid (the geoid without topography). It can be explicitly expressed:

$$d{g}_I=2\frac{{\gamma }_0}{R}d{N}_I^T$$

(4)

where $d{N}_I^T$ is given in Equation (5).

Finally, after Stokes’ integration, the effect of restoration of the topography is executed by adding the primary indirect topographic effect (PITE):

$$d{N}_I^T=\frac{V_g^T}{{\gamma }_0},$$

(5)

where $V_g^T$ is the potential of the topography on the geoid (denoted by subscript g).

The whole computational procedure can therefore be expressed as:

$$N_1=\frac{R}{4\pi {\gamma }_0}{\displaystyle \underset{\sigma }{\iint }S\left(\psi \right)}\Delta g_{red}d\sigma +d{N}_I^T$$

(6a)

where:

$$\Delta g_{red}=g-B+F-{\gamma }_0+d{g}_I=\Delta g^B+d{g}_I$$

(6b)

At this point it should be stated that the above corrections (DITE, SITE, and PITE) can be reduced by compensating/removing from each of them the effect of a suitable topographic model located on (e.g., a Helmert coating layer) or below sea level. Although this is a fruitful numerical trick to obtain smaller numbers in the numerical application of each correction, it does not change anything in theory, as the result is theoretically independent of the choice of the compensation model.

It was shown in [13] that Equation (6a) suffers from a geoid error dN of the order of the difference between the height anomaly $\left(\zeta \right)$ and the geoid height:

$$dN\approx \zeta -N,$$

(7)

which is negligible for low elevation topography, but ranges to a few metres in the highest mountains (e.g., [14]).

The method in Equation (6) is frequently included as part of the so-called remove-compute-restore (RCR) technique, but it usually also includes an Earth gravitational model (EGM) that cares that the long-wavelength part of the gravity anomaly is removed prior to Stokes integration and added back as a contribution to the final geoid estimator. Some general remarks on these types of applications of the RCR method are provided in [15].

3. The RCR Technique with Surface Gravity

By utilizing the surface gravity anomaly in the RCR technique, the above problem resulting in the error expressed by Equation (7), can be avoided. This approach will be presented here.

First, considering the boundary condition of physical geodesy (([2], Sections 2–14), ([16], Section 8.3)):

$$\Delta g=-\frac{\partial T}{\partial r}+\frac{\partial \gamma }{\partial h}\frac{T}{\gamma }\approx -\frac{\partial T}{\partial r}-2\frac{T}{r},$$

(8)

where the approximation is the result of the spherical approximation of the derivative of $\gamma$ with respect to the geodetic height (h). If the disturbing potential T is decomposed into the no-topography (NT) and topographic potentials $T^{NT}$ and $V^T$, respectively:

$$T=T^{NT}+V^T,$$

(9)

and inserted into Equation (8), it follows (see also [17,18]) that the gravity anomaly can be decomposed into a NT gravity anomaly ($\Delta g^{NT}$) and a topographic gravity anomaly ($\Delta g^T$):

$$\Delta g=\Delta g^{NT}+\Delta g^T,$$

(10a)

where:

$$\Delta g^{NT}=-\frac{\partial T^{NT}}{\partial r}+\frac{\partial \gamma }{\partial h}\frac{T^{NT}}{\gamma }\approx -\frac{\partial T^{NT}}{\partial r}-2\frac{T^{NT}}{r},$$

(10b)

and:

$$\Delta g^T=-\frac{\partial V^T}{\partial r}+\frac{\partial \gamma }{\partial h}\frac{V^T}{\gamma }\approx -\frac{\partial V^T}{\partial r}-2\frac{V^T}{r}.$$

(10c)

Equation (10a–c) show that only $\Delta g^{NT}$, and not $\Delta g^B$, includes the total topographic effect $\left(-\Delta g^T\right)$, which differs from the Bouguer gravity effect –B of Equation (2). Therefore, the NT anomaly, but not the Bouguer anomaly, can be correctly analytically downward continued to the sphere, from the surface to the geoid as a harmonic function, denoted by ${\left(\Delta g^{NT}\right)}^{*}$, yielding Stokes’ integral:

$$N_2=\frac{R}{4\pi \gamma }{\displaystyle \underset{\sigma }{\iint }S\left(\psi \right)}{\left(\Delta g^{NT}\right)}^{*}d\sigma +d{N}_I^T,$$

(11)

The first term on the right hand-side is the NT geoid height:

$$\frac{R}{4\pi {\gamma }_0}{\displaystyle \underset{\sigma }{\iint }S\left(\psi \right)}{\left(\Delta g^{NT}\right)}^{*}d\sigma =\frac{T_g^{NT}}{{\gamma }_0}=N^{NT},$$

(12)

So that Equation (11) equates to:

$$N_2=N^{NT}+d{N}_I^T=\frac{T_g^{NT}+V_g^T}{{\gamma }_0}=\frac{T_g}{{\gamma }_0}=N,$$

(13)

i.e., N₂ is the correct geoid height.

The main problem in this version of the RCR technique is the procedure to downward continue the surface gravity anomaly to the sphere. As one example, this method (with topographic compensation by Helmert condensation) is used in the University of New Brunswick (UNB) geoid model [5,19], where the dwc of the gravity anomaly is carried out by solving the discrete Poisson’s integral equation for $\Delta g^{NT}$, which is likely to become numerically singular for high-resolution gravity anomalies in high mountains (in his pioneer approach to determine the quasigeoid by reducing gravity anomalies to an internal sphere, A. Bjerhammar [20,21] also employed Poisson’s integral equation, but his method implies an analytical continuation into the topographic masses, formally a strictly singular problem).

4. The KTH Approach

Equation (10a–c) are also the basis in the KTH geoid modelling technique ([9,10] and ([16], Section 6.2.2 and Chapter 5)). However, this approach differs in a number of ways from the RCR technique. Generally, the solution starts from a preliminary geoid estimator obtained by a direct integration of the surface anomaly in Stokes’ formula, a long-wavelength contribution is added from an EGM and the two contributors are spectrally weighted by a least squares procedure. The topographic corrections DITE and PITE are added as a combined topographic effect for the geoid height, and the dwc effect is directly computed as another correction to the geoid height. In this way the final formulas for the corrections become simpler, and the dwc effect becomes more stable ([16], Chapter 5). In the case of Stokes’ integral solution, the combined topographic correction of the geoid height $\left(N^T\right)$ is the Bouguer shell correction (the negative of the topographic bias of [22,23,24,25]; see also Section 4.1 below):

$$N^T=-2\pi G\rho \left[H^2+\frac{2H^3}{3R}\right],$$

(14)

where G is the gravitational constant and $\rho$ is the topographic density. In [22,25] Equation (14) is also generalized for a variable topographic density with depth. That is, Stokes’ formula may be expressed (cf. Equation (11) for the RCR technique):

$$N_{KTH}=\frac{R}{4\pi \gamma }{\displaystyle \underset{\sigma }{\iint }S\left(\psi \right)}{\left(\Delta g\right)}^{*}d\sigma +N^T,$$

(15)

where the downward continued gravity anomaly ${\left(\Delta g\right)}^{*}$ can be determined, e.g., by Taylor expansion. This is why there is no need for adding a terrain correction in the KTH approach.

As the KTH method implies a least squares solution, an internal mean square error estimate of the geoid height is also provided. For applications, see ([16], Section 6.6.1 and its reference list).

Finally we mention that Equation (14) can also be used to estimate the error $\left({\epsilon }_N\right)$ in the geoid height for an estimated error $d\mu =Gd\rho$ in the topographic density (([16], Section 5.2.8), [26]),

$${\epsilon }_N\left(d\mu \right)=2\pi d\mu \left[H^2+\frac{2H^3}{3R}\right].$$

(16)

4.1. More on the Topographic Bias

The topographic bias in geoid determination implies that the external topographic potential $V^T$ is analytically continued as a harmonic function down to seal level, denoted ${\left(V^T\right)}^{*}$. As $V^T$ is not harmonic inside the topographic masses, the geoid height will be biased by:

$$bias=\frac{{\left(V^T\right)}^{*}-V_g^T}{{\gamma }_0}=-N^T,$$

(17)

where $N^T$ is the correction presented in Equation (14) in the case of a constant topographic density. This bias was first demonstrated in [22] and was later proved in different ways in [22,23,24,25] (a refinement for a vertical variation of density was also shown in [22]). As the resulting bias is nothing but the correction for analytical continuation through a spherical Bouguer plate, there are still many colleagues who believe that Equation (17) is an approximation, and that a terrain correction is missing. Here we will show that this is not the case.

Outside topographic masses the topographic potential is harmonic and, therefore, obeying the Laplace differential equation, $\Delta V^T=0$, where $\Delta$ is the Laplace operator, and there is no topographic bias. On the contrary, inside the masses Poisson’s differential equation holds:

$$\Delta V^T=-4\pi \mu ,$$

(18)

suggesting that there is a topographic bias.

The topographic bias in Equation (17) can be decomposed into the biases of the Bouguer shell and the terrain (with obvious notations):

$$bias=bia{s}^{shell}+bia{s}^{terr}.$$

(19)

Applying the Laplace operator on all terms one obtains:

$$\Delta bias=\Delta bia{s}^{shell}+\Delta bia{s}^{terr},$$

(20a)

or, as both the total bias and that of the shell obey Equation (18), it follows that:

$$\frac{-4\pi \mu }{{\gamma }_0}=\frac{-4\pi \mu }{{\gamma }_0}+\Delta bia{s}^{terr},$$

(20b)

yielding:

$$\Delta bia{s}^{terr}=\frac{\Delta {\left(V^{terr}\right)}^{*}-\Delta V_g^{terr}}{{\gamma }_0}=-\frac{\Delta V_g^{terr}}{{\gamma }_0}=0.$$

(21)

Equation (21) shows that the terrain potential is harmonic at the geoid, which by definition implies that there is no topographic bias. Hence, the topographic correction in geoid determination can be determined without a terrain correction [25], which is the case in the KTH method. This is because the terrain effect is already considered in the analytical continuation of the potential, which is utilized in the KTH method.

5. Geoid Heights by an EGM

Frequently the geoid height is directly computed by an EGM:

$$N=\frac{R}{{\gamma }_0}{{\displaystyle \sum _{n=0}^{n_0}\left(\frac{a}{r_e}\right)}}^{n+1}{\displaystyle \sum _{m=-n}^n{C}_{nm}^{*}{Y}_{nm}}\left(\theta ,\lambda \right),$$

(22)

where a is the radius of definition of the EGM, $r_e$ is the radius of the reference ellipsoid, $\left(\theta ,\lambda \right)$ = (colatitude, longitude), $Y_{nm}$ is the spherical harmonic of degree n and order m, and $C_{nm}^{*}$ is the corresponding EGM potential coefficient minus the coefficient of the normal gravity field. The resolution of the estimate is governed by the degree of truncation $\left(n_0\right)$ of the series.

However, in continental areas Equation (22) must also be corrected for topography, as the series is analytically continued inside the topography down to sea-level. This correction can be approximated by:

$$dN=-2\pi G\rho \left[H_{n_0}^2+\frac{2H_{n_0}^3}{3R}\right]$$

(23a)

where:

$$H_{n_0}^i={\displaystyle \sum _{n=0}^{n_0}{\displaystyle \sum _{m=-n}^n{H}_{nm}^i}}{Y}_{nm};i=2,3.$$

(23b)

Here $H_{nm}^i$ are the spherical harmonics of the orthometric height.

In case Equation (22) converges down to sea level, Equation (23a) holds strictly. This was demonstrated for a homogeneous ellipsoid in [27]. In the general case one can assume that the series does not converge down to sea level, and a more precise formula for the topographic bias correction was given in [28] and ([16], pp. 160–163).

6. RCR and KTH Methods Combined with an EGM

Both RCR methods [Equations (6a) and (12)], as well as the KTH method (Equation (15), can be combined with an EGM ([16], Sections 6.2.2.1 and 6.2.2.2). In the case of the RCR approach this combination is performed by removing the low-degree gravity anomaly from the observed gravity anomaly prior to Stokes’ integration and adding back the EGM in the form of a low-degree geoid height after integration. This approach is sometimes called “geoid determination by a higher order reference field” (e.g., [29]), but this terminology may be misleading, as for both approaches topographic corrections are also needed for the EGM parts (which is frequently overlooked).

Theoretically, the RCR and KTH methods should give the same result, but this would only happen if all corrections are based on the same direct and indirect effects (corresponding to the combined/additive effect in the KTH method). However, as the KTH method includes a least squares spectral combination of the EGM and Stokes’ contributions, it is (at least theoretically) superior. In addition, the latter method yields also an internal error estimate of the computed geoid solution based on the least squares residuals.

7. GNSS-Levelling

Frequently the geoid height is determined as the difference between geodetic heights determined by satellite positioning (usually GNSS) and orthometric heights:

$$N=h-H.$$

(24)

Similar to the gravimetric geoid height one error source to this (geometric) geoid height solution is the uncertainty in the topographic density distribution, say $d\mu$, which affects the orthometric height by the error $dH(d\mu )$ (but not the geodetic height given by the space positioning). Hence, one obtains from Equation (24) the geoid error:

$$dN(d\mu )=-dH(d\mu ).$$

(25)

In [30] it was shown that:

$$dH(d\mu )=-{\epsilon }_N\left(d\mu \right),$$

(26)

where ${\epsilon }_N\left(d\mu \right)$ is the error of $N$ in Stokes’ formula for an incorrect topographic density as given in Equation (16). Hence, in view of Equations (16), (25), and (26) one can see that the errors due to the uncertainty in topographic density are the same in geoid determinations by gravimetric (e.g., by Stokes’ formula) and geometric (e.g., by GNSS/levelling) methods (see also [30] for an application).

8. Concluding Remarks

First we conclude that the traditional Bouguer gravity anomaly, used in Equation (6a,b), is not sufficiently accurate for today’s demands, say, for the 1-cm geoid geoid determination, as its error is most significant in mountainous areas. Only the no-topography (NT) gravity anomaly, which also corrects for the SITE at the point level, fulfils this requirement.

The RCR and KTH methods based on the NT gravity anomaly with basically the same corrections for both methods should provide the same geoid estimates, at least theoretically. However, the KTH method includes a least squares spectral combination of the Stokes and EGM parts of the combined estimator, which is normally not the case in the RCR application. Other comparisons of the KTH and RCR methods can be found in ([16], Section 6.4).

GNSS-levelling is a typical tool to verify a gravimetric geoid model. However, as discussed in the previous section, a suitable choice of the density distribution for the topographic corrections cannot be validated in this way.