Geoid determination in the mountains using ultra-high

resolution spherical harmonic models – The Auvergne case

Rene Forsberg
ational Space Institute, Technical University of Denmark,
Juliane Maries Vej 30, DK2100 Copenhagen Ø, Denmark.
[email protected]

Preface
This paper illustrates the use of GRAVSOFT for geoid determination in the mountains,
with comparisons to high-quality GPS-leveling data in an international test data case. It
represents a typical application of the GRAVSOFT set of programs, which prof. Dimitris
Arabelos have worked with and influenced through many years, a.o. through frequent dis-
cussions and requests at his extended research visits to Denmark (especially at the former
Geodetic Institute in the late 1980’s, where many of the GRAVSOFT modules where de-
veloped), as well as research discussions in Thessaloniki and elsewhere.

Abstract
In this paper a geoid computation using the remove-restore technique and spherical FFT
is used to derive a geoid of the Auvergne region, central France. Long and medium wave-
lengths are removed by the EGM08 spherical harmonic model, while terrain reduction is
handled by the RTM method, implemented by prism integration and Fourier techniques.
The results are compared to accurate 1st order GPS-leveling data, showing an r.m.s. accu-
racy of 29 mm in the best case, which is at the best reported level for the Auvergne test data
set. The geoid comparison results vary depending on how the topography is handled rela-
tive to the inferred EGM08 topography. Due to some inherent problems in implementing
the RTM method for highly varying reference topography, it appears that best results are
obtained for a relatively low-resolution (30ʹ) reference heights, irrespective of whether
EGM08 is used at a corresponding resolution (360) or to full resolution.

1. Introduction
A major goal for physical geodesy has been the determination of the geoid at an
accuracy of 1 cm, matching the accuracy of GPS height determination. Although
the cm-geoid has been demonstrated in numerous cases in lowland areas with
dense gravity data coverage, combining gravity and (sparse) GPS-levelling data to
make an operational “GPS geoid”, so far no convincing case of attaining a cm-
geoid in mountainous regions has been reported. This is likely a consequence of

Geoid determination in the mountains using ultra-high resolution spherical harmonic models 101
– The Auvergne case
both insufficient gravity data coverage, theoretical shortcomings, and – especially –
insufficient quality of the leveling data, used to compute “ground thruth” geoid (or
quasi-geoid) values by
 = h – H or ζ = h – H* (1)
Here N is the classical geoid height, ζ the height anomaly (“geoid at the surface of
the terrain”), h the ellipsoidal height of GPS, and H and H* the orthometric and
normal heights, respectively.
A “best” comparison data set for geoid determination methods in the mountains
was distributed in 2008 by the late H. Duquenne, IGN, France, on behalf of the
International Geoid Service (IGeS), as a “ground thruth” example for precise geoid
determination method. A number of international groups have already worked with
this data set and reported results using different methods (for a summary see, e,g,
Bazarghi and Sanso, 2009).
The Auvergne data set consist of about 240.000 gravity data points from the
Bureau Gravimetrique (provided by the French geological agency BRGM), cover-
ing a 6° × 8° area including most of France; a DEM based on 3” SRTM height data
covering a somewhat larger 8° × 10° area; and a set of 75 GPS-levelling points in
the “Massif Centrale” area, all with 1st order leveling connections, and a quoted
GPS ellipsoidal height accuracy of 2-3 cm (RBF points) or “slightly better” (NI-
VAG points), cf. Fig. 1 (H. Duquenne, pers.comm.). The elevations of the GPS-
leveling points range from 206 to 1235 m, and the highest mountain in the central
area is 1886 m; it is therefore a relatively moderate mountainous area.

Fig. 1: Left: data coverage for geoid computation. Right: distribution of GPS-leveling
data in the central area.

102 Rene Forsberg

2. Geoid computation by the remove-restore method: Methodol-
ogy
The methodology for geoid construction is based on remove-restore techniques.
The anomalous gravity potential T is split into three parts
T = TEGM 08 + TRTM + Tres (2)

where TEGM08 is the contribution of the EGM08 spherical harmonic expansion

(Pavlis and Kenyon, 2008), of form
 n n
