Creation of the Gravimetric Network and Determination of the GEOÏD in the Mediterranean Basin (Region-Tunisia)

Marcel Ghannem¹, Atef Belhaj Ali²

¹Department de Geology, Faculty of Sciences University of El Manar, Tunis, Tunisie.
²Department of Geography and GIS, Faculty of Arts and Humanities, King Abdulaziz University, Jeddah, KSA.

DOI: 10.4236/jgis.2020.123010                

Abstract

The determination of geoid models with great precision (centimeter) was always at the center of interest in geodesy research [1] [2]. One of the most used methods to calculate the geoid model is the method called Remove-Compute-Restore (R-C-R). This method applies the stokes’ integral formula by the use of short, medium and long wavelength information via the three main stages R-C-R. The GravSoft software implemented by [3] was used for this study. Geodetic heights, a digital terrain model (SRTM) and leveled GPS points were used as inputs. The geoid modeling was carried out on the North region of Tunisia (Grand Tunis and Bizerte) on an area of 83 × 83 km. The accuracy of the quasi-geoid provisional reached 3.1 cm.

Keywords

Gravimetric Network, Geoid, Geodesy, Grand Tunis, Bizerte, Stokes Integral, R-C-R Method

1. Introduction

The most accurate approximation of the actual shape of the earth is the geoid. It is so called equipotential surface of gravity field [4]. It is affected by hollows and bumps, it is therefore completely irregular. We can imagine it as being the average level of the oceans and its imaginary extension under the continents. The geoid, whose mathematical definition is relatively complex, is not very easy to use. Its use was mainly reserved for research on vertical references and mean sea level [5] [6]. Since the appearance of spatial positioning techniques, and more particularly with the rapid development of the GPS system, the situation has radically changed: the geoid has become an essential tool for converting heights of differences in heights from GPS into altitudes or differences in altitudes [7].

Given the importance of knowing the height of the geoids with respect to the reference ellipsoid in the geodetic works, especially for the reduction of distances, the upgrading work of the Tunisian geodesy and the establishment of the spatial network, we have proceed to its determination according to the needs and the dispensability of the data:

· The 80s, we calculate the heights of the geoids on a certain set of points of the Tunisian Geodesic Primordial network by the method of astro-geodesic leveling which is based on the determination of the difference of the height of the geoids between two points compared to the ellipsoid Clarke 1880. The precision is estimated at 17.5 cm/100 km.

· In 2004, establishment of a local geoids model over the Tunis area by adjusting the global model EGM96 to GPS-derived leveling points, this model had a standard deviation of 10 cm.

· -The calculation is made by assimilating the geoids to a plane passing closer to the points of support and by calculating the distance between the searched point and the plan modeling the geoids. 28 points of the place on all Tunisia are used as points of support what gave a maximum error of 45 cm on the point of Gsar Ghilene and a standard deviation of 22 cm on the 28 points.

2. Study Area

For the calculation of the quasi-geoid, the zone is chosen according to the availability of gravimetric data, Situated in the North of Tunisia (Figure 1), it is between latitudes from 36.61˚ to 37.36˚ and longitudes 9.64˚ to 10.39˚. The study area covers an area of approximately 9400 km²: The marine part (southwest corner) will be excluded from the operation of the grid because it is devoid of gravity points.

Red stars represent gravimetric points with a precision of 0.2 mGal.

Figure 1. Study area (background image: GoogleEarth images).

3. Data and Methodology

3.1. Data and Software

3.1.1. Data

The national gravity network consists of three networks types:

An absolute network:

An absolute zero order gravimetric network of Tunisia: made up of 09 points distributed throughout the Tunisian territory. It was produced in cooperation with the Strasburg Globe Physics Institute in 2009. This network is subdivided into three sub-networks:

1^st Order network:

A first order network composed by 66 points distributed homogeneously throughout the national territory is based on the absolute network with a mesh network between 200 and 300 km.

2^nd Order network:

A secondary order network is based on the first two networks with a mesh network between 70 and 100 Km.

3^rd Order network:

Based on previous orders it is a mesh network of 5 to 10 km. for example, a map at 1/50,000 scale should have about 130 gravity points. The number of gravity points needed to cover the Tunisian territory is therefore of the order of 32,000 points. The number of gravity points made since the beginning of the project engaged by the Tunisian national topographic and cartographic office in 2009 is 1200 points covering the northern, north-eastern and central parts of the country.

Gravity network

To establish the gravity network, we have three relative gravimeters CG5 AUTOGRAV (Figure 2) which have the specifications below:

Figure 2. CG5 Gravimeter (CG-5 Manual, SCINTREX Limited: https://scintrexltd.com/).

· Field instrument: portability, autonomy, 12 V power supply, robustness, thermal insulation.

· “Geodetic” instrument (measurement anywhere in the world): measurement range up to 7000 - 8000 mGal.

· Automated digital instrument: microprocessor, memory.

· Means of N measurements (1 Hz), Std Dev, releases.

· Storage and pre-processing of measurements.

· Multi-parameter measurement sensor: gravity + internal temperature and tilt of the sensor, GPS (option).

· “On-line” digital corrections (Earth tide, Temp & Tilts instrument, seismic noise filtering, etc.)

· Resolution: 1 μGal/Standard deviation: 5 μGal.

· 2 modes of use: “Field mode” and “Cycling mode” (continuous recording).

· Fieldwork consists of two phases:

· Construction:

The gravimetric points are constructed of concrete after having chosen their location, generally on a scale of 1/50,000 to guarantee a more or less homogeneous distribution, sometimes existing geodetic points are exploited in the gravity network.

· Observations:

The observations are made by triangulation, for each triangle is assigned a gravimetric closure which will be used for the calculation and compensation of the network.

To determine the local geoid model of Grand Tunis and Bizerte, we used 12 leveling marks whose ellipsoidal heights are known (used for the validation of the accuracy of the conversion grid we will use) and a Digital Terrain Model (DTM) showing the distribution of the altitude in the study area, which is obtained from USGS (https://gdex.cr.usgs.gov/gdex/). The DTM is with a grid of 30 m (Figure 3).

3.1.2. Software

The software used for calculation is GravSoft [3]. GravSoft is a software suite that has been developed and supplemented since the 1970s in order to perform physical geodesy calculations, in particular the determination of the gravimetric geoid.

This software has been used by many organizations in several countries with success and it contains many programs with several choices of methods that can be used. The software is written in FORTRAN and is composed of independent programs; in 2008 a graphical interface was created by the University of Copenhagen (Figure 4). This software was provided by the international geoid service (official service of the International Association of Geodesy).

3.2. Methodology

The measurement of the gravity field on a global scale has become possible due

Figure 3. DTM of the study area.

Figure 4. main interface of GravSoft Software.

to the evolution of spatial gravimetric measurement techniques. These models contain a large part of the content of the gravimetric signal, which motivated the fact that the method of Withdrawal-Calculation-Restoration uses, instead of complete anomalies, the residual anomalies of the gravitational field compared to the global models. The principle of the Remove-Compute-Restore method, called R-C-R in the rest of this paper, is to combine the information which covers three different wavelength ranges [8]:

· the low frequencies given by global models;

· the high frequencies which come mainly from the topography;

· The medium frequencies accessible by gravimetric measurements.

First, the residual anomalies are kept in the gravity anomalies by removing the anomalies from the global model and the corrections from the terrain. These anomalies are then interpolated to a regular grid. The Stokes integral in the next step calculates the residual altitude anomalies. Finally, the altitude anomalies from the topography and the global model are added to those of Stokes in order to have the heights of the quasi-geoid. This method is illustrated by Figure 5: Let ∆g be the anomalies in the open air at the measurement points.

These anomalies can be written using Equation (1), depending on the disturbing potential:

${\Delta }_g=-\frac{\partial T}{\partial r}-\frac{2}{r}T$ (1)

The signal is considered to be composed by three wavelength ranges:

$T=T_M+T_{RT}+T_{\mathrm{Re}s}$ (2)

T_M: is the long-wavelength part determined by spatial gravimetry (M for model).

T_RT: is the short wavelength part that comes from the topography.

T_Res: is the medium wavelength part.

In the first step we calculate the residual anomalies such as:

${\Delta }_{g_{Res}}={\Delta }_g-{\Delta }_{g_M}-{\Delta }_{g_{RT}}$ (3)

In the second step, the Stokes integral is applied to the values of interpolated ${\Delta }_{g_{Res}}$ and we obtain ${\xi }_{Res}$, the residual altitude anomaly value can be written as

Figure 5. general workflow.

${\xi }_{Res}=\frac{T}{\gamma }=\frac{R}{4\pi \gamma }{\displaystyle {\int }_{sigma}S\left(\psi \right){\Delta }_{g_{Res}}}$ (4)

At the last stage, the altitude anomalies can be written as:

$\xi ={\xi }_M+{\xi }_{Res}+{\xi }_{RT}$ (5)

“Remove” Calculation

The first step in the R-C-R method consists in calculating the residual anomalies ∆_gRes at the points of the measurements:

${\Delta }_{g_{Res}}={\Delta }_g-{\Delta }_{g_M}-{\Delta }_{g_{RT}}$ (6)

- ∆_g: the anomaly in the open air;

- ∆_gM: the anomaly of the global model which contains the long wavelength part. This part is obtained by taking the fundamental equation of geodesy and applying it to the decomposition into spherical harmonics of the field:

$\begin{array}{l}{\Delta }_{g_M}=\left(r,\theta ,\lambda \right)\\ =\frac{GM}{r2}{\displaystyle {\sum }_{n=2}^{n_{\max }}\left(n-1\right){\left(\frac{a}{r}\right)}^n}{\displaystyle {\sum }_{m=0}^n{P}_{n,m}\cos \theta }\left(\Delta C_{n,m}\cos m\lambda +\text{Δ}{S}_{n,m}\sin m\lambda \right)\end{array}$ (7)

where: G is the universal gravitational constant, M the mass of the Earth, r is the radius, n and m are the degree and order of the geopotential model, a is a scale factor of the field model, n_max is the maximum degree of development in spherical harmonics, P_nm is the function of Legendre, ∆C_n, , ∆S_n,m are the coefficients of development;

- ∆_gRT: contains the short wavelengths of the disturbing field created by the topography. To calculate this component, its long wavelength effects must be removed from the ground, since they are already taken into account in ∆_gM from the field model. The terrain is filtered in order to keep only the short wavelengths of the relief:

${\Delta }_{g_{RT}}=-\frac{\partial T_{RT}}{\partial R}-\frac{2}{R}{T}_{RT}$ (8)

With

$-\frac{\partial T_{RT}}{\partial R}=G{\displaystyle {\iiint }_{RT}\frac{\rho d\tau }{l^3}}$ (9)

$T_{RT}=G{\displaystyle {\iiint }_{RT}\frac{\rho d\tau }{l}}$ (10)

Before being integrated by the Stokes formula, the residual point anomalies must be transformed into a regular grid by interpolation.

The residual altitude anomalies on the grid points are then calculated from the Stokes formula (Equation (11)) applied to the residual anomalies at these points. It is first considered that the area around the gravity station P is flat and horizontal and that the density between the geoid and the earth’s surface is constant ρ = 2.67 g/cm³. The Bouguer plateau exerts an attraction (2πGρh). All the masses above the geoid are removed (the simple Bouguer reduction) [9]. Simple Bouguer anomalies are calculated by the Equation (11) [10]:

${\Delta }_g=g_p-\left(g_o-0.3086h+2\pi G\rho h\right)$ (11)

To calculate the terrain correction we perform the direct summation of the contributions of each model element (prism, tesseroides etc.).

4. Case Study: Region of GRAND TUNIS and BIZERTE-Tunisia

4.1. Geoid Determination

As mentioned in the general process of the R-C-R method (Figure 5), the gravimetric anomalies used in the calculation are the free air anomalies, their determination is ensured by a translation of the Somigliana and free air formulas described above which give the value of g₀ (gravity on the ellipsoid) and the value of the anomalies. The following figure (Figure 6) shows the free air anomalies.

“Remove” calculation:

This step is the most important step in the geoid calculation; in fact it requires an optimized manipulation of the different types of data as well as a good configuration of the software.

Calculation of the reduction due to the global model:

This step allows subtracting anomalies due to the global model from anomalies in the free air.

Calculation of terrain corrections at gravimetric stations:

Terrain correction (TC) makes it possible to calculate the effects of the terrain on the potential and its derivatives. We detail here the parameters in the case of

Figure 6. Free air anomalies.

application of the residual ground. The effects of the terrain are represented by straight prisms taking into account the local curvature of the Earth. The computation of the terrain corrections is performed by numerical integration using the formulas of the gravitational attraction of the right prisms to represent the terrain in the near area and approximate formulas in the far area.

Terrain corrections are obtained by a simple summation of the contribution of these elements. This allows us to calculate residual anomalies at gravity stations which will be interpolated over the entire study area in order to transform them into residual ripples through the Stokes formula.

The following figure (Figure 7) shows the residual anomalies:

Integration:

The Stokes module makes it possible to calculate residual altitude anomalies from the grid of residual gravity anomalies by a discretized convolution with subsampling in the vicinity of zero.

It is important here to point out that the choice of the integration radius is very important; in fact it directly influences the precision of the geoid.

The choice here was based on a publication by Z. Ismail and O. Jamet [11] which introduces the choice of the radius of integration according to the nature of the terrain (mountainous, semi-mountainous or plain) and the used global geopotential model (Figure 8).

Calculation of the terrain effect:

The terrain effect is calculated using the RTM method; the principle is to calculate this effect from the two DEMs, filtered DEM and reference DEM, over the entire study area.

Figure 7. Residual anomalies.

The result of this calculation is a grid of corrections (Figure 9) to restore for the determination of the final geoid.

The radii R1 and R2 are fixed at 0 and 999 km to optimize the corrections.

Calculation of the undulations of the global model:

These undulations are obtained from the file of coefficients in spherical harmonics by using the GEOEGM module of GravSoft software. Figure 10 shows the undulations of the global model.

Restoration

Restoration is the last step in the R-C-R method. It consists in restoring the anomalies of altitudes calculated by combining the three grids, that of the anomalies of residual altitudes, the corrections of the ground and the anomalies of

Figure 8. Integration radius.

Figure 9. Terrain corrections.

Figure 10. Undulations of the global model EGM96.

altitudes due to the global model.

This step is used to determine the provisional geoid which will later be adjusted to leveling marks. The following figure illustrates the provisional geoid (Figure 11).

4.2. The Geoid Accuracy

The main factors that determine the accuracy of a geoid model are mainly:

· The resolution of the digital terrain model, the higher it is, the more the terrain corrections are punctual and precise.

· The accuracy of gravimetric measurements.

· The density of the gravimetric points and their distributions on the study area: an optimized density and a homogeneous distribution lead to a precise interpolation.

For our case study two factors can have an influence on the precision of the calculated model:

· The absence of a high resolution digital terrain model (terrain resolution 30 m × 30 m).

· The non-homogeneous distribution of gravity points:

- Points spaced 1.5 km apart, others 20 km apart.

- Areas almost devoid of points (Tunis area).

To assess the accuracy of the calculated quasi-geoid, a comparison with other models calculated in different countries makes it possible to rank in a precision scale, for this, a research made by Z. Ismail within the framework of a PhD thesis [12] allowed to know the order of precision of certain models calculated in different countries the table below (Table 1) summarizes the results of this research of which I checked the bibliographical references and we selected the models having conditions similar to our case study.

Verification of the accuracy of the model is ensured by calculating the undulation of the geoid in leveling marks with undulations considered to be known (surveyed by GPS) distributed almost throughout the study area (Figure 12).

Figure 11. provisional geoid.

Table 1. Accuracy of geoid models.

Figure 12. Leveling marks used for verification.

The following table (Table 2) summarizes the difference between the values of the undulations.

The bias is the average of the deviations from the model to the data. It identifies a gap between the model and the data.

$biais=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{e}_i}$ (11)

With:

n: the number of control points, e_i: deviation from the benchmarks. We have a bias of −0.0005 cm, that we consider almost negligible.

The Root Mean Square RMS provides an indication of the accuracy of the model. It is defined by:

Table 2. Verification with leveling marks.

$RMS=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{e}_i^2}}$ (12)

