Engineering, 2010, 2, 160-165
doi:10.4236/eng.2010.23022 lished Online March 2010 (http://www.SciRP.org/journal/eng/)
Pub
Determination of a Gravimetric Geoid Solution for
Andalusia (South Spain)
Francisco Manzano1, Victor Corchete1, Mimoun Chourak2,3, Gil Manzano1
1Higher Polytechnic School, University of Almeria, Almeria, Spain
2Faculté Polidisciplinaire d'Errachidia, University of Moulay Ismaïl, B.P. Boutalamine, Morocco
3NASG (North Africa Seismological Group), Gelderland, The Netherlands
Received November 24, 2009; revised December 20, 2009; accepted December 22, 2009
Abstract
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 con-
sequence, 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/lev-
elling points, distributed within the study area. The new geoid shows an improvement in accuracy and reli-
ability, fitting the geoidal heights determined for these GPS-levelling points better than any previous geoid.
Keywords: Gravity, Geoid, FFT, GPS/Levelling, Andalusia, South Spain
1. 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 rela-
tion 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 ge-
odesists 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 ob-
tained as a data grid in the GRS80 reference system, dis-
tributed 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 bathym-
etry 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 computa-
tion of IGG2005 is based on the EIGEN-CG01C geopo-
tential 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 ac-
curate 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 influ-
ence of a bathymetry with a poorer resolution than to-
pography, 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

F. Manzano ET AL. 161
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 corre-
ction would be applied to obtain the gridded reduced gra-
vity anomalies. The corresponding indirect effect would
be taken into account. After the computation of this An-
dalusian Regional Geoid of Spain (ARGOS), it would be
compared to the IGG2005 geoid, to demonstrate the im-
provement in accuracy and reliability achieved by the
new geoid.
2. Data Preparation and Pre-Processing
A complete data set is necessary for the gravimetric ge-
oid computation. This data set consists of: 1) free-air
gravity anomalies; 2) a geopotential model for the com-
puting 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-proc-
essing 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 Cen-
ter (NGDC). This data set consists of 18621 point free-
air 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 cor-
rection 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 GEOdynam-
ics 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 geo-
potential 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 condensa-
tion [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 com-
pute the Residual Terrain Model reduction (RTM reduc-
tion) 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 mini-
mize the loss of accuracy associated to the low resolution
of the ETOPO2 bathymetry, the data window was se-
lected 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 (ob-
tained by means of high-precision levelling), the geoid
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
C opyright © 2010 SciRes. ENG

F. Manzano ET AL. 162
undulations can be obtained. They are called GPS/Lev-
elling geoid undulations. This procedure provides a dis-
crete 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 reli-
able validation of the obtained geoid gravimetric model
is possible, as the GPS and levelling techniques are high-
precision techniques. Therefore, the geoid undulation
obtained by means of these joint techniques, jointly, will
lead to a very precise determination. In Spain, there ex-
ists 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 distribu-
tion of the GPS/Levelling points, over the study area.
3. 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 in-
terpolation, it is highly necessary to remove the short-
wavelength and the long-wavelength effects applying the
well-known relationship (the RTM correction)
pts
GM
ptspts
ref
pts
free
pts
red gchhkgg  )(2

(1)
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,
free
g

is the density of the topography (2.67 g/cm3) for the
RTM correction on land or the density of the topography
minus seawater density (2.67 – 1.03 = 1.64 g/cm3) for
marine RTM, h is the elevation, denotes the eleva-
tion of the reference surface (this reference surface is
ref
h
Figure 2. Geographical distribution of the GPS/Levelling
points over the study area (262 points).
obtained by applying of a 2D low-pass filter with a reso-
lution of 60’, to the elevations field), c is the terrain cor-
rection computed for land and marine points, and GM
g

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 dis-
tributed over the study area from 36 to 39 degrees of lati-
tude 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 grid-
ded anomalies to obtain the true free-air anomalies rela-
tive to EIGEN-GL04C. This RTM effect can be restored
by
gridgrid
ref
grid
red
grid
free chhkgg )(2

(2)
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 regu-
lar grid of points.
free
gred
g
)(2 ref hhk  
4. 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 pro-
posed methodology will be included in this paper.
Geoid computation. The present geoid has been com-
puted applying the classical remove-restore technique.
Following this method, the geoid model is obtained by
the sum of three terms
N = N1 + N2 + N3 (3)
where N1 is the long-wavelength contribution, N2 is the
indirect effect (obtained from the elevations model above
described) and N3 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 N1 (shown in Fig-
ure 3(a)) can be computed from a spherical harmonic
expansion by

0
0
2
1
sincos) (sin
max
NmKmJP
r
a
r
kM
N
n
m
nmnmnm
n
n
n














(4)
where (  Jnm,  Knm) are the fully normalized geopotential

F. Manzano ET AL. 163
coefficients of the anomalous potential, a is the semi-
major axis of the reference ellipsoid, kM is the Earth’s
geocentric gravitational constant, Pnm are the fully nor-
malized Legendre functions, (r,  ,  ) are the geocentric
position of the points over the geoid surface,  is the
normal gravity over the reference ellipsoid, nmax is the
maximum degree of the considered geopotential model
(EIGEN-GL04C model) and N0 is the zero degree term.
The second term of the Equation (3) is the indirect ef-
fect of Helmert’s second method of condensation on the
geoid. This term N2 consists of two terms in the planar
approximation

dxdy
yyxx
yxhyxh
k
yxh
k
N
App
pp
pp
 



2/3
22
33
2
2
)()(
),(),(
6
),(




(5)
where k is the Newton’s gravitational constant,  is the
density of the topography (assumed constant with a value
of 2.67 g/cm3), (xp, yp) are the coordinates of the compu-
tation point, (x, y) are the coordinates of the integration
points,  is the normal gravity over the reference ellip-
soid, 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 N3 in Equation (3) is the con-
tribution of the residual gravity obtained by means of



dSg
R
N)( ),(
4
3 (6)
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 re-
duced anomalies in agreement to Helmert’s second met-
hod of condensation, obtained by
g = + c +
 g (7)
free
g
where is the gravity 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 Equa-
tion (7) was computed by a means of a linear approxima-
tion. The indirect effect on gravity  g was computed
from the indirect effect N2 by means of
free
g
g = 0.3086 N2 (in mgal for N2 in meters)
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
(a)
(b)
(c)
Figure 3. (a) The EIGEN-GL04C geoid model computed for
the study area. The contour interval is 1.0 m; (b) The indi-
rect effect on the geoid (plotted positive). The grey scale is
non-linear; (c) The residual geoid undulation.
C opyright © 2010 SciRes. ENG

F. Manzano ET AL. 164
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 cen-
timetre 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 measure-
ments 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 measure-
ments made with modern precise gravimeters. In spite of
this, ARGOS has better accuracy than any previous ge-
oid solutions obtained in the South Spain. This new
model will be useful for orthometric height determina-
tion by GPS as it allows the orthometric height determi-
nation in the mountains and remote areas, where level-
ling has many logistic problems.
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.
Geoid validation. To check our geoid, we have pro-
ceeded 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 ac-
curacy and reliability, fitting the geoidal heights deter-
mined for the GPS-levelling points, better than the pre-
vious geoid (IGG2005).
5. Conclusions
The new data available for the study area and the com-
putation methods based on the FFT analysis, have al-
lowed 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 Sto-
kes 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,
6. 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, Natio-
nal 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.
Table 1. Statistics of the differences between the geoid heights predicted by the available geoids and the geoid heights of the
GPS/levelling points (used as a control data set).
Geoid (name) Mean (m) St. D. (m) Maximum (m) Minimum (m) Range (m) Stations (number)
IGG2005 –0.098 0.278 0.996 0.051 0.945 262
EIGEN 0.172 0.297 0.918 0.013 0.905 262
ARGOS –0.002 0.251 0.852 0.004 0.848 262

F. Manzano ET AL.
C opyright © 2010 SciRes. ENG
165
7. References
[1] V. Corchete, M. Chourak and D. Khattach, “The
high-resolution gravimetric geoid of Iberia: IGG2005,”
Geophysical Journal International, Vol. 162, pp. 676-684,
2005.
[2] C. Wichiencharoen, “FORTRAN programs for comput-
ing geoid undulations from potential coefficients and
gravity anomalies,” Internal Report, Department of Geo-
detic Science and Surveying Report, Ohio State Univer-
sity, Columbus, 1982.
[3] I. Kuroishi, “Precise gravimetric determination of geoid
in the vicinity of Japan,” Bulletin of Geographical Survey
Institute, Vol. 41, pp. 1-94, 1995.
[4] C. Förste, F. Flechtner, R. Schmidt, R. König, U. Meyer,
R. Stubenvoll, M. Rothacher, F. Barthelmes, H. Neu-
mayer, R. Biancale, S. Bruinsma, J.-M. Lemoine and S.
Loyer, “A mean global gravity field model from the com-
bination of satellite mission and altimetry/gravimetry
surface data: EIGEN-GL04C,” Geophysical Research
Abstracts, Vol. 8, 03462, 2006.
[5] C. Förste, F. Flechtner, R. Schmidt, U. Meyer, R.
Stubenvoll, F. Barthelmes, R. König, K. H. Neumayer, M.
Rothacher, Ch. Reigber, R. Biancale, S. Bruinsma, J.-M.
Lemoine and J. C. Raimondo, “A new high-resolution
global gravity field model derived from combination of
GRACE and CHAMP mission and altimetry/gravimetry
surface gravity data,” Poster presented at EGU General
Assembly, Vienna, Austria, pp. 24-29, April 2005.
[6] W. A. Heiskanen and H. Moritz, “Physical geodesy,” W.
H. Freeman, San Francisco, 1967.