Computing Local Geoid Model Using DTM and GPS Geodetic Points. Case Study: Mejez El Bab-Tunisia

Noamen Rebaï¹, Olfa Zenned², Hichem Trabelsi³, Hammadi Achour³

¹Département de Génie Civil, Université Tunis El Manar/Ecole Nationale d’Ingénieur de Tunis, Tunis, Tunisia.
²Schlumberger Limited, Houston, USA.
³Laboratoire des Ressources Sylvo-Pastorales, ISPT—University of Jandouba, I.S.P., Tabarka, Tunisia.

DOI: 10.4236/ijg.2018.93011                

Abstract

Different methods have been deployed to compute the geoid, the altimetry reference for surveying applications. One of their main goals is to allow the use of GPS (Global Positioning System) or GNSS heights, which are related to an ellipsoid and therefore must be corrected. Some of these methods are accurate but quite heavy as developed by [1], but one of them is easy to use while giving very good results in a local system: some mm for a 10 × 10 km² area developed by [2] [3]. In our study, we have used software called “Géoide Program”, previously used at the CERN in Switzerland and set up by [4], which they complete this software allowing a parameterization of general data to provide results in a general system. Then, tests have shown the way to optimize computations without any loss of accuracy. For our computations we use gridded of geodetic heights, from Lambert or WGS 84 datum’s, DTM (Digital Terrain Model) and leveled GPS points. To obtain these results, components of the vertical deflection are computed for every point on the grid, deduced from the attraction exerted by the mass Model. Then, geodetic heights are computed by an incremental way from an arbitrary reference. Once the calculation is performed, the geodetic height of any point located in the modelled area can be interpolated. The variations of parameters (mainly size and increments of the DTM and of the modeled area, and ground density) have shown that they do not play a significant role although DTM must be large enough to take into account an important area around a selected zone. However, the choice of the levelled GPS points is primordial. We have performed tests with real data concerning Mejez El Bab zone, in north of Tunisia. Nevertheless, for a few hundreds of square kilometers area, and just by using a DTM and a few levelled GPS points, this method provides results that look extremely promising, at least for surveying activities, as it shows a good possibility to use GPS for coarse precision levelling, and as DTM are now widely available in many countries.

Keywords

Geoid Model, “Géoide Soft”, GPS, Digital Terrain Model, Height, Vertical Deflection, Mejez El Bab-Tunisia

1. Introduction

Knowledge about the shape of the Earth and the variations of its field of gravity is now crucial for any precise work in the spatial domain of artificial satellites trajectory, machines location and astronomy of position. It also allows the elaboration of geophysical hypotheses concerning the inside of the globe [5] [6] [7] .

Since the geoid of a good approximation was a reference of the heights until the end of the 80’s, the geoid has remained an object of scientific study. The fast development of the GPS system and the necessity of converting the ellipsoidal elevations into heights have made the geoid an important stake [3] [8] .

Knowledge about the geoid is essential for the levelling because it is considered as the reference surface [9] [10] [11] . The shape of the geoid depends on the distribution of the masses inside the globe. Indeed, the fact that the Earth is not a homogeneous entity makes the geoid present undulations with regard to the ellipsoid. These undulations reflect the different aspects of heterogeneous densities [12] .

The determination of a Tunisian geoid of centimetric accuracy turns out to be more and more essential. However, the gravimetric data, necessary for its determination, are very scattered. Their homogeneity, density and quality are unknown [13] [14] .

Thus, the main objective of this work is the determination of local geoid model of Mejez El Bab using the DTM (Digital Terrain Model) of the sector and levelled GPS points. It is, in fact, about a quasi-geoid model. The purpose of this work is, also, to master the methodology of the determination of a geoïd in order to spread this calculation all over Tunisia.

2. Geographical and Geological Context of the Study Area

The choice of the study area was made according to several criteria in particular: its chaotic relief, the knowledge of its geology and a precise DTM obtained from the interpolation of the contour line of topographic map (scale 1:50,000).

The region of Mejez El Bab is situated in northern Tunisia (Figure 1). The studied sector is part of the “zone of domes” [15] [16] [17] [18] . This sector is of interest to the structures that become integrated into the front-country of the alpine chain, situated on the borders of the valley of Mejerda [19] [20] [21] . In this