The RMS is about of 0.031 m (3.1 cm)

The maximum deviation is defined by: Maximal deviation = Max |ei| = 0.069 m.

In terms of justification of the value of the standard deviation (3.1 cm) we can note that:

o The area of the nearest leveling mark FF13 is almost devoid of gravity point (distance of almost 20 Km), which has produced a wide interpolation residual anomalies.

Gravimetric points are not compensated for a block compensation so as to homogenize their precision (points of different order).

o The altimetric network has not yet undergone block compensation, in fact the calculation and compensation are done by tracking and not by mesh.

o The ellipsoidal heights of the leveling marks are observed in RTK mode which produced an error of the order of 2 cm on the undulation at these points.

So, according to the statistics of the deviations, we can assimilate the accuracy of the provisional geoid to: 3.1 cm. Figure 13 the deviation between exact and calculated undulations, this deviation is superimposed with the spatial distribution of the control network.

After having assessed the accuracy of the provisional quasi-geoid, an adaptation of this surface to the leveling marks is performed.

The following table (Table 3) shows the exact undulations of the leveling marks as well as the calculated undulations after adjusting the provisional geoid to the leveling marks.

Figure 13. Deviation between exact and calculated undulations

Table 3. Exact and calculated undulations of the leveling marks after adjusting provisional geoid.

