The Gravity Effect of Topography: A Comparison among Three Different Methods

Abstract

In this paper, three different methods for computing the terrain correction have been compared. The terrain effect has been accounted for by using the standard right parallelepiped closed formula, the spherical tesseroid and the flat tesseroid formulas. Particularly, the flat tesseroid approximation is obtained by flattening the top and the bottom sides of the spherical tesseroid. Its gravitational effect can be computed as the gravitational effect of a polyhedron, i.e. a three-dimensional body with flat polygonal faces, straight edges and sharp corners or vertices. These three methods have been applied in the context of a Bouguer reduction scheme. Two tests were devised in the Alpine area in order to quantify possible discrepancies. In the first test, the terrain correction has been evaluated on a grid of points on the DTM. In the second test, Bouguer gravity anomalies were computed on sparse observed gravity data points. The results prove that the three methods are practically equivalent even in an area of rough topography though, in the second test, the Bouguer anomalies obtained by using the tesseroid and the flat tesseroid formulas have slightly smaller RMSs than the one obtained by applying the standard right parallelepiped formula.

1. Introduction

The gravity effect of topography has been always intensively analyzed and modeled. In all the classical books of Geodesy and Geophysics there are sections devoted to this topic. Among the standard methods for modeling the gravity effect of topography one can enumerate the Bouguer reduction, the Helmert reduction and the isostasy reduction according to Airy-Heiskanen and Pratt-Hayford models (see e.g. [1]). In more recent years, the Residual Terrain Correction (RTC) has been devised as a method to be used in geodetic applications for gravity field and geoid estimation (see e.g. [2]).

In all these approaches, the gravity effect of topography is computed by using a Digital Terrain Model (DTM) at a given resolution, which is assumed to represent the actual shape of the Earth surface. Thus, the topography is discretized at the DTM resolution and the gravitational effect of each topography element is computed. In doing so, different formulas can be applied. Usually, the single terrain element is modeled as a right parallelepiped (see e.g. [2]) or as a spherical or ellipsoidal tesseroid [3]. In this study, we compare the two aforementioned approach with that introduced by Tsoulis [4] in which the bottom and the top sides of the tesseroid are flat surfaces (flat tesseroid). This is done for the computation of the terrain correction in the framework of the Bouguer reduction.

The formulas giving the gravitational effect of a right parallelepiped, a spherical tesseroid and a flat tesseroid are given in Section 2. In Section 3, the Bouguer reduction is computed using these three different approaches in an area of the Alps, both on a grid of points and on a set of observed gravity values, and comments on the obtained results are given. Conclusions are then stated in Section 4.

2. The terrain correction and its computation

The terrain correction is commonly applied in the computation of the Bouguer gravity anomalies. The gravity observation in a point P on the Earth surface is strongly influenced by gravitational effect of the topographic masses. In order to use gravity in geophysical analyses, e.g. for estimating the Moho depth or intra-crustal mass anomalies, the topographic gravity signal is removed from the observed data. In this context, the most frequently used reduction is the Bouguer reduction (see e.g. [1]).

From the observed gravity value g(P) one removes the gravity effect of a plate of height H_P equal to the height of point P, the so-called Bouguer plate. This plate is usually considered as an infinite horizontal slab of known constant density ρ even though spherical plate models have been proposed [5]. If the infinite horizontal slab model is assumed, the gravitational effect A_B is given by the well-known formula (see e.g. [1]):

$$A_B=2\mathit{\pi G\rho }{H}_P$$E1

Where G = 6.67 × 10⁻¹¹ m³/kg s² is the universal gravitational constant. The topography signal reduction is then refined through the so-called Terrain Correction (TC). This is computed accounting for the masses that are above or below the plate of height H_P (see Figure 1). This is a quantity always positive that must be summed to the observed gravity reduced by the effect of the plate. Thus, one has:

Figure 1.

The plate and the terrain correction (H_P = orthometric height, h_P = ellipsoidal height).

$$g\left(P\right)-A_B+\mathit{TC}\left(P\right)$$E2

It must be underlined that both the plate and the TC effects are usually computed at constant density ρ set at 2670 kg/m³. This means that residual effects coming from masses having density different from the standard value given above may still affect the reduced gravity values. However, it is assumed that, in this way, the topography effect is substantially removed from the observed gravity values g(P).