(a) (b)

Figure 1. Study area: (a) DTM and hydrographic network of Mejez El Bab area [34] modified; (b) Geometry of Mourra Mountain. Partial geological Mejez el Bab map (sheet 27 of Mejez El Bab zone, scale 1:50,000, ONM).

zone the tectonic structure is characterized by Triassicoutcrops aligned according to a NE-SW direction [22] [23] [24] .

On the paleogeographic plan, the “zone of domes”, which is of about 50 km, corresponds to one bet of the qualified deep marine domain known as the “Tunisian groove”. This domain presents a complex paleogeographic and structural evolution, bound to the multiple periods of opening and closing, accompanied with Triassic halokinetic movements [25] [26] .

The studied structures in Mejez El Bab, Soughia sector correspond to a mega-anticline Plio-Quaternary resuming previous anticline and synclinal structures.

Indeed, the massive mountains of Mejez El Bab are part of the Southern sub-zone of the zone of Diapirs and constitute a mega-anticlinal which is directed NE-SW, lined in the SE by the plain of Goubellat and in the NW by the bed of Oued Mejerda [27] [28] [29] [30] [31] . We name: Jebel BouRahal, Jebel Mourra-BouMous, Jebel Kechtilou the syncline of Oued el Hamra, Jebel Jebs and Jebel Rihane.

2. Data and Methodology

2.1. Data

The Vertical reference system adopted by the Tunisian Topographical and Cadastral Office (OTC) is the orthometric height. It is a point height calculated according to the vertical line above the geoid. However, this system presents an inconvenience. The orthometric height of a point to the ground depends on values of g along the vertical line. These values are generally undefined and inaccessible in reality [7] [32] [11] . In consequence, the corrections and the orthometric heights are not physically determined. Thus, they can be calculated from the moderate raw, measuring the differences of raw level, only by means of theoretically approached formula giving the gravity g in each point of the globe. The obtained values are only approximate orthometric corrections. Resulting in the fundamental mark, the origin of the heights of the network is the mark of “Port de France”, which was built 60 years ago, and having a new height of 7 m. The 10 cm correction with regard to the former network was necessary for following the information provided by the tide gauge of la Goulette [33] . However, it is not an obvious solution because this mark is situated in a geologically unstable ground that is far away from the site of the tide gauge of la Goulette.

The Lambert conical conform projection used in the geodetic system. In order to avoid serious distortions, Tunisia has been divided into two parts: Lambert North of Tunisia (11 g, 40 g) and Lambert South of Tunisia (11 g, 37 g). The extension of the geodesic network southward requires a 3rd zone Lambert that can be characterized by (11 g, 34 g).

In fact, it is important to draw the attention on the fact that we have two types of coordinates system in Lambert projection used by the Topographic Service of Tunisia, called “STT” coordinates are as follows: X has a positive sense northward and Y has a positive sense westward.

With X, Y counted from the origin of the projection axes in which X₀ = 0 and Y₀ = 0.

The formulas allowing the shift from “IGN” coordinate system to “STT” (shift of the origin in meters, 300,000, 500,000) are:

$X_{STT}=Y_{IGN}-300000$

$Y_{STT}=500000-X_{IGN}$

The parameters of transformation of Clarke 1880 IGN (Clarke 1880 National Geographic Institute, Datum Carthage) towards GRS 80 (Datum WGS 84) are:

$DX=-266.112\text{m},DY=7.984\text{m},DZ=443.699\text{m}$

$Da=-112.2\text{m},Df=-{\mathrm{54.750714.10}}^{-}{}^6$

Among the objectives of the upgrade of the Tunisian geodesy is setting up a three-dimensional geodesic network of spatial reference, meeting the requirements of GPS levelling points. Thus, the use of the GPS observations thus requires a better determination of the height [3] and the undulations of the geoid with regard to the topographic reference surfaces of the Tunisian territory [26] .

To determine the local geoid model of Mejez El Bab, we used 12 GPS points and a Digital Terrain Model (DTM) showing the distribution of the altitude in Mejez el Bab, which is obtained from the Topographic Map (sheet n˚ 27, scale 1:50,000. OTC, 1922). The DTM is with a grid of 50 m. It presents the entire cover zone in Lambert Nord projection. The south West area of Mejez El Bab presents important reliefs which can affect the geoid model (Figure 2).

