Determination of a Gravimetric Geoid Solution for Andalusia ( South Spain )

The orthometric heights can be obtained without levelling by means of the ellipsoidal and geoidal heights. For engineering purposes, these orthometric heights must be determined with high accuracy. For this reason, the determination of a high-resolution geoid is necessary. In Andalusia (South Spain) a new geopotential model (EIGEN-GL04C) has been available since the publication of a more recent regional geoid. As a consequence, these new data bring about improvements that ought to be included in a new regional geoid of Andalusia. With this aim in mind, a new gravimetric geoid determination has been carried out, in which these new data have been included. Thus, a new geoid is provided as a data grid distributed for the South Spain area from 36 to 39 degrees of latitude and –7 to –1 degrees of longitude (extending to 3 × 6 degrees), in a 120 × 240 regular grid with a mesh size of 1.5’ × 1.5’ and 28800 points in the GRS80 reference system. This calculated geoid and previous geoids are compared to the geoid undulations obtained for 262 GPS/levelling points, distributed within the study area. The new geoid shows an improvement in accuracy and reliability, fitting the geoidal heights determined for these GPS-levelling points better than any previous geoid.


Introduction
In order to convert the high accuracy ellipsoidal heights (usually obtained from the modern GPS technology) to orthometric heights, as required for engineering purposes, the determination of the geoid is necessary, as the relation between ellipsoidal and orthometric heights involves the geoidal height.For this reason, the determination of the geoid has been one of the main objectives of the geodesists in the last few decades.
In the southern part of Iberia, a previous study had as main objective the geoid computation [1].In this study, the IGG2005 geoid was calculated.This geoid was obtained as a data grid in the GRS80 reference system, distributed for the Iberian area from 35 to 44 degrees of latitude and -10 to 4 degrees of longitude and provided as a 361 × 561 regular grid with a mesh size of 1.5' × 1.5'.Nevertheless, the negative influence of a bathymetry with much less resolution than the topography and the scarcity of gravity data in some areas of Iberia (like Portugal), resulted in a geoid of poor accuracy (27.8 cm of standard deviation).On the other hand, the computation of IGG2005 is based on the EIGEN-CG01C geopotential model, and as a result of the publishing of IGG2005, a new geopotential model (EIGEN-GL04C) has been available since.Without any doubt, this new model (EIGEN-GL04C) means an improvement that ought to be included in any new geoid to be computed from now on.
Thus, it would be very desirable to obtain a more accurate geoid solution for the southern part of Spain based on this new model, considering study windows fitted to the areas with higher density of gravity data, to reduce as much as possible the computation errors associated to the scarcity of gravity data.Furthermore, the negative influence of a bathymetry with a poorer resolution than topography, would be reduced considering the smallest marine area possible.Thus, the new geoid would give a complete picture of the geod for the South Spain area, with more accuracy than the previous geoid IGG2005.
Thereby, a new geoid should be computed as a 120 × 240 regular data grid in the GRS80 reference system, with a mesh size of 1.5' × 1.5', extended from 36 to 39 degrees of latitude and -7 to -1 degrees of longitude.This new geoid would be computed using the Stokes integral in convolution form.The necessary terrain correction would be applied to obtain the gridded reduced gravity anomalies.The corresponding indirect effect would be taken into account.After the computation of this Andalusian Regional Geoid of Spain (ARGOS), it would be compared to the IGG2005 geoid, to demonstrate the improvement in accuracy and reliability achieved by the new geoid.

Data Preparation and Pre-Processing
A complete data set is necessary for the gravimetric geoid computation.This data set consists of: 1) free-air gravity anomalies; 2) a geopotential model for the computing the long-wave length contribution to the geoid and the gravity anomaly; 3) a high-resolution DTM for the computation of the terrain correction and the indirect effect; and 4) GPS/Levelling geoid undulations to test the geoid model obtained.The data set used for the computation of the Andalusian geoid has been described below.We also described the preparation and pre-processing required for each kind of data used in this study.
Land and marine gravity data bank.The land and marine gravity data used in this study has been provided in two CD-ROM by the National Geophysical Data Center (NGDC).This data set consists of 18621 point freeair gravity anomalies (10180 in land and 8441 at sea), distributed in the study area (36 to 39 degrees on latitude and -7 to -1 degrees on longitude) as shown in Figure 1.
The accuracy of all these data ranges from 0.1 to 0.2 mgal.The compiling gravity data have been checked to remove repeated points.All the data have been converted in the GRS80 reference system and the atmospheric correction has been taken into account [2,3].