The reduced gravity data are then moved to the geoid by considering the gravity gradient, which is approximated by the free-air normal gravity gradient ∂γ/∂h:

$$g\left(P\right)-A_B+\mathit{TC}\left(P\right)-\frac{\partial g}{\partial H}{H}_P\cong g\left(P\right)-A_B+\mathit{TC}\left(P\right)-\frac{\partial \gamma }{\partial h}{H}_P$$E3

Finally, the normal gravity at point Q (see Figure 1) on the ellipsoid is subtracted and the standard Bouguer anomaly is obtained as:

$$\mathrm{\Delta }{g}_B=g\left(P\right)-A_B+\mathit{TC}\left(P\right)-\frac{\partial \gamma }{\partial h}{H}_P-\gamma \left(Q\right)$$E4

While the Bouguer plate effect can be easily computed via e.g. Eq. (1), the computation of the terrain effect is more complex. In planar approximation, this effect is given by the integral formula:

$$\mathit{TC}\left(P\right)=\mathit{G\rho }{\iint }_{-\infty }^{+\infty }\underset{z_P}{\overset{z}{\int }}\frac{\zeta -z_P}{l^3}\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta$$E5

where $=\sqrt{{\left(x_P-\xi \right)}^2+{\left(y_P-\eta \right)}^2+{\left(z_P-\zeta \right)}^2}$, $\left(x_P,y_P,z_P\right)$ are the coordinates of the computational point $P$ and $\left(\xi ,\eta ,\zeta \right)$ are the coordinates of the integration point. The integral is numerically evaluated by using a Digital Terrain Model (DTM). This can be done using FFT methods in planar approximation as described in e.g. [6, 7].

Alternatively, the TC is computed by quadrature of the integral formula using the DTM in the area $S$ having a radius that depends on the topography roughness (in high mountain range this radius can be 200 km). The TC effect can be thus evaluated as:

$$\mathit{TC}\left(P\right)\cong \sum _{k=1}^n{\mathit{TC}}_k$$E6

where ${\mathit{TC}}_k$ is the volume integral giving the gravitational effect of the $k$-th DTM element in $S$ and $n$ is the total number of DTM elements (see Figure 1). Different mathematical models for the computation of ${\mathit{TC}}_k$ have been proposed. In this paper three of them have been considered and compared, namely the formula of the gravitational effect of a right parallelepiped, a spherical tesseroid and a flat tesseroid.

Given the gravitational potential $V$ of a body $B$ of constant density ρ as the integral:

$$V\left(P\right)=\mathit{G\rho }\underset{B}{\iiint }\frac{1}{l}\mathrm{d}v\left(Q\right)$$E7

$$l=\left|\underset{¯}{r}\left(P\right)-\underset{¯}{r}\left(Q\right)\right|$$

$\underset{¯}{r}\left(P\right)$ = position vector of the computation point;

$\underset{¯}{r}\left(Q\right)$ = position vector of the integration point;

one can compute the gravitational effect of a right parallelepiped assuming the computational point $P$ at the origin of a Cartesian reference system as the closed formula [2]:

$$\frac{\partial V}{\partial z}=\mathit{G\rho }{\left|\Vert xln\left(y+r\right)+yln\left(x+r\right)-zarctan\left(\frac{\mathit{xy}}{\mathit{zr}}\right)\right|}_{x_2}^{x_1}{\left|{}_{y_2}^{y_1}\right|}_{z_2}^{z_1}$$E8

where $r=\sqrt{x^2+y^2+z^2}$ and $\left(x_1,x_2\right),\left(y_1,y_2\right),\left(z_1,z_2\right)$ are the coordinates of the parallelepiped edges.

This formula is implemented in e.g. the TC program of the GRAVSOFT package [8] that will be used in the computations described in the next section.

The gravitational effect of spherical tesseroid has been studied in several papers. Relevant studies on this topic have been carried out by [3, 9, 10, 11, 12]. In this study the approach presented in [11], hereinafter referred to as UNIPOL, is used.

The UNIPOL approach grounds on the result that a closed formula of the volume Newtonian integral (Eq. (8)) exists when the observation point $P$ is located along the polar axis (see e.g. [13]) and the tesseroids coincide with sectors of a spherical zonal band of a spherical cap.