2.2. Methodology

For the calculation of the Mejez El Bab Local geoid model, we principally used two methods. The first principle method consists to the integration of the vertical deflection (ξ and η˚). The second method is based on the calculation of the potential engendered by the DTM of the area.

Our aim in this study is to elaborate “Géoide PROGRAM”. It allows us to calculate geodetic heights in points of a grid defined by the user, using DTM defined in “Lambert IGN (Clark 1880)” or in “WGS 84” of the concerned zone, and GPS Geodetic pointslevelled on this zone. These results allow then the determination of normal height at any point of the grid in the “WGS 84” system. This program uses a DTM to calculate at each points of the grid the attraction of the surrounding masses. It also allows the calculations of the vertical deflection at each point of the grid (Figure 3), defined beforehand [4] [10] [35] [36] [37] [38] [39] . Then, these results, providing only the relief of the geoid, and after

(a) (b)

Figure 2. Digital Terrain Model showing the distribution of the altitude in the area [34] , (a), and location of 12 levelled GPS points (b).

Figure 3. Calculation of the zone’s prisms attraction: [4] .

corrected by georeferencing, we obtained geoid model. This georeferencing is realized bylevelled GPS points of which the geodetic height is known by difference of the measured heights [37] [40] [41] ).

Figure xx : Flow chart of the potential’s method

$N=h-H$

where:

h indicates the ellipsoidal height obtained by GPS, expressed in WG S84 datum,

H refers to the height expressed in an altimetry system and・ N thus expresses the altitude of the height 0 of the altimetry system on the ellipsoid of geodesic reference system.

The attraction of a prism at a point of the grid is calculated by way of two formulas:

The first (F1) consists of a direct application of Newton’s law:

$F=G\cdot m_1\cdot m_2/r^2$ (F1)

$t_x=\frac{G\rho Vx}{r^3}$

where r is the distance between two bodies (m₁, m₂), V the volume of the prism supposed constant and ρ is the density.

The second which is more rigorous, is obtained by integrating the first one on the volume of the prism:

$t_x={\displaystyle \iiint \left(G\cdot \rho x/r^3\right)}\text{d}x\cdot \text{d}y\cdot \text{d}z$ (F2)

the integration gives (F3):

$t_x=-G\rho \cdot \left({\varphi }_{x2y2z2}-{\varphi }_{x2y2z1}-{\varphi }_{{}_{x2y1z2}}+{\varphi }_{{}_{x2y1z1}}-{\varphi }_{x1y2z2}+{\varphi }_{x1y2z1}+{\varphi }_{{}_{{}_{x1y1z2}}}-{\varphi }_{{}_{x1y1z1}}\right)$

where:

${\varphi }_{xyz}=y\ln \left(z+r\right)+z\ln \left(y+r\right)+x\arctan \left(\frac{xr}{yz}\right)$

The formula F1 is used until a critical distance beyond which differences of results between both formulae are unimportant, and where the formula F3 is thus applied. After several attempts this critical distance was fixed to 50,000 m because when it was considered superior, the calculation time increased strikingly without significant variation of the result.

If the zone has a more heterogeneous ground, the choice of the critical distance could have more consequences. Since the attraction is being inversely proportional to the distance, it is the closest prisms which have a dominating place in the calculation of the attraction and thus the vertical deflection. Yet, differences between the attraction’s calculation, according to the rigorous formula or not, are more important as the relief is pronounced.

2.2.1. Calculation of the Components of the Vertical Deflection x and h˚

In the studied area, the attraction t of every prism is calculated on the point.

Dg is the sum of the attractions t: $Dg={\displaystyle \sum t}$

The vertical deflection according to the X axis is given by:

$V_x=\frac{D{g}_x}{\gamma h}$

With normal gravity at the point of the grid is calculated by the formula of the IUGG 1967:

${\gamma }_0=978.03185\left(1+0.005278895{\sin }^2\phi +0.000023462{\sin }^4\phi \right)gal$

(with φ = latitude)

${\gamma }_h={\gamma }_0\frac{\left(1-2h\right)}{R}$

ξ and η are obtained from V_x and V_y by a rotation of angle α = azimuth of the Y axis.

$\left(\begin{array}{c}\xi \\ \eta \end{array}\right)=\left(\begin{array}{cc}\cos & \sin \\ -\sin & \cos \end{array}\right)=\left(\begin{array}{c}{V}_Y\\ V_X\end{array}\right)$

However, noting that this formula is approximate

Indeed, we normally have:

$\sin \theta =\frac{Dg}{g}$

With θ representing here the vertical deflection

Consequently, we can say that:

$\sin \xi =\frac{D{g}_x}{g},\sin \eta =\frac{D{g}_y}{g}$

Yet, since the acceleration of the gravity g is hardly measurable, Comet and Maillard [4] replaced it by γ (The normal acceleration).

2.2.2. Calculation of the Geodetic Height

The difference of height between two successive points of the grid is obtained through multiplying the distance between both points by the value of the vertical deflection in the direction obtained from both points. Indeed, in a point the angle between the normal and the vertical line is equal to the one between the tangent in the ellipsoid and the tangent in the geoid [4] .

In the consideration of the geoid as a line between two consecutive points of the grid the geodetic heights can be so calculated step by step. This method is applied for general leveling and not only in local leveling coordinate system, the case of “Géoïde” Soft developed by CERN (Figure 4). The modelling was made possible from a DTM with rectangular grid, in “Lambert IGN” or in WGS 84.

Once the components of the vertical deflection (ξ, η) are calculated, a wedging is made to add the undulations which are not caused by the relief. It consists in a translation and a rotation, and is made by leveled GPS points, for which the geodetic heights are known by subtraction [42] .

The principle of the wedging is to overlay the geodetic height points calculated by the model and leveled GPS points in the same geodetic system.

For this, a fundamental point will be used; it is going to be made a vertical translation and a rotation considering it as chosen in the middle of the Calculated geoid model in order to avoid losing precision too important for boundaries zone.

Indeed, an imprecision in the determination of the rotation parameters has a proportional effect on the height of the point to the distance between it and to the fundamental one. If this point is taken in the middle, the maximal distance is reduced and therefore allows minimizing errors [4] .

For wedging the geoid Model, we calculate the value of the vertical translation dN₀ and two rotation angles according to the direction of ξ and η˚: dξ₀ et dη₀. Thus, we calculate at first step as described previously, the geodetic heights of all the points of the grid and the wedging points.

During the calculating, we have chosen N₀. The latter is the approached value: It is N of wedging point that is the closest to the fundamental point. Other parameters dξ₀ and dη₀ are fixed at first to 0.

Then, the least mean square method is applied with three unknown parameters (dN₀, dξ₀, dη₀) and a number of relations of observations equal to those of the wedging points. The relations of observations are explained in [4] .

2.2.3. The Potential Calculation Method

We have based our calculation of the gravitational potential on a formula used by [42] . This method consists in the calculation of the potential of gravity. Indeed, a model made up of n prisms can result in a potential T. This allows the conclusion that a constant volume generates a constant potential.

Figure 4. Flow Chart of “GEOIDE PROGRAM”.

The contribution of this potential in the calculation of the geodetic height is based on the following formula of Browns: [43]

$N=-\frac{T}{\gamma }$

where: N is the geodetic height, T represents the Potential and $\gamma$ denotes the normal gravity.

According to this approach, the undulations of short wavelengths of the geoid ( ${\gamma }_N\le 300\text{km}$ ) are determined from a 3D model of the lithosphere. The Coordinates X and Y of the components of the vertical deflection are:

$\xi =-\frac{1}{\gamma }{T}_x\eta =-\frac{1}{\gamma }{T}_y$

The gravimetric anomaly engendered by this model is:

${\Delta }_g\equiv -T_z$

It is about calculating the potential generated by the whole DTM at each point of the grid to be modeled.

3. Case Study: Mejez El Bab Zone, North of Tunisia

The present work leads to the elaboration of a “GEOIDE Soft” presented in (Figure 4). This program enables the design of a first prototype of local geoid model in Tunisia, and precisely in Mejez El Bab.

This program allows the modelling of a local surface of height 0 at a grid by use of a DTM, and the wedging of this latter by means of levelled GPS points. It was so expressed in the Tunisian parameters. Besides, we added to this program the second method based on the calculation of the potential of gravity [43] [44] [45] .

Every method has its own specificities and scopes. For the case of Tunisia which has a system of orthometric height, we were able to make the wedging only for the method of integration of the components of the vertical deflection. Indeed, the method of the Potential is based on normal heights, but it has allowed us to validate the model of geoid obtained and to make a comparison between two different methods (Figure 5). In both cases, we obtain a geoid of a centimetric precision. This one reproduces globally the relief of Mejez El Bab.

In all what follows, the used data are plotted using the curve line functions of “Surfer Software” to modelling the relief.

What is represented in the (Figure 6) is the entire zone covered by the DTM in “Lambert North IGN”. What is worth noting is that Southwest half is made up of a rather strong relief that has an influence on the geoid.

(a) (b)

Figure 5. Comparison of the geoid model (color scale is in meters), obtained from 1-Integration of the vertical deflection of ξ and η, 2-Potential method.

Figure 6. Undulations of Mejez El Bab Geoid Model with respect to the DTM.

In order to prop up the model, we used 12 levelled GPS points. They were supplied to us by the Cadastral and Topographical Tunisian Office (named OTC), with the following information: codes, heights, longitudes and ellipsoidal heights. The Cartesian coordinates, (X,Y,Z) are presented in the systems CLARKE 1880 and WGS 84 [32] .

4. Results and Discussion

The obtained model represents a geoid that reproduces globally the relief, as shown in the (Figure 7). Indeed, we discern in its South West part a bump of direction NE-SW where the geodetic height reaches 0.03 meters, bounded by two hollows.

This bump would correspond to Jebel Bou Mous, Jebel Mourra, Jebel Jebs and Jebel Kechtilou (Figure 1), whereas the two hollows on both sides would represent the plain of Mejerda to the North (between 0.05 m and 0.03 m) and that of Goubellat in the South (between −0.03 m and 0.01 m). Also, in the southeast corner, we notice a light rise of the geoid reaching 0.02 m corresponding to Jebel Khamfir and Jebel Ed Deridjah. Whereas in the extreme North-West, the geoid rises again but achieves only −0.01 m in 0 m. On the map this rise corresponds to Jebel Jedidi and Jebel Chaabane.

The model obtained by the potential method (Figure 8) is a geoid which slightly reproduces the relief but by being too much inflated. Hence, we realized

Figure 7. Calculated Geoid by the Potential method.

Figure 8. View of Geoid Model obtained by the potential method with a grid space of: 1600 m.

that is necessary to subtract the average value from the potential generated by the whole DTM in the middle point. We obtain then the geoid of the Figure 5 with a minimum of 0.07 m and a maximum of 0.06, corresponding at 13 cm of undulations.

The tests cover the influence of the various parameters (size and stitch of the DTM and the zone). In every test, a parameter is adapted to a standard, and the differences on the obtained heights are then studied.

The objective of these experiences is to compare the various models with what is in priori the best. The comparison is a quadratic gap average, and must be examined with regard to the scattering, which is the gap quadratic average of the geoid in its average.

For all the parameters, we note that their choice is bound, for the majority of them, to be bound to the relief of the considered region. A precise analysis of the reliefs surrounding the zone must thus be made. The size of the DTM is fundamental as far as a bigger DTM is needed, but it is not enough. It is useful that the limits of the DTM are fixed by taking into account the relief with a maximum integration of the stressed parts [4] . The stitch of the DTM does not seem, as far as it is concerned, to be the most sensitive parameter. However, the realized tests converge on the choice of 200 m (Table 1; Table 2).

For the stitch of the zone to be modelled, better results are obtained with a reduced stitch ξ and η. Seen the limit of the calculation, the choice of 500 m seems suitable.

Regarding the points of wedging, their choice is determining for the final precision (Table 3). Seen their limited number, we are not able to test their influence on the modelling of the geoid. However, the more this number decreases, the more the position in the zone and the quality of points become more important during the wedging of the modeled surface. A compromise between the minimization of the number of points and the precision of the result must thus be reached.

(a) (b)

Table 1. Influence MNT’s size on the modelling of the geoid. 1) Grid of MNT used: 200*200 m. 2) Size of the modelling area: 5*5 km. 3) Grid of the modelling area: 500*500 m. (a) (Dispersion of total MNT = 1.417 × 10⁻² m); (b) (Dispersion of total MNT = 2.58 × 10⁻² m).

(a) (b)

Table 2. Influence of Grid’s size of MNT on the modelling geoid. 1) Grid of MNT used: total MNT = 33*22 km. 2) Size of the modelling area: 5* 5 km. 3) Grid of the modelling area: 500*500 m. (a) (Dispersion of the Grid 200 = 1.417 × 10⁻² m); (b) (Dispersion of the Gird 200 = 2.58 × 10⁻² m).

(a) (b)

Table 3. Influence areas size on the modelling of geoid. 1) Size of MNT used: 33*22 km. 2) Grid of MNT used: 200*200 m. 3) Grid of the modelling area: 500*500 m. (a) (Dispersion of the area 5*5 km = 1.417 × 10⁻² m); (b) (Dispersion of the area 5*5 km = 2.58 × 10⁻² m).

5. Conclusions

The present work led to the elaboration of a program allowing the design of a first local of geoid model of the region of Mejez El Bab which is located in the North of Tunisia.

A prospective improvement of the model would be possible by certain techniques. First, varying the density on the DTM while considering a chaotic terrain has important implications. Second, having precise DTM is helpful. Third, the size of the DTM plays a crucial role. Therefore, it could be interesting to increase its size without too much increasing the mass of the calculations by using a much bigger stitch in periphery. The “Géoïde Soft”, being for the moment limited in capacity, can’t handle a DTM of more than 178*178 knots. These limits could be exceeded by finding the way to develop more the programs without amplifying the calculation. Also, it would be necessary to arrange more leveled GPS points with the maximum of accuracy, which requires a precise leveling.

Through this example, by using the data in the best way, it appears that rms on the obtained geodetic heights is better than 3 cm, which is close to GPS altimetry precision in z. Therefore, the total rms might be perhaps better, which could be confirmed with more accurate GPS measurements.

Nevertheless, for a few hundreds of square kilometers area, and just by using a DTM and a few levelled GPS points, this method provides results that look extremely promising, at least for surveying activities, as it shows a good possibility to use GPS for coarse precision levelling, and as DTM is now widely available in many countries.

The Terrestrial geoids are so difficult to be interpreted as they represent the effect resulting from the heterogeneity existing in the various internal layers of the Earth having complex structure and dynamics, which interact with each other through mechanical or thermal couplings, etc. However, by setting more data and ways, an alliance geoid-gravimetric would be interesting [41] . It could help us to emit various hypotheses and to interpret the heterogeneity of the interior of the earth.

Indeed, we consider in this study that the attraction value or the integration of vertical deflection components is directly proportional to the density. Therefore, an increase in density accentuates the relief. In our case, we have considered that the density is constant and equal to the continental crust (2.67).

In order to have a better definition of the geological structures, it is necessary to refine the gravimetric models according to the Regional scale. These models will contribute, in a significant way, to the realization of a Tunisian geoid of high accuracy. Indeed, it would be necessary to integrate various methods to reach the required precision (gravimetric collocation) [44] [45] .

Today, the essential results are known by the global models of the geopotential. Such models integrate all the information collected from artificial satellites, measures of gravity and, for several years, even altimetry spatial measures.

Acknowledgements

Authors are thanking fully the Laboratorymembers of geodetic research in IGN (France) in particular Olivier JAMET and Michel KASSER for accepting to offer us the necessary data and to help us to run the “Géoïde Soft”. Also, we appreciate the collaboration of the Tunisian Office of Cadastral and Topography for providing us with GPS measurements.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Torgue, W. (1991) Geodesy. Second Edition, Walter de Gruyter, Berlin. https://doi.org/10.1515/9783111542683
[2]	BELL, B. (1984) A Simulation of the Gravity Field around LEP. CERN, Geneva, 100.
[3]	Gervaise, J., Mayoud, M., Beutler, G. and Gurtner, W. (1985) Test of Global Positioning System (NAVSTAR SATELLITES) on the CERN-LEP Control Network. Joint Meeting FIG-Study Groups 5B and 5C on Inertial, Doppler and GPS Measurements for National and Engineering Surveys, Munich, 1-3 July 1985.
[4]	Comet, G. and Maillard, D. (1994) Détermination d'un géoide local par utilisation d'un MNT, Mémoire de Diplome d'Ingénieur. Ecole Supérieur des Géomètres et Topographes-France, 82.
[5]	Lacombe, H. and Costabel, P. (1990) La Figure de la Terre du XVIIIème siècle à l'ère spatial. Académie des Sciences.
[6]	Levallois, J.-J. (1970) Géodésie générale, tome III: Le champ de pesanteur-Eyrolles.
[7]	Heiskanen, W.A. and Moritz, H. (1967) Physical Geodesy. Institute of Physical Geodesy, Technical University, Graz.
[8]	DMA—Department of Defense Administration (1987) Geodetic Systeme 1984 Works. Its Definition and Relationships with Local Geodetic System—DMA TR B350.2, USA.
[9]	Gurtner, W. (1978) Das Geoid in der Schweizerland Institute für Geodasie und Photogrammetrie. Mitteilungen Nr 20, Zurich.
[10]	Delomenie, M. (1987) Calcul d'un géoide gravimétrique sur la France. BGI, IGN.
[11]	Duquenne, H., Zhiheng, J. and Sechaud, N. (1994) Nivellement par GPS, vers un géoide local centimétrique. Géomètre, 42-47.
[12]	Iang, Z. and Duquenne, H. (1995) Fast Integration for the Integrals of Stokes, Potential and Terrain Correction in Geoid Determination. 20th General Assembly of Europe Geophysical Society, Hamburg, April 1995, 23-25.
[13]	Jallouli, C. and Mickus, K. (2000) Regional Gravity Analysis of the Crustal Structure of Tunisia. Journal of African Earth Sciences, 30, 63-78. https://doi.org/10.1016/S0899-5362(00)00008-7
[14]	Jallouli, C., Rebai, N., Turki, M.M. and Mickus, K. (2000) L’importance de la Gravimétrie en Géodésie; Géo-Top. Revue scientifique et technique; Office de la Topographie et de la Cartographie, 72.
[15]	Bolze, J. (1954) Ascension et percée des diapirs au Crétacé moyen dans les monts de Téboursouk. C.R. Somm. Soc. Géol. France, II, 303, 831.
[16]	Perthuisot, V. (1978) Dynamique et pétrogenèse des extrusions triasiques en Tunisie Septentrionale. Thèse de doctorat es-Sciences, Trav. Labo. Géol. (Ecole Normale Supérieure); Paris, 360 p.
[17]	Vila, J.M. and Charriere, A. (1993) Découverte d’Albien calcaire et de Trias resédimenté au Jebel Bou Jaber Ac. Sci. Paris, II, 316.
[18]	Hammami, M. (1999) Tectonique, Halocinèse et mise en place des minéralisations dans la zone des diapirs. Thèse de doctorat, Univ. Tunis II, FST, 183 p.
[19]	Zargouni, F. (1975) Etude géologique de la chaine de Lansarine (région Tébourba, Atlas tunisien). Thèse 3ème cycle, Univ. Paris and M. CIRIE., 86 p.
[20]	Ben Ayed, N. (1994) Les décrochements-chevauchements E-W et N-S convergents de la Tunisie septentrionale: Géométrie et essai de reconstitution des conditions de déformations. Actes des IVème Journées de l’Exploration pétrolière en Tunisie, 25-37.
[21]	El Ouardi, H., Rami, A. and Turki, M.M. (1999) Géodynamique des bassins crétacés d’une partie de l’Atlas tunisien septentrional. Géologie Méditerranéenne, 26, 203-216.
[22]	Chikhaoui, M. and Turki, M.M. (1995) Roles et importance de la fracturation méridienne dans les déformations crétacées et alpines de la zone des diapirs. African Earth Sciences, 21, 271-281. https://doi.org/10.1016/0899-5362(95)00066-3
[23]	Perthuisot, V. (1981) Diapirism in Northern Tunisia. Journal of Structural Geology, 3, 231-235. https://doi.org/10.1016/0191-8141(81)90019-5
[24]	Aoudjehane, M., Bouzenoune, A., Rouvier, H. and Thieberaz, J. (1992) Halocinèse et dispositifs d’extrusion du Trias de l’Atlas Saharien Oriental. Géologie Méditerranéenne, 4, 273-280.
[25]	Castany, G. (1956) Etude géologique de l’Atlas tunisien oriental. Thèse de doctorat Sciences, Ann. Mines et Géol. Tunisie; 8, 632p.
[26]	Chikhaoui, M. (1988) Succession Distension-Compression dans le sillon tunisien. Secteur de Nabeul-El Kef, Tunisie centre-nord. Role des extrusions triasiques précoces lord des serrages alpins. Thèse de doctorat, Univ. De Nice, 143 p.
[27]	Adil, S. (1993) Dynamique de Trias dans le nord de la Tunisie: Bassins en relais multiples de décrochement, magmatisme et implications minières. Thèse de doctorat, Univ. De Tunis II, FST, 284 p.
[28]	El Ouardi, H. and Turki, M.M. (1995) Tectonique salifère polyphasée dans la région de Mejez El Bab-Testour (Zone des domes, Tunisie septentrionale): Controle de la sédimentation méso-cénozoique. Géologie Méditerranéenne, 22, 73-84.
[29]	El Ouardi, H. (1996) Halocinèse et role des décrochements dans l’évolution géodynamique de la partie médiane de la zone des domes. Thèse de doctorat, Univ. Tunis II, FST, 242 p.
[30]	Chikhaoui, M., Hammami, M., Rabhi, M. and Turki, M.M. (1999) Tectonique synsédimentaire distensive dans le Jurassique du sillon tunisien. Exemple du Jurassique de Mjez El Bab. Notes du Service Géol., Tunisie, 43-48.
[31]	Chikhaoui, M. (2002) La zone des diapirs en Tunisie; cadre structural et évolution géodynamique de la sédimentation méso-cénosoique et géométrie des corps triasiques. Thèse de doctorat, Université Tunis El Manar, FST, 323.
[32]	Boucher, C. (1981) Définition des systèmes géodésiques utilisés en France-IGN, rapport technique SGNM/DEC, RT/G No. 7.
[33]	Zid, J. (1990) Le réseau général de nivellement de la Tunisie; Office de la Topographie et du Cadastre, direction de la Géodésie et des Travaux spéciaux, Tunis-Tunisie.
[34]	Slama, T. (2002) DEA Etude méthodologique pour la détermination d’un modèle numérique de terrain. Application à la région de Mejez El Bab, FST.
[35]	Achilli, V. and Baldi, P. (1982) Computations of Local Anomalies of the Vertical Deflections in Geodetic Network. Survey Review, 26, 327-335.
[36]	Tscherning, C.C. (1985) Geoid Modelling using Collocation in Scandinavia and Greenland. Manuscripta Geodetica, 10, 136-149.
[37]	Sechaud, N. (1993) Détermination d’un géoide local par géodésie intégrée-Diplome d’ingénieur; IGN, Décembre 1993.
[38]	Duquenne, H. and Sarrailh, M. (1996) Improvement of Gravemetric Geoid Determination in the French Alps. Finnish Geodetic Institute Report, No. 96, 71-76.
[39]	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
[40]	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
[41]	Nagy, D., Papp, G. and Benedek, J. (2000) The Gravitational Potential and Its Derivatives for the Prism. Journal of Geodesy, 74, 552-560. https://doi.org/10.1007/s001900000116
[42]	Jiang, T. and Wang, Y.M. (2016) On the Spectral Combination of Satellite Gravity Model, Terrestrial and Airborne Gravity Data for Local Gravimetric Geoid Computation. Journal of Geodesy, 90, 1405-1418. https://doi.org/10.1007/s00190-016-0932-7
[43]	Mishra, U.N. and Ghosh, J.K. (2016) Development of a Gravimetric Geoid Model and a Comparative Study. Journal Geodesy and Cartography, 42, 75-84.
[44]	Nagy, D. and Fury, R.J. (1990) Local Geoid Computation from Gravity using the Fast Fourier Transform Technique. Bulletin Géodésique, 64, 283-294. https://doi.org/10.1007/BF02519181
[45]	Arana, D., Camargo, P.O., Guimaraes and Gabriel, N. (2017) Hybrid Geoid Model: Theory and Application in Brazil. Anais da Academia Brasileira de Ciências, 89, 1943-1959. https://doi.org/10.1590/0001-3765201720160802