**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 *HP* 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 *AB* is given by the well-known formula (see e.g. [1]):

$$A\_B = 2\pi G\rho H\_P \tag{1}$$

It must be underlined that both the plate and the TC effects are usually com-

The reduced gravity data are then moved to the geoid by considering the gravity

*<sup>∂</sup><sup>H</sup> HP* ffi *g P*ð Þ� *AB* <sup>þ</sup> *TC P*ð Þ� *<sup>∂</sup><sup>γ</sup>*

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

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

gradient, which is approximated by the free-air normal gravity gradient *∂γ/∂h*:

Finally, the normal gravity at point *Q* (see **Figure 1**) on the ellipsoid is

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

> ððþ<sup>∞</sup> �∞

ð*z*

*ζ* � *zP l*

, *xP*, *yP*, *zP*

*zP*

the computational point *P* and ð Þ *ξ*, *η*, *ζ* 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

*TC P*ð Þffi <sup>X</sup>*<sup>n</sup>*

*k*¼1

where *TCk* 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 *TCk* 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

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

*B* 1 *l*

*l* ¼ j j *r P*ð Þ� *r Q*ð Þ

*V P*ð Þ¼ *Gρ*∭

*r P*ð Þ = position vector of the computation point; *r Q*ð Þ = position vector of the integration point;

<sup>Δ</sup>*gB* <sup>¼</sup> *g P*ð Þ� *AB* <sup>þ</sup> *TC P*ð Þ� *<sup>∂</sup><sup>γ</sup>*

. This means that residual effects

*<sup>∂</sup><sup>h</sup> HP* (3)

*<sup>∂</sup><sup>h</sup> HP* � *<sup>γ</sup>*ð Þ *<sup>Q</sup>* (4)

<sup>3</sup> d*ξ*d*η*d*ζ* (5)

� � are the coordinates of

*TCk* (6)

d*v Q*ð Þ (7)

puted at constant density *ρ* set at 2670 kg/m<sup>3</sup>

*DOI: http://dx.doi.org/10.5772/intechopen.97718*

*g P*ð Þ� *AB* <sup>þ</sup> *TC P*ð Þ� *<sup>∂</sup><sup>g</sup>*

effect is given by the integral formula:

subtracted and the standard Bouguer anomaly is obtained as:

*TC P*ð Þ¼ *Gρ*

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *xP* � *<sup>ξ</sup>* <sup>2</sup> <sup>þ</sup> *yP* � *<sup>η</sup>* � �<sup>2</sup> <sup>þ</sup> ð Þ *zP* � *<sup>ζ</sup>* <sup>2</sup>

values *g(P)*.

where ¼

evaluated as:

tesseroid.

**45**

q

Where G = 6.67 � <sup>10</sup>�<sup>11</sup> <sup>m</sup><sup>3</sup> /kg s2 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 *HP* (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:

$$\text{g}(P) - A\_B + TC(P) \tag{2}$$

**Figure 1.** *The plate and the terrain correction (HP = orthometric height, hP = ellipsoidal height).*

*The Gravity Effect of Topography: A Comparison among Three Different Methods DOI: http://dx.doi.org/10.5772/intechopen.97718*

It must be underlined that both the plate and the TC effects are usually computed at constant density *ρ* set at 2670 kg/m<sup>3</sup> . 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*:

$$\text{g}(P) - A\_B + TC(P) - \frac{\partial \text{g}}{\partial H} H\_P \cong \text{g}(P) - A\_B + TC(P) - \frac{\partial \text{g}}{\partial h} H\_P \tag{3}$$

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

$$
\Delta \mathbf{g}\_B = \mathbf{g}(P) - A\_B + TC(P) - \frac{\partial \chi}{\partial h} H\_P - \chi(\mathbf{Q}) \tag{4}
$$

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:

$$TC(P) = G\rho \left\| \int\_{-\infty}^{+\infty} \int\_{x\_P}^{x} \frac{\zeta - x\_P}{l^3} \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta \right\|\,\tag{5}$$

where ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *xP* � *<sup>ξ</sup>* <sup>2</sup> <sup>þ</sup> *yP* � *<sup>η</sup>* � �<sup>2</sup> <sup>þ</sup> ð Þ *zP* � *<sup>ζ</sup>* <sup>2</sup> q , *xP*, *yP*, *zP* � � are the coordinates of the computational point *P* and ð Þ *ξ*, *η*, *ζ* 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:

$$TC(P) \cong \sum\_{k=1}^{n} TC\_k \tag{6}$$

where *TCk* 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 *TCk* 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(P) = G\rho \iint\_{B} \frac{1}{l} \mathrm{d}\nu(Q) \tag{7}$$

$$l = |\underline{r}(P) - \underline{r}(Q)|$$

*r P*ð Þ = position vector of the computation point; *r Q*ð Þ = position vector of the integration point;

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 formulas giving the gravitational effect of a right parallelepiped, a spherical

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 *HP* 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 *AB* is given by the well-known formula

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 *HP* (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:

*The plate and the terrain correction (HP = orthometric height, hP = ellipsoidal height).*

*AB* ¼ 2*πGρHP* (1)

/kg s2 is the universal gravitational constant. The

*g P*ð Þ� *AB* þ *TC P*ð Þ (2)

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**

(see e.g. [1]):

**Figure 1.**

**44**

Where G = 6.67 � <sup>10</sup>�<sup>11</sup> <sup>m</sup><sup>3</sup>

the terrain correction in the framework of the Bouguer reduction.

*Geodetic Sciences - Theory, Applications and Recent Developments*

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 \mathbf{z}} = G\rho \Big| \left| \mathbf{x} \cdot \ln \left( y + r \right) + y \cdot \ln \left( \mathbf{x} + r \right) - z \cdot \arctan \left( \frac{\mathbf{x} \mathbf{y}}{\mathbf{z} \mathbf{r}} \right) \Big|\_{\mathbf{x}\_2}^{\mathbf{x}\_1} \Big|\_{\mathbf{y}\_2}^{\mathbf{x}\_1} \tag{8}$$

and polar axis coinciding with the line connecting *O* to the observation point *P* (ECPTesseroid), for which it is possible to use the exact solution in Eq. (9) (**Figure 2**). Mapping can be done by means of two procedures, depending on the spherical

*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; φ*

*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*

The first procedure (ST procedure) involves a local second-order decomposition of each tesseroid into a number *NS* 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 *NS* 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 *NS* ¼ 4, in

distance of the tesseroid from the observation point (**Figure 3a**).

*DOI: http://dx.doi.org/10.5772/intechopen.97718*

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

agreement with [11].

**Figure 2.**

**Figure 3.**

**47**

*Marotta and Barzaghi [11].*

*and λ stand for colatitude and longitude, respectively.*

where *<sup>r</sup>* <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *<sup>x</sup>*<sup>2</sup> <sup>þ</sup> *<sup>y</sup>*<sup>2</sup> <sup>þ</sup> *<sup>z</sup>*<sup>2</sup> <sup>p</sup> and ð Þ *<sup>x</sup>*1, *<sup>x</sup>*<sup>2</sup> , *<sup>y</sup>*1, *<sup>y</sup>*<sup>2</sup> � �,ð Þ *<sup>z</sup>*1, *<sup>z</sup>*<sup>2</sup> 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–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{G\rho}{r^2} \Delta \lambda (I\_1 + I\_2 + I\_3 + I\_4 + I\_5) \tag{9}$$

with

$$\begin{aligned} I\_1 &= \frac{1}{3} \left[ (R+h)^2 + r^2 \left( 3\cos^2 \varphi\_2 - 2 \right) + r(R+h) \cos \varphi\_2 \right] \\ &\times \sqrt{(R+h)^2 + r^2 - 2r(R+h) \cos \varphi\_2} \\\\ I\_2 &= \frac{1}{3} \left[ R^2 + r^2 \left( 3\cos^2 \varphi\_2 - 2 \right) + rR \cos \varphi\_2 \right] \times \sqrt{R^2 + r^2 - 2r(R+h) \cos \varphi\_2} \\\\ I\_3 &= -\frac{1}{3} \left[ (R+h)^2 + r^2 \left( 3\cos^2 \varphi\_1 - 2 \right) + r(R+h) \cos \varphi\_1 \right] \\\\ &\times \sqrt{(R+h)^2 + r^2 - 2r(R+h) \cos \varphi\_1} \\\\ I\_4 &= \frac{1}{3} \left[ R^2 + r^2 \left( 3\cos^2 \varphi\_1 - 2 \right) + rR \cos \varphi\_1 \right] \times \sqrt{R^2 + r^2 - 2r(R+h) \cos \varphi\_1} \\\\ I\_5 &= -\cos \varphi\_2 \sin^2 \varphi\_2 r^3 \\\\ &\times \ln \left( \frac{\sqrt{(R+h)^2 + r^2 - 2r(R+h) \cos \varphi\_2} + (R+h) - r \cos \varphi\_2}{\sqrt{R^2 + r^2 - 2r(R+h) \cos \varphi\_2} + R - r \cos \varphi\_2} \right) \end{aligned}$$

In Eq. (9), *R* is the radius of the Earth and *φ* and *λ* 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 (ECRTesseroid) 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*
