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

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.


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.

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 2 : 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.
• Storage and pre-processing of measurements.
• Multi-parameter measurement sensor: gravity + internal temperature and tilt of the sensor, GPS (option).
• 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).

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).

Methodology
The measurement of the gravity field on a global scale has become possible due Journal of Geographic Information System  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: The signal is considered to be composed by three wavelength ranges: 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: In the second step, the Stokes integral is applied to the values of interpolated At the last stage, the altitude anomalies can be written as: "Remove" Calculation The first step in the R-C-R method consists in calculating the residual anomalies ∆ gRes at the points of the measurements: 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: 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 3 . The Bouguer plateau exerts an attraction (2πGρh). All the masses above the geoid are removed (the simple Bouguer reduction) [9]. Simple Bougu-Journal of Geographic Information System er anomalies are calculated by the Equation (11) [10]: To calculate the terrain correction we perform the direct summation of the contributions of each model element (prism, tesseroides etc.).

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 0 (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 M. Ghannem, A. Belhaj Ali 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. 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   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).

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. To assess the accuracy of the calculated quasi-geoid, a comparison with other Journal of Geographic Information System 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.  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.   (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 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.

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.