In this case the gravitational contribution of spherical tesseroid of height $h$ can be expressed as:

$$\frac{\partial V}{\partial r}=\frac{\mathit{G\rho }}{r^2}\mathrm{\Delta }\lambda \left(I_1+I_2+I_3+I_4+I_5\right)$$E9

with

$$\begin{array}{cc}{I}_1& =\frac{1}{3}\left[{\left(R+h\right)}^2+r^2\left(3{cos}^2{\phi }_2-2\right)+r\left(R+h\right)cos{\phi }_2\right]\\ & \times \sqrt{{\left(R+h\right)}^2+r^2-2r\left(R+h\right)cos{\phi }_2}\\ I_2& =\frac{1}{3}\left[R^2+r^2\left(3{cos}^2{\phi }_2-2\right)+\mathit{rR}cos{\phi }_2\right]\times \sqrt{R^2+r^2-2r\left(R+h\right)cos{\phi }_2}\\ I_3& =-\frac{1}{3}\left[{\left(R+h\right)}^2+r^2\left(3{cos}^2{\phi }_1-2\right)+r\left(R+h\right)cos{\phi }_1\right]\\ & \times \sqrt{{\left(R+h\right)}^2+r^2-2r\left(R+h\right)cos{\phi }_1}\\ I_4& =\frac{1}{3}\left[R^2+r^2\left(3{cos}^2{\phi }_1-2\right)+\mathit{rR}cos{\phi }_1\right]\times \sqrt{R^2+r^2-2r\left(R+h\right)cos{\phi }_1}\\ I_5& =-cos{\phi }_2{sin}^2{\phi }_2{r}^3\\ & \times ln\left(\frac{\sqrt{{\left(R+h\right)}^2+r^2-2r\left(R+h\right)cos{\phi }_2}+\left(R+h\right)-rcos{\phi }_2}{\sqrt{R^2+r^2-2r\left(R+h\right)cos{\phi }_2}+R-rcos{\phi }_2}\right)\end{array}$$

In Eq. (9), $R$ is the radius of the Earth and $\phi$ and $\lambda$ colatitude and longitude, respectively. When the observation point $P$ is not located along the polar axis, the UNIPOL approach maps each tesseroid defined in the Earth-Centered Rotational reference frame (ECR_Tesseroid) into a sector of a spherical zonal band of a spherical cap in the Earth-Centered $P$-Rotational reference frame, having the same origin $O$ and polar axis coinciding with the line connecting $O$ to the observation point $P$ (ECP_Tesseroid), for which it is possible to use the exact solution in Eq. (9) (Figure 2).

Figure 2.

Geometry used to calculate the contribution of a spherical tesseroid (green rectangle) to the gravitational acceleration at a point P outside a spherically symmetric earth. a) Representation in the earth-Centered rotational reference frame (ECR). b) Representation in the earth-Centered P-rotational reference frame (ECP), defined with respect to the observation point $P$. $R$ stands for the mean radius of the spherical earth; $\phi$ and $\lambda$ stand for colatitude and longitude, respectively.

Mapping can be done by means of two procedures, depending on the spherical distance of the tesseroid from the observation point (Figure 3a).

Figure 3.

Scheme of the UNIPOL approach. a) Set of tesseroids in the earth-centred rotational (ECR) reference frame. The yellow circlet indicates the observation point. The blue dashed circle around the observation point marks the area where tesseroids are at a distance ≤0.1° from the observation point and the ST procedure is required. Cyan and pink colors are used to indicate the tesseroids that require the application of the ST (cyan) and the RT (pink) procedure. b) Tesseroids mapped in the earth-Centered P-rotational (ECP) reference frame ECP and decomposition of some of them into sectors of a spherical band (light blue-colored tesseroids) in the local ECP reference frame. c) Tesseroids mapped in the earth-Centered P-rotational (ECP) reference frame ECP (pink areas) and re-orientation of some of them in the local ECP reference frame (red sectors). Modified from Marotta and Barzaghi [11].