Figure 14. Adjusted geoid.

Figure 15. Undulations of the geoid after adjustments.

Figure 14 and Figure 15 illustrate the adjusted geoid and the undulations of the geoid after adjustments.

5. Conclusions

This modeling is promising if we have a gravimetric network and an altimetric network compensated in block, the existence of a high resolution DTM will also make it possible to calculate a high precision geoid model. Determining the exact undulation at the leveling marks by GPS survey in static mode will allow a millimeter adjustment of the altimetric conversion grid.

In any case, this local quasi-geoid model has given encouraging results given the novelty of this kind of treatment on a national scale. This work can be considered as a basis for a generalization of this prototype model especially using the R-C-R method, the manipulation of the GRAVSOFT software and the selection of the data used.

This conversion grid makes it possible to make precise topographic surveys with GPS (standard deviation = 0.015 m) especially on the area where there is a good resolution of the gravimetric points, the precision can reach a few millimeters.

Finally, some recommendations are to cite for optimal accuracy of a national geoid, hence a great performance of GPS for altimetric surveys and geodetic objectives:

· Calculation of gravity observations and compensation in a single block of the completed part of the network, final compensation will be carried out upon completion.

· The compensation of all the meshes of the leveling of precision in a single block by the theory of least squares.

· Carry out a national program with short, medium and long-term objectives throughout the territory in this area.