GM Ê Rˆ
TEGM 08 =
r
 ÁË r ˜¯  (C 'nm cos mλ + Snm sin mλ)Pnm (sin φ) (3)
n=2 m=0

EGM08 is here used to maximal degree either 360 or 2190. TRTM are the terrain
effects, and Tres the residual gravity field. T is treated as a spatial function, and
with the French national height system being normal heights, I will in the sequel
work consistently with the quasigeoid as the fundamental “geoid”, i.e.
T (φ, λ, H *)
ζ= (4)
γ0

where γ0 is the normal gravity, and ϕ and λ the geographical latitude and longi-
tude. Another geoid-like surface to be used is the quasi-geoid at sea level
T (φ, λ, 0)
ζ*= (5)
γ0 (φ)

ζ* is the “geoid” corresponding to the harmonically downward continued anoma-

lous potential T, neglecting any mass associated with the topography. It is there-
fore not the same as the classical geoid , which corresponds to the actual equipo-
tential surface inside the mass. The separation of the geoid and the quasigeoid is to
first order expressed by the Bouguer anomaly ΔgB
Δg B
ζ - ª- H (6)
γ0

whereas the separation between ζ* and the quasigeoid ζ to first order involves
the free-air anomaly (or, as a slightly better approximation, the gravity disturbance
δg)
δg
ζ - ζ*ª - H (7)
γ0

The above equations follows easily from the basic Chapter 8 on Molodensky’s
theory in Heiskanen and Moritz (1967). Using equation (6) or (7) is therefore

Geoid determination in the mountains using ultra-high resolution spherical harmonic models 103
– The Auvergne case
straightforward in geoid determination computations to either work with ζ or ;
it is a simple task to convert between the two quantities, as long as the first-order
approximation can be accepted.
For the terrain-reduction, the terrain effects are reduced relative to a mean ele-
vation surface, as shown schematically in Fig. 2. The terrain potential is subtracted
from the observations using a prism integration, i.e. representing the mass between
the actual topography and the mean elevation surface (m.e.s.) as mass prisms of
either positive or negative density, nominally 2.67 g/cm3. In practice the geoid “re-
store” signal is computed by Fourier methods, for details see Forsberg (1984,
1985).
In the RTM method the resolution of the m.e.s. is controlled by the user,
through a suitable low-pass filter. Ideally such a filter should ideally correspond to
the equivalent resolution of the highest spherical harmonic expansion degree N in
the reference potential. This has worked well with e.g. reference fields to degree
and order 360 (such as EGM96).
The prism implementation of the RTM method has an inherent problem: the
method leaves a point P above the m.e.s. in the mass-free domain, whereas a point
Q below the m.e.s. after the reduction would correspond to the value inside the
reference topography mass, cf. Fig. 2.

Fig. 2: Principle of RTM reduction.

The mean elevation surface is computed by low-pass filtering of a DEM.

Since all geodetic gravity field modeling methods require observations derived
from a harmonic function, i.e. in a mass-free environment above the geoid, a cor-
rection must be made after the prism computations. This correction – the harmonic
correction – to first order only applies to gravity anomalies, and only for points
below the m.e.s. It has the magnitude
Δg hc = -4πGρ(hQ - href ) (8)

where G is the gravitational constant, and ρ the density; for details see Forsberg
(1984). The equation (8), equivalent to the Prey reduction of sub-surface gravim-
etry, relies on a Bouguer plate assumption for the reference topography, and there-

104 Rene Forsberg

fore breaks down for a too quickly varying reference topography (as in the case of
the EGM08 resolution of 5ʹ). What is needed is a series expansion of the harmonic
correction, taking into account the non-level surface of the m.e.s., and being able to
compute values anywhere in the space between the m.e.s. and the zero level. This
problem has to my knowledge not been solved, in spite of the similarity to the gen-
eral “Molodensky” downward continuation problem, and the theoretical advantage
of using reference topography, not real topography.
To overcome the harmonic correction problem, alternative formulations for the
RTM method in the frequency domain may be used; such methods have e.g. been
used in the spectral combination method of the European geoid (Denker and Torge,
1998), and also in the derivation of EGM08 itself, which is based on an RTM-like
interpolation of gravity field effects from a global 15’ grid to a 5’ grid over large
parts of the continents (as a consequence of lack of data or classification of the
available 5’ gravity grid data). An alternative RTM approach, based on multipole
mass expansions, have been proposed by Vermeer and Forsberg (1992). None of
these alternative methods have been used here, so the RTM-reductions used in the
Auvergne tests are based on the “simple” prism method with the harmonic correc-
tion (8).
With the EGM08 and RTM terrain effects removed, the geoid contribution is
found by Fourier implementation of Stokes/Molodensky’s formula
R
ÚÚ ( Δgres + g1 )S (ψ )dσ
c
ζ res = (9)
4πγ
σ

where S is Stokes’ function, modified as outlined below.

When the RTM reduction is used, the Molodensky g1c-term will generally be
insignificant (Forsberg and Sideris, 1989), and the formula converts into the con-
ventional Stokes’ formula for the geoid, in principle applied to gravity field quanti-
ties at the geoid (i.e, ζ* and Δg*). The evaluation of the Stokes formula is done
using spherical Fourier techniques, for a details of the methods see Schwarz et al
(1990) and Forsberg and Sideris (1993). In the method the Stokes integral is evalu-
ated by a series of convolutions, each accurate around a certain reference latitude ϕ

( )
zres = Sref ( Dj , Dl )*[ Dg res ( j , l )sin j ˘˚ = F -1 ÈÎ F S ref F ( Dg res sin j )] (10)

where F[⋅] is the two-dimensional Fourier transform.

To keep the inherently highly accurate GRACE gravity field information in
EGM08 to be “overruled” by the influence from terrestrial gravity data, modified
Stokes functions are used. The modified Wong-Gore formulation used here is of
form
2
2n + 1
Smod (ψ ) = S (ψ ) - Â α (n) n - 1 Pn cos (ψ ) (11)
n=2

Geoid determination in the mountains using ultra-high resolution spherical harmonic models 105
– The Auvergne case
where the “GRACE-transition” coefficient α(n) increase linearly from 0 to 1 be-
tween degrees 1 and 2
Ï 1 for 2 £ n £ 1
ÔÔ  - n
α(n) = Ì 2 for 1 £ n £ 2 , n = 2,…,  (12)
Ô 2 - 1
ÔÓ 0 for 2 £ n

The estimation of 1 and 2 can only be done empirically, but values in the range
60-90 would be expected, based on the error characteristics of GRACE.

3. The Auvergne gravimetric geoid computation

The GRAVSOFT system (Tscherning et al, 2002) implements the different pro-
cedures outlined in section 2. Before the computations a slightly thinned gravity
data of 59097 points were selected (pixel binning to approx. 1 km resolution), and
at the same time the DEM was thinned to 6ʺ resolution in latitude and longitude for
ease of handling.
The EGM08 was computed to either degree 360 or 2190 in sandwich grid form
(Fig. 3), so that the gravity and potential (i.e., Δg and ζ) could be evaluated any-
where in space by three-dimensional grid interpolation; the values computed in this
way are thus harmonic continued values inside the terrain mass, not the physical
values.

Fig. 3: Principle of EGM08 sandwich grid interpolation. From the dense grids the full
three-dimensional field in the relevant height ranges are computed by linear interpolation.

The RTM terrain effects were computed either with a m.e.s. of resolution 5ʹ or
30ʹ, with the m.e.s. simply constructed by a Gaussian filter with a corresponding
half-width resolution. The mean elevation surfaces are shown in Fig. 4. Terrain
effects were subsequently subtracted from the gravity data, with statistical results
shown in Table 1. It is seen that the available gravity data show a very good fit to
EGM08, with little bias. The larger standard deviation seen for the degree 360 use
of EGM08 (and the corresponding lower resolution m.e.s.) is in principle showing
the magnitude of the non-terrain related gravity field signal in the spherical har-
monic band between 360 and 2190.

106 Rene Forsberg

Table 1: Statistics of the gravity data reductions (mgal)
EGM degree 2190 / EGM degree 360 /
No of gravity points:
RTM 5ʹ RTM 30ʺ
59097
Mean Std. dev. Mean Std. dev.
Original data 4.2 22.8 4.2 22.8
Free-air – EGM08 –2.5 11.6 –3.5 17.9
Free-air – EGM08 – RTM 0.3 3.3 0.6 8.2

Fig. 4: RTM reference elevations: Left: 5ʹ, right 30ʹ resolution. Colour scale 0 to 1500 m.

The reduced gravity data were subsequently gridded in the data coverage region
(Fig.1) to a grid of 0.01° × 0.0125° resolution by least squares collocation. The
spherical FFT conversion was done using zero-padding in a zero-padded grid of
dimension 1200 x 1280 data points, using 3 reference latitude parallels. The RTM
geoid terrain effects were subsequently restored using ζ-terrain effects computed
by FFT in a similar grid, and finally the gravimetric geoid was obtained by adding
the geoid EGM08 effects, computed at the level of the topographic surface. Fig 5
shows the magnitude of the geoid “restore” effects for the “full” EGM08 solution
(with 5’ resolution RTM) from the terrain and the residual gravity. The r.m.s. val-
ues of the terrain effects were 19 cm and 3 cm r.m.s., for the 30’ and 5’ RTM re-
ductions, respectively, illustrating the different frequency content.
The outcome of the above processing scheme is in principle the quasigeoid ζ,
when the Stokes function (9) is evaluated directly with the reduced gravity anoma-
lies Δgres. In principle a more rigorous estimate of the quasigeoid would involve a
downward continuation of Δgres to the ground level by Δgres* ≈ Δgres – Tzz*h,
with the second order gradient Tzz estimated Δgres by Fourier methods, and the
subsequent conversion of ζ* to ζ by (7); this Molodensky-style principle is im-

Geoid determination in the mountains using ultra-high resolution spherical harmonic models 107
– The Auvergne case
plemented in one of the GRAVSOFT programs (geofour) but tests did not show
any significant improvement.

Fig. 5: The geoid “restore” signals for EGM08 to degree 2190. Left: RTM terrain effect
(colour scale -20 to 40 cm). Right: Geoid effect from residual gravity data (colour
scale -20 to 20 cm).

4. Comparison to GPS-levelling
The GPS-levelling data were used first to estimate an “optimal” Wong-Gore
modification of Stokes’ function. This was done for the EGM08 field used to de-
gree 360, and 30ʹ RTM. The results are shown in Table 2. It is seen that the data
supports an optimal degree band of degree 80-90, in good agreements with the es-
timated accuracy of GRACE.

Table 2: Statistics of fit of the geoid to GPS-leveling geoid heights for different
kernel modification parameters.
Modification Mean Std.dev.
degree band (m) (m)
None (normal Stokes) –0.176 0.034
Degree 40-50 –0.084 0.035
Degree 60-70 –0.100 0.036
Degree 80-90 –0.138 0.029
Degree 100-110 –0.149 0.041
Degree 120-130 –0.147 0.057

108 Rene Forsberg

 Table 3 shows the geoid fits for the 75 GPS-levelling points, for both EGM08
alone (to degree 2190), the EGM08 (degree 360) gravimetric geoid solution, the
EGM08 quasigeoid fit, as well as a geoid solution using the full EGM08 to degree
2190, but a 30’-resolution RTM surface (i.e., taking the shorter-wavelength topog-
raphy into account twice). It is seen that the latter geoid model actually gives the
best r.m.s. fit of 29 mm, indicating that the “double accounting” of the topography
does not matter in practice (which it should not in principle, since the remove-
restore principle will account for this). That the RTM 5’ model seems to fit poorer
is likely due to the approximation errors in the harmonic correction.

Table 3: Comparison of different geoid models to GPS-leveling

Geoid model fit, 75 points Std.dev.
Mean
(unit: m)
EGM08 to degree 360 -.109 0.172
EGM08 to degree 2190 -.109 0.036
EGM08 to 360, RTM 30ʹ -.128 0.030
EGM08 to 2190, RTM 5ʹ -.114 0.036
EGM08 to 2190, RTM 30ʹ -.138 0.029

Fig. 6: Left: Differences between GPS-leveling and geoid in the central region in cm, with
the geoid model (EGM08-2190, colour range 46 to 53 m) in background. Right:
Distribution of gravity points (crosses) and GPS-leveling points (circles) in central
region; colours show DEM heights (colour scale 0 to 1800 m)

To judge the errors of the GPS-leveling comparison in details, Figure 6 shows

the geoid outlier values in cm, as well as the underlying gravity data coverage in
the central area. This coverage is relatively good, although parts of the area have

Geoid determination in the mountains using ultra-high resolution spherical harmonic models 109
– The Auvergne case
gaps of 5 km extent or more, a possible source of some of the discrepancy between
geoid and GPS.
To further check for systematic errors, Fig. 7 shows the correlation of the dif-
ferences with height; the lack of such correlation confirms a high-quality geoid
solution and been obtained, and has been compared to high-quality GPS-levelling
data. Overall, however, the agreement between GPS-leveling and geoid is probably
as good as one could expect, given the earlier quoted 2-3 cm rms accuracy of the
GPS-leveling. The experiment thus confirms that a very precise geoid computation
and validation can be done in the Auvergne area by GRAVSOFT.

0
Difference (m)

-0.05

-0.10

-0.15

-0.20
0 250 500 750 1000 1250

Height

Fig. 7: Distribution of outliers for the gravimetric geoid (using EGM08 to degree
2190) as a function of height

5. Conclusions
The computations in the Auvergne test area confirms that it is possible to carry
out a geoid determination fitting at the 3 cm r.m.s. level to high-quality GPS-
leveling, and that the residuals are uncorrelated with height. In reality it is not
really possible to estimate the actual accuracy of the geoid, since the comparison
error is close to the estimated errors of the GPS-leveling.
It is also shown that the EGM08 in itself provides an excellent fit to the geoid of
the Auvergne region, illustrating that when the underlying data are good, EGM08
also provides excellent results. Combining EGM08 with the RTM method, it is also
seen that some theoretical errors prevent the fully consistent use of residual topog-
raphy, when implemented in space domain as prism computations. It should there-
fore be avoided to use too high-resolution reference elevation surfaces in RTM
with GRAVSOFT.
For a comparison of the results of this paper with other results in the region, as
given in Barzaghi and Sanso (2009), the 29 mm r.m.s. fit of the best results com-
pare very well with the results of other methods; the 6 other groups reporting r.m.s.

110 Rene Forsberg

comparing results in the range of 29 to 65 mm, with a median result of 35 mm
r.m.s. Some possible theoretical improvements to the GRAVSOFT implementa-
tion, e.g. using improved higher-order atmospheric corrections, series expansions
of the harmonic correction, and a rigorous use of Tzz-gradients in Molodensky har-
monic continuation, could make results better. However, the ultimate 1 cm geoid in
the mountains can probably not be attained until more gravity data is available, and
the terrain reduction is done using improved density estimates.

References
Barzaghi, R. and F. Sanso: Report of the International Geoid Service. Travaux de
l’Association Internationale de Geodesie 2007-9, available at www.iges.polimi.it, 2009.
Denker, H and W. Torge: The European Gravimetric Quasigeoid EGG97 - An IAG sup-
ported continental enterprise. In: R. Forsberg, M. Feissel, R. Dietrich (eds.): Geodesy
on the Move - Gravity, Geoid, Geodynamics and Antarctica. IAG Symp., Vol. 119,
249-254, Springer Verlag, Berlin, Heidelberg, New York, 1998
Forsberg, R.: A Study of Terrain Reductions, Density Anomalies and Geophysical Inversion
Methods in Gravity Field Modelling. Reports of the Department of Geodetic Science
and Surveying, No. 355, The Ohio State University, Columbus, Ohio, 1984.
Forsberg, R.: Gravity field terrain effect computations by FFT. Bull. Geod., 59, pp. 342-
360, 1985.
Forsberg, R. and M.G. Sideris: On topographic effects in gravity field approximation. In:
Kejlsø, Poder and Tscherning (Eds.): Festschrift to Torben Krarup, Geodætisk Institut
Meddelelse (Reports of the Geodetic Institute - Denmark) no. 58, pp. 129-148, 1989.
Forsberg, R., M. G. Sideris: Geoid computations by the multi-band spherical FFT ap-
proach. Manuscripta Geodaetica, 18, pp. 82-90, 1993.
Heiskanen, W. A., H. Moritz: Physical Geodesy. Freeman, San Francisco, 1967.
Pavlis, N, S. Holmes, S. Kenyon and J. Factor: An Earth Gravitational Model to Degree
2160: EGM2008. EGU General Assembly 2008, Vienna, Austria, April, 2008.
Schwarz, K.-P., M. G. Sideris, R. Forsberg: Use of FFT methods in Physical Geodesy.
Geophysical Journal International, vol. 100, pp. 485-514, 1990.
Tscherning, C. C., R. Forsberg and P. Knudsen: The GRAVSOFT package for geoid deter-
mination. Proc. 1st continental workshop on the geoid in Europe, Prague 1992, pp. 327-
334.
Vermeer, M. and R. Forsberg: Filtered terrain effects: A frequency domain approach to
terrain effect evaluation. Manuscripta Geodaetica, vol. 17, pp. 215-226, 1992.

Geoid determination in the mountains using ultra-high resolution spherical harmonic models 111
– The Auvergne case