The first procedure (ST procedure) involves a local second-order decomposition of each tesseroid into a number $N_S$ of small equal-area sectors of a spherical band, which develop along the local ECP meridians (dashed black lines) and parallels (dashed blue lines in Figure 3b). Sensitivity tests [11] show that the optimum value of $N_S$ for which the ST procedure converges varies only with the latitude and decreases from the equator to the North Pole. In the present study we use $N_S=$ 4, in agreement with [11].

The second procedure (RT procedure) is based on the rotation of an ECR_Tesseroid (pink areas in Figure 3c) around its center and on its resizing to form a sector of a spherical band (red lines in Figure 3c) that develops along the local ECP-meridian and has the ECP-longitudinal dimension such that the ECP_Tesseroid maintains the same area of the original ECR_Tesseroid.

In agreement with [11], in the present study we use the ST procedure for angular distance between the observational point and the center of the tesseroid less than 0.1°. For angular distance greater than 0.1° we use the RT procedure.

Finally, the Flat Tesseroid (FT) approximation is considered. In this case, the bottom and the top of the tesseroid are flattened thus obtaining a particular polyhedron. A polyhedron is by definition a three-dimensional body with flat polygonal faces, straight edges and sharp corners or vertices and, through proper determination of its vertices, it can realize volumes approximating both right and spherical prisms.

By assuming the FT approximation, one can apply the linear integral algorithm for the computation of the gravitational effect implied by the topographic masses as proved in [4]. The implementation of such approach is based on the double application of the divergence theorem of Gauss to volume integrals of the general expressions for the gravity field resulting from a polyhedral constant density mass and evaluated at an arbitrary point in space. Given the potential of the body B in Eq. (7), the three components of the attraction of this body in a Cartesian orthogonal frame (x, y, z) can be expressed as

$$V_{x_i}\left(P\right)=\mathit{G\rho }\underset{B}{\iiint }\frac{\partial }{\partial x_i}\left(\frac{1}{l}\right)\mathit{dv}\left(Q\right)i=1,2,3$$E10

where $V_{x_i}\left(P\right)$ is the partial derivative of the potential $V\left(P\right)$ along the $x_i$ direction.

The double application of the divergence theorem transforms these expressions first into an equal set of surface integrals (of the same number of faces of the polyhedral source) and subsequently each of them into a set of line integrals defined for each individual segment belonging on that face of the polyhedron. The solution of each line integral produces the final analytical formulas that, in terms of first-derivatives of the potential, are given by [4]:

$$V_{x_i}\left(P\right)=\mathit{G\rho }\sum _{p=1}^{n}cos\left({\mathit{N}}_p,{\mathit{e}}_i\right)\left[\sum _{q=1}^m{\sigma }_{\mathit{pq}}{h}_{\mathit{pq}}L{N}_{\mathit{pq}}+h_p\sum _{q=1}^m{\sigma }_{\mathit{pq}}A{N}_{\mathit{pq}}+sin\left(g_{A_p}\right)\right]$$E11

In Eq. (11), p defines one of the polygonal surfaces $S_p$, q defines one of the segments delimitating the $p$ polygonal surface, ${\mathit{N}}_p$ is the outer unit normal of the polyhedral plane $p$, ${\mathit{e}}_i$ is the unit vector situated in the computation point $P$, ${\sigma }_{\mathit{pq}}$ is equal to −1 when the normal of segment $q$ lying on the plane of polygon $S_p$ points to the half-plane that contains the projection $P’$ of $P$ on $S_p$, ${\sigma }_{\mathit{pq}}$ is equal to +1 otherwise, $h_{\mathit{pq}}$ is the distance between $P’$ and the segment $q$, $h_p$ is the distance between $P$ and the plane $p$. $L{N}_{\mathit{pq}}$ and $A{N}_{\mathit{pq}}$ are transcendental expressions:

$${\mathit{LN}}_{\mathit{pq}}=ln\frac{s_{2\mathit{pq}}+l_{2\mathit{pq}}}{s_{1\mathit{pq}}+l_{1\mathit{pq}}}$$E12

$${\mathit{\text{AN}}}_{\mathit{pq}}=arctan\frac{h_p{s}_{2\mathit{pq}}}{h_{\mathit{pq}}{l}_{2\mathit{pq}}}-arctan\frac{h_p{s}_{1\mathit{pq}}}{h_{\mathit{pq}}{l}_{1\mathit{pq}}}$$E13

The $l_{1\mathit{pq}}$ and $l_{2\mathit{pq}}$ terms are the three-dimensional distances between P and the end points of the segment $\mathit{pq}$, while $s_{1\mathit{pq}}$ and $s_{2\mathit{pq}}$ are the one-dimensional distances between the origin $P”$ of a 1D local coordinate system defined on the segment $\mathit{pq}$ and its two end points. Last term of Eq. (11), $\mathit{sin}\left(g_{\mathit{Ap}}\right)$, is the singularity term that appears for specific locations of $P’$ with respect to the closed polygonal line $G_p$. It expresses the analytical solutions of the corresponding limiting values of the line integrals that are obtained from the partial application of the divergence theorem for a small circle containing the singularity point when its radius tends toward zero. This singularity terms yield finally the values $-2\pi h_p$ when $P’$ lies inside $S_p,$$-\pi h_p$ when $P’$ is located on $G_p$ but not at any of its vertices, $-\vartheta h_p$ when $P’$ is located at one of the vertices of $G_p$ and $\vartheta$ is the angle defined by two subsequent segments that meet at the corresponding vertex, and $0$ when $P’$ is located outside $S_p$. In Eq. (11), the two coordinate transformations occurring in every face of the polyhedron are given by two nested summations, firstly over faces and secondly over segments, of the same transcendental expressions depending on the vertices of the polyhedron. This means properly managing the relative position of points $P$, $P’$ and $P”$ with respect to every surface $S_p$ that can be done algorithmically by standard tools of vector calculus. Figure 4 presents a sketch of their relative positions.

Figure 4.

Sketch of the relative position of points $\mathrm{P}$, $\mathrm{P}’$ and $\mathrm{P}”$ and of the local reference system defined in the computation point $\mathrm{P}$.

Once clarified that the linear integral approach requires the control of the relative position of the vertices of a polyhedral mass with respect to the computation point, in the following a brief description of the implementation of such approach in the framework of residual terrain correction algorithms assessment is provided. Let us consider a set of n points in space P_i (i = 1 to n) on which the gravitational effect of the terrain must be computed and a set of m gridded points $Q_j$ ($j=1$ to $m$) with grid spacings $\Delta {\phi }_j$ in latitude and $\Delta {\lambda }_j$ in longitude describing the DTM considered to compute this effect. The coordinates of $P_i$ and $Q_j$ are expressed in terms of spherical coordinates $\left(\phi ,\lambda ,r\right)$ referring to the same reference system, as well as the angular grid spacings $\Delta {\phi }_j$ and $\Delta {\lambda }_j$. On the basis of the spherical coordinates of $Q_j\left({\phi }_j,{\lambda }_j,r_j\right)$, the grid spacings $\Delta {\phi }_j$, $\Delta {\lambda }_j$ and the radius of the computation point $r_i$, a six-facet polyhedron can be defined. In particular, the eight vertices of this polyhedron will be: $A\left({\phi }_j,{\lambda }_j,r_i\right)$, $B\left({\phi }_j,{\lambda }_j+\Delta {\lambda }_j,r_i\right)$, $C\left({\phi }_j+\Delta {\phi }_j,{\lambda }_j+\Delta {\lambda }_j,r_i\right)$, $D\left({\phi }_j+\Delta {\phi }_j,{\lambda }_j,r_i\right)$, $E\left({\phi }_j,{\lambda }_j,r_j\right)$, $F\left({\phi }_j,{\lambda }_j+\Delta {\lambda }_j,r_j\right)$, $G\left({\phi }_j+\Delta {\phi }_j,{\lambda }_j+\Delta {\lambda }_j,r_j\right)$, $H\left({\phi }_j+\Delta {\phi }_j,{\lambda }_j,r_j\right)$. Note that the two planar surfaces identified by the closed polygons $\mathit{\text{ABCD}}$ and $\mathit{\text{EFGH}}$ have their vertices whose radiuses depend in the first case on the radius of the computation point $P_i$ while in the second case on the radius of the DTM point $Q_j$. Figure 5 illustrates in graphical form how the single polyhedron is built.

Figure 5.

Sketch of the polyhedron vertices building procedure on the basis of the DTM point ${\mathrm{Q}}_{\mathrm{j}}$ and the computation point ${\mathrm{P}}_{\mathrm{i}}$.

This procedure defines the spherical coordinates of the vertices of the polyhedron that are at the height of the terrain and at the height of the computation point, i.e. the level at which the Bouguer plate is computed. The computation of the first-order derivative along the direction of $r_i$ of the gravitational effect of such polyhedron corresponds to the terrain correction to be applied at $P_i$as contribute of $Q_j$. Such value is obtained by running the code polyhedron.f made available by the author [4] implementing the linear integral approach. As input, the relative cartesian coordinates of the polyhedron vertices with respect to the computation point and their topological relationships are required. These are obtained applying a change of reference system to the polyhedron vertices. In particular, their coordinates were roto-translated into a local reference system having the origin at the computation point $P_i$, the $z$ axis pointing up along the direction of $r_i$ and the $x$ and $y$ axes parallel to the local East and North directions, respectively. Regarding the topological relationships defining the outer normal direction of the six planes of the polyhedron, they are defined by a topology matrix containing the counterclockwise sequence of the vertices as seen from outside. As output, the absolute value of the computed $V_r$ is taken. This procedure contemplates two nested loops over all the $m$ DTM points $Q_j$ and the $n$ computation points $P_i$. Within the loop over the computation points, different local reference systems are defined. This leads to slight changes in the directions of the x and y axes but not on the z axis, always normal and pointing outside the reference sphere defined on $P_i$, then maintaining consistency between the different $V_r$ computed by the same DTM point $Q_j$ with respect to the different computation points $P_i$.

The three different models adopted in computing the $\mathit{TC}$ account for a different geometry of the single terrain element. The tesseroid formula (Eq. (9)) considers the radial convergence component of the vertical edges of each terrain element and, besides that, gives a spherical approximation of the top and the bottom of this element. The flat tesseroid has the same geometry in the radial component but top and bottom are planar surfaces. Thus, its geometry is, in principle, less accurate than the one of the tesseroid. Finally, the computation based on the right parallelepiped disregards also the radial convergence of the vertical edges. Hence, this model describes the topography geometry of each element in a way that does not adhere properly the two main features of the given DTM element. However, since in computing the $\mathit{TC}$ effect only the differences between the height of the computation point and the heights of the DTM needs to be considered, the above-mentioned differences between the geometries of the elements used in computing the $\mathit{TC}$ should have a limited impact on the computed values. To prove this, we have devised a test in the Alpine area where the largest discrepancies are expected in the $\mathit{TC}$ effect evaluated with the three different geometries of each single DTM element.

3. The numerical tests on the three proposed approaches

The three different mathematical models presented in Section 2 have been applied in two $\mathit{TC}$ computation tests. Both tests have been carried out in the Alpine area.

In the first test, the SRTM3 DTM (see [14]) have been selected in the area (called AREA_1):

$$46°\le \mathrm{\phi }\le 47°11°\le \mathrm{\lambda }\le 12°$$

The statistics of the height data in this area (see Figure 6) are given in Table 1.

Figure 6.

The DTM (AREA_1) and the points for TC computation.

Number of points	μ [m]	σ [m]	Min [m]	Max [m]
1442401	1502.0	644.1	160.0	3460.0

Table 1.

The statistics of the DTM data in AREA_1.

Points for the TC computation have been selected on a 3′ × 3′ regular grid in the inner area:

$$46.25°\le \mathrm{\phi }\le 46.75°11.25°\le \mathrm{\lambda }\le 11.75°$$

The computation points are thus in an inner area which have 0.25° from the extents of the outer one containing the DTM. So, the terrain correction in points close to the DTM boundaries is not computed according to the standard. However, in the relative comparison between methods, this should not affect the results.

The heights of these 121 prediction points have been assumed coincident with those of the SRTM3 DTM and their statistics are listed in Table 2.

Number of points	μ [m]	σ [m]	Min [m]	Max [m]
121	1471.3	611.7	212.0	2577.0

Table 2.

The statistics of the heights of the computation points in AREA_1.

Given the geometry of tesseroids and flat tesseroids, SRTM3 orthometric heights were transformed into ellipsoidal heights via the EGM96 geoid undulation [15]. Based on the ellipsoidal coordinates $\left({\phi }_{\mathit{\text{ellipsoidal}}},\lambda ,h\right)$ of both DTM and prediction points were converted into spherical coordinates $\left(\phi ,\lambda ,r\right)$ and used in the terrain effect computation with tesseroids and flat tesseroids.

The statistics of the differences among the three methods are given in Table 3.

	μ [mGal]	σ [mGal]	Min [mGal]	Max [mGal]
TC (GRAVSOFT)	13.686	5.194	5.185	35.526
TESSEROID (UNIPOL)	13.553	5.111	5.313	35.065
FLAT TESSEROID (FT)	13.472	5.077	5.275	34.883
TC (GRAVSOFT) - UNIPOL	0.133	0.175	−0.246	0.782
TC (GRAVSOFT) - FT	0.214	0.203	−0.156	1.050
UNIPOL - FT	0.081	0.056	0.006	0.306

Table 3.

The statistics of the TC computations in AREA_1.

As a first overall comment, it can be stated that the results are in good agreement even in such a rough mountain area with sharp height variations (see Tables 1 and 2). By inspecting in more detail the statistics, one can see that values computed by the TC-GRAVSOFT and the UNIPOL approaches are in better agreement than those computed by the TC-GRAVSOFT and the Flat Tesseroid (FT) approaches. The mean of the differences in the first case is nearly 60% of that of the second comparison and the standard deviation is the 86%. This is quite an unexpected result as the geometry of the FT tesseroid is closer to that of the TC prism than to the geometry of the UNIPOL tesseroid. Further investigations will be performed in order to understand this behavior.

On the other hands, the terrain correction values based on the UNIPOL tesseroid procedure and the ones obtained with the FT approach agree very well.

The mean and the standard deviation of the differences between these terrain corrections are of the order of some hundredth of mGal and the maximum difference is of the order of one third of mGal. Although in principle this is quite foreseen, it is important to quantify the differences in view of practical applications.

Larger differences can be seen when comparing these two methods with the one based on prism effect. In these cases, the maximum differences are of the order of 1 mGal. Even though this value is high if compared with the precision of the gravity observations (which can reach few μGals), one has to consider that other error sources in the topography reduction process can have a larger impact. As an example, the discrepancy between the heights of the point associated with the gravity observations as compared with those obtained by the DTM in the same points can amount to ten meters (or even more) in mountain areas. Given that the absolute value of the free-air gradient is $0.30877\text{mGal}/\mathrm{m}$, this implies 3 mGal in 10 m due to this mismatch. Also, biases can occur due to the assumption of constant density. In view of that, even the maximum difference between the GRAVSOFT terrain correction and the spherical tesseroid/flat tesseroid values are not so significant.

A second test was then devised. Observed gravity data were selected in the area (AREA_2):

$$46°\le \mathrm{\phi }\le 47°11°\le \mathrm{\lambda }\le 12°$$

Gravity point coordinates were surveyed with GNSS and framed to ITRF94. Statistics of the ellipsoidal heights of these gravity points are listed in Table 4.

Number of points	μ [m]	σ [m]	Min [m]	Max [m]
116	1161.74	384.22	312.67	2217.08

Table 4.

The statistics of the heights of the computation points in AREA_2.

Gravity values have been measured with a Lacoste&Romberg G-367 relative gravimeter. The standard deviation of the observed values is of the order of 0.02 mGal. Gravity data are referred to IGSN71 and their statistics are summarized in Table 5.

Number of points	μ [mGal]	σ [mGal]	Min [mGal]	Max [mGal]
116	980404.459	71.605	980208.060	980545.806

Table 5.

The statistics of the observed gravity values.

For the computation of the terrain component, the SRTM3 DTM have been selected in the 3° × 3° area centered on the one containing the gravity data area (AREA_2)

$$45°\le \mathrm{\phi }\le 48°10°\le \mathrm{\lambda }\le 13°$$

The statistics of the SRTM3 in AREA_2 are described in Table 6.

Number of points	μ [m]	σ [m]	Min [m]	Max [m]
12967201	1054.4	853.2	−35.0	3865.0

Table 6.

The statistics of the DTM data in AREA_2.

Figure 7 shows the DTM features of AREA_2 and the position of the gravity points.

Figure 7.

The DTM (AREA_2) and the points for TC computation.

Similarly to what has been done in the first test, SRTM3 and gravity point coordinates were transformed into spherical coordinates for the computation of the terrain correction with the UNIPOL and FT approaches. On the other hands, ellipsoidal heights of gravity points have been converted into orthometric heights via the EGM96 geoid undulations when computing the terrain correction with the TC software of the GRAVSOFT package.

Given the three different terrain corrections, by applying Eq. (4), three different sets of Bouguer anomalies have been derived.

In Eq. (4), as previously stated, the orthometric height $H$ is derived via the EGM96 geoid undulation and the Bouguer plate is accounted for by using Eq. (1). We further assumed that

$$\frac{\partial \gamma }{\partial h}=-0.30877\mathit{\text{mGal}}/m$$E14

and the normal gravity in a point $Q$ of latitude ${\phi }_Q$ on the ellipsoid is given by the GRS80 normal gravity formula [16]:

$$\gamma \left(Q\right)=978032.7\left(1+0.0053024{\mathit{sin}}^2{\phi }_Q-0.0000058{\mathit{sin}}^{2}2{\phi }_Q\right)\mathit{\text{mGal}}$$E15

Although this formula has an accuracy of 0.1 mGal (see [16]), it can be used in the context of this relative comparison among different terrain correction computation methods.

Table 7 summarizes the statistics of the Bouguer anomalies obtained with the three terrain correction methods.

$∆g_B$	μ [mGal]	σ [mGal]	RMS [mGal]	Min [mGal]	Max [mGal]
TC (GRAVSOFT)	−137.846	16.131	138.786	−168.873	−90.511
TESSEROID (UNIPOL)	−137.074	16.205	138.028	−166.337	−90.299
FLAT TESSEROID (FT)	−136.894	16.221	137.852	−166.262	−90.069

Table 7.

The statistics of the Bouguer anomalies.

Comments similar to those given on Table 3 hold for the Bouguer values in Table 7. The Bouguer anomalies obtained by applying the three methods have quite similar statistics. Those computed via TC-GRAVSOFT software have the smallest standard deviation and the highest mean while those obtained with the other two methods have smaller mean and higher standard deviations. If the RMSs are considered, one can see that the Bouguer anomalies based on the flat tesseroid have the smallest value. However, as pointed out before, even the largest difference among the maximum values (around 2.6 mGal) is not so significant if compared with other biases occurring in the Bouguer reduction.

Thus, the statistics of Tables 3 and 7 prove the substantial equivalence of the three approaches used for the $\mathit{TC}$ computation.

4. Conclusions

Three different methods for terrain correction have been compared in two areas over the Alps. The standard computation given by the TC-GRAVSOFT program has been compared with the terrain corrections evaluated via spherical tesseroid and flat tesseroid formulas. In the first test, the SRTM DTM was clipped in a 1° × 1° window and $\mathit{TC}$effect was computed in a set of gridded points in the same area. In the second test, observed gravity values in a 1° × 1° area have been used in the computation of Bouguer anomalies considering the 3° × 3° SRTM DTM values centered on the area containing the gravity data. Despite the fact that the topography in the two selected DTM windows is quite rough, no significant differences among the methods have been revealed. The statistics of the values obtained by modeling in different ways the shape of the discretized topography elements are practically equivalent. Differences among $\mathit{TC}$ effects and Bouguer anomalies computed with parallelepiped, spherical tesseroid and flat tesseroid amount to maximum values that are around 1 and 3 mGal respectively. As a matter of fact, there are other error sources (e.g., density heterogeneities, DTM and gravity point heights mismatch) that can have impacts on the terrain correction computation larger than 3 mGal. However, if in the second test on Bouguer computation we consider the values per se, spherical tesseroid and flat tesseroid models perform slightly better when RMS values are compared, i.e. the spherical tesseroid and flat tesseroid based Bouguer anomalies are smoother.

Finally, we remark that the concept applied in the flat tesseroid modeling can be adapted to the terrain effect computation when shaping the topography according to the Triangulated Irregular Network model [17]. In this way, a more detailed terrain effect evaluation will be possible, particularly in the neighbor of the computational points, by better modeling the terrain slopes.

Acknowledgments

We thank DTU Space for providing us with the TC program of the GRAVSOFT package.

Conflict of interest

The authors declare no conflict of interest.