· The establishment of a specific DTM throughout Tunisian national territory.

Having compared the different methods and using new technologies, we can see that the combinations of the appropriate modes lead to high-performance precision which we recommend to generalize in order to promote a sustainable geodesic and topographic sector.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Rebaï, N., Zenned, O., Trabelsi, H. and Achour, H. (2018) Computing Local Geoid Model Using DTM and GPS Geodetic Points. Case Study: Mejez El Bab-Tunisia. International Journal of Geosciences, 9, 161-178. https://doi.org/10.4236/ijg.2018.93011
[2]	Sabril, L.M., Helianil, L.S., Sunantyol, T.A. and Widjajanti, N. (2019) Geoid Determination with Hotine’s Integral Based on Terrestrial Gravity Data in Semarang City. Journal of Physics: Conference Series, 1127, Article ID: 012047. https://doi.org/10.1088/1742-6596/1127/1/012047
[3]	Forsberg, R. and Tscherning, C. (2008) An Overview Manual for the GRAVSOFT Geodetic Gravity Field Modelling Programs. Report, 2nd Edition, National Space Institute (DTU-Space), Denmark. https://www.scribd.com/document/207911750/An-Overview-Manual-for-the-GRAVSOFT
[4]	Gurtner, W. (1978) Das Geoid in der Schweizerland Institute fur Geodesie une photogrammetrie. Mitteilungen Nr 20, Zurich.
[5]	Li, X. and Hans, J.G. (2001) Ellipsoid, Geoid, Gravity, Geodesy, and Geophysics. Geophysics, 66, 1660-1668. https://doi.org/10.1190/1.1487109
[6]	Bingham, R.J., Haines, K. and Hughes, C.W. (2008) Calculating the Ocean’s Mean Dynamic Topography from a Mean Sea Surface and a Geoid. Journal of Atmospheric and Oceanic Technology, 25, 1808-1822. https://doi.org/10.1175/2008JTECHO568.1
[7]	Nahavandchi, H. and Soltanpour, A. (2006) Improved Determination of Heights Using a Conversion Surface by Combining Gravimetric Quasi-Geoid/Geoid and GPS-Levelling Height Differences. Studia Geophysica et Geodaetica, 50, 165-180. https://doi.org/10.1007/s11200-006-0010-3
[8]	Tscherning, C.C., Forsberg, R. and Knudsen, P. (1992) The GRAVSOFT Package for Geoid Determination. Continental Workshop on the Geoid in Europe, Prague, May 1992, 327-334.
[9]	LaFehr, T.R. (1991) Standardization in Gravity Reduction. Geophysics, 56, 1140-1295. https://doi.org/10.1190/1.1443137
[10]	Heiskanen, W.A. and Moritz, H. (1967) Physical Geodesy. Freeman and Company, San Francisco. https://doi.org/10.1007/BF02525647
[11]	Ismail, Z. and Jamet, O. (2015) Accuracy of Unmodified Stokes’ Integration in the R-C-R Procedure for Geoid Computation. Journal of Applied Geodesy, 9, 112-122. https://doi.org/10.1515/jag-2014-0026
[12]	Ismail, Z. (2016) Détermination de l’exactitude d’un géoïde gravimétrique. Thèse de Doctorat de l’Université de recherche Paris Sciences et Lettres, Observatoire de Paris, 156 p.
[13]	Corchete, V., Chourak, M. and Khattach, D. (2005) The High-Resolution Gravimetric Geoid of Iberia: IGG2005. Geophysical Journal International, 162, 676-684. https://doi.org/10.1111/j.1365-246X.2005.02690.x
[14]	Véronneau, M. and Huang, J. (2007) The Canadian Gravimetric Geoid Model 2005 (CGG2005). Report, Division des levés géodésiques, Secteur des sciences de la Terre, Ressources naturelles Canada, Ottawa.
[15]	Corchete, V. (2010) The High-Resolution Gravimetric Geoid of Italy: ITG2009. Journal of African Earth Sciences, 58, 580-584. https://doi.org/10.1016/j.jafrearsci.2010.05.010
[16]	Abdalla, A. and Tenzer, R. (2011) The Evaluation of the New Zealand’s Geoid Model Using the KTH Method. Geodesy and Cartography, 37, 5-14. https://doi.org/10.3846/13921541.2011.558326
[17]	Corchete, V. (2013) The First High-Resolution Gravimetric Geoid for Ukraine: UGG2013. http://www.isgeoid.polimi.it/Geoid/Europe/Ukraine/Ukraine.pdf
[18]	Corchete, V. (2013) The First High-Precision Gravimetric Geoid of Hungary: HGG2013. http://www.isgeoid.polimi.it/Geoid/Europe/Hungary/Hungary.pdf
[19]	Sjöberg, L., Gidudu, A. and Ssengendo, R. (2015) The Uganda Gravimetric Geoid Model 2014 Computed by the KTH Method. Journal of Geodetic Science, 5, 35-46. https://doi.org/10.1515/jogs-2015-0007