Geopotential model.
The EIGEN-GL04C model [4] is an upgrade of the EIGEN-CG03C model [5], which is a combination of the GRACE (Gravity Recovery and Climate Experiment) and LAGEOS (LAser GEOdynamics Satellite) mission solution adding a 0.5 × 0.5 degrees gravimetry and altimetry surface data.The surface data are identical to EIGEN-CG03C set except of the geoid undulations over the oceans.The EIGEN-GL04C geopotential model means a major advance in the modelling of the Earth's gravity and geoid.Thereby, this model is the geopotential model that ought to be used for the computation of the long-wave contribution to the geoid and the gravity anomaly, to obtain a high-precision geoid in the study area.
Digital terrain model (DTM).Any gravimetric geoid computation based on the Stokes's integral must use anomalies that have been reduced to the geoid, usually by means of the Helmert's second method of condensation [6].This involves the computation of the terrain correction and the indirect effect on the geoid, which are computed from a DTM, which is also necessary to compute the Residual Terrain Model reduction (RTM reduction) for the point anomalies in order to obtain smooth gravity anomalies, which are more easily gridded.In this study a new elevation model for the whole study area, with a 3'' × 3'' spacing, has been obtained from the Shuttle Radar Topography Mission (SRTM) elevation data and the ETOPO2 bathymetry data, following the process described by Corchete et al. [1].In order to minimize the loss of accuracy associated to the low resolution of the ETOPO2 bathymetry, the data window was selected so as to include the minor marine data as possible.
GPS/Levelling control data.The Global Positioning System (GPS) yields the ellipsoidal height of a terrestrial point, therefore, subtracting the orthometric height (obtained by means of high-precision levelling), the geoid undulations be obtained.They are called GPS/Levelling geoid undulations.This procedure provides a discrete geometrical estimation of the geoid height, in the study region, which can be compared with the geoid height predicted by the gravimetric model.Thus, a reliable validation of the obtained geoid gravimetric model is possible, as the GPS and levelling techniques are highprecision techniques.Therefore, the geoid undulation obtained by means of these joint techniques, jointly, will lead to a very precise determination.In Spain, there exists a high precision geodetic network that covers all the Spanish territory.The Instituto Geográfico Nacional (Madrid, Spain) developed the REGENTE cover all the Spanish territory with a high precision three-dimensional geodetic network.The number of stations covering all the national territory is of 1200 (including the Canary and Balearic Islands).As a result of this project there is a high precision geodetic network in Spain which allows to perform geodetic studies.Figure 2 shows the distribution of the GPS/Levelling points, over the study area.

Gravity Data Gridding
As the gravity data set consists of point data anomalies distributed randomly, an interpolation process should be applied to obtain a regular data grid [1].Before the interpolation, it is highly necessary to remove the shortwavelength and the long-wavelength effects applying the well-known relationship (the RTM correction) where the superscript pts denotes each point randomly distributed over the study area, is the free-air gravity anomaly, k is the Newton's gravitational constant, obtained by applying of a 2D low-pass filter with a resolution of 60', to the elevations field), c is the terrain correction computed for land and marine points, and is the gravity anomaly computed from the geopotential model EIGEN-GL04C.
After the smoothing procedure given by (1), a regular grid has been calculated by using Kriging-based routines, which are a part of OriginLab software package (© 1991-2003 OriginLab Corporation).The gridded data are distributed over the study area from 36 to 39 degrees of latitude and -7 to -1 degrees of longitude (extending to 3 × 6 degrees), in a 120 × 240 regular grid with a mesh size of 1.5' × 1.5'.Finally, RTM must be restored in the gridded anomalies to obtain the true free-air anomalies relative to EIGEN-GL04C.This RTM effect can be restored by where the superscript grid denotes each point of the regular grid considered (241 × 321 = 77361 points), is the free-air gravity anomaly, is the gravity anomaly reduced by (1) and gridded.It should be noted that and c are computed in the same way as described above, but this time over a regular grid of points.

Geoid Computation and Validation
The computation method adopted for the calculation of a gravimetric geoid was the method described in detail by Corchete et al. [1].So only a brief description of the proposed methodology will be included in this paper.
Geoid computation.The present geoid has been computed applying the classical remove-restore technique.Following this method, the geoid model is obtained by the sum of three terms where N 1 is the long-wavelength contribution, N 2 is the indirect effect (obtained from the elevations model above described) and N 3 is the short-wavelength contribution obtained from the gravity anomalies over the study area reduced to the geoid (following the Helmert's second method of condensation).The long-wavelength contribution N 1 (shown in where (J nm , K nm ) are the fully normalized geopotential coefficients of the anomalous potential, a is the semimajor axis of the reference ellipsoid, kM is the Earth's geocentric gravitational constant, P nm are the fully normalized Legendre functions, (r, , ) are the geocentric position of the points over the geoid surface,  is the normal gravity over the reference ellipsoid, n max is the maximum degree of the considered geopotential model (EIGEN-GL04C model) and N 0 is the zero degree term.
The second term of the Equation ( 3) is the indirect effect of Helmert's second method of condensation on the geoid.This term N 2 consists of two terms in the planar approximation where k is the Newton's gravitational constant,  is the density of the topography (assumed constant with a value of 2.67 g/cm 3 ), (x p , y p ) are the coordinates of the computation point, (x, y) are the coordinates of the integration points,  is the normal gravity over the reference ellipsoid, A is the study area in which the data area provided and h(x,y) are the elevations model (considering only positive elevations, i.e. masses above the geoid).The planar approximation expressed by the Equation ( 5) can be easily written in convolution form and computed by the FFT [1].The results of this computation are shown in Figure 3(b).
Finally, the third term N 3 in Equation ( 3) is the contribution of the residual gravity obtained by means of where  is the normal gravity over the reference ellipsoid, R is the mean Earth radius, S() is the Stokes function,  is the unit sphere of integration and g are the fully reduced anomalies in agreement to Helmert's second method of condensation, obtained by When we have computed all terms in the Formula (7), we can proceed to integrate g by means of Equation ( 6).
This integration can be expressed in convolution form by using of 1D FFT (Haagmans et al., 1993).The results of th d of 28800 points in the GRS80 rence system. is computation are shown in Figure 3(c).Thus, the gravimetric geoid solution (ARGOS geoid) is reached by the sum of all previously computed terms according to Equation (3).This geoid with a mesh size of 1.5' × 1.5' (extended 3 × 6 degrees over the study area) is shown in Figure 4.This new model is provided as a data grid distributed for the South Spain area from 36 to 39 degrees of latitude and -7 to -1 degrees of longitude, in a 120 × 240 regular gri refe with a mesh size of 1.5' × 1.5'.The new geoid shows a minor discrepancy with the observed geoid undulations obtained for 262 GPS/Levelling points with respect to the other available geoid.Thus, a new high-resolution geoid has been achieved for Andalusia (South Spain).Nevertheless, this geoid is only a first step towards the obtaining of a geoid with centimetric precision.The centimetre accuracy should be reached when new gravity data (in some zones with scarce gravity data coverage) and a high-resolution bathymetry are available for South Spain and the Gibraltar Strait area.New gravity data are needed as some of the existing gravity data are very old (a lot of these points were measured a long time ago, from 1950 onwards).For this reason, the computation of a gravimetric geoid with a centimetre precision is not possible with the present gravity data.This centimetre precision in the geoid model will be obtained, when new gravity data are available to replace the older measurements of the gravity data compilations for Iberia.An updating of these international compilations is greatly needed to supply new gravity data measured with the latest technology, to replace the older measurements, which obviously are not as accurate as recent measurements made with modern precise gravimeters.In spite of this, ARGOS has better accuracy than any previous geoid solutions obtained in the South Spain.This new model will be useful for orthometric height determination by GPS as it allows the orthometric height determination in the mountains and remote areas, where levelling has many logistic problems.Geoid validation.To check our geoid, we have proceeded to compare the geoid undulations predicted by the other previous geoid existing for this area (IGG2005), with the undulations predicted by our geoid.The newly calculated geoid (ARGOS) and the previous geoid (IGG 2005) are compared to the geoid undulations obtained for the 262 GPS/levelling points, used as the control data set.The statistics of these differences are shown in Table 1.The new geoid (ARGOS) shows an improvement in accuracy and reliability, fitting the geoidal heights determined for the GPS-levelling points, better than the previous geoid (IGG2005).

Conclusions
The new data available for the study area and the computation methods based on the FFT analysis, have allowed the calculation of a new gravimetric geoid for the Andalusian region, which is a major advance in the modelling of the geoid for this region.The good quality of the gravity data considered joint to a high-resolution DTM and the consideration of EIGEN-GL04C, as the reference global model, have provided the high-quality data set needed to compute a new high-resolution geoid.By using this high-quality data set, we have carried out the gravimetric geoid determination by means of the Stokes integral in convolution form.This method is shown as an efficient method to compute a high-resolution geoid, obtaining a regular gridded geoid of 120 × 240 points distributed for the study area (South Spain), from 36 to 39 degrees of latitude and -7 to -1 degrees of longitude,

Acknowledgements
The authors are grateful to the National Geophysical Data Center (NGDC) who has provided the gravity data used in this study.To David Dater, Dan Metzger and Allen Hittelman (U.S.Department of Commerce, National Oceanic and Atmospheric Administration) who have compiled this gravity data set.To NGDC and the United States Geological Survey (USGS) who have supplied the elevation data required to compute the necessary terrain corrections, through the databases: ETOPO2 and SRTM 90M (available by FTP internet protocol).The authors are also grateful to the Instituto Geográfico Nacional (Madrid, Spain) who has provided GPS/Levelling data used for validation of the geoid obtained.

Figure 1 .
Figure 1.Geographical distribution of the observed data over the study area (18621 free air gravity anomalies: 10180 in land nd 8441 at sea, from the National Geophysical Data Center).a

Figure 2 .
Figure 2. Geographical distribution of the GPS/Levelling points over the study area (262 points).
anomalies computed by Equation (2), c is the terrain correction (only considering masses above the geoid) and g is the indirect effect on gravity.The terrain correction c that appears in the Equation (7) was computed by a means of a linear approximation.The indirect effect on gravity g was computed from the indirect effect N 2 by means of free g  g = 0.3086 N 2 (in mgal for N 2 in meters)

Figure 3 .
Figure 3. (a) The EIGEN-GL04C geoid model computed for the study area.The contour interval is 1.0 m; (b) The indirect effect on the geoid (plotted positive).The grey scale is non-linear; (c) The residual geoid undulation.

Figure 4 .
Figure 4.The geoid obtained for Andalusia (South Spain) called ARGOS, as a sum of the terms given by the Equation (3).The contour interval is 1.0 m.