Physical Geodesy
Physical Geodesy
Physical Geodesy
Mass excess
Mass
deficit
−N
Martin Vermeer1
1
[email protected]
Preface
This book aims to present an overview over the current state of the study
of the Earth’s gravity field and those parts of geophysics closely related
to it, such as especially geodynamics, the study of the changing Earth.
It grew out of two decades of teaching in Helsinki’s two universities:
Helsinki University of Technology — today absorbed into Aalto University
— and the University of Helsinki. As such, it presents a somewhat
Fennoscandian perspective on a very global subject. Also the author’s
own research, on gravimetric geoid determination, helped shape the
presentation. While there exist excellent textbooks on all parts of what is
presented here, he may still hope that this text will find a niche to fill.
Martin Vermeer
Acknowledgements
Thanks are due to the many students and colleagues, both in academia
and at the Finnish Geodetic Institute, giving useful comments and cor-
rections over the course of many years of lecturing this course both at
the University of Helsinki and at the Helsinki University of Technology,
today Aalto University.
Special thanks are due to the foreign students at Aalto University, who
forced me during recent years to provide an English version of this text.
The translation work prompted a basic revision also to the Finnish text,
which was long overdue as parts were written in the 1990s before the
author had had the benefit of pedagogical training. Thanks are thus also
due to the Aalto University’s pedagogical training programme.
– iii –
Contents
Chapters
Preface iii
–v–
vi Contents
Abbreviations xix
Óîá .
Contents vii
Óîá .
viii Contents
5. Geophysical reductions 95
5.1 General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.2 Bouguer anomalies . . . . . . . . . . . . . . . . . . . . . . . 96
5.3 Terrain effect and terrain correction . . . . . . . . . . . . . 99
5.4 Spherical Bouguer anomalies . . . . . . . . . . . . . . . . . 103
5.5 Helmert condensation . . . . . . . . . . . . . . . . . . . . . 104
5.6 Isostasy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.7 Isostatic reductions . . . . . . . . . . . . . . . . . . . . . . . 112
5.8 The “isostatic geoid” . . . . . . . . . . . . . . . . . . . . . . 113
Self-test questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Exercise 5 – 1: Gravity anomaly . . . . . . . . . . . . . . . . . . . . 117
Exercise 5 – 2: Bouguer reduction . . . . . . . . . . . . . . . . . . . 117
Exercise 5 – 3: Terrain correction, Bouguer . . . . . . . . . . . . . 118
Exercise 5 – 4: Isostasy . . . . . . . . . . . . . . . . . . . . . . . . . 118
Óîá .
Contents ix
Óîá .
x Contents
Óîá .
Contents xi
Óîá .
xii List of Tables
Bibliography 331
Index 339
List of Tables
Óîá .
List of Figures xiii
List of Figures
Óîá .
xiv List of Figures
Óîá .
List of Figures xv
Óîá .
xvi List of Figures
Óîá .
Abbreviations
– xvii –
xviii Abbreviations
GRAVSOFT
Geopotential determination software, mainly developed in
Denmark
GRS67 Geodetic Reference System 1967
GRS80 Geodetic Reference System 1980
HUT Helsinki University of Technology, today a part of Aalto Uni-
versity
IAG International Association of Geodesy
ICET International Center for Earth Tides
ICGEM International Center for Global Earth Models
IDEMS International Digital Elevation Model Service
IGeC International Geoid Commission (obsolete)
IGeS International Geoid Service (obsolete)
IGFS International Gravity Field Service
ISG International Service for the Geoid
IUGG International Union of Geodesy and Geophysics
JILA Joint Institute for Laboratory Astrophysics (Boulder, Col-
orado, USA)
KKJ National Grid Co-ordinate System (Finland)
Lego™ “Leg Godt”, en. “Play Well”, Danish toy brand
LLR Lunar laser ranging
LSC least-squares collocation
Mf Moon, fortnightly tide
moho Mohorovičić discontinuity
N60 height system (Finland)
N2000 height system (Finland)
NAP Normaal Amsterdams Peil, height system (Western Europe)
NASA National Aeronautics and Space Administration (USA)
NAVD88 North American Vertical Datum 1988
NC normal correction
NGA National Geospatial-Intelligence Agency (USA, formerly
NIMA)
NIMA National Imagery and Mapping Agency (USA, formerly DMA)
NKG Nordiska Kommissionen för Geodesi, Nordic Geodetic Com-
mission
NKG2004 geoid model (Nordic area)
NKG2015 geoid model (Nordic area)
NOAA National Oceanic and Atmospheric Administration (USA)
OSU Ohio State University, Columbus, Ohio, USA
OC orthometric correction
ppm parts per million
ppb parts per billion
RTM residual terrain modelling
SI Système International d’Unités
Óîá.
Abbreviations xix
Óîá.
D Fundamentals of the theory of gravitation
1
D 1.1 General
In this chapter we discuss the foundations of Newton’s theory of gravi-
tation. Intuitively, the theory of gravitation is easiest to understand as
“action at a distance”, latin actio ad distans, where the force between two
masses is proportional to the masses themselves, and inversely propor-
tional to the square of the distance between them. This is the form of
Newton’s general law of gravitation familiar to all.
There exists an alternative but equivalent presentation, field theory,
which describes gravitation as a phenomenon propagating through space,
a field. The propagation is described in the field equations. The field
approach isn’t quite as intuitive, but is a powerful theoretical tool1 .
In this chapter we acquaint ourselves with the central concept of field
theory, the gravitational potential. We also get to know the potential fields
of the theoretically interesting single and double mass density layers. Of
the practical and theoretical applications of these may be mentioned the
Bouguer plate and Helmert condensation. In the following we will discuss
their properties in detail. Mass density layers are also used in deriving
the theorems of Green. We will learn about important integral theorems
like the theorems of Gauss and Green, with the aid of which we may
infer the whole potential field in space from field values given only on a
certain surface. Other similar examples are the Chasles theorem and the
solution to Dirichlet’s problem.
In chapter 2 we apply these fundamentals of potential theory to derive
a spectral representation of the Earth’s gravitational field, a so-called
1
There is also a philosophical side to this. To many, e.g., to Leibniz, the idea
of a force that jumps from object to object was an abomination. Many tried to
explain gravitation — and also electromagnetism etc. — by a “world aether”.
It wasn’t until the advent of relativity theory, that the understanding gained
ground that a physical theory doesn’t have to satisfy our preconceptions about
what constitutes a “sensible explanation” — as long as it describes the physical
phenomena correctly.
–1–
2 Fundamentals of the theory of gravitation
spherical-harmonic expansion.
Here, in the beginning of the text, we derive an extensive set of math-
ematical equations, such as well-known integral equations. This is an
unfortunately necessary preliminary work. These equations, however,
are no end in themselves and they are not worth committing to memory.
Try rather to understand their logic, and how historically these vari-
ous results have been arrived at. Try as well to acquire an intuition, a
fingertip feeling, about the nature of this theory.
Ó »îá .
Gravitation between two masses 3
Let us call the small body, the test mass, e.g., a satellite, m, and the
large mass, the planet or the Sun, M . Then, m 1 = M may be called the
attracting mass, and m 2 = m the attracted mass, and we obtain
mM
F =G .
ℓ2
According to Newton’s law of motion
F = ma,
M
a=G .
ℓ2
From this equation, the quantity m = m 2 has vanished. This is the famous
observation by Galileo, that all bodies fall equally fast4 , irrespective of
their mass. This is also known as Einstein5 ’s principle of equivalence.
2
Sir Isaac Newton PRS (1642 – 1727) was an English universal genius who
derived the mathematical underpinning of astronomy, and much of geophysics, in
his main work “Philosophiæ Naturalis Principia Mathematica”, “Mathematical
Foundations of Physics”.
3
Henry Cavendish FRS (1731 – 1810) was a British natural scientist from a
wealthy nobility background. He did also pioneering work in chemistry. He was
extremely shy, and the renowned neurologist Oliver Sacks retrodiagnosed him
as afflicted with the Asperger syndrome (Sacks, 2001).
4
At least in vacuum. The Apollo astronauts showed impressively, how on the
Moon a feather and a hammer fall equally fast! https://www.youtube.com/watch?
feature=player_embedded&v=KDp1tiUsZw8#!.
5
Albert Einstein (1879 – 1955) was a theoretical physicist of Jewish German
descent, who created both the special and general theories of relativity, applying
the latter to cosmology, and did fundamental work in quantum theory.
Ó »îá .
4 Fundamentals of the theory of gravitation
Both the force F and the acceleration a have the same direction as the
line connecting the bodies. For this reason one often writes equation 1.1
as a vector equation, which is more expressive:
r−R
a = −GM , (1.2)
ℓ3
where the three-dimensional vectors of place of both the attracting and
attracted masses are defined as follows in rectangular co-ordinates6 :
r = xi + yj + zk,
R = X i + Y j + Z k,
Ó »îá .
Potential of a spherical shell 5
bd θ
p b
P
r
O
θ
ℓ
b
This is an integral over mass elements dm, where every such mass
element is located at place R. The potential V is evaluated at location r,
and the distance ℓ = ∥r − R∥.
We now derive the formula for the potential of a thin spherical shell, see
figure 1.2, where we have placed the centre of the sphere at the origin O .
Because the circumference of a narrow ring, width b · d θ , is 2π b sin θ , its
surface area is
(2π b sin θ ) ( b · d θ ) .
Ó »îá .
6 Fundamentals of the theory of gravitation
Let the thickness of the shell be p (small) and its density ρ . We obtain
for the total mass of the ring
2π pρ b2 sin θ d θ .
Because every point of the ring is at the same distance ℓ from point P,
we may write for the potential at point P :
2πG pρ b2 sin θ d θ
VP = .
ℓ
With the cosine rule,
ℓ2 = r 2 + b2 − 2 rb cos θ , (1.6)
we obtain, using equation 1.5, for the potential of the whole shell
ˆ
2 sin θ d θ
VP = 2πG ρ pb p .
r 2 + b2 − 2 rb cos θ
In order to evaluate this integral, we must replace the integration vari-
able θ by ℓ. Differentiating equation 1.6 yields
ℓ d ℓ = br sin θ d θ ,
p
and remembering that ℓ = r 2 + b2 − 2 rb cos θ we obtain
ˆ ℓ2
2 dℓ
VP = 2πG ρ pb .
ℓ1 br
In the case that point P is outside the shell, the integration bounds of ℓ
are ℓ1 = r − b and ℓ2 = r + b, and we obtain for the potential of point P
]ℓ=r+b
ℓ 4πG ρ pb2
[
2
VP = 2πG ρ pb = .
br ℓ= r − b r
8
Here, the ∇ (nabla) operator is used, to be explained in section 1.5.
Ó »îá .
Computing the attraction from the potential 7
4πG ρ pb
Acceleration
2
4πG ρ p br2
4πG ρ pb br
Potential
0
0 b r
Figure 1.3. Dependence of potential and attraction on distance r from the centre
D of a spherical shell.
Ó »îá .
8 Fundamentals of the theory of gravitation
−
→ −−→ ∂V ∂V ∂V
a = ∇ V = gradV = i +j +k . (1.7)
∂x ∂y ∂z
∂V ∂V ∂ℓ 1 x− X x− X
= = GM · − = −GM 3 .
∂x ∂ℓ ∂ x ℓ 2 ℓ ℓ
Similarly we compute the y and z components:
∂V y − Y ∂V z−Z
= −GM ; = −GM 3 .
∂y ℓ 3 ∂z ℓ
E pot = −V m.
In practice one calls the vector of gravitational acceleration the “gravita-
tional vector”.
9
From the equation
√ ] 12
(x − X )2 + (y − Y )2 + (z − Z)2 = (x − X )2 + (y − Y )2 + (z − Z)2
[
ℓ=
Ó »îá .
Potential of a solid body 9
dm ( x, y, z)
ρ ( x, y, z) = ,
d V ( x, y, z)
As we showed already above for mass points, also the first derivative with
respect to place or gradient of the potential V of a solid body,
−−→ −
→
gradV = ∇ V = a, (1.9)
Ó »îá .
10 Fundamentals of the theory of gravitation
and thus
1
→ 0 when ∥r∥ → ∞.
ℓ
For the acceleration of gravitation the same applies for all three compo-
nents, and thus also for the length of this vectorial quantity:
−→
∥r∥ → ∞ =⇒ ∇ V → 0.
This result can still be sharpened: if ∥r∥ → ∞, then, again by the triangle
inequality,
ℓ = ∥r − R∥ ≤ ∥r∥ + ∥R∥ ≤ ∥r∥ + ϵ,
and thus
1 1 1 1 1 1 1 1
≤ ≤ =⇒ ≤ ≤ .
∥r∥ + ϵ ℓ ∥r∥ − ϵ ∥r∥ 1 + /∥r∥ ℓ ∥r∥ 1 − ϵ/∥r∥
ϵ
1 1
r → ∞ =⇒ → .
ℓ r
When we substitute this into the above integral 1.8, it follows that for
large distances r → ∞:
˚ ˚
ρ G GM
V =G dV ≈ ρ dV = ,
body ℓ r body r
where M, the integral of density over the volume of the body, is precisely
its total mass. From this we see, that at great distance, the field of a
finite sized body M is nearly identical with the field generated by a point
mass the total mass of which is equal to the total mass of the body. This
important observation was already made by Newton. As a result of this
phenomenon we can treat, in the celestial mechanics of the Solar system,
the Sun and planets10 as mass points, although we know that they are
not.
10
The only important exception is formed by the forces between a planet and its
moons, both due to the flattening of the planet and due to tidal effects.
Ó »îá .
Example: The potential of a line of mass 11
Now we can expand this into a Taylor series in H around the point H = 0:
the first derivative of equation 1.10 is
∂V G G
=√ =
∂H 2 2
( X − x ) + (Y − y ) + ( H − z )2 ℓ
√
in which ℓ ( H ) = ( X − x)2 + (Y − y)2 + ( H − z)2 . The second derivative is
(chain rule)
∂2 V ∂ 1 H−z
= ℓ ( H )−1 = − ℓ−3 · 2 ( H − z) = −G 3 .
∂H 2 ∂H 2 ℓ
The third derivative, in the same way:
∂3 V ∂ 3 ( H − z)2 ℓ2 3 ( H − z)2 − ℓ2
( )
H−z
= − 3 = − = G ,
∂H 3 ∂H ℓ ℓ5 ℓ5 ℓ5
1 1 z 1 3 z2 − ℓ20 3
V= 0 +G H + G 3 H2 + G H +... (1.11)
ˆ ℓ0 2 ℓ0 6 ℓ50
H =0
1
dZ
0 ℓ
√
in which ℓ0 = ( X − x)2 + (Y − y)2 + z2 .
Ó »îá .
12 Fundamentals of the theory of gravitation
Question: how can we exploit this result for computing the gravitational
potential of a complete, realistic topography?
1 1 z
Answer: in this expansion, the coefficients , , . . . like ℓ0 , depend
ℓ0 2 ℓ30
only on the differences ∆ x = X − x and ∆ y = Y − y between the co-
ordinates of the location of the mass line ( X , Y ) and those of the
evaluation location ( x, y) — and of the height z of the evaluation
point. If the topography is given in the form of a grid, we may
evaluate the above expansion 1.11 for the given z value and for all
possible value pairs (∆ x, ∆ y). Then, if the grid is, e.g., N × N in size,
we need only N 2 operations for calculating every coefficient. The
brute-force evaluation of the Taylor expansion itself for the whole
topography, i.e., for every point of the terrain grid and every point of
the evaluation grid, requires after that N 4 operations, but those are
much simpler: the coefficients themselves have been precomputed.
And brute force isn’t even the best approach: as we shall see, the
above convolution can be computed much faster using the Fast
Fourier Transform.
We shall return to this subject more extensively in connection with
terrain effect evaluation, sections 5.3 and 8.7.
where
def
⟨−
→ −
→⟩ ∂2 ∂2 ∂2
∆= ∇·∇ = + +
∂x 2 ∂ y2 ∂ z2
is a well known symbol called the Laplace11 operator.
In equation 1.4 for the potential of a point mass we may show, by per-
forming all partial derivations 1.12, that
∆V = 0, (1.13)
11
Pierre-Simon, marquis de Laplace (1749 – 1827) was a French universal genius
who contributed to mathematics and natural sciences. He is one of the 72 French
scientists, engineers and mathematicians whose names were inscribed on the
Eiffel Tower, https://en.wikipedia.org/wiki/List_of_the_72_names_on_the_Eiffel_
Tower.
Ó »îá .
Gauge invariance 13
the well known Laplace equation. This equation applies outside a point
mass, and more generally everywhere in empty space: all masses can
in the limit be considered to consist of point-like mass elements. Or, in
equation 1.8 we may directly differentiate inside the triple integral sign,
exploiting the circumstance that it is allowed to interchange integration
and partial differentiation, if both are defined.
A potential field for which the Laplace equation 1.13 applies, is called a
harmonic field or function.
In the case where the mass density doesn’t vanish everywhere, we have a
different equation, with ρ the mass density:
∆V = −4πG ρ . (1.14)
is known as the field equations of the gravitational field. They play the
same role as Maxwell’s13 field equations in electromagnetism. Unlike in
Maxwell’s equations, however, in the above there is no time co-ordinate.
Because of this, it is not possible to derive a formula for the propagation
in space of gravitational waves, like the one for electromagnetic waves in
Maxwell’s theory.
We know today that these “Newton field equations” are only approximate,
and that a more precise theory is Einstein’s general theory of relativity.
Nevertheless, in physical geodesy Newton’s theory is generally precise
enough and we shall use it exclusively.
12
Siméon Denis Poisson (1781 – 1840) was a French mathematician, physicist
and geodesist, one of the 72 names inscribed on the Eiffel Tower.
13
James Clerk Maxwell FRS FRSE (1831 – 1879) was a Scottish physicist, the
discoverer of the field equations of electromagnetism. He found a wave-like
solution to the equations, and, based on propagation speed, identified light as
such.
Ó »îá .
14 Fundamentals of the theory of gravitation
Ó »îá .
Double mass density layer 15
integral 1.15 we may prove (in a complicated way) that the exterior
potential is the same as it would be if all of the mass of the body were
concentrated in the sphere’s centre. Earlier, in section 1.4, we proved
that the potential interior to the sphere is constant.
Thus, the exterior attraction (ℓ > R ) is
( )2
M κ · 4π R 2 R
a e (ℓ ) = G 2 = G = 4πG κ .
ℓ ℓ2 ℓ
a i (ℓ) = 0.
This means that on the surface of the sphere, the attraction is discontin-
uous:
a e (R ) − a i (R ) = 4πG κ.
∂V
a = ∥a∥ = (1.16)
∂n
where the differentiation variable n represents the normal direction,
i.e., the direction perpendicular to the surface S . If the surface S is an
equipotential surface of the potential V , equation 1.16 applies generally;
then, the attraction vector — more precisely, the acceleration vector — is
perpendicular to the surface S , and its magnitude is equal to that of the
normal derivative of the potential.
dM
µ= ,
dS
in which d M is a “dipole layer element”; this layer may be seen as made
up of two single layers. If we have a positive layer at density κ and a
negative layer at density −κ, and the distance between them is δ, we get
for small values of δ an approximate correspondence:
µ ≈ δκ.
Ó »îá .
16 Fundamentals of the theory of gravitation
P n
ℓ
δ κ
−κ
The combined potential of the two single mass density layers computed
as explained in the previous section is
¨ ( )
1 1
V =G κ − dS .
surface ℓ1 ℓ2
Ó »îá .
The Gauss integral theorem 17
by computing the surface integral using the sphere’s centre as the evalu-
ation point, and using the earlier established circumstance that inside a
sphere covered by a single mass density layer the potential is constant.
Now in the limit ℓ → R the result is different for the exterior and interior
potentials. The difference is
Ve (R ) − Vi (R ) = 4πG µ.
The Gauss14 integral theorem, famous from physics, looks in its vector
form like this: ˚ ¨
div a d V = 〈a · n〉 dS, (1.18)
V ∂V
in which n is the exterior normal to surface S , now as a vector: the length
of the vector is assumed ∥n∥ = 1. ∂V is the surface of body V .
This theorem applies to all differentiable vector fields a and all “well
behaved” bodies V on whose surface ∂V , everywhere a normal direction
n exists. In other words, this is not a special property of the gravitational
acceleration vector, though it applies also to this vector field.
14
Johann Carl Friedrich Gauss (1777 – 1855) was a German mathematician and
universal genius. “Princeps mathematicorum”.
15
But the charges for gravitation, i.e., masses, are always positive.
Ó »îá .
18 Fundamentals of the theory of gravitation
Field line
.
Flux
.
.
.
. .
Body
Sources surface
.
.
. .
.
.
.
Figure 1.5. A graphical explanation of the Gauss integral theorem. The concept
D of field line was Michael Faraday’s invention.
just like a liquid flow — from the space inside the surface S to the outside
through S .
The Gauss integral theorem states the two amounts to be equal: it is
in a way a book-keeping statement demanding that everything which is
produced inside a surface — div a — has also to come out through the
surface — 〈a · n〉.
In figure 1.5 it is graphically explained, that the sum of “sources” over
∑
the inner space of the body, i.e., (+ + + . . .), has to be the same as the
∑
sum of “flux” (↑↑↑ . . .) over the whole boundary surface delimiting this
inner space.
Let us write the Gauss equation a bit differently, using the potential
instead of the gravitational vector:
˚ ¨
∂V
∆V d V = dS, (1.19)
V ∂V ∂n
where we have done the above substitutions. We also here see the popular
notation ∂V for the surface of the body V . The presentational forms 1.19
and 1.18 are connected by the equations 1.12 and 1.9, connecting V and
a.
Ó »îá .
The Gauss integral theorem 19
a+z
a+y
∆z a+
x
a−
x
a−y
a−z
∆y
∆x
Here a+x is the value of component a x on the one face in the x direction
and a−x its value on the other face, etc. E.g., a+z is the value of a z on the
box’s upper and a−z on its lower face. A box has of course six faces, in each
of three co-ordinate directions both “up- and downstream”.
Then
∂a x
a+x − a−x ≈ ∆ x,
∂x
∂a y
a+y − a−y ≈ ∆ y,
∂y
∂a z
a+z − a−z ≈ ∆ z,
∂z
and by substitution we see that
¨
∂a x ∂a y ∂a z
〈a · n〉 dS ≈ ∆x · ∆ y ∆z + ∆ y · ∆x ∆z + ∆z · ∆x ∆ y =
∂V ∂x ∂y ∂z
∂a x ∂a y ∂a z
( )
= + + ∆ x ∆ y ∆ z,
∂x ∂y ∂z
Ó »îá .
20 Fundamentals of the theory of gravitation
the same equation as 1.20. So, in this simple case the Gauss equation
applies.
Obviously the equation works also, if we build out of these “Lego™ bricks”
a larger body, because the faces of the bricks touching each other are
oppositely oriented and cancel from the surface integral of the whole body.
It is a bit harder to prove that the equation also applies to bodies having
inclined surfaces.
∆V = −4πG ρ .
16
Paul Adrien Maurice Dirac (1902 – 1984) was an English quantum physicist
who found the relativistic wave equation for the electron, and theoretically
predicted the existence of antimatter. Physics Nobel laureate 1933, shared with
Erwin Schrödinger. He is also believed to have been on the autism spectrum
(https://en.wikipedia.org/wiki/The_Strangest_Man).
Ó »îá .
Green’s theorems 21
+1
z
0
GM, @ (0, 0, 0)
x
−1 y
0 −1
0
y +1
x
integration. The surface integral is six times that over the top face
ˆ +1
[ˆ
+1
]
1
−GM ( )3/2 dx d y.
−1 −1 x2 + y2 + 1
Ó »îá .
22 Fundamentals of the theory of gravitation
17
George Green (1793 – 1841) was a British mathematical physicist, an au-
todidact, working as a miller near Nottingham. He also invented the word
“potential”.
http://www-history.mcs.st-and.ac.uk/Biographies/Green.html,
http://www.greensmill.org.uk/,
http://en.wikipedia.org/wiki/An_Essay_on_the_Application_of_Mathematical_
Analysis_to_the_Theories_of_Electricity_and_Magnetism.
Ó »îá .
Green’s theorems 23
Distance ℓ
Surface element d V
Surface
normal
n
Body V
Surface S = ∂V
Figure 1.8. Geometry for deriving Green’s third theorem if point P is outside
D surface ∂V .
Ó »îá .
24 Fundamentals of the theory of gravitation
Surface ∂V , part 1
Surface ∂V , part 2
Point P Volume V
Figure 1.9. Geometry for deriving Green’s third theorem if point P is inside
D surface ∂V .
Ó »îá .
The Chasles theorem 25
Point P
Integration space V
Boundary ∂V , part 2
Matter
Boundary ∂V , part 1
Boundary ∂V , part 3
(Limit)
D Figure 1.10. Green’s third theorem for the space external to a body.
From the Green theorem 1.24 derived above, we may derive for a har-
monic function V (so, ∆V = 0) in the exterior space:
¨ ¨
1 1 ∂V 1 ∂ (1)
VP = − dS + V dS . (1.25)
4π ∂V ℓ ∂n 4π ∂V ∂n ℓ
Ó »îá .
26 Fundamentals of the theory of gravitation
the right-hand side integral vanishes based on the Gauss integral theo-
rem. This is because the function 1/ℓ (with ℓ the distance from point P ) is
harmonic inside the Earth’s body of which ∂V is the surface. This is the
Chasles theorem18 , also called the Green equivalent-layer theorem.
The theorem is used in Molodensky’s19 theory. Also the representation of
the Earth’s gravity field by underground point-mass layers, e.g., Vermeer
(1984), might be justified with this theorem.
The case where ∂V is an equipotential surface is realized if the body
is fluid and seeks by itself an external form equal to an equipotential
surface. For planet Earth, this applies for the ocean surface. Also in
electrostatic theory, for a conductor in which the electrons can move
freely, the physical surface will become one equipotential surface. This is
why it is often stated that the electric charges are on the surface of the
conductor. This isn’t necessarily so, but from a practical viewpoint the
result is the same.
Equation 1.25 simplifies in this case as follows:
¨ ¨
1 1 ∂V κ
VP = − dS = G dS . (1.26)
4π ∂V ℓ ∂n ∂V ℓ
The equation tells us that we can compute the whole potential exterior to
the Earth, if only on the surface of the Earth (the shape of which we also
assume given in order to compute 1 /ℓ!) is given the gradient, in the normal,
i.e., vertical, direction, of the potential, ∂V /∂n. This gradient is precisely
the gravitational acceleration, a quantity derivable from measurement.
All of gravimetric geopotential determination (“geoid determination”),
ever since G. G. Stokes, is based on this.
18
Michel Chasles (1793 – 1880) was a French mathematician and geometrician,
one of the 72 whose names are inscribed on the Eiffel Tower.
19
Mikhail Sergeevich Molodensky (1909 – 1991) was an illustrious Russian phys-
ical geodesist.
20
Peter Gustav Lejeune Dirichlet (1805 – 1859) was a German mathematician
also known for his contributions to number theory.
Ó »îá .
What the boundary-value problem cannot compute 27
L {V } ,
Ó »îá .
28 Fundamentals of the theory of gravitation
D Self-test questions
Ó »îá .
Exercise 1 – 2: Atmosphere 29
g
g0
Σ1
d
∞ ∞
Σ2
3. What is the attraction g at the centre of the core? What can you say
in general about the geopotential in this point (don’t try to calculate
it)?
∂g
4. Derive the equation for the gravity gradient on the surface of a
∂r
homogeneous-density sphere of density ρ .
D Exercise 1 – 2: Atmosphere
1. The mean pressure of the atmosphere at sea level is 103 hPa (the
unit Pascal: Pa = Nm−2 .) On the Earth surface gravity is 10 m /s2 .
Calculate the mean surface density as a thin layer κ in units of
kg /m2 .
2. Calculate the total mass of the atmosphere using the spherical shell
approximation. You may take as its radius 6378 km.
3. Calculate the attraction generated by the atmosphere outside it,
both as acceleration and as a fraction of the total Earth attraction.
4. How much is the attraction from the atmosphere inside the atmo-
sphere?
There is a deposit (body) of iron ore inside the Earth, which generates (in
the flat Earth approximation!) an attraction on the Earth surface, which
has been drawn as the g curve. See figure 1.11.
The true attraction curve is approximated by a simple function
{
g 0 if r≤d
g=
0 if r>d
Ó »îá .
30 Fundamentals of the theory of gravitation
(red dashed line), where r is the distance from the point on the Earth
surface straight above the ore deposit. So, the area where g ̸= 0 for a disk
of radius d on the surface of the Earth.
= − 4πGMbody ,
calculate GMbody .
3. Assuming that the deposit is a sphere at depth d , calculate GM
using Newton’s law of gravitation from the value g 0 straight above
the deposit at the Earth surface.
4. Compare results 2. and 3. and draw conclusions. Is the function g
given above a good approximation?
21
Be careful with the algebraic signs!
Ó »îá .
D The Laplace equation and its solutions
2
D 2.1 General
An equation central to the study of the Earth’s gravitational field is the
Laplace equation,
∂2 ∂2 ∂2
( )
∆V = + + V = 0.
∂ x2 ∂ y2 ∂ z2
We call the symbol ∆ the Laplace operator. In many texts, the alternative
notation ∇2 is used.
If we study gravitation as a field, then the Laplace equation is more
natural than Newton’s formalism. Newton’s equations are used when
the mass distribution is known: it yields directly the gravitational force
caused by the masses.
The Laplace equation on the other hand is a partial differential equation.
Its solution gives the potential of the gravitational field thoughout space
or a part of space. From this potential one may then calculate the effect
of the field on a body moving in space at the location where the body is.
This is a two-phased process. The conceptual difference is, that a certain
property, a field, is attributed to empty space, and we no longer talk about
action at a distance directly between two bodies.
Solving the Laplace equation in the general case may be difficult. The
approach generally taken is, that we choose some co-ordinate frame —
a rectangular frame (as above), spherical co-ordinates, cylindrical co-
ordinates, toroidal co-ordinates, or whatever — which fits best with
the geometry of the problem at hand. Then, we transform the Laplace
equation to those co-ordinates, we find special solutions of a certain form,
and finally we compose a general — or not-so-general — solution as a
linear combination of those special solutions, i.e., a series expansion.
Fortunately the theory of linear partial differential equations is well
developed. Similar theoretical problems are encountered in the theory
of the electromagnetic field (Maxwell theory) and quantum mechanics
– 31 –
32 The Laplace equation and its solutions
V = αV1 + βV2 , α, β ∈ R
V ( x, y, z) = X ( x) · Y ( y) · Z ( z) .
1
Erwin Rudolf Josef Alexander Schrödinger (1887 – 1961) was a German physi-
cist and quantum theorist, the inventor of the wave equation of matter named
after him which earned him the 1933 physics Nobel (shared with Paul Dirac),
and of the eponymous unobserved cat, which finds itself in a superposition state
of being both alive and dead.
Ó »îá .
The Laplace equation in rectangular co-ordinates 33
∂2 ∂2 ∂2
YZ X +XZ Y + X Y Z = 0.
∂x 2 ∂ y2 ∂ z2
Because this has to be true for all values x, y, z, it follows that each term
must be a constant. If we take for the first and second term − k21 and
− k22 , we get in conclusion for the third constant k21 + k22 . By writing this
definition and result out and moving the denominator to the other side,
we obtain
∂2
X ( x) = − k21 X ( x) ,
∂ x2
(why the minus sign? We shall presently see. . . )
∂2
Y ( y) = − k22 Y ( y) ,
∂ y2
and
∂2
Z ( z) = k21 + k22 Z ( z) .
( )
∂z 2
Now, the solution is readily found at least to the first two equations: they
harmoninen are harmonic oscillators, and their basis solutions2 are
värähtelijä
X ( x) = exp (± ik 1 x) ,
Y ( y) = exp (± ik 2 y) .
2
Alternative basis solutions are X (x) = sin k 1 x, X (x) = cos k 1 x etc. They
are equivalent with those presented because exp (ik 1 x) = cos k 1 x + i sin k 1 x,
exp (− ik 1 x) = cos k 1 x − i sin k 1 x.
Ó »îá .
34 The Laplace equation and its solutions
The general solution is obtained by summing the terms Vk1 k2 with varying
coefficients, for different values of k 1 , k 2 .
We cannot choose the values k 1 , k 2 entirely freely. Which values are
allowed, will depend on the boundary conditions given.
Let us assume that both in the x and in the y direction the size of
our world is L (“shoebox world3 ”). Let us make things a little simpler
by assuming that, on the borders of our shoebox world, we have the
boundary conditions
V (0, y) = V (L, y) = V ( x, 0) = V ( x, L) = 0.
It then follows that the only pairs ( k 1 , k 2 ) yielding a solution that fits the
box are
πj πk
k1 = , k2 = , j, k ∈ Z,
L L
and the only suitable functions are sine functions. Thus we obtain as a
solution:
x) ( y) )z)
( ( √(
V jk ( x, y, z) = sin π j sin π k exp ±π j 2 + k2 .
L L L
This particular solution may now be generalized by multiplying it with
suitable coefficients, and summing it over different index values j =
0, ±1, ±2, . . . ; k = 0, ±1, ±2, . . .. We may however remark, that the terms
for which j = 0 or k = 0 will always vanish, and the terms that contain
j = + n and j = − n, or k = + n and k = − n, n ∈ N, are (apart from their
algebraic signs) identical. Therefore in practice we sum over the values
j = 1, 2, . . . ; k = 1, 2, . . ..
Different boundary conditions will give slightly different general solu-
tions. Their general form is however always similar.
The zero-level z = 0 expansion resulting from the general solution is the
familiar Fourier4 sine expansion:
∞ ∑
∑ ∞
V ( x, y, 0) = v jk V jk ( x, y) =
j =1 k=1
∞ ∑ ∞ [ ( [ ]) ( [ ])]
∑ jx ky
= v jk sin π sin π , (2.1)
L L
j =1 k=1
3
. . . though real-world shoeboxes are rarely square.
4
Joseph Fourier (1768 – 1830) was a French mathematician and physicist — and
some would say, climatologist —, one of the Eiffel Tower’s 72 names.
Ó »îá .
The Laplace equation in rectangular co-ordinates 35
sin 13Lπ x
Sea level
sin 5Lπ x
Ó »îá .
36 The Laplace equation and its solutions
D Table 2.1. Vertically shifting the potential field V , in the space and frequency
domains.
V (x, y, 0) F =⇒ v jk
↓ (hard) ⇓(× (easy) )
√
V (x, y, z) ⇐= F −1
v jk exp −π j 2 + k2 Lz
∂2 V 1 ∂V 1 ∂2 V
∆V = + + 2 = 0.
∂r2 r ∂ r r ∂α2
Question:
Perform on this the same kind of separation of variables as in
section 2.2, i.e., write first
V ( r, α) = R ( r ) A (α)
and then split the above equation into two equations, one for the
right-hand side function R ( r ) and one for the function A (α).
What form does the A (α) function of the general solution have?
Answer:
Substitution yields
∂2 R ( r ) A (α) ∂R ( r ) R ( r ) ∂2 A (α)
A (α) + + 2 = 0.
∂r2 r ∂r r ∂α2
r2
Multiply by the expression :
A (α) R ( r )
∂2 A(α)
r 2 ∂2 R ( r ) r ∂R ( r )
( )
∂α2
+ + = 0.
R (r) ∂r 2 R (r) ∂r A (α)
Ó »îá .
Example: The Laplace equation in polar co-ordinates 37
5 ∂V2
In fact, lim V2 → ∞ but lim = 0.
r →∞ r →∞ ∂r
Ó »îá .
38 The Laplace equation and its solutions
Pole Z
r cos φ P
Equator
r r sin φ
φ Y
Greenwich meridian
X = r cos φ cos λ,
Y = r cos φ sin λ, (2.3)
Z = r sin φ.
Ó »îá .
Spherical, geodetic, ellipsoidal co-ordinates 39
Z Ellipsoidal
P normal
h
1 − e2 N + h sin ϕ
[( ) ]
ϕ
X,Y
O
Reference (N + h) cos ϕ
ellipsoid
much used. On the Earth’s surface on the other hand, most often geodetic
— or geographical — co-ordinates ϕ, λ, h are used:
X = ( N + h) cos ϕ cos λ,
Y = ( N + h) cos ϕ sin λ, (2.4)
2
( )
Z = N + h − e N sin ϕ,
where
( ) a a2
N ϕ =√ =√ . (2.5)
1 − e2 sin2 ϕ a2 cos2 ϕ + b2 sin2 ϕ
The quantity N defined in equation 2.5 is the East-West direction, or
transversal, radius of curvature of the reference ellipsoid. In the equation,
2 2
a is the equatorial radius of the reference ellipsoid used, e2 = a a−2b is the
square of the so-called first eccentricity6 , and in equations 2.4, h is the
height of the point above the reference ellipsoid, see figure 2.3.
Converting rectangular co-ordinates into geodetic ones is easiest to do
iteratively, although the literature also offers closed formulas.
Spherical co-ordinates and geodetic, i.e., geographical, co-ordinates are
considerably different. In latitude, the difference is up to 11 minutes of
arc, or almost 20 km. This maximum is attained for latitudes ±45◦ .
In theoretical work one also uses ellipsoidal co-ordinates u and β. The
redukoitu co-ordinate β is called the reduced latitude. The relationship with rectan-
leveysaste
6
The parameter is connected to the Earth’s flattening f through the equation
e2 = 2 f − f 2 .
Ó »îá .
40 The Laplace equation and its solutions
gular co-ordinates is
√
X= u2 + E 2 cos β cos λ,
√
Y= u2 + E 2 cos β sin λ, (2.6)
Z = u sin β.
If the semi-major axis of the Earth ellipsoid is a and its semi-minor axis
b, it follows — exercise! — that E 2 = a2 − b2 .
∂2 V 2 ∂V 1 ∂2 V tan φ ∂V 1 ∂2 V
∆V = + + 2 − + = 0, (2.7)
∂r2 r ∂ r r ∂φ2 r 2 ∂φ r 2 sin2 φ ∂λ2
in which φ is the (geocentric) latitude, λ is the longitude, and r is the
distance from the origin or centre of the Earth.
We shall here not derive the solution of this equation by separation of
variables, as it is pretty complicated and can be found in ready form in the
literature (Heiskanen and Moritz, 1967, section 1 – 9). What is significant
is, that the solution looks somewhat similar to the above solution in
rectangular co-ordinates. The basis solutions of the Laplace equation are
( )
( ) n
( ) ( ) Y n φ, λ
Vn,1 φ, λ, r = r Yn φ, λ , Vn,2 φ, λ, r = , (2.8)
r n+1
where the first is again nonphysical in outer space, because, unlike the
true geopotential, these expressions grow to infinity for r → ∞.
In the above equations, the functions Yn are called surface spherical pintapallofunktiot
harmonics, whereas the functions Vn are solid spherical harmonics. The avaruus-
latter are harmonic functions everywhere in space except at the origin pallofunktiot
(2.8, rightmost equation) or at infinity (leftmost, physically unrealistic
equation).
The functions Yn are
n
∑
( ) ( )
Y n φ, λ = P nm sin φ (a nm cos mλ + b nm sin mλ) . (2.9)
m=0
Ó »îá .
Dependence on height 41
n, and so on.
7
Adrien-Marie Legendre (1752 – 1833) was a French mathematician known for
his work on number theory, statistics — he invented independently from Gauss
the method of least squares — and on elliptical functions. His name is inscribed
on the Eiffel Tower.
Ó »îá .
42 The Laplace equation and its solutions
Similar equations exist also for the functions P nm , m > 0. There are
even alternatives to choose from, though most equations are complicated.
One should be careful that in their computation, the factorials don’t go
overboard. Already 30! (factorial of 30) is a larger number than most
computers can handle as integers. . . not to mention 360!. Heiskanen
and Moritz (1967) equation 1 – 62, unlike is stated there, is not directy
suitable for computer use!
The first Legendre polynomials are listed in table 2.2.
Higher polynomials than this are rarely needed in manual computation.
Note that the even polynomials are mirror symmetric through the ori-
( ( )) ( )
gin, P n (− t) = P n ( t) — or equivalently P n sin −φ = P n sin φ — and
( ( ))
the odd ones are antisymmetric, i.e., P n (− t) = −P n ( t) or P n sin −φ =
( )
−P n sin φ .
For comparison: also the Fourier basis functions (like, in a somewhat
more complicated way, also sines and cosines!)
( x)
F j ( x) = exp 2π i j ,
L
F j+1 ( x) = F j ( x) · F1 ( x) .
Ó »îá .
Legendre’s functions 43
P6
−0.5 P4 P3
P2
P1
−1
−1 −0.5 0 −→ t 0.5 1
Ó »îá .
44 The Laplace equation and its solutions
5
P22
P21
P31 P11
−5 P32
P11
P21
−10 P22
P33 P31
P32
P33
−15
−1 −0.5 0 −→ t 0.5 1
2πR R
=π ,
2m m
where R is again the radius of the Earth.
( )
A similar formula applies also for functions P nm sin φ : as the function
passes through zero n − m times on the interval — from pole to pole —
Ó »îá .
Legendre’s functions 45
Figure 2.6. The algebraic signs of spherical harmonics on the Earth surface.
Grey means positive, white negative. The functions “wave” in a sine
D or cosine function like fashion.
Figure 2.7. Surface spherical harmonics as maps. Horizontal axis λ ∈ [0, 360◦ ),
vertical axis φ ∈ [−90◦ , 90◦ ]. Functions depicted are
( ) ( ) ( )
P50 sin φ P66 sin φ cos 6λ P11,6 sin φ cos 6λ
( ) ( ) ( ) .
D P40 sin φ P65 sin φ cos 5λ P10,6 sin φ cos 6λ
Ó »îá .
46 The Laplace equation and its solutions
D Table 2.4. Semi-wavelengths for different degrees and orders of spherical har-
monics.
10 2000 18◦
40 500 4◦ .5
180 111 1◦
360 55 0◦ .5 = 30′
1800 11 0◦ .1 = 6′
10,800 1.85 1′
or (antisymmetric case)
( ) ( ( ))
P nm sin φ = −P nm sin −φ .
P nm ( t) = P nm (− t)
or (antisymmetric case)
P nm ( t) = −P nm (− t) .
Ó »îá .
Orthogonality, orthonormality of Legendre functions 47
This dependence works though the “Fourier basis functions” cos mλ and
sin mλ. The interesting property here is rotational symmetry: does the
spherical-harmonic expansion 2.10 change when we change λ?
We see immediately that, for m ̸= 0, there will be dependence on λ if
any coefficient a nm , b nm is non-zero. So, in order to obtain rotational
symmetry, all coefficients a nm , b nm for values m > 0 must be suppressed:
a 11 = b 11 = a 21 = b 21 = a 22 = b 22 = · · · = 0.
Of the remaining coefficients, we can say that for m = 0, sin mλ = 0
identically, so the coefficients b 00 , b 10 , b 20 , . . . simply don’t matter. They
may be any value, including zero. The coefficients a 00 , a 10 , a 20 however
do matter, as for m = 0, cos mλ = 1 identically. So we obtain as the
rotationally symmetric expansion
∞
( ) ( ) ∑ 1 ( )
V φ, λ, r = V φ, r = a P sin φ ,
n+1 n n
r
n=0
in which P n = P n0 and a n = a n0 .
Ó »îá .
48 The Laplace equation and its solutions
in which ψ is the angular distance from some point on the surface of the
sphere. Equation 2.14 tells us, that Legendre polynomials are mutually
orthogonal if the vectorial product is defined as an integral over the
surface of the unit sphere σ. Alternatively we may define also fully
normalized Legendre polynomials
( ) p ( )
P n cos ψ = 2 n + 1P n cos ψ , (2.15)
Ó »îá .
Splitting a function into degree constituents 49
in which case the now modified scalar product — the mean square value
over the unit sphere — is
¨
n = n′ ,
{
1 ( ) ( ) 1 if
P n cos ψ P n′ cos ψ d σ =
4π σ 0 if n ̸= n′ ,
¨
2n + 1
f φ′ , λ′ P n cos ψ d σ′ ,
( ) ( ) ( )
fn φ, λ = (2.16)
4π σ
8
And also ˆ +1
n = n′ ,
{
1 1 if
P n (t) P n′ (t) dt =
2 −1 0 if n ̸= n′ ,
again the mean square value.
Ó »îá .
50 The Laplace equation and its solutions
and substituting this into the degree constituent equation 2.16, we obtain,
by exploiting the orthogonality of the Legendre functions, on the right-
hand side of the degree constituent equation:
¨
2n + 1
f φ′ , λ′ P n cos ψ d σ′ =
( ) ( )
IR =
4π σ
¨
2n + 1
P n2 cos ψ d σ′ =
( )
= a n0
4π σ
ˆ 2π [ ˆ 2π ]
2n + 1
P n2 ′
( )
= an cos ψ · sin ψ dλ dψ
4π 0 0
+1 ˆ [ ˆ 2π ][ ]
2n + 1 1
= an P n2 ( t) · sin ψ dλ ′
dt =
4π −1 0 sin ψ
2n + 1 2
= · 2π a n · = an,
4π 2n + 1
def
where we used the notation a n = a n0 as well as equation 2.13.
On the left-hand side of the degree constituent equation we obtain, be-
cause, on the assumed north pole, φ = 90◦ and thus sin φ = 1:
n
∑
( ) def ( ◦
)
I L = f n φ, λ = Yn 90 , λ = P nm (1) (a nm cos mλ + b nm sin mλ)
m=0
= P n0 (1) a n0 = a n ,
by using
P n0 (1) = 1,
P nm (1) = 0 if m ̸= 0.
( )
As this applies for every point φ, λ — and note that the values a n depend
on this choice! — it follows that the degree constituent equation 2.16 is
generally true.
Ó »îá .
Spectral representations of various quantities 51
where the spectral constituents Vn are the familiar Yn from equation 2.9:
n
∑
aVnm cos mλ + bVnm sin mλ =
( ) ( )( )
Vn φ, λ = P nm sin φ
m=0
∑n
( )
= vnm Ynm φ, λ ,
m=− n
and coefficients
aVnm if m ≥ 0,
{
vnm = (2.19)
bVn|m| if m < 0.
On the Earth’s surface ( r = R ) we obtain
∞
∑ n
∞ ∑
∑
( ) ( ) ( )
V φ, λ, R = Vn φ, λ = vnm Ynm φ, λ .
n=0 n=0 m=− n
Ó »îá .
52 The Laplace equation and its solutions
D 2.11.2 Gravitation
n+1
g nm = − vnm , (2.24)
R
see equation 2.19. This is an interesting result worth thinking about:
9
Carl Gottfried Neumann (1832 – 1925) was a German mathematician.
Ó »îá .
Low-degree spherical harmonics 53
1. Firstly, note how simple the connection 2.24 between potential vnm
and gravitation g nm has become in the frequency domain!
2. Secondly, if there are available, over the whole surface area of the
Earth, measurement values of gravitational acceleration g, we may
derive from these the spherical-harmonic coefficients g nm and the
( )
degree constituent functions g n φ, λ using the method explained
earlier. In this way we can then obtain the solution by means of
equation 2.23 for the whole exterior geopotential field! This is the
basic idea of geopotential — or geoid — determination, from the
spectral perspective.
a 11 i + b 11 j + a 10 k = G d,
Ó »îá .
54 The Laplace equation and its solutions
Figure 2.8. Monopole, dipole, and quadrupole, at the Earth’s centre, and their
D effects on the geoid.
In that case we may compute the dipole moment of the whole Earth by
integration:
˚ ˚
d⊕ = r dm = ρ r d V = M⊕ · rcentre of mass ,
Earth Earth
10
Note that here is used a, the equatorial radius of the Earth’s reference ellipsoid,
( )
not R, and φ, the geocentric latitude. The co-ordinates r, φ, λ form a spherical
Ó »îá .
Often used spherical-harmonic expansions 55
360 ( ) ∑n
( )
GM⊕ ∑ a n ( )[ ]
V= 1+ P nm sin φ C nm cos mλ + S nm sin mλ . (2.25)
r r
n=2 m=0
def
Then we use the notation Jn = Jn0 , and J2 is the most important, the
parameter of the Earth’s gravity field describing the flattening of the
Earth. The relationship with the parameters C, S is
p
{ } { }
Jn0 C n0
= − 2n + 1 ,
K n0 S n0
co-ordinate system.
Ó »îá .
56 The Laplace equation and its solutions
{ } √ { }
Jnm ( n − m)! C nm
=− 2 (2 n + 1) , m ̸= 0.
K nm ( n + m)! S nm
Ó »îá .
Ellipsoidal harmonics 57
in which Q nm (z) are the so-called Legendre functions of the second kind,
sampled in table 2.7. Though the general argument z is complex, equation
2.27 gives a real result for real-valued coefficients aenm , benm .
Those interested in the derivation of the above equation can find it in
Heiskanen and Moritz (1967) or other textbooks on potential theory.
Heiskanen and Moritz give a slightly different form to the equation, the
auxiliary equations needed for the normalization used here can be found
on their pages 66 – 67.
Assume ae10 = ae11 = be11 = 0, i.e., the vanishing of the dipole moment12 .
We can also show that in the expansion 2.27 the first coefficient is
e GM⊕ E
a 00 = arctan
E b
and the expansion specialized for a rotationally symmetric field becomes
∞ ∞
Q n i Eu e
( )
( ) ∑ ( ) ∑ ( )
V u, β = Vn u, β = ( b ) a n0 P n sin β ,
n=0 n=0
Qn i E
yields now
Q 0 i Eu GM⊕
( )
E
V0 ( u) = ( b) arctan ,
Q0 i E E b
the gravitational potential of the field constituent of ellipsoidal degree
zero.
With the substitutions (Heiskanen and Moritz, 1967, page 66)
( u) ( )
E b E
Q0 i = − i arctan , Q 0 i = − i arctan (2.28)
E u E b
we obtain
GM⊕ E
V0 ( u) = arctan .
E u
11
This expansion for the ellipsoid of revolution differs from the expansion into
Lamé functions found for the triaxial ellipsoid.
12
Strictly speaking this works only in case of a spherical-harmonic expansion, or
in the limit E → 0.
Ó »îá .
58 The Laplace equation and its solutions
( ) GM⊕ E
V u, β, λ = arctan ·
[ E∞ n u
( e )]
∑ ∑ arctan E Q nm i u
( )
b C nm cos m λ+
( Eb ) P nm sin β
( )
· 1+ e , (2.29)
arctan E
u Q nm i + S nm sin m λ
n=2 m=0 E
Q nm i Eu
( ) ( a )n+1
lim ( b
)=
E →0 Q nm iE r
Ó »îá .
Self-test questions 59
1 z+1
Q 0 (z) = ln (n + 1)Q n+1 (z) − (2n + 1) zQ n (z)+
2 z−1
z z+1 + nQ n−1 (z) = 0
Q 1 (z) = ln −1
2 z−1
3z2 − 1 z + 1 3z
Q 2 (z) = ln − )m/2 d m
Q mn (z) = 1 − z2
(
4 z−1 2 Q n (z)
3
5z − 3x z + 1 5z2 2 dzm
Q 3 (z) = ln − +
4 z−1 2 3
D Self-test questions
Ó »îá .
60 The Laplace equation and its solutions
If
n
∞ ( ) n+1 ∑
( ) ∑ R ( )
V φ, λ, r = P nm sin φ (a nm sin mλ + b nm cos mλ) =
r
n=0 m=0
∞ ( )n+1
∑ R ( )
= Vn φ, λ, R , (2.30)
r
n=0
Ó »îá .
Exercise 2 – 2: Symmetries of spherical harmonics 61
3. For what n value will they be less than 10−4 × of what they are on
the Earth’s surface?
1 d n [( 2 )n ]
P n ( t) = n n
t −1 .
2 n! dt
1. Note first, that
(a) Differentiating a symmetric function of t will produce an anti-
symmetric function, and vice versa.
(b) The function t2 − 1 and its powers are symmetric.
( )
Ó »îá .
62 The Laplace equation and its solutions
( )
(Hint: see the example formulas and graphs for P nm sin φ in the
lecture notes and try to guess a general rule. Then, verify.)
3. The same question if the potential is rotationally symmetric about
( ) ( )
the Earth rotation axis, i.e., V φ, λ, r = V φ, r .
Ó »îá .
D The normal gravity field
3
−
→
γ
−
→
γ
−
→
γ Equipotential surfaces of
the normal gravity field
Reference ellipsoid
(flattening exaggerated)
– 63 –
64 The normal gravity field
∂U
γ ( x, y, z) = − ,
∂n
U = Ψ + Φ,
p = X i + Y j,
Ó »îá .
The centrifugal force and its potential 65
p
Centrifugal force
k
Gravitation
Gravity
i j
Y
X
We may also derive from the centrifugal potential Φ the following equa-
tion by differentiating it twice:
∂ ∂
∆Φ = ∇2 Φ = ∇f = ω2⊕ X + ω2⊕ Y + 0 = 2ω2⊕ , (3.1)
∂X ∂Y
1
Gaspard-Gustave Coriolis (1792 – 1843) was a French mathematician, physicist
and mechanical engineer. His name is inscribed on the Eiffel Tower.
Ó »îá .
66 The normal gravity field
W ( x, y, z) = W0 = const.
Ó »îá .
Level surfaces and plumb lines 67
tangent plane
P
W = WP − δW
x0 W = WP
ρ1
Level surfaces and gravity vectors, or plumb lines, are always
perpendicular to each other.
Given in point P a plane that in P has the same direction as the level
surface, i.e., its tangent plane. If the local curvature of the level surface
in the x direction is ρ 1 , and the x co-ordinate of point P is x0 , we may
develop the distance between the surfaces in a Taylor series:
1
ϵ≈ ( x − x0 )2 .
2ρ 1
Ó »îá .
68 The normal gravity field
def 1 Wxx
κ1 = =− , (3.3)
ρ1 g
def 1 Wyy
κ2 = =− , (3.4)
ρ2 g
we obtain
−2 gJ + Wzz = −4πG ρ + 2ω2⊕ .
By using
∂g ∂g
Wzz = − =−
∂z ∂H
— in which H is the height co-ordinate — we obtain for the vertical
gradient of gravity (Heiskanen and Moritz, 1967, equation 2 – 20):
∂g
= −2 gJ + 4πG ρ − 2ω2⊕ ,
∂H
an equation found by Ernst Heinrich Bruns.
3
Marie-Sophie Germain (1776 – 1831) was a brilliant French mathematician,
number theorist and student of elasticity. She corresponded with Gauss, among
others, on number theory, and did foundational work toward a proof of Fermat’s
last theorem. Her name is missing from the Eiffel Tower.
Ó »îá .
Natural co-ordinates 69
Astronomical
co-ordinates Φ, Λ
Plumb line
Greenwich
O Φ
where the integral is taken along the plumb line of point P . ∂∂H = ∂∂n is
the derivative in the direction of the plumb line, i.e., the local normal
to the level surfaces. g is the acceleration of gravity along the plumb
line as a function of place — or of geopotential level. In this case of
orthometric heights, g is the true gravity inside the rock, which is a
non-linear function of place and will also depend on rock density. This
trickiness of their determination is a problem specific to orthometric
heights. We will return to this later on (Heiskanen and Moritz, 1967
chapter 4).
Also the co-ordinates Φ, Λ and H form a natural co-ordinate system.
Ó »îá .
70 The normal gravity field
Now
U u, β = Ψ u, β + Φ u, β .
( ) ( ) ( )
( )
On the reference ellipsoid u = b we have as a requirement U b, β = U0 ,
which is possible only if
1
U0 = A 0 + ω2⊕ b2 + E 2 ,
( )
3
0 = A1,
1
0 = A 2 − ω2⊕ b2 + E 2 ,
( )
3
0 = A n , n = 3, 4, 5, . . .
GM⊕ E 1
U0 = arctan + ω2⊕ a2 .
E b 3
1 GM⊕ E
A 0 = U0 − ω2⊕ a2 = arctan .
3 E b
The normal gravity potential U is obtained as follows (remember the
Ó »îá .
Normal gravity on the reference ellipsoid 71
U ( u, β) = Ψ u, β + Φ u, β =
( ) ( )
Q 2 i Eu ( 3
( )
GM⊕ E 1 1)
= arctan + ω2⊕ a2 ( b) sin2 β − +
E u 3 Q 2 i E 2 2
Ψ0 (u) A2 P2 (sin β)
1
+ ω2⊕ u2 + E 2 cos2 β =
( )
2
Φ( u,β)
= C 1 ( u) sin2 β + C 2 ( u) cos2 β,
and
Z
sin ϕ (1− e2 ) N 1 Z a2
tan ϕ = = p 2 2 = p = tan φ,
cos ϕ X +Y 1 − e2 X 2 + Y 2 b2
N
in which φ is the geocentric latitude, see equation 2.3. From this follows
directly
b
tan β = tan ϕ,
a
Ó »îá .
72 The normal gravity field
.P
b
β
φ ϕ . X /Y
a O
D Figure 3.5. The geometry of the meridian ellipse and various types of latitude.
4
Carlo Somigliana (1860 – 1955) was an Italian mathematician and physicist.
Paolo Pizzetti (1860 – 1918) was an Italian geodesist.
Ó »îá .
Numerical values and formulas 73
a = 6,378,137 m,
1
= 298.257,222,101,
f
ω⊕ = 7,292,115 · 10−11 s−1 ,
GM = 3,986,005 · 108 m3 s−2 .
U = 62,636,860.8500 +
−9.780,326,77 − 0.051,630,75 sin2 ϕ −
[ ]
+ h+
−0.000,227,61 sin4 ϕ − 0.000,001,23 sin6 ϕ
0.015,438,99 · 10−4 − 0.000,021,95 · 10−4 sin2 ϕ −
[ ]
+ h2 −
− 0.000,000,10 · 10−4 sin4 ϕ
− −0.000,024,22 · 10−8 + 0.000,000,07 · 10−8 sin2 ϕ h3 ,
[ ]
and normal gravity (note the minus sign, U is positive and diminishes
going upward):
∂U
γ= − =
∂h
= + 9.780,326,77 + 0.051,630,75 sin2 ϕ +
+ 0.000,227,61 sin4 ϕ + 0.000,001,23 sin6 ϕ −
0.030,877,98 · 10−4 − 0.000,043,90 · 10−4 sin2 ϕ −
[ ]
− h−
− 0.000,000,200 · 10−4 sin4 ϕ
− −0.000,072,65 · 10−8 + 0.000,000,21 · 10−8 sin2 ϕ h2 .
[ ]
(3.7)
Here, the unit of potential is m2 /s2 , and the unit of gravity, m /s2 . More
precise equations can be found from Heikkinen (1981). In these equa-
tions, the coefficient 9.780,32 . . . m /s2 is equatorial gravity, and the value
0.030,87 . . . s−2 is the vertical gradient of gravity on the equator. ϕ is
geodetic latitude, h (in metres) is the height above the reference ellipsoid.
Ó »îá .
74 The normal gravity field
Other gravity formulas and reference ellipsoids still in legacy use (and
slowly vanishing) are Helmert’s 1906 ellipsoid, the Krasovsky ellipsoid or
SK-42 in the countries of Eastern Europe, the International or Hayford
ellipsoid (1924) and its gravity formula, and Geodetic Reference System
1967.
According to the above equation, the normal potential over the equator is
Questions:
Answers:
1. See figure 3.6. The minima are at heights 3000 km and 2000 km,
approximately.
2. Not very physical: the stationary point for potential U (the normal
potential in a co-rotating reference system) should be located at
approx. 36,000 km height, at the geostationary orbit.
This tells us that polynomial approximation cannot be extrapolated
very far. In this case the interval of extrapolation is of the same
order as the radius of the Earth, and that won’t work any more.
then we may also write the normal gravitational potential, Ψ, into the
form [ ]
∞ ( a )2n
GM⊕ ∑
Ψ=
( )
1− J2n P2n sin φ ,
r r
n=1
Ó »îá .
The normal potential as a spherical harmonic expansion 75
.
140,000,000
Cubic
120,000,000 Quadratic
Linear
100,000,000 Realistic
80,000,000
60,000,000
40,000,000
20,000,000
Figure 3.6. The normal field’s potential curve over the equator. Heights in
D kilometres.
which contains only even coefficients J2n = J2n,0 as the normal field is
symmetric about the equatorial plane.
The coefficients for the GRS80 normal gravitational potential are
found5 in table 3.1. Higher terms are usually not needed. The
relationship between fully normalized and non-normalized coefficients is
p
J n = J n 2 n + 1.
Note for comparison that in the expansion of the same field into ellipsoidal
harmonics only the degree zero and degree two coefficients are non-zero!
5
They can also be calculated from equation (2 – 97) given in Heiskanen and
Moritz (1967) :
( )n
3 e2
( )
n+1 J2
J2n = (−1) 1 − n + 5n 2 ,
(2n + 1) (2n + 3) e
starting from the values J2 and e2 . The results are the same as in the table’s
left column.
Ó »îá .
76 The normal gravity field
This is the main reason why these functions are used at all.
Instead of using an ellipsoidal model, we may use as a normal gravity
potential formula also the first two, three terms of the spherical-harmonic
expansion of the real geopotential. Then we obtain, taking the centrifugal
potential along:
Y0 Y2 (φ, λ) 1 2 2
U= + + ω⊕ ( X + Y 2 ),
r r3 2
with the corresponding equipotential surface being the “Bruns spheroid”,
or
Y0 Y2 (φ, λ) Y4 (φ, λ) 1 2 2
U= + + + ω⊕ ( X + Y 2 ),
r r3 r5 2
def ( ) ( )
the “Helmert spheroid”. Here, Y0 = GM and Y2 φ, λ , Y4 φ, λ are taken
from the true geopotential.
These equations are easy to compute, but their equipotential or level
surfaces are not ellipsoids of revolution, and in fact not even rotationally
symmetric. The are, in fact, quite complicated surfaces (Heiskanen and
Moritz, 1967, 2 – 12)!
However, in geometric geodesy we always use a reference ellipsoid, so
this is also a wise thing to do in physical geodesy.
W = V + Φ,
in which Φ is the centrifugal potential (see above), and the normal poten-
tial
U = Ψ + Φ.
The difference between them is
def
T = W − U = V − Ψ,
W = V +Φ =
n
{ ∞ ( ) ∑
}
GM⊕ ∑ a n
= Φ+
( )
1− P nm sin φ [ Jnm cos mλ + K nm sin mλ] ,
r r
n=2 m=0
Ó »îá .
The disturbing potential 77
T = W −U =
n
{∞ }
GM⊕ ∑ ( a )n ∑ ( )
=− P nm sin φ [δ Jnm cos mλ + δK nm sin mλ] ,
r r
n=2 m=0
in which {
δ Jn0 = Jn0 − Jn∗ if n even,
δ Jnm = Jnm if n odd,
and δK nm = K nm .
where, in every term, the degree constituent T n has the same dimension
as T , and
n
( ) GM⊕ ∑ ( )
Tn φ, λ = − P nm sin φ [δ Jnm cos mλ + δK nm sin mλ] .
a
m=0
1. the total mass of the Earth GM⊕ assumed by the normal field is
realistic, and
2. the origin of the co-ordinate reference system is assumed to be at
the centre of mass of the Earth.
6
Earlier on we have used for this reference radius (in spherical approximation)
also the symbol R.
Ó »îá .
78 The normal gravity field
D Self-test questions
1
Φ = ω 2 ( X 2 + Y 2 ),
2
derive the centrifugal acceleration as a vector. X , Y , Z are rectan-
gular co-ordinates of a frame rotating at angular rate ω around the
Z axis.
4. Explain the idea of natural co-ordinates.
5. What relationhip was there between Mr. Leblanc and C.F. Gauss?
Use Google.
6. What makes the Somigliana–Pizzetti equation (3.6) valuable?
7. What are the defining parameters of the Geodetic Reference System
1980?
8. Why does the spherical-harmonic expansion of the normal potential
contain only a small number of terms and coefficients?
9. Why doesn’t the spherical-harmonic expansion of the normal poten-
tial contain any terms with order m ̸= 0?
10. Why does the spherical-harmonic expansion of the normal potential
contain only terms with even degree numbers n?
Given is the rotation rate of the Earth, radians per second: ω⊕ = 7292115 ·
10−11 s−1 .
Ó »îá .
Exercise 3 – 2: The centrifugal force 79
Ó »îá .
D Anomalous quantities of the gravity field
4
All other anomalous quantities are various functions of this, like the
geoid height N and the plumb-line deflections ξ, η. They are generally
obtained by subtracting from each other
ξ = Φ − ϕ,
η = (Λ − λ) cos ϕ,
in which (Φ, Λ) are astronomical latitude and longitude, i.e., the direction
( )
of the local plumb line, and the geodetic latitude and longitude ϕ, λ
form correspondingly the direction of the surface normal on the reference
ellipsoid. See figure 4.1.
The geoid height or geoid undulation is
N = H − h,
– 81 –
82 Anomalous quantities of the gravity field
Plumb-line
deflection ξ, η Topography
Reference ellipsoid
Geoid Geoid height N
1
This is not self-evident! In a local vertical datum the potential of the zero point
could well differ by even as much as a metre from the normal potential of a
global reference ellipsoid.
Ó »îá .
Disturbing potential, geoid height. . . 83
Figure 4.2. A geoid model for Finland from 1984. Deflections of the plumb line
D from observations in red (Vermeer, 1984).
of the height system used. It follows that also the first term vanishes (it
would in any case always be small, except in the mountains). So
ˆ UP
1 1 TP T
NP = − dU ≈ (WP − UP ) = or N= . (4.2)
WP γ γP γP γ
2
Ernst Heinrich Bruns (1848 – 1919) was an eminent German mathematician
and mathematical geodesist.
Ó »îá .
84 Anomalous quantities of the gravity field
−
→
γ = gradU
g = gradW
P
WP Geoid
UP
N Ellipsoid
UQ (= WP )
Q
Figure 4.3. Equipotential surfaces of the gravity field (W) and the normal
D gravity field (U).
Figure 4.3 depicts the situation still better. In this figure, the gradient
−−→ −−→
vectors g = grad W and −→
γ = grad U have lengths ∂∂W ∂U
H and ∂ H , from which it
follows, with equation T = W − U , that the separation between “matching”
surfaces W = WP and U = UQ , when WP = UQ , is
UQ − UP WP − UP T
N= = = .
γ γ γ
∂W ∂U
( )
δg = − − ;
∂H ∂h
where differentiation takes place along the plumb line for W , and along
the ellipsoidal normal for U . The directions of plumb line and ellipsoidal
surface normal are actually very close to each other.
In spherical approximation we have
∂W ∂U ∂T
( )
δg = − − =− .
∂r ∂r ∂r
Ó »îá .
Gravity anomalies 85
using GPS, but traditionally it has been impossible. For this reason
gravity disturbances are little used. One rather uses gravity anomalies,
about which more below.
3
Note that, if the normal field that defines T has a reference ellipsoid centred
on the Earth’s centre of mass, and the implied total mass GM of the normal
( )
field equals that of the true Earth, then the first two constituents T0 φ, λ =
( ) ∑∞
T1 φ, λ = 0, and the sum may be taken as n=2 . This was assumed here. See
section 2.12.
Ó »îá .
86 Anomalous quantities of the gravity field
Topography
Telluroid P = observation point
ζ
Mean sea Q
level (geoid)
hQ HP
Ellipsoid
N
Figure 4.4. Reference ellipsoid, mean sea level (geoid), and gravity measure-
D ment.
however expects a height above the reference ellipsoid! This special trait
of the definition of gravity anomalies may be called a“free boundary-value
problem”.
According to this we calculate gravity anomalies as follows:
∆ g P = g P − γQ = g P − γP + γP − γQ =
( ) ( )
∂W ⏐⏐ ∂U ⏐⏐ ∂U ⏐⏐ ∂U ⏐⏐
( ⏐ ⏐ ) ( ⏐ ⏐ )
=− − − −
∂ H ⏐P ∂ h ⏐P ∂ h ⏐P ∂ h ⏐Q
∂ (W − U ) ⏐⏐ ) ∂γ ⏐
⏐ ⏐
(
≈− ⏐ + h P − hQ ∂ H ⏐ =
⏐
∂H P P
∂T ⏐⏐ ∂γ ⏐⏐
⏐ ⏐
=− + (h P − HP ) =
∂ H ⏐P ∂ H ⏐P
∂T ⏐⏐ ∂γ ⏐⏐
⏐ ⏐
=− + NP =
∂ H ⏐P ∂ H ⏐P
∂T T ∂γ ⏐⏐
( )⏐
− + ,
∂H γ ∂H ⏐ P
∂T 2
∆g = − − T, (4.5)
∂r r
in which r = R + H is the distance from the Earth’s centre.
Ó »îá .
Units used for gravity anomalies 87
def n−1
∆gn = Tn. (4.7)
R
The presence of the factor n − 1 shows that gravity anomalies cannot
contain n = 1 constituents, even if T would contain them. It is always
wise to choose the origin of the co-ordinate system to be in the centre
of mass of the Earth, but if it is not, at least gravity anomalies do not
change.
Equation 4.7 applies only on a spherical Earth of radius R . In the
space outside the Earth we obtain, using equations 4.3 and 4.5, the
corresponding equation
∞ ( )n+1
1∑ R
∆g = [( n + 1) − 2] Tn =
r r
n=2
∞ ( )n+1
1∑ R
= ( n − 1) Tn =
r r
n=2
∞ ( ) n+2
∑ R
= ∆ gn. (4.8)
r
n=2
Ó »îá .
88 Anomalous quantities of the gravity field
∂T 2
∆g = − − T. (4.9)
∂n R
4
The third or mixed boundary-value problem is associated with Victor Gustave
Robin (1855 – 1897), a French mathematician. Then, the Dirichlet problem could
be called the first, the Neumann problem the second boundary-value problem.
Ó »îá .
The boundary-value problem of physical geodesy 89
5
In fact, the atmosphere complicates this matter.
6
Note, however, that GRS80 has an equatorial radius of 6,378,137.0 m, while
the newer models like EGM2008 give a smaller value of 6,378,136.3 m as the
location of global mean sea level. Uncertainty continues to be in the decimetres.
) def
7
We write ∆ g φ′ , λ′ = ∆ g φ′ , λ′ , R , the integral being evaluated on the Earth
( ( )
sphere r = R.
Ó »îá .
90 Anomalous quantities of the gravity field
In section 7.1 we will give a closed formula for this function (for the case
r = R ) and a graph.
◦ the true potential field of the Earth is the normal potential, and
◦ the mathematical mean sea surface or geoid, the reference surface
for height measurement, coincides with the reference ellipsoid.
In other words, the telluroid is a model for the Earth’s topographic surface
that is obtained by taking levelled heights — more precisely, geopotential
numbers obtained from levelling — as if they represented differences in
normal potential with that of the reference ellipsoid.
8
Sir George Gabriel Stokes PRS (1819 – 1903) was an Irish-born, gifted mathe-
matician and physicist making his career in Cambridge.
Ó »îá .
Free-air anomalies 91
9
For greatest precision one should consider that also the latitude Φ may not be
a latitude on a geocentric reference ellipsoid, but, e.g., astronomical latitude, or
latitude in some old national co-ordinate system computed on a non-geocentric
ellipsoid, like in Finland KKJ, the National Grid Co-ordinate System, which was
computed on the Hayford ellipsoid. The error caused by this is however of order
a thousand times smaller than the effect caused by H − h.
Ó »îá .
92 Anomalous quantities of the gravity field
Questions:
Answers:
105 × 9.8
1. At −0.3 mGal /m, it takes m = 3267 km to go to zero.
0. 3
2. Not very. The gravity gradient itself drops quickly from the value of
−0.3 mGal /m going up, so this linear extrapolation is simply wrong.
D Self-test questions
Ó »îá .
Exercise 4 – 1: The spectrum of gravity anomalies 93
18 20 22 24 26 28 30 32
64 64
62 62
60 60
18 20 22 24 26 28 30 32
Figure 4.5. Free-air gravity anomalies for Southern Finland, computed from
the EGM2008 spherical-harmonic expansion. Data © Bureau
Gravimétrique International (BGI) / International Association of
Geodesy. Extraction address:
http://bgi.omp.obs-mip.fr/data-products/Toolbox/
D EGM2008-anomaly-maps-visualization.
Use equation 4.7. If we assume that the mean magnitude of the spectral
components ∆ g n of gravity anomalies
√ ¨
def 1
∆gn = ∆ g2n φ, λ d σ
( )
4π σ
does not depend on the chosen degree number n, how then does the
similar T n depend on n?
In other words: which degree numbers of the gravity field are relatively
strongest in the disturbing potential, and which in the gravity anomalies?
Ó »îá .
94 Anomalous quantities of the gravity field
T
N= ,
γ
Ó »îá .
D
5 Geophysical reductions
D 5.1 General
We see that integral equations, like Green’s third theorem 1.24, offer a
possibility to calculate the whole exterior potential of the Earth — as
well as all quantities that may be calculated from the potential, like the
acceleration of gravity, etc. — from observed values V and ∂∂Vn on the
boundary surface only. Green’s third theorem is but one example out of
many: every integral theorem is the solution of some boundary-value
problem.
There are three alternatives concerning the choice of boundary surface:
1. Choose the topographic surface of the Earth.
2. Choose mean sea level, more precisely, a equipotential surface close
to mean sea level called the geoid.
3. Choose the reference ellipsoid.
1
This influence is called the “indirect effect”.
– 95 –
96 Geophysical reductions
D 5.2.1 Calculation
ℓ
ds dz = d β dz,
cos β
2
Pierre Bouguer (1698 – 1758) was a French professor in hydrography, who par-
ticipated in the public discussion on the figure of the Earth, and in 1735 – 1743
led an expedition of the French Academy of Sciences doing a grade measurement
in Peru, South America, at the same time when De Maupertuis executed a
similar grade measurement in Lapland. In addition to geodesy, he was also
active in astronomy.
Ó »îá .
Bouguer anomalies 97
z
P
dz
ds ℓ
β cos β d β
H ℓ
d
s dm x
Ó »îá .
98 Geophysical reductions
Evaluation point P
II
Bouguer plate
I I
d=H
Topography
If we define √
def
ℓ ( z) = ( H − z)2 + r 2
we obtain for the whole integral
2πG ρ [ d + ℓ ( d ) − ℓ (0)] .
where we assume for the density ρ of the plate an often used value for
the average density of the Earth’s crust, ρ = 2670 kg /m³. By substituting
into this equation 4.12 we obtain
Ó »îá .
Terrain effect and terrain correction 99
Topography
Free-air anomaly
Geoid
Bouguer anomaly
D 5.2.2 Properties
Bouguer anomalies are, unlike free-air anomalies that vary on both sides
of zero, especially in the mountains strongly negative. E.g., if the mean
elevation of a mountain range is H = 1000 m, the Bouguer anomalies
systematiikka will, as a consequence of this, contain a bias of 1000 × (−0.1119 mGal) =
−112 mGal, about −100 mGal for every kilometre of elevation.
The advantage of Bouguer anomalies is their smaller variation with place.
For this reason they are suited especially for interpolation and predic-
tion of gravity anomalies, in situations where the available gravimetric
material is sparse. However one then has to have access to topographic
heights.
Both errors work in the same direction! Because volumes I are below
the point of evaluation, their attraction — which the simple Bouguer
reduction considers present, and removes — would act downward. And
because volumes II are above the point of evaluation, their attraction —
which in the simple Bouguer reduction is not corrected for — acts upward.
The error made is in the same direction as in the previous case.
The terrain correction is always positive!
Ó »îá .
100 Geophysical reductions
18 20 22 24 26 28 30 32
64 64
62 62
60 60
18 20 22 24 26 28 30 32
We write
∆ g′B = ∆ g B + TC,
where TC — the “terrain correction” — is positive. ∆ g′B is called the
terrain corrected Bouguer anomaly.
The terrain correction is calculated by numerical integration. Figure
5.5 shows the prism integration method, and how both prisms, I and II,
lead to a positive correction, because prism I is computationally added
and prism II removed when applying the terrain correction. One needs
a digital terrain model, DTM, which must be, especially around the
evaluation point, extremely dense: according to experience, 500 m is
the maximum inter-point separation in a country like Finland, in the
mountains one needs even 50 m. The systematic nature of the terrain
correction makes a too sparse terrain model cause, possibly serious, biases
in the insufficiently corrected gravity anomalies.
Ó »îá .
Terrain effect and terrain correction 101
( )
H x′ , y′
Topography
θ II
H (x, y) P
I
( )
H x ′ , y′
Geoid
x, y x′ , y′
D Figure 5.5. Calculating the classical terrain correction by the prism method.
For computing the terrain correction with the prism method we use the
following equation, assuming a constant crustal density ρ and a flat
Earth, in rectangular map co-ordinates x, y:
ˆ +Dˆ +D [
1 ]2
H x′ , y′ − H ( x, y) ℓ−3 dx′ d y′ ,
( )
TC ( x, y) = G ρ
2 −D −D
where √
{1 }2
ℓ= ( x′ − x)2 + ( y′ − y) +2
[ H ( x′ , y′ ) − H ( x, y)]
2
[ ]T
is the distance between the evaluation point x y H ( x, y) and the
[ 1
[ ( )] ]T
centre point of the prism x y 2 H ( x, y) + H x′ , y . Of course ′
Given the special terrain shape rendered in quasi 3D in figure 5.7. Here,
the height differences are PQ ′ = 300 m and QQ ′ = 200 m. Rock density is
the standard crustal density, 2670 kg /m³.
Questions:
Ó »îá .
102 Geophysical reductions
Given: gravity g
on terrain Terrain Bouguer
correction plate
correction
Free-air
reduction Subtract normal gravity
to sea at sea level:
level −γ0 (ϕ)
Figure 5.6. The steps in calculating the Bouguer anomaly. Note that the reduc-
tion to sea level uses the standard free-air vertical gravity gradient,
D −0.3084 mGal /m, i.e., the vertical gradient of normal gravity.
Answers:
Q′ Sea level
Figure 5.7. A special terrain shape. The vertical rock wall at PQ is also straight
D on a map and extends to infinity in both directions.
Ó »îá .
Spherical Bouguer anomalies 103
TC = 5.595 mGal,
∆ g FA +TC −Plate ∆ gB
∆ g B +Plate −TC ∆ g FA
3
One can do so, and often does, also in connection with the Bouguer plate
correction.
Ó »îá .
104 Geophysical reductions
Just for fun, we compute the net mass effect of doing the complete spher-
ical Bouguer reduction globally. The mean height of the land topogra-
phy is 800 m, land occupying 29% of the globe. The mean ocean depth
is 3700 m, corresponding to an equivalent rock depth to be “filled in”
of 3700 × 2. 67 −1. 03
2. 67 m = 2272 m, assuming a density for crustal rock of
2670 kg /m3 and a sea-water density of 1030 kg /m3 , and ocean occupying 71%
of the globe. The sum weighted by area is thus
4
Friedrich Robert Helmert (1843 – 1917) was an eminent German geodesist
known for his work on mathematical and statistical geodesy.
Ó »îá .
Helmert condensation 105
Equipotential surface
Topography
g
g′
Condensation layer
D Figure 5.8. Helmert condensation and the changes it causes in the gravity field.
Ó »îá .
106 Geophysical reductions
densation over Bouguer reduction is, that no mass is being removed. The
Bouguer reduction amounts to the computational removal of topographic
masses on a large scale. Therefore, unlike with Bouguer reduction, in
Helmert condensation gravity anomalies will not change systematically.
In appendix D we derive series expansions in spherical geometry, which
describe both the external and the internal potential as functions of the
( )
“degree constituents” of the various powers of the topography H φ, λ .
The extensively presented derivation in the appendix is much used in
gravity field theory to model the gravity effect of the topography. In this
theory, issues of convergence are difficult, though we gloss over those
here.
D 5.6 Isostasy
D 5.6.1 Classical hypotheses
Already in the 18th and 19th centuries, e.g., thanks to Bouguer’s work
in South America, as well as that of British geodesists in the Indian
Himalayas, it was understood that mountain ranges weren’t just piles
of rock on top of the Earth’s crust. The gravity field surrounding the
mountains, specifically the plumb line deflections, could only be explained
by assuming that under every mountain range there was also a “root”
made from lighter rock species. The origin of this root was speculated to
be the almost hydrostatic behaviour of the Earth’s crust over geological
time scales. This assumption of hydrostatic equilibrium was called the
hypothesis of isostasy, also isostatic compensation.
Back then, unlike now, it was not yet possible to get a precise or even cor-
rect picture using physical methods (seismology) of how these mountain
roots were really shaped. That’s why simplified working hypotheses were
formulated.
One older isostatic hypothesis is the Pratt–Hayford hypothesis. This was
proposed by J. H. Pratt5 in the middle of the 19th century (Pratt, 1855,
1858, 1864), and J.F. Hayford6 developed the mathematical tools needed
for computation. According to this hypothesis, the density of the “root”
under a mountain would vary with the height of the mountain, so that
under the highest mountains would be the lightest material, and the
boundary between this light root material and the denser Earth mantle
5
John Henry Pratt (1809 – 1871) was a British clergyman and mathematician
who worked as the archdeacon of Kolkata, India.
https://en.wikipedia.org/wiki/John_Pratt_(Archdeacon_of_Calcutta).
6
John Fillmore Hayford (1868 – 1925) was a United States geodesist who studied
isostasy and the figure of the Earth.
Ó »îá .
Isostasy 107
Plumb-line deflections
Mountain Geoid
Earth’s crust
“Root”
Earth’s mantle
D Figure 5.10. Isostasy and the bending of plumb lines towards the mountain.
7
George Biddell Airy PRS (1801 – 1892) was an English mathematician and
astronomer, “Astronomer Royal” 1835 – 1881.
8
Veikko Aleksanteri Heiskanen (1895 – 1971), ““the great Heiskanen” (http:
//en.wikipedia.org/wiki/Beyond_Sleep) was an eminent Finnish geodesist who
also worked in Ohio, USA He is known for his work on isostasy and global geoid
modelling. See Kakkuri (2008).
Mountains
Sea
Compen-
sation Earth’s crust
depth
Compen-
sation
level
Mantle
Ó »îá .
108 Geophysical reductions
Mountains
Sea σw
Earth’s crust
t0
σc
Anti-root
Mountain root
σm
Mantle
mass density of the “root” is fixed, and that the isostatic compensation
is realized by varying the depth to which the root extends down into the
Earth’s mantle. In our current understanding this corresponds better to
what is really happening inside the Earth. This hypothesis is illustrated
in figure 5.12.
Airy’s isostatic hypothesis assumes that in every place the total mass of a
column of matter is the same. So, let the density of the Earth crust be ρ c ,
the density of the mantle ρ m , and the density of sea water ρ w ; sea depth
d , crustal thickness t and topographic height H . We have
( )
d ρm − ρw + c
tρ c + d ρ w − ( t + d ) ρ m = c =⇒ t = −
ρm − ρc
on the sea, and
Hρm − c
tρ c − ( t − H ) ρ m = c =⇒ t =
ρm − ρc
9
Its dimension is pressure: according to Archimedes’ law, the pressure of the
crustal (plus sea-water) column minus the pressure of the column of displaced
mantle material.
Ó »îá .
Isostasy 109
t r
t0 t
r
Anti-root
Root
Note that the constant c is arbitrary and expresses the fact that the
level from which one computes the depth of the root — less precisely, the
“average thickness of the crust” — can be chosen arbitrarily.
Another approach: instead of c, use the “zero topography compensation
level” t 0 , to be computed from the above equations by setting H = d = 0:
( )
t 0 ρ c − ρ m = c.
In other words,
∑
(deviation × density contrast) = const.
interfaces
Ó »îá .
110 Geophysical reductions
1. What is the depth of the root under the Hardanger plateau, relative
to the compensation depth t 0 ?
2. What is the negative depth of the anti-root under the Norwegian
Sea, relative to the same compensation depth?
3. What is the relative depth of the root of the Hardanger plateau,
compared to the nearby Norwegian Sea?
Answers:
Ó »îá .
Isostasy 111
Mid-Atlantic ridge
Plate motion
Deep-sea trench Conrad
Earth’s crust discontinuity
Sea X
Mohorovičić
discontinuity
Lithosphere
Convection Bottom of
Subduction
lithosphere
Asthenosphere
Figure 5.14. The modern understanding of isostasy and plate tectonics. Deep-
D sea trenches are known to be in isostatic disequilibrium.
During the last glacial maximum, some 20,000 years ago, Fennoscandia
was covered by a continental ice sheet of thickness up to 3 km.
Questions:
1. How much was the Earth’s surface depressed by this load, assuming
isostatic equilibrium?
2. Currently the land is rising in central Fennoscandia, where the ice
thickness was maximal, at a rate of 10 mm /a. How long would it take
at this rate for the depression to disappear?
Answers:
1. We assume for the ice density a value of 920 kg /m3 . Then, with an
upper mantle density of 3370 kg /m3 — note that it is Earth’s mantle
material that is being displaced by the ice, the Earth’s crust just
transmits the load! See figure 11.1a — we find for the depression:
920 kg /m3
∆ H = 3000 m × = 819 m.
3370 kg /m3
2. At the rate of 10 mm /a it will take 819 m/0. 01 m /a = 81,900 years total.
Part of this uplift has already taken place since the last deglacia-
tion.
In reality, of course, the rate has decreased substantially, and will
continue to decrease, over time.
Ó »îá .
112 Geophysical reductions
10
Felix Andries Vening Meinesz (1887 – 1966) was a Dutch geophysicist, geode-
sist and gravimetrist. He wrote together with V. A. Heiskanen the textbook The
Earth and its Gravity Field (1958).
Ó »îá .
The “isostatic geoid” 113
18 20 22 24 26 28 30 32
64 64
62 62
60 60
18 20 22 24 26 28 30 32
Ó »îá .
114 Geophysical reductions
We obtain for the potential field mass density layer at sea level, when
also the evaluation point is placed at at sea level, H = 0 =⇒ r = R :
¨ ∞
∑ ( )
Ttop = GR κ P n cos ψ d σ
σ n=0
11
Of course Bouguer reduction is even worse! The indirect effect can be hundreds
of metres.
Ó »îá .
The “isostatic geoid” 115
Sea level
Compensation depth
12
This works on dry land and on the ocean. Lakes, glaciers and areas like the
Dead Sea are more complicated.
Ó »îá .
116 Geophysical reductions
so that
∞ [ ( )n ]
∑ 2 R−D
δTiso = − R 1− 2πG κn =
2n + 1 R
n=1
∞ [ ( )n ]
∑ 2 R−D
=− R 1− [ A B ]n .
2n + 1 R
n=1
and
D AB δTiso
δ Niso = . ≈− (5.8)
γ γ
This is the indirect effect of isostatic reduction.
Let’s substitute realistic values. Let the Mohorovičić14 discontinuity’s
depth be on average ∼ 20 km15 .
13
The contribution from degree numbers n > R/D is
∞
∑ 2R
δTiso ≈ − [A B ]n ,
2n + 1
n= N +1
where the terms are small and rapidly falling to zero. In this degree range also
the mass density layer approximation for the topography breaks down.
14
Andrija Mohorovičić (1857 – 1936) was a Croatian meteorologist and a pioneer
of modern seismology.
15
Under the continents 35 km, under the oceans 7 km below the sea floor, accord-
ing to Encyclopædia Brittannica. Using these values, we find δ Niso = −3.2 m on
land, +2.8 m on the ocean.
Ó »îá .
Self-test questions 117
D Self-test questions
Ó »îá .
118 Geophysical reductions
2. See section 5.2: Bouguer anomalies. Derive the equations 5.2 and
5.3 anew, assuming that the mean density of the Earth crust were
ρ = 3370 kg /m3 .
Q′ Sea level
The vertical rock wall PQ is also straight on a map and extends in both
directions (“into” and “out of” the paper) to infinity.
Height differences: PQ ′ = 600 m,QQ ′ = 300 m.
1. Compute in point P the terrain correction (hint: use the formula
for the attraction of a Bouguer plate. We have here a half Bouguer
plate, with only half the attraction of a full one.)
2. Compute in point Q the terrain correction. Algebraic sign?
3. If in point P is given that the free-air anomaly is 60 mGal, how
much is then the Bouguer anomaly in the point? (Use the complete
Bouguer reduction.)
4. If it is given in point Q that the Bouguer anomaly is 10 mGal, how
much is the point’s free-air anomaly?
D Exercise 5 – 4: Isostasy
Ó »îá .
D Vertical reference systems
6
10
1001 11
Level
back
11001100
00
00
11
fore
20
View
∆ H = back − fore
– 119 –
120 Vertical reference systems
⃝ ∆H
∑
around a closed path is generally not zero.
Geometric height is not a conservative field.
This is why, in precise levelling, the height differences are always con-
verted to potential differences: ∆W = −∆ H · g, in which g is the local
gravity, which is either measured or — e.g., in Finland — interpolated
from an existing gravity survey data base. The sum of potential differ-
ences around a closed loop is always zero: ⃝ ∆W = 0.
∑
P
∑
C P = − (WP − W0 ) = (∆ H · g ) ,
sea level
positive above sea level, is called the geopotential number of point P . geopotentiaaliluku
1
However, the value engraved in the pillar is the reference height of the still
older system NN, not of N60. The correct reference value for N60 for this pillar,
30.51376 m, is given in the publication Kääriäinen (1966).
2
The district Pointe-au-Père of the city of Rimouski was named after the Jesuit
priest Father Henri Nouvel (1621? – 1701?), who served forty years with the na-
tive population of New France, today’s Quebec. Pointe-au-Père is also notorious
as the location of the RMS Empress of Ireland shipwreck in 1914, in which over
a thousand passengers perished.
Ó »îá .
Orthometric heights 121
Ó »îá .
122 Vertical reference systems
g
P
g WP
∆ H3
∆ H2 H ∆ H2′
∆ H1 ∆ H1′
W0
Geoid
Figure 6.3. Levelled heights and geopotential numbers. The height obtained by
∑3
summing levelled height differences, i=1 ∆ H i , is not the “correct”
∑3
height above the geoid, i.e., i=1 ∆ H ′i computed along the plumb
line.
Note how the equipotential or level surfaces of the geopotential are
not parallel: because of this, a journey along the Earth’s surface
may well go “upward”, to increasing heights above the geoid, al-
though the geopotential number decreases. Thus, water may flow
“upward”.
The gravity vector g is everywhere perpendicular to the level sur-
faces, and its length is inversely proportional to the distance sepa-
D rating the surfaces.
Ó »îá .
Normal heights 123
North South
gS
Lake Päijänne: C = − (W − W0 ) = 76.9 GPU
gN
HS
Lake Päijänne
HN
Geoid: W = W0
and z is the measured distance along the plumb line. Because the formula
for g already itself contains H , we obtain the solution iteratively, using
initially a crude estimate for H . The iteration converges fast.
We shall see that determining very precise orthometric heights is chal-
lenging, especially in the mountains.
def C
H∗ = ,
γ0H
where γ0H is the average normal gravity computed between the zero level
(reference ellipsoid) and H ∗ along the ellipsoidal normal. So, the same
way of computing as in the case of orthometric heights, but using the
normal gravity field instead of the true gravity field.
Ó »îá .
124 Vertical reference systems
Heights “above sea level” are for practical reasons given in metres. For
large, continental networks we want to give heights above a computa-
tional reference ellipsoid in metres, and thus also heights above “sea
level” have to be in metres.
Molodensky proposed also that instead of the geoid, height anomalies
would be used, the definition of which is
def T
ζ= , (6.1)
γHh
where now γHh is the average normal gravity at terrain level, more
precisely: between the heights H ∗ (but reckoned from the ellipsoid) and
h. T is the disturbing potential.
Based on these assumptions, he showed that
H ∗ + ζ = h,
Ó »îá .
Normal heights 125
numbers, and that also would be compatible with similarly defined, so-
called height anomalies, and with geometric heights h reckoned from the
reference ellipsoid.
The geometric height h from the reference ellipsoid may be connected to
the potential U of the normal gravity field indirectly, though the following
integral equation:
ˆ h
U = U0 − γ ( z) dz .
0
Here, U is the normal potential and γ normal gravity. One level surface
of U , U = U0 , is also the reference ellipsoid. The variable z is the distance
from the ellipsoid along the local normal.
By defining
ˆ h
1
def
γ0h = γ ( z) dz
h 0
we obtain
U − U0
h=− .
γ0h
By using W = U + T and dividing by γ0h we obtain:
W − W0 T
= −h
γ0h γ0h
? T
N + = h − H+ =
γ0h
as the corresponding new geoid height type. It has however the esthetic
flaw, that we divide here by the average normal gravity computed between
the levels 0 and h. This quantity is not operational without a means of
determining the ellipsoidal height h.
This suggests the following improvement based on the circumstance that
γ ( z) is a nearly linear function. This means that the vertical derivative
∂γ
∂r
is nearly constant in the height interval considered.
We define
ˆ h ˆ H+ ˆ h
1
def 1
def def1
γ0h = γ ( z) dz, γ0H = + γ ( z) dz, γHh = + γ ( z) dz .
h 0 H 0 N H+
Now
dγ N+
( )
1
γ0H ≈ γ0h − N + ≈ γ0h 1 − (6.2)
2 dr R
Ó »îá .
126 Vertical reference systems
dγ 2γ
(R is the Earth’s radius, and dr ≈ R is assumed constant), and
dγ H+
( )
1
γHh ≈ γ0h + H + ≈ γ0h 1 + . (6.3)
2 dr R
N+ H+
Next, we also exploit that both R ≪ 1 and R ≪ 1, so
)−1 )−1
N+ N+ H+ H+
( ( ) ( ( )
1− ≈ 1+ , 1− ≈ 1+ ,
R R R R
W − W0 γ0h N+ N + H+
( )
∗ W − W0
def
H =− =− · ≈ H+ 1 + = H+ + ,
γ0H γ0h γ0H R R
T γ0h H+ N + H+
( )
def T +
ζ= = · ≈N 1− = N+ − .
γHh γ0h γHh R R
N + H+
Because the, already small, correction terms R cancel, we finally
obtain
H ∗ + ζ = H + + N + = h. (6.4)
The quantity γ0H , and thus also normal height H ∗ , can be, unlike γ0h ,
computed using only information obtained by (spirit or trigonometric) lev-
elling, without having to know the height h above the reference ellipsoid,
which would require again knowledge of the local geoid.
This was Molodensky’s realization (Molodensky et al., 1962) already in
1945, long before the global positioning system GPS, or a global, geo-
centric reference ellipsoid, existed. Back then, continental triangulation
networks were computed on their own, regional reference ellipsoids.
+ +
The size of the correction term N RH is, taking for the heights of the global
geoid up to 110 m, 17 mm for each kilometre of terrain height. The errors
remaining after applying this term are microscopically small, because
normal gravity is, unlike true gravity, extremely linear along the plumb
line — as equations 6.2 and 6.3 already assumed.
Figure 6.6 attempts to visualize the derivation.
Normal height:
C W − W0
H∗ = =− , (6.5)
γ γ
Ó »îá .
Normal heights 127
N + H+
R
ζ N+
z
h
H+ H∗
N + H+
R
γ (z)
Figure 6.6. A graphic cartoon of Molodensky’s proof. The blue and red areas,
which are equal, represent the correction terms which convert N +
to ζ and H + to H ∗ respectively. The red and blue arrows describe
the conversion process. The balls represent midpoints of averaging
D intervals for the function γ (z).
Height anomaly:
W −U
ζ=
γHh
where ˆ h
1
γHh = γ ( z) dz.
ζ H∗
The height anomaly ζ, which otherwise is a quantity similar to the
geoid height N , however is located at the level of the topography,
not at sea level. The surface formed by points which are a distance
H ∗ above the reference ellipsoid (and thus a distance ζ below the
topography), is called the telluroid. It is a mapping of sorts of the
topographic surface: the set of points Q , whose normal potential
UQ is the same as the true potential WP of the true topography’s
Topography
ζ
Telluroid
H∗ H h
Figure 6.7. Geoid, quasi-geoid, telluroid and topography. Note the correlation
D between quasi-geoid and topography.
Ó »îá .
128 Vertical reference systems
U − U0
h= ,
γ0h
where ˆ h
1
γ0h = γ ( z) dz.
h 0
h = H ∗ + ζ.
In all three cases, the quantity is defined by dividing the potential differ-
ence by some sort of “average normal gravity”, suitably computed along
a segment of the local plumb line. In the case of the height anomaly ζ,
a piece of plumb line is used high up, close to the topographic surface,
between level H ∗ (telluroid) and level h (topography).
Ó »îá .
Difference between geoid height and height anomaly 129
TH
and by using the Bruns equation twice, ζ = (height anomaly
γ
T0
or quasi-geoid height) and NFA = (”free-air geoid height”, FA =
γ
Free Air) we obtain3
∆ g FA H
ζ − NFA ≈ − . (6.7)
γ
2πGH ρ .
3
Here we made the approximation, that γ is the same on topography level as on
sea level.
Ó »îá .
130 Vertical reference systems
See also Heiskanen and Moritz (1967, pages 8 – 13). As in the mountains
the Bouguer anomaly is strongly negative, it follows that the quasi-geoid
is there always above the geoid: approximately, using equation 5.2:
h = H ∗ + ζ,
Ó »îá .
Calculating orthometric heights precisely 131
The method is recursive: H appears both on the left and on the right side.
This is not a problem: both H and g are obtained iteratively. Convergence
is fast.
In practice one calculates orthometic height using an approximate for-
mula. In Finland, Helmert orthometric heights have long been used,
where gravity measured on the Earth’s surface, g ( H ), is extrapolated
downward by using the estimated vertical gravity gradient interior to
the rock. It is assumed that its standard value outside the rock, the
value, −0.3086 mGal /m (the free-air gradient), changes to a value that is
+0.2238 mGal /m greater (double Bouguer plate effect): the end result is the
total inside-rock gravity gradient, −0.0848 mGal /m.
This is called the Prey4 reduction. The end result are the following
equations (the coefficient is half the gradient, i.e., the mean gravity along
the plumb line is the same as gravity at the midpoint of the plumb line):
( 1 )
g = g ( H ) − 0.0848 mGal /m − H = g ( H ) + 0.0424 mGal /m · H, thus
2
C
H= , (6.10)
g ( H ) + 0.0424 mGal /m · H
in which C is the geopotential number (potential difference with mean
sea level) and g ( H ) is gravity at the Earth’s surface. See also Heiskanen
and Moritz (1967) pages 163 – 167. Note that the term 0.0424 mGal /m · H is
typically much smaller than g ( H ) , which is about 9.8 m /s² = 980,000 mGal!
So, iteration, where the above denominator is first calculated using a
crude H value, converges really fast.
The use of Helmert heights as an approximation to orthometric heights
is imprecise for the following reasons:
◦ The assumption that gravity changes linearly along the plumb line.
This is not the case, especially not because of the terrain correction.
In the precise computation of orthometric heights, one ought to
compute the terrain correction separately for every point on the
plumb line.
◦ The assumption that the free-air vertical gravity gradient is a
constant, −0.3086 mGal /m. This is not the case, the gradient can
easily vary by ±10%.
◦ The assumption that rock density is ρ = 2.67 g /cm³. The true density
value may easily vary by ±10% or more around this assumed value.
The first approximation, neglecting the terrain effect, can be corrected by
using Niethammer5 ’s method (see Heiskanen and Moritz (1967) page 167).
4
Adalbert Prey (1873 – 1949) was an Austrian astronomer and geodesist and an
author of textbooks.
5
Theodor Niethammer (1876 – 1947) was a Swiss astronomer and geodesist who
was the first to map the gravity field of the Swiss Alps.
Ó »îá .
132 Vertical reference systems
Ó »îá .
Calculation example for heights 133
Questions:
1. Calculate the free-air anomaly ∆ g FA of point P .
2. Calculate the orthometric height of the point.
3. Calculate the Bouguer anomaly (without terrain correction)
∆ g B of the point.
4. Calculate the normal height of point P .
5. If the geoid height at point P is N = 25.000 m, how much is
then the height anomaly (“quasi-geoid height”) ζ?
Answers:
1. The free-air anomaly is
ζ = N − (−0.055 m) = 25.055 m.
Ó »îá .
134 Vertical reference systems
The fact that the difference in orthometric heights between two points A
and B is not equal to the sum of the levelled height differences is due to
gravity not being the same everywhere.
With C A , C B , ∆C the geopotential numbers at A and B, and the geopo-
∑B
tential differences along the levelling line, we have C B − C A − A ∆C = 0
because of the conservative nature of the geopotential. Dividing by a
constant γ0 yields
B
C B C A ∑ ∆C
− − = 0.
γ0 γ0 γ0
A
On the other hand we have
B B
∑ C B C A ∑ ∆C
OC AB = HB − H A − ∆H = − − ,
gB g A g
A A
in which
γ0 − g B C B γ0 − g B
( ) ( )
CB CB
− = = HB ,
gB γ0 γ0 gB γ0
γ0 − g A
( )
CA CA
− = HA,
gA γ0 γ0
∆C ∆C γ0 − g
( )
− = ∆ H,
g γ0 γ0
Ó »îá .
A vision for the future: relativistic levelling 135
B (
g − γ0 g A − γ0 gB − γ0
∑ ) ( ) ( )
OC AB = ∆H + HA − HB , (6.13)
γ0 γ0 γ0
A
B B
∑ C B C A ∑ ∆C
∗
NC AB = HB − H ∗A − ∆H = − − ,
γB γ A g
A A
B (
g − γ0 γ A − γ0 γB − γ0
∑ ) ( ) ( )
NC AB = ∆H + H ∗A − ∗
HB . (6.14)
γ0 γ0 γ0
A
Note that the identical first term in both equation 6.13 and equation 6.14
derives from the term
B B
∑ ∆C ∑
= ∆ H,
g
A A
B
∑
∗
HB = H ∗A + ∆ H + NC AB .
A
Ó »îá .
136 Vertical reference systems
6
Karl Schwarzschild (1873 – 1916) was a German physicist who was the first to
derive, in 1915 while serving on the Russian front, a closed spherically symmetric
solution to the field equation of Albert Einstein’s general theory of relativity, the
Schwarzschild metric.
Ó »îá .
Self-test questions 137
Braunschweig
Garching
Figure 6.8. An optical lattice clock: the ultra-precise atomic clock of the future
operates at optical wavelengths. To the right, the trajectory of the
D Predehl et al. experiment.
some 100 km, which must be replaced by modified ones (Predehl et al.,
2012). In this way both the traditional precise levelling networks and the
height systems based on GNSS technology and geoid determination may
be replaced by this hi-tech (and hi-science!) solution.
D Self-test questions
Ó »îá .
138 Vertical reference systems
In Finland, in connection with the N60 system, heights used were so-
called Helmert heights, a good approximation of orthometric heights.
The potential difference with sea level at point P , − (W − W0 ), equals
1000 m2 s−2 . Gravity in the point is g P = 9.820000 m s−2 . Calculate the
precise orthometric height of the point.
− (W − W0 ) = 5000 m2 /s2 .
Ó »îá .
D
7 The Stokes equation and other
integral equations
(in which, for the degree numbers n = 0, 1, the ∆ g n are assumed to vanish,
as n = 0 means a different total mass for the normal field than for the
Earth, and n = 1 an offset of the co-ordinate origin from the Earth’s centre
of mass, see section 2.12). This is now the Stokes equation’s spectral form.
Using the degree constituent equation 2.16 one obtains the integral
equation
∞ ¨
R ∑ 2n + 1
∆ gP n cos ψ d σ =
( )
T=
4π n−1 σ
n=2
¨ [∑∞
]
R 2n + 1
∆ gd σ =
( )
= P n cos ψ
4π σ n−1
n=2
¨
R
S ψ ∆ gd σ,
( )
=
4π σ
in which
∞
( ) ∑ 2n + 1 ( )
S ψ = P n cos ψ ,
n−1
n=2
the Stokes kernel function. The angle ψ is the geocentric angular distance
between evaluation point and moving observation point, see figure 7.2.
The equation above allows the calculation, from global gravimetric data
and for every point on the Earth surface, of the disturbing potential T ,
and from that the geoid height N using the Bruns equation 4.2, N = Tγ ,
with the result
– 139 –
140 The Stokes equation and other integral equations
Mass excess
Mass
deficit
−N
¨
R
S ψ ∆ g φ′ , λ′ d σ′ ,
( ) ( ) ( )
N φ, λ = (7.1)
4πγ σ
( ) ( )
in which φ, λ and φ′ , λ′ are the evaluation point and the moving
point (“observation point”) respectively, and the distance between them
is ψ. Equation 7.1 is the classical Stokes equation of gravimetric geoid
( )
N φ, λ S(ψ)
Evaluation
point
∆ g φ′ , λ′ d σ
( )
ψ
Moving data or
integration point
Earth’s centre
Ó »îá .
The Stokes equation and the Stokes integral kernel 141
25
( )
S ψ
1
20
sin (ψ/2)
ψ
−6 sin 2 + 1 − 5 cos ψ
ψ ψ)
+ sin2
(
15 −3 cos ψ ln sin 2 2
10
5
S(ψ) →
ψ→
−5
0 0.5 1 1.5 2 2.5 3 3.5
( )
Figure 7.3. The Stokes kernel function S ψ . The argument ψ is in radians
[0, π). Also plotted are the three parts of the analytical expression
D 7.2 with their different asymptotic behaviours.
determination.
The above illustrates the correspondence between integral equations
and spectral equations. There are other examples of this. Earlier we
presented the spectral representation of the function 1ℓ , Heiskanen and
Moritz (1967) equation 1 – 81. Of course 1ℓ is also the kernel function of an
integral equation, the one yielding the potential V if given is the single
layer mass density κ.
Also a version of the Stokes equation for the exterior space exists. We
gave it earlier, equation 4.10. The spectral form of its kernel function, see
equation 4.11, is
∞ ( ) n+1
( ) ∑ R 2n + 1 ( )
S r, ψ, R = P n cos ψ .
r n−1
n=2
1 ψ ( ψ ψ)
S (ψ) = − 6 sin + 1 − 5 cos ψ − 3 cos ψ ln sin + sin2 . (7.2)
sin (ψ/2) 2 2 2
Ó »îá .
142 The Stokes equation and other integral equations
ψ
when ψ → 0. The next three terms, −6 sin 2 + 1 − 5 cos ψ, are all bounded
on the whole interval [0, π) and the limit for ψ → 0 is −4. The last, com-
( ψ ψ)
plicated term −3 cos ψ ln sin 2 + sin2 2 goes also to infinity — positive
infinity! — for ψ → 0, but much more slowly, because of the logarithm.
These equations were derived for the first time by the Dutch geophysicist
F. A. Vening Meinesz. The angle α is the azimuth or direction angle
( )
between the calculation or evaluation point φ, λ and the moving inte-
( )
gration or observation point φ′ , λ′ . These equations are much harder to
write in spectral form, as the kernel functions are now also functions of
the azimuth direction α, in other words, they are anisotropic.
The disturbing potential, the gravity disturbance, and the gravity anom-
aly, are all so-called isotropic quantities: they do not depend on the
azimuth, and therefore in the spectral representation the transformations
between them are only functions of harmonic degree n.
1
We may also write the function ℓ
as the following expansion (for proof,
1
Carl Gustav Jacob Jacobi (1804 – 1851) was a German mathematician.
Ó »îá .
The Poisson integral equation 143
z
P
ℓ r
y
Q ψ
R
Figure 7.4. The geometry of the generating function of the Legendre polynomi-
D als.
in which r = ∥r∥ and R = ∥R∥ are the distances of points P and Q from the
origin, the centre of the Earth. The function 7.4 is called the generating
function of the Legendre polynomials.
Differentiating equation 7.4 with respect to r yields
∞ ( )n+1
r − R cos ψ 1 ∑ n+1 R
− = − P n (cos ψ).
ℓ3 R r r
n=0
This we multiply by 2 r :
∞ ( )n+1
2 r 2 − 2 rR cos ψ 1∑ R
− = − (2 n + 2) P n (cos ψ).
ℓ3 R r
n=0
−2 r 2 + 2 rR cos ψ + ℓ2 R 2 − r 2
= ,
ℓ3 ℓ3
and the end result is, by multiplying with −R ,
∞ ( )n+1
R r2 − R 2
( )
∑ R ( )
= (2 n + 1) P n cos ψ . (7.5)
ℓ3 r
n=0
Ó »îá .
144 The Stokes equation and other integral equations
Applying now the degree constituent equation 2.16 to the harmonic po-
tential field V on the spherical Earth’s surface, radius R :
¨
2n + 1
V φ′ , λ′ , R P n cos ψ d σ′ ,
( ) ( ) ( )
Vn φ, λ =
4π σ
we obtain
∞ ( )n+1 ¨
1 ∑ R
V φ′ , λ′ , R P n cos ψ d σ′ =
( ) ( ) ( )
V φ, λ, r = (2 n + 1)
4π r σ
n=0
¨ [∞ ( )n+1 ]
1 ′ ′
) ∑ R
d σ′ =
( ( )
= V φ ,λ ,R (2 n + 1) P n cos ψ
4π σ r
n=0
¨ (
r 2 − R 2 V φ′ , λ′ , R
) ( )
R
= d σ′
4π σ ℓ3
by substituting the expression in square brackets directly into equation
7.5.
Thus we have obtained the Poisson equation for computing a harmonic
field V from values given on the Earth’s surface:
¨ (
r2 − R 2
)
R
VP = VQ d σQ , (7.6)
4π σ ℓ3
in which ℓ is again the straight distance between evaluation point P
(where VP is being computed) and moving data point Q (on the surface of
the sphere, VQ under the integral sign). In this equation we have given
the points symbolic names: the co-ordinates of evaluation point P are
( ) ( )
φ, λ, r , the co-ordinates of the data point Q are φ′ , λ′ , R .
Still a third way to write the same equation, useful when the function
or field V isn’t actually defined between the topographic Earth’s surface
and sea level, is
¨ ( 2
r − R2 ∗
)
R
V= V d σ,
4π σ ℓ3
in which V ∗ denotes the value of a harmonically downward continued
function V — downward continued into the topography, all the way down
to sea level, or, in spherical approximation, to the surface of the sphere
r = R . This is a function that above the topography is identical to V ,
that is harmonic, and that also exists between topography and sea level.
The existence of such a function has been a classical theoretical nut to
crack. . .
Equation 7.6 solves the so-called Dirichlet boundary-value problem, find-
ing a function in an area of space when the value of the function on its
boundary has been given.
Ó »îá .
Gravity anomalies in the exterior space 145
i.e.,
∞ ( ) n+1 ∞ ( ) n+1
∑ R ∑ R
r∆ g = ( n − 1) T n = Sn,
r r
n=2 n=2
( ) def ( )
in which S n φ, λ = ( n − 1) T n φ, λ is a perfectly legal surface spherical
( )
harmonic just like T n φ, λ itself. Also, the dependence on the radius
( )n+1
r , the factor Rr , is the same as for the (harmonic) potential. So,
Poisson’s integral equation 7.6 applies to function r ∆ g:
¨ (
r2 − R 2 R ∆ g φ′ , λ′ , R
)[ ( )]
R
r ∆ g φ, λ, r d σ′ ,
[ ( )]
=
4π σ ℓ3
or
¨ (
r2 − R 2
)
R2
∆ g φ, λ, r ∆ g φ′ , λ′ , R d σ′ .
( ) ( )
= (7.8)
4π r σ ℓ3
An alternative notation:
¨
R2 r2 − R 2
∆g = ∆ g ∗ d σ,
4π r σ ℓ 3
Ó »îá .
146 The Stokes equation and other integral equations
in which
∞ ( ) n+2
( ) def ∑ R ( )
K r, ψ, R = (2 n + 1) P n cos ψ
r
n=2
is a (modified) Poisson kernel for gravity anomalies. Its closed form can
be lifted from equation 7.8:
( ) R 2 r2 − R 2
K r, ψ, R = .
r ℓ3
Compared to the Stokes kernel, the Poisson kernel drops off fast to zero
for growing ℓ values. In other words, the evaluation of the integral
equation may be restricted to a very local area, e.g., a cap of radius kalotti
1◦ . See figure 7.5. The main use of Poisson’s kernel is the harmonic
continuation, upward or downward, i.e., the shifting of gravity anomalies
measured and computed at various levels to the same reference level.
In the limit r → R (sea level becomes the level of evaluation) this kernel
function goes asymptotically to the Dirac δ function.
Ó »îá .
The vertical gradient of the gravity anomaly 147
0.2
Distance (km) −→
0 1 2 3 4 5 6 7
Figure 7.5. The Poisson kernel function for gravity anomalies as well as the
kernel for the anomalous vertical gravity gradient, both at various
height levels. These are kernel values to be used when evaluating
D the surface integral in map co-ordinates in km.
i.e.,
∞ ( )n+3
∂∆ g φ, λ, r
( )
1 ∑ R
=− (2 n + 1) ( n + 2) ·
∂r 4π R r
n=2
¨
∆ g φ′ , λ′ , R P n cos ψ d σ′ =
( ) ( )
·
σ
¨
1
K ′ r, ψ, R ∆ g φ′ , λ′ , R d σ′ ,
( ) ( )
= (7.10)
4πR σ
2
Hint: use symbolic algebra software.
Ó »îá .
148 The Stokes equation and other integral equations
¨ (
r 2 − R 2 ∆ g φ′ , λ′ , R
) ( )
1 R2
− d σ′ =
r 4π r ℓ3
⎡σ
3 r2 − R 2
( ) ⎤
¨ 2− −
R2 1 ⎢ 2 r 2
⎥ ∆ g φ , λ , R dσ −
⎥ ( ′ ′ ) ′
= ⎢
4π ℓ 3⎣ 2 2 2 2 ⎦
( )( )
σ 3 r − R r − R
−
2 r 2 ℓ2
1
− ∆ g φ, λ, r =
( )
r
¨ [ )2 ]
3 r2 − R 2
(
R2 1
∆ g φ′ , λ′ , R d σ′ −
( )
= 2−
4π σℓ
3 2 r 2 ℓ2
[1
3]
∆ g φ, λ, r =
( )
− +
r 2 r[
¨ )2 ]
3 r2 − R 2
(
R2 1 ( ′ ′ ) ′ 5
∆ ∆
( )
= 2 − g φ , λ , R d σ − g φ , λ, r .
4π σ ℓ 3 2 r 2 ℓ2 2r
(7.11)
In the final right-hand side, the rightmost term is very small, typically
under one part in a thousand, compared to the leftmost term. Both terms
inside the square brackets are of the same order of magnitude.
Finally, note that, if we are integrating over the surface of the Earth
sphere, radius R , rather than the unit sphere σ, radius 1, the coefficient
R 2 drops out from both equation 7.8 and equation 7.11.
∆ g(0) φ′ , λ′ , R ≈ ∆ g φ′ , λ′ , r = ∆ g φ′ , λ′ , R + H ,
( ) ( ) ( )
( ) ( )
where H = H φ′ , λ′ is the topographic height at point φ′ , λ′ . When the
first, crude anomalous gradient has been calculated, e.g., using equation
7.11, we may perform a real reduction to sea level, initially linearly:
⏐(0)
(1) ∂∆ g ⏐⏐
∆g φ′ , λ′ , R ≈ ∆ g φ′ , λ′ , R + H −
( ) ( )
H,
∂r ⏐
and so forth.
3
In the derivation there it is assumed that ∆ g is harmonic. Not so: r ∆ g is
harmonic. The error made is small.
Ó »îá .
Gravity reductions in geoid determination 149
Imagine that, conceptually, the topographic masses are shifted into the
geoid, to below sea level, in a way that does not change the exterior field.
This is materially the same as determining the geoid associated with the
harmonically downward continued exterior field.
The problem here is, that such a mass distribution below sea level, which
produces the harmonically downward continued external potential in the
space between topographic surface and geoid, doesn’t always precisely
exist. Or, that a suitable mass distribution will contain extremely large
positive and negative masses close to each other, which are physically
unrealistic.
huonosti asetettu One expresses this by saying that the problem is “ill posed”. In such cases
one uses regularization: one changes a little — as little as possible —
the exterior field, so that it corresponds precisely to some sensible field
that is harmonically continued below the topographic surface, and some
sensible mass distribution interior to the geoid that produces it.
Ó »îá .
150 The Stokes equation and other integral equations
One can start, e.g., by filtering out the short-wave parts caused by the
topography using a high resolution digital terrain model. This is called
the RTM (residual terrain modelling) method.
In this method, we do not actually move all topographic masses to be-
low the geoid. Instead, only masses close to the topographic surface are
either removed or filled in, in a way which creates a smooth replace-
ment topography that is long-wavelength only. The exterior field of this
smoothed topography, unlike that of the original topography, lacks the
shortest wavelengths. It may thus be downward continued to the geoid
with sufficient precision.
First we computationally remove from the topography only the short
wavelengths (under 30 km) by moving the masses of the peaks into the
valleys, i.e., we do a low-pass filtering. The effect of this on the free-air
gravity anomalies ∆ g calculated from measurements is evaluated and
taken into account: the “Remove” step.
In detail:
1. In each point P we apply the terrain correction as described in
section 5.3 to the gravity anomalies.
2. Next, we remove the attraction of a Bouguer plate, thickness H −
HRTM , where H stands for the terrain height of point P , and HRTM
for the height of the smoothed, or low-pass filtered, terrain at the
location of P . This effect is, according to equation 5.2 on page 98,
equal to
2πG ρ ( H − HRTM ) ,
in which ρ is the rock density assumed in the calculation.
3. After this, the location of the gravity anomaly is moved down4 —
“downward continuation” — from the original terrain level H to the
surface of the new, smoothed terrain, HRTM . For this, the Poisson
equation may be used as described in section 7.8.
If the vertical gravity gradient of the terrain-reduced field equals
the normal vertical gradient of gravity — according to section 4.4,
0.3 mGal /m —, this will leave the anomaly unchanged. Typically,
there will be a small change: one may show — exercise 1 – 1 item
4 — that on the surface of a buried sphere of anomalous density
∆ρ , there will be a radial anomalous gravity gradient of 83 πG ∆ρ .
For ∆ρ ≪ ρ , this will be negligible compared to the Bouguer-plate
coefficient in item 2.
4. Rigorously speaking, an inverse terrain correction for the shapes
of the smoothed terrain should be applied, to arrive at gravity
anomalies realistic for this new replacement topography. Often this
step is left out as the effect is small.
4
Or up!
Ó »îá .
Gravity reductions in geoid determination 151
P −
− − P′ +
+ +
−
Bouguer plate, down- Inverse
Terrain correction ward continuation terrain correction
Figure 7.6. Residual terrain modelling (RTM). One removes from the terrain
computationally the short wavelengths, i.e., the differences from
the red dashed line: the masses rising above it are removed, the
valleys below it are filled. After reduction, the red line, smoother
than the original terrain, is the new terrain surface. The exterior
potential of the new mass distribution will differ only little from
the original one, but may be harmonically downward continued to
sea level.
Left, terrain correction for point P; middle, Bouguer plate and
gradient reduction to the level of the smoothed terrain point P ′ ;
D and right, the inverse terrain correction for point P ′ .
Because the mass shifts in the RTM method are so small, take place over
such small distances, and are so short wavelength in nature, the indirect
effect or “Restore” step — the change in geopotential due to the mass
shifts that has to be applied in reverse to arrive at the final geopotential
or geoid solution — is so small as to be often negligible. For the same
reason the effect of unknown topographic density will remain small.
Finally we note that, because the RTM method removes the effect of the
short-wavelength topography, it is also a suitable method for interpolating
gravity anomalies. See Märdla (2017).
Ó »îá .
152 The Stokes equation and other integral equations
After this, we apply, at sea level, the Stokes equation, and obtain the
disturbing potential at sea level T ∗ . After this, the disturbing potential
is “unreduced” back to terrain level, to the evaluation point, with the
equation
∂T
T = T∗ + H.
∂H
In these equations T , its vertical derivative, and ∆ g and its vertical
derivative always belong to the exterior harmonic gravity field, and the
connection between them is the fundamental equation of physical geodesy,
equation 4.5, in spherical geometry:
∂T 2
∆g = − − T,
∂H r
in which r = R + H . Here, we need firstly the vertical derivative of the
disturbing potential. This is easy: we have
∂T 2
= −∆ g − T,
∂H r
where the first term on the right is directly measured, and the second
term’s T is obtained iteratively from the main product of the solution
process.
Calculating the vertical gradient of gravity anomalies, i.e., the anomalous
vertical gradient of gravity, is harder. For this, we have the integral
equation 7.10:
¨
∂∆ g φ, λ, r
( )
1
K ′ r, ψ, R ∆ g φ′ , λ′ , R d σ′ ,
( ) ( )
=
∂r 4π R σ
∞ ( ) n+3
′
( ) ∑ R ( )
K r, ψ, R = − (2 n + 1) ( n + 2) P n cos ψ .
r
n=2
In the above equation 7.12 we used as the reference level the sea surface.
This is arbitrary: we may use whatever reference level, e.g., H0 , in which
case
¨ (
) ∂∆ g ( ′
)
R + H0
∆ g φ′ , λ′ , H ′ − H − H 0 S ψ d σ′ +
( ) ( ) ( )
T φ, λ, H =
4π σ ∂H
∂T
+ ( H − H0 ) .
∂H
If we now choose H0 = H , the last term drops off, and we obtain
¨ (
∂∆ g (
)
R+H
∆g φ ,λ , H −′ ′ ′ ′
S ψ d σ′ .
( ) ( ) ) ( )
T φ, λ, H = H −H
4π σ ∂h
Ó »îá .
The Remove-Restore method 153
In this case the reduction takes place from the height of the ∆ g measure-
ment point to the height of the T evaluation point, probably a shorter
distance than from sea level to evaluation height, especially in the imme-
diate surrounding of the evaluation point. This means that the lineariza-
tion error will remain smaller. What is bad, on the other hand, is that
the expression in parentheses is now different for each evaluation point.
This complicates the use of FFT based computation techniques, on which
more later.
Here, we were all the time discussing the determination of the disturbing
( )
potential T φ, λ, H ; this is in practice the same as determining the
height anomaly
( ) ( )
( ) T φ, λ, H T φ, λ, H
ζ φ, λ, H = ≈ ) ,
γHh
(
γ φ, H
1. From the observed gravity values, first the effect of the a global
gravity field model is removed. This model is generally given in the
form of a spherical-harmonic expansion. Thus, a residual gravity
field is obtained
◦ that has numerically smaller values — easier to work with —,
and
◦ that is more local: the long “wavelengths”, the patterns ex-
tending over large areas, have been removed from the residual
field, only the local details remain.
2. From the observed gravity, the effects are removed of all masses
that are outside the geoid. The purpose of this is to obtain a residual
gravity field
◦ on which the Stokes equation may be applied, because no
masses are left outside the boundary surface, and
◦ from which especially the very small “wavelengths” — details
the size of which is of order a few kilometres — caused by the
topography, are gone. After this, prediction of gravity values
from sparse measurement values will work better.
Ó »îá .
154 The Stokes equation and other integral equations
“Remove” “Restore”
∆g (−→ “Brute force” −→) N
⇓− L99 Global grav. field model 99K +⇑
∆ g loc Nloc
⇓− L99 Exterior masses 99K +⇑
∆ g red =⇒ Stokes =⇒ Nred
In this diagram the thick arrows denote calculations that are recom-
mended, because they are easy and accurate. The arrows in parentheses
refer to direct computation, which again is troublesome and compute
intensive.
and
∞
∑
∆ g red φ, λ = ∆ g n φ, λ ,
( ) ( )
n=L+1
Ó »îá .
Kernel modification in the Remove-Restore method 155
assuming that L is the largest degree number that is still along in the
global spherical-harmonic expansion, or gravity model, that was sub-
tracted from the data5 .
Now, because ∆ g n is a certain linear combination of the surface spherical
harmonics
{ ( )
( ) P n|m| cos ψ sin | m| α m = − n, . . . , −1,
Ynm ψ, α = ( )
P nm cos ψ cos mα m = 0, . . . , n,
i.e.,
n
∑
∆ g n ψ, α = ∆ g nm Ynm ψ, α ,
( ) ( )
m=− n
and also
( ) ( )
Yn0 ψ, α = P n cos ψ ,
it follows from the orthogonality of the Y functions, that
¨ ¨
( )
P n cos ψ Yn′ m d σ = Yn0 Yn′ m d σ = 0
σ σ
S ψ ∆ g red φ, λ = S ψ ∆ g red ψ, α =
( ) ( ) ( ) ( )
n
[ ∞
] [ ∞ ]
∑ 2n + 1 ∑ ∑
∆ g nm Ynm ψ, α
( ) ( )
= P n cos ψ × =
n−1
n=2 n=L+1 m=− n
n
[ ∞ ] [ ∞ ]
∑ 2n + 1 ( ∑ ∑
∆ g nm Ynm ψ, α
) ( )
= P n cos ψ × =
n−1
n=L+1 n=L+1 m=− n
L
ψ ∆ g red φ, λ ,
( ) ( )
=S
in which
∞
L
( ) ∑ 2n + 1 ( )
S ψ = P n cos ψ
n−1
n=L+1
is a so-called modified Stokes kernel function. The degree number L
is called the modification degree. The size of the evaluation area σ0 is
chosen to be compatible with this.
The modification method described here, restricting the Legendre expan-
sion of the S function to higher degree numbers, is called the Wong-Gore6
modification (Wong and Gore, 1969). A desirable property of the new
kernel function S L is, that it would be — at least compared to the original
5
. . . and that the model is accurate!
6
L. Wong and R.C. Gore worked at the Aerospace Corporation, a space technology
research institution in California. http://en.wikipedia.org/wiki/The_Aerospace_
Corporation.
Ó »îá .
156 The Stokes equation and other integral equations
25
S(ψ)
S 3 (ψ)
20 S 4 (ψ)
S 5 (ψ)
15 S 6 (ψ)
10
Figure 7.7. Modified Stokes kernel functions. Note how the kernel value outside
D the local area goes to zero for higher L values.
function S — small outside the cap area σ0 . In that case restricting the
integral to the cap instead of the whole unit sphere (equation 7.13) does
not do much damage. It is clear that S L is much narrower than S , as in
it, only the higher harmonic degrees are represented. This can be verified
by plotting a graph of both curves. It doesn’t go totally to zero outside the
cap, however, but oscillates somewhat.
The reason for this oscillation is, that in the frequency (i.e., degree num-
ber) domain, the cut-off of the modified kernel is quite sharp. Trans-
forming such a sharp edge between space and frequency domains will
invariably produce an oscillation, which is related to the so-called Gibbs7
phenomenon. Gibbsin ilmiö
∞ L
L
( ) ∑ 2n + 1 ( ) ∑ 2n + 1 ( )
S ψ = P n cos ψ + (1 − s n ) P n cos ψ =
n−1 n−1
n=L+1 n=2
L
( ) ∑ 2n + 1 ( )
=S ψ − sn P n cos ψ , (7.14)
n−1
n=2
7
Josiah Willard Gibbs (1839 – 1903) was an American physicist, chemist, mathe-
matician and engineer.
Ó »îá .
Advanced kernel modifications 157
over the area outside a local cap, σ − σ0 . Let us multiply this expression
( )
with each of the Legendre polynomials P n cos ψ , n = 2, . . . , L in turn,
and integrate over the area σ − σ0 outside the local cap:
ˆ
( ) ( )
S ψ P n cos ψ d σ −
σ−σ0
L ˆ
∑ 2 n′ + 1 ( ) ( )
− s n′ ′ P n′ cos ψ P n cos ψ d σ = 0, n = 2, . . . , L,
n −1 σ−σ0
n′ =2
with ˆ
( ) ( )
bn = S ψ P n cos ψ d σ
σ−σ0
and ˆ
2 n′ + 1 ( ) ( )
A nn′ = ′ P n′ cos ψ P n cos ψ d σ.
n −1 σ−σ0
From this we can solve the s n for every degree number n from 2 to L.
This solution sets to zero the expressions
ˆ
S L ψ P n cos ψ d σ,
( ) ( )
(7.15)
σ−σ0
8
The choice s n = 1 again gives the simply modified Stokes kernel from which the
low degrees have been completely removed.
Ó »îá .
158 The Stokes equation and other integral equations
9
Thomas Simpson FRS (1710 – 1761) was an English mathematician and text-
book writer. Actually Simpson’s rule was used already a century earlier by
Johannes Kepler.
Ó »îá .
The effect of the local zone 159
4
1 36 1
36 36
k=1
16
4 36 4
36 36 k=0
1 1
k = −1
36 4 36
36
j= -1 0 1
where ∆λ and ∆ϕ are the block sizes in the latitude and longitude direc-
tions, and w−1 = w1 = 16 , w0 = 46 are the weights.
Ó »îá .
160 The Stokes equation and other integral equations
and substitute:
ˆ ψ0 ˆ 2π
∂∆ g ∂∆ g
( )[ ( )]
1 2
ξint ≈ − 2 ∆ g 0 + R sin ψ cos α + sin α ·
4πγ 0 0 ψ ∂x ∂y
· cos α sin ψ d α d ψ,
ˆ ψ0 ˆ 2π
∂∆ g ∂∆ g
( )[ ( )]
1 2
η int ≈ − 2 ∆ g 0 + R sin ψ cos α + sin α ·
4πγ 0 0 ψ ∂x ∂y
· sin α sin ψ d α d ψ .
´ 2π
Here, the terms in ∆ g 0 drop out in α integration (because 0 sin α = 0),
as do the mixed terms in sin α cos α. What remains is
ˆ ψ0 ˆ 2π
1 2 ∂∆ g
ξint ≈ − R sin ψ cos α cos α · sin ψ d α d ψ ≈
4πγ 0 0 ψ2 ∂x
ˆ ψ0 ˆ 2π ˆ ψ0
R ∂∆ g 2 R ∂∆ g
≈− cos α · d α d ψ = − d ψ,
2πγ 0 0 ∂x 2γ 0 ∂x
ˆ ψ0 ˆ 2π
1 2 ∂∆ g
η int ≈ − R sin ψ sin α sin α · sin ψ d α d ψ ≈
4πγ 0 0 ψ2 ∂y
ˆ ψ0 ˆ 2π ˆ ψ0
R ∂∆ g 2 R ∂∆ g
≈− sin α · d α d ψ = − dψ .
2πγ 0 0 ∂y 2γ 0 ∂y
d ∂∆ g d ∂∆ g
ξint = − , η int = − .
2γ ∂ x 2γ ∂ y
Sometimes these equations are useful, e.g., when estimating the errors of
grid based methods. Let the mesh size of a grid be ∆ x, then we may set
in the above equations d ≈ ∆2x , and for ∆ g we take
∆ gobs − ∆ ggrid ,
in which ∆ ggrid is the gravity anomaly value interpolated from the grid
file at the evaluation point. In this way one obtains a rough estimate of
how much error is due to the mesh size of the grid.
D Self-test questions
1. What do the Stokes equation and its spectral form look like?
( )
2. What does the Stokes kernel function S ψ look like when expanded
in Legendre polynomials?
3. What is a suitable aaproximation of the Stokes kernel when ψ is
small?
4. What is an isotropic, what an anisotropic quantity on the Earth’s
surface? Give an example of the latter.
Ó »îá .
Exercise 7 – 1: The Stokes equation 161
Ó »îá .
D Spectral techniques, FFT
8
( )
in which φ′ , λ′ is the location of the moving point — i.e., the integration
( )
or observation point — on the Earth’s surface, and φ, λ the location
of the evaluation point, it too on the Earth’s surface. Generally the
( )
locations of both points are given in spherical co-ordinates φ, λ , and
correspondingly the integration is executed over the surface of the unit
sphere σ: a surface element is d σ = cos φ d φ d λ, in which the factor cos φ
represents the determinant of Jacobi, for these spherical co-ordinates
( )
φ, λ .
However locally, in a sufficiently small area, one may write the point
co-ordinates also in rectangular form, and then express also the integral
in rectangular co-ordinates. Suitable rectangular co-ordinates are, e.g.,
map projection co-ordinates, see figure 8.1.
A simple example of rectangular co-ordinates in the tangent plane would
be
x = ψR sin α,
(8.1)
y = ψR cos α,
where α is the azimuth of the connection between evaluation point and
moving data point. The centre of this projection is the point where
the tangent plane touches the sphere. The locations of other points are
measured by the angle at the Earth’s centre, ψ, i.e., the geocentric angular
distance, and by the direction angle in the tangent plane or azimuth α.
A more realistic example uses a popular conformal map projection, the
stereographic projection:
(ψ)
x = 2 sin 2 R sin α,
(ψ)
y = 2 sin 2 R cos α.
– 163 –
164 Spectral techniques, FFT
y
Data point
α
x
Evaluation point
R ψ
Earth’s centre
In the limit for small values of ψ this is the corresponds to equations 8.1.
Taking the squares of equations 8.1, summing them, and dividing them
by R 2 yields
x2 + y2
ψ2 ≈ .
R2
More generally ψ is the angular distance between the two points ( x, y)
( )
(evaluation point) and x′ , y′ (data, integration or moving point) seen
from the Earth’s centre, approximately
)2 ( )2
x − x′ y − y′
(
2
ψ ≈ + .
R R
d σ = R −2 dx d y
a two-dimensional convolution1 .
Convolutions have nice properties in Fourier theory. If we designate the
Fourier transform with the symbol F , and convolution with the symbol
⊛, we may abbreviate the above equation as follows:
1
T= S ⊛ ∆ g,
4π R
1
The integration extends from minus to plus infinity in both co-ordinates x and y.
This can only be kept physically realistic on a curved Earth if it is assumed that
the kernel S is of bounded support, i.e., it differs from zero only in a bounded
area. This is the case for the modified kernels discussed in section 7.8.
Ó »îá .
Integration by FFT 165
∆ g i j = ∆ g xi , y j ,
( )
N N
x i = i δ x, i=− , . . . , − 1,
2 2
N N
y j = j δ y, j = − , . . . , − 1,
2 2
2
There exist alternatives to this. E.g., one could calculate for every grid point
the average over a square cell surrounding the point.
Ó »îá .
166 Spectral techniques, FFT
S ψ = S x − x′ , y − y′ = S (∆ x, ∆ y) ,
( ) ( )
i.e., we write
S i j = S ∆xi , ∆ y j ,
( )
3
Without this measure the result of the calculation would be correct, but in the
wrong place. . .
4
In fact, for both ∆ g and T we could choose the simpler grid geometry where the
subscript sequences are i, j = 0, · · · , N − 1; however, for S it is mandatory to have
the origin in the middle of the grid.
Ó »îá .
Solution in rectangular co-ordinates 167
This method is good for computing the disturbing potential T — and sim-
ilarly the geoid height N = Tγ — from gravity anomalies using the Stokes
equation. It is just as good for evaluating other quantities, like, e.g., the
vertical gradient of the gravity anomaly by the Poisson equation. The
only requirement is, that the equation can be expressed as a convolution.
Also the inversion calculation is easy, as we shall see: in the frequency
domain it is just a simple division.
Using the discrete Fourier transform requires that the input data in a
field to be integrated — in the example, gravity anomalies — is given
on a regular grid covering the area of computation, or can be converted
to one. The result — e.g., the disturbing potential — is obtained on a
regular grid in the same geometry. Values can then be interpolated to
chosen locations.
kommutoiva kaavio The FFT method may again be depicted as a commutative diagram:
Appendix C offers a short explanation of why FFT works and what makes
it as efficient as it is.
Ó »îá .
168 Spectral techniques, FFT
( )
φ′ , λ′ . The angular distance may be written as follows (cosine rule on
the sphere):
Substitute
( )
( ′
) 2 λ′ − λ
cos λ − λ = 1 − 2 sin ,
2
ψ
cos ψ = 1 − 2 sin2
,
2
φ′ − φ
cos φ′ − φ = 1 − 2 sin2
( )
,
2
and obtain the half-angle cosine rule:
( ′
) ′ 2 λ′ − λ
cos ψ = cos φ − φ − 2 cos φ cos φ sin =⇒
2
ψ φ ′
− φ λ′
− λ
sin2 = sin2 + cos φ′ cos φ sin2 .
2 2 2
Here we may use the following approximation:
def
cos φ′ ≈ cos φ = cos φ0 ,
ψ φ′ − φ λ′ − λ
sin2 ≈ sin2 + cos2 φ0 sin2 , (8.3)
2 2 2
def def
which depends only on the differences ∆φ = φ′ − φ and ∆λ = λ′ − λ, a
requirement for convolution.
( )5
After this, the FFT method may be applied by using co-ordinates φ, λ
and the re-written Stokes kernel
√
∆φ ∆λ
( )
) def 2
∗
∆φ, ∆λ = S 2 arcsin + cos2 φ0 sin2
(
S sin .
2 2
5
In practice one uses the geodetic or geographic latitude ϕ instead of φ without
significant error.
6
Govert L. Strang van Hees (1932 – 2012) was a Dutch gravimetric geodesist.
Ó »îá .
Solution in rectangular co-ordinates 169
where ∆φ = φ − φ′ , ∆λ = λ − λ′ and
S i ∆φ, ∆λ = S φ − φ′ , λ − λ′ ; φ i ,
( ) ( )
φ i+1 − φ φ − φi
{[ ] [ ] }
( R)
N φ, λ = Ii + I i+1 (8.5)
4πγ φ i+1 − φ i φ i+1 − φ i
with
¨
S i ∆φ, ∆λ ∆ g φ′ , λ′ cos φ′ d φ′ d λ′ ,
( )[ ( ) ]
Ii =
¨
S i+1 ∆φ, ∆λ ∆ g φ′ , λ′ cos φ′ d φ′ d λ′ .
( )[ ( ) ]
I i+1 =
ψ φ′ − φ λ′ − λ
sin2 = sin2 + cos φ′ cos φ sin2
=
2 2 2
∆φ ∆λ
= sin2 + cos φ − ∆φ cos φ sin2
( )
.
2 2
Here again, we calculate S i and S i+1 for the values φ = φ i and φ = φ i+1 ,
we evaluate the integrals with the aid of the convolution theorem, and
( )
interpolate N φ, λ according to equation 8.5 when φ i ≤ φ < φ i+1 . After
this, the solution isn’t entirely exact, because inside every band we still
use linear interpolation. However by making the bands narrower, we can
keep the error arbitrarily small.
This somewhat more complicated but also more versatile approach ex-
pands the Stokes kernel into a Taylor expansion with respect to latitude
Ó »îá .
170 Spectral techniques, FFT
C = C ϕ, ϕ′ , ∆λ = C ∆ϕ, ∆λ, ϕ .
( ) ( )
Linearize:
C = C 0 ∆ϕ, ∆λ + ϕ − ϕ0 C ϕ ∆ϕ, ∆λ + . . .
( ) ( ) ( )
) def ∂
⏐
∆ϕ, ∆λ = C ∆ϕ, ∆λ, ϕ ⏐⏐
( )⏐ (
Cϕ .
∂ϕ ϕ=ϕ0
Note that this expansion into two terms will work only for a limited
range on ∆ϕ, and that the kernel function C is assumed to be of bounded
support. In this case, the integrals may be calculated over a limited area
instead of over the whole Earth.
Substitution yields
¨
C ∆ϕ, ∆λ, ϕ · m ϕ′ , λ′ cos ϕ′ d ϕ′ d λ′ =
( ) ( ) ( )
ℓ ϕ, λ =
¨
C 0 + ϕ − ϕ0 C ϕ · m cos ϕ′ d ϕ′ d λ′ =
[ ( ) ]
=
¨
= C 0 · m cos ϕ′ d ϕ′ d λ′ +
¨
C ϕ · m cos ϕ′ d ϕ′ d λ′ .
( )
+ ϕ − ϕ0 (8.6)
7
In the literature the method has been generalized by expanding the kernel also
with respect to height.
Ó »îá .
Solution in rectangular co-ordinates 171
Important here is now, that the integrals in the first and second terms,
¨
C 0 ∆ϕ, ∆λ m ϕ′ , λ′ cos ϕ′ d ϕ′ d λ′ = C 0 ⊛ m cos ϕ ,
( )[ ( ) ] [ ]
¨
C ϕ ∆ϕ, ∆λ m ϕ′ , λ′ cos ϕ′ d ϕ′ d λ′ = C ϕ ⊛ m cos ϕ ,
( )[ ( ) ] [ ]
are both convolutions: both C functions depend only on ∆ϕ and ∆λ. Both
integrals can be calculated only if the corresponding ∆ϕ = ϕ − ϕ′ and
∆λ = λ − λ′ , and the corresponding coefficient grids C 0 , C ϕ , are calculated
first in preparation. After this (in principle expensive, but, thanks to
FFT and the convolution theorem, a lot cheaper) integration, computing
the compound 8.6 is cheap: one multiplication and one addition for each
( )
evaluation point ϕ, λ .
Example: let the evaluation area at latitude 60◦ be 10◦ × 20◦ in size. If
the grid mesh size is 5′ × 10′ , the number of cells is 120 × 120. Let
us choose, e.g., a 256 × 256 grid (i.e., N = 256) and fill the missing
values with extrapolated values.
Also the values of the kernel functions C 0 and C ϕ are calculated
on a 256 × 256 size ∆ϕ, ∆λ grid. The number of these is thus
( )
[ ]
also 65,536. Calculating the convolutions C 0 ⊛ m cos ϕ and C ϕ ⊛
[ ]
m cos ϕ by means of FFT — i.e.,
¨
C 0 ∆ϕ, ∆λ m ϕ′ , λ′ cos ϕ′ d ϕ′ d λ′ = C 0 ⊛ m cos ϕ =
( ) ( ) [ ]
{ { } { }}
= F −1 F C 0 F m cos ϕ ,
¨
C ϕ ∆ϕ, ∆λ m ϕ′ , λ′ cos ϕ′ d ϕ′ d λ′ = C ϕ ⊛ m cos ϕ =
( ) ( ) [ ]
{ { } { }}
= F −1 F C ϕ F m cos ϕ ,
C 0 = C ∆ϕ, ∆λ, ϕ0 ,
( )
C +1 − C −1
Cϕ ≈ .
ϕ+1 − ϕ−1
Ó »îá .
172 Spectral techniques, FFT
F ℓ − ℓ(0)
{ { }}
](1) ](0)
+ F −1
[ [
m cos ϕ = m cos ϕ ,
F C0
{ }
and so on, iteratively. Two, three interations are usually enough. This
method has been used to compute underground mass points from gravity
anomalies to represent the exterior gravity field of the Earth. More is
explained in Forsberg and Vermeer (1992).
D 8.3.4 “1D-FFT”
This is a limiting case of the previous ones, in which FFT is used only
in the longitude direction. In other words, a zones method where the
zones have a width of only a single grid row. This method is exact if all
longitudes (0◦ − 360◦ ) are along in the calculation. It requires a bit more
computing time compared to the previous methods. In fact, it is identical
to a Fourier transform in variable λ, longitude. Details are found in
Haagmans et al. (1993).
1. the data on the opposite side of an edge must be so far away to have
no noticeable influence on the result of the calculation, and
2. the data is continuous across the edges.
Therefore, always when using FFT with the convolution theorem 8.2, two
measures need to be taken.
8
Topologically the area with the edges thus connected is equivalent to a torus,
and the data is presupposed to be continuous on the surface of the torus.
Ó »îá .
Bordering and tapering of the data area 173
Ó »îá .
174 Spectral techniques, FFT
Figure 8.3. Example images for FFT transform without (above) and with
(below) tapering. Used on-line FFT http://www.ejectamenta.com/
D Imaging-Experiments/fourierimagefiltering.html.
9
Carl Christian Tscherning (1942 – 2014) was a Danish physical geodesist well
known for his research into the gravity field of the Earth. He did groundbreaking
work on statistical computation methods for modelling the Earth’s gravity field
from many different measurement types.
Ó »îá .
Computing a geoid model with FFT 175
20
70˚ 70˚
19
68˚ 68˚
24 25
22
21
23
31
28
29
30
20
26
27
66˚ 19 66˚
17
18
18
64˚ 64˚
18
5
23 24 2
20 212
2
19
62˚ 62˚
18
17
19 16
60˚ 60˚
15
16
20˚ 32˚
24˚ 28˚
Figure 8.4. The Finnish FIN2000 geoid. Data source: Finnish Geodetic Insti-
D tute.
Currently two geoid models are in use in Finland: FIN2000 (figure 8.4)
and FIN2005N00 (Bilker-Koivula and Ollikainen, 2009). The first model
is a reference surface for the N60 height system: using it together with
GNSS positioning allows determination of the N60 heights of points.
Ó »îá .
176 Spectral techniques, FFT
The model gives geoid heights above the GRS80 reference ellipsoid. The
second model is similarly a reference surface for the new N2000 height
system. It too gives heights from the GRS80 reference ellipsoid.
The precisions (mean errors) of FIN2000 and FIN2005N00 are on the
level of ±2 − 3 cm.
The Danish researchers Per Knudsen and Ole Balthasar Andersen have
computed a gravity map of the world ocean by starting from satellite
altimetry derived “geoid heights” and inverting them to gravity anomalies
(Andersen et al., 2010). A pioneer of the method has been David Sandwell
from the Scripps Institute of Oceanography in California (e.g., Garcia
et al., 2014). The short-wavelength features in the map can tell us about
the sea-floor topography.
Also the data from satellite gravity missions (like CHAMP, GRACE and
GOCE) can be regionally processed using the FFT method: in the case
of GOCE, the inversion of gradiometric measurements, i.e., calculating
geoid heights on the Earth’s surface from measurements made at satellite
level. Also airborne gravity measurements are processed in this way. The
problem is called “harmonic downward continuation” and is in principle
unstable.
Airborne gravimetry is a practical method for gravimetric mapping of
large areas. In the pioneering days, the gravity field over Greenland
was mapped, as well as many areas around the Arctic and Antarctic.
Later, areas were measured like the Brazilian Amazonas, Mongolia and
Ethiopia (Bedada, 2010), where no full-coverage terrestrial gravimetric
data existed. The advantage of this method is that one measures rapidly
large areas in a homogeneous way. Also for the processing of airborne
gravimetry data, FFT is suitable.
Ó »îá .
Computing terrain corrections with FFT 177
is the cosine of the angle θ between the force vector — assumed coming
from the midpoint of the rock column — and the vertical direction. This
is the so-called prism method.
We will make a linear approximation, wherein ℓ, the slant distance
( )
between the evaluation point ( x, y) and the moving data point x′ , y′ , is
also the horizontal distance:
)2 ( )2
ℓ2 ≈ x ′ − x + y′ − y
(
.
Ó »îá .
178 Spectral techniques, FFT
Then, the terms in the above sum are large numbers that almost cancel
each other, giving a nearly correct result. Numerically this is however an
unpleasant situation.
If ℓ is defined according to equation 8.9, then the Fourier transform of
kernel ℓ−3 is (Harrison and Dickinson, 1989):
2π 2π 4π2 q2 δ2
{ } [ ]
−3
F ℓ = exp (−2πδ q) = 1 − 2πδ q + −... ,
δ δ 1·2
def p
in which q = u2 + v2 , and u, v are wave numbers (i.e., “frequencies”)
in the x and y directions in the ( x, y) plane. If we substitute this into
equation 8.8, we notice that the terms containing δ−1 sum to zero, and of
course also the terms containing positive powers of δ vanish when δ → 0.
As follows (Harrison and Dickinson, 1989):
{ } 1 { }[ 2π ]
2
F TC ≈ Gρ H F 1 (1 − 2πδ q) −
2 δ
{ }[ 2π ]
′
− Gρ HF H (1 − 2πδ q) +
δ
2 [ 2π
{( ) }
1 ]
+ GρF H′ (1 − 2πδ q)
2 δ
2πG ρ [ 1 2 1 ( ′ )2 ]
TC = H − H′ H + H +
δ 2 { {2 } }
+G ρ H F −1 F H ′ · 4π2 q −
1 { {( ) }
2
}
− G ρ F −1 F H′ · 4π 2 q .
2
In the first term
1 2 1 ( ′ )2 1 ( ′ )2
H − H′ H + H = H −H =0
2 2 2
Ó »îá .
Self-test questions 179
D Self-test questions
Ó »îá .
D Statistical methods
9
∂
⏐
⏐
f→
↦ f ( x)⏐⏐ .
∂x x = x0
– 181 –
182 Statistical methods
f ↦→ f ( x0 ) .
and so on.
We may write symbolically
∂ ⏐⏐ ∂ f ⏐⏐
⏐ ⏐
L= , i.e., L { f } = .
∂ x ⏐ x= x0 ∂ x ⏐ x= x0
Note that all partial derivatives, as also the Laplace operator ∆, are
linear.
In physical geodesy, all interesting functionals are functionals of the
( ) ( )⏐
function T φ, λ, R = T φ, λ, r ⏐r=R , i.e., the disturbing potential on the
Earth’s surface. The theory thus uses the spherical approximation1 ,
and the surface of the sphere, radius R , corresponds to mean sea level.
def ( )
E.g., the disturbing potential T P = T φP , λP , R in a point P at sea-level
( )
location φP , λP is such a functional:
( )
T (·, ·, R ) ↦→ T φ, λ, R .
Even if the quantity is not the disturbing potential, but, e.g., the gravity
anomaly or the deflection of the plumb line:
T (·, ·, R ) ↦→ ∆ g φ, λ, r ,
( )
( )
T (·, ·, R ) ↦→ ξ φ, λ, r ,
( )
T (·, ·, R ) ↦→ η φ, λ, r .
1
This is not mandatory, but the error of approximation is small.
Ó »îá .
Statistics on the Earth’s surface 183
T (·, ·, R ) ↦−→ a nm ,
T (·, ·, R ) ↦−→ b nm .
we may also do so for a stochastic process. The only difference is, that by
doing so we obtain functions.
Let, e.g., the stochastic process x ( t) be a function of time. Then we may
compute its variance function as follows:
{ }
C xx ( t) = Var x ( t) .
For a stochastic process however, much more can be computed: e.g., the
covariance of values of the same function taken at different points in
time, the so-called autocovariance:
{ }
A xx ( t 1 , t 2 ) = Cov x ( t 1 ) , x ( t 2 ) =
{[ { }][ { }]}
=E x ( t1 ) − E x ( t1 ) x ( t2 ) − E x ( t2 ) .
Ó »îá .
184 Statistical methods
and again, and one can practice the art of statistics on the results of
the throws. Another classic example is measurement. Measurement of
the same quantity can be repeated, and is repeated, in order to improve
precision.
For a stochastic process defined on the Earth’s surface, the situation is
different.
We have only one Earth.
For this reason, statistics must be done in a somewhat different fashion.
( )
Given a stochastic process on the surface of the Earth, x φ, λ , we define
a quantity similar to the statistical expectancy E {·}, the geographic mean
¨ ˆ 2πˆ +π/2
def 1 ( ) 1 ( )
M { x} = x φ, λ d σ = x φ, λ cos φ d φ d λ. (9.1)
4π σ 4π 0 −π/2
( )
Here x φ, λ is the one and only realization of process x that we have
available on this Earth.
Clearly this definition makes sense only in the case where the statisti-
( )
cal behaviour of the process x φ, λ is the same everywhere on Earth,
( )
independently of the value of φ, λ . This is called the assumption of
homogeneity. It is in fact the assumption that the spherical symmetry of homogeenisuus
the Earth extends to the statistical behaviour of her gravity field.
Similarly to the statistical variance based on expectancy, we may define
the geographic variance:
) def )} def
= M [ x − M { x}]2 .
( { ( { }
C xx φ, λ = Var x φ, λ (9.2)
The global average of gravity anomalies ∆ g φ, λ vanishes based on their
( )
definition:
M {∆ g} = 0.
In that case, equation 9.2 is simplified as follows:
¨
2 1 )]2
C ∆ g∆ g φ, λ = Var ∆ g φ, λ = M ∆g ∆ g φ, λ
( ) { ( )} { } [ (
= d σ.
4π σ
The definition given here of the geographic mean M {·} is based on inte-
gration of the one and only realization over the surface of the Earth. As
has been seen, in statistics the mean is defined slightly differently, as
the expectancy of a stochastic process. For gravity anomalies this means
E ∆ g , where ∆ g is the anomaly considered as a stochastic process, i.e.,
{ }
Ó »îá .
The covariance function of the gravity field 185
α
P
Q
∆ g P = ∆ g φP , λP ,
( )
∆ g Q = ∆ g φQ , λQ .
( )
angular distance ψPQ and the azimuth angle αPQ . See figure 9.1.
Now we may write
2
This is called the geodetic forward problem on the sphere.
Ó »îá .
186 Statistical methods
ˆ 2π
1
′
∆ g P ∆ g Q(P) M ∆ g P ∆ g Q(P) d αPQ =
( ) { } { }
C ∆ g∆ g ψPQ = M =
2π 0
ˆ 2πˆ +π/2ˆ 2π
1
= ∆ g P ∆ g Q(P) d λP cos φP d φP d αPQ . (9.3)
8π2 0 −π/2 0
Remark. The true gravity field of the Earth isn’t terribly homogeneous
or isotropic, but in spite of this, both hypotheses are widely used.
Ó »îá .
Least-squares collocation 187
variance matrix):
⎡ ⎤
C ( t1 , t1 ) C ( t2 , t1 ) C ( t1 , t N )
···
⎢ .. ⎥
{ } ⎢ C ( t1 , t2 ) C ( t2 , t2 ) · · · . ⎥
Var s i =⎢ .. .. .. ..
⎥.
.
⎢ ⎥
⎣ . . . ⎦
C ( t1 , t N ) C ( t2 , t N ) · · · C ( t N , t N )
( )
We use for this the symbol C i j , both for one element C i j = C t i , t j of
[ ( ) ]
the matrix, and for the whole matrix, C i j = C t i , t j , i, j = 1, . . . , N . The
symbol s i again denotes a vector s ( t i ) , i = 1, . . . , N consisting of process
values — or one of its elements s ( t i ).
Note that, if the function C t i , t j or C (∆ t) is known, then the whole
( )
matrix and all of its elements can be calculated if only all argument
values (observation times) t i are known.
Let the shape of the problem now be, that one should estimate, i.e., predict,
the value of process s at the moment in time T , i.e., s (T ), based on our
knowledge of the above described observations s ( t i ) , i = 1, . . . , N .
In the same way as we calculated above the covariances between s ( t i )
( )
and s t j (elements of the signal variance matrix C i j ), we also compute
the covariances between s (T ) and all s ( t i ) , i = 1, . . . , N . We obtain
{ } [ ]
Cov s (T ) , s ( t i ) = C (T, t 1 ) C (T, t 2 ) · · · C (T, t N ) .
For this we may again use the notation C T j . It is assumed here, that there
is only one point in time T for which estimation is done. Generalization
to the case where there are several T p , p = 1, . . . , M , is straightforward.
In that case, the signal covariance matrix will be of size M × N :
⎡ ⎤
C (T1 , t 1 ) C (T1 , t 2 ) ··· C (T1 , t N )
{ ( ) } ⎢ C (T2 , t 1 ) C (T2 , t 2 ) ··· C (T2 , t N ) ⎥
Cov s T p , s ( t i ) = ⎢ ⎥.
⎢ ⎥
.. .. .. ..
⎣ . . . . ⎦
C (T M , t 1 ) C (T M , t 2 ) · · · C (T M , t N )
ℓi = s (t i ) + ni . (9.4)
Ó »îá .
188 Statistical methods
ŝ T p − s T p = Λ pi ℓ i − s T p = Λ pi s ( t i ) + n i − s T p .
( ) ( ) ( ) [ ] ( )
∑
Here, for the sake of writing convenience, we left the summation sign
away (Einstein summation convention): We always sum over adjacent,
identical indices, in this case i .
Study the variance of this difference, the so-called variance of prediction:
def
Σ p p = Var ŝ T p − s T p
{ ( ) ( )}
.
and
Σ pq = Cov
{[ ( ) ( )] [ ( ) ( )]}
ŝ T p − s T p , ŝ T q − s T q =
= Λ pi Cov ΛTjq + Cov s T p , s T q
{[ ] [ ( ) ]} { ( ) ( )}
s (t i ) + ni , s t j + n j −
− Λ pi Cov s ( t i ) , s T q ΛTjq =
{ ( )} { ( ) ( )}
− Cov s T p , s t j
= Λ pi C i j + D i j ΛTjq + C pq − Λ pi C iq − C p j ΛTjq .
( )
(9.5)
Ó »îá .
Least-squares collocation 189
Here we show that the optimal estimator is indeed the one producing the
minimum possible variances.
Choose
def )−1
Λ p j = C pi C i j + D i j
(
.
Then, from equation 9.5, and exploiting the symmetry of the C and D
matrices, we obtain
)−1
Σ p p = C pi C i j + D i j
(
C j p + C p p−
( )−1 ( )−1
− C pi C i j + D i j C j p − C pi C i j + D i j C jp =
( )−1
= C p p − C pi C i j + D i j C j p. (9.6)
+ δΛ pi C i j + D i j δΛTj p − δΛ pi C i p − C p j δΛTj p =
( )
)−1
C j p + δΛ pi C i p + C p j δΛTj p −
(
= C p p − C pi C i j + D i j
− δΛ pi C i p − C p j δΛTj p +δΛ pi C i j + D i j δΛTj p =
( )
)−1
C j p + δΛ pi C i j + D i j δΛTj p .
( ( )
= C p p − C pi C i j + D i j
Here, the last term — the only difference with result 9.6 — is positive,
because the matrices C i j and D i j are positive definite: Σ′p p > Σ p p , except
when δΛ pi = 0. In other words, the solution given above,
)−1 )−1
Λ p j = C pi C i j + D i j
( ( ) (
=⇒ ŝ T p = C pi C i j + D i j ℓ j,
Ó »îá .
190 Statistical methods
{ }
Cov ∆ g P , ∆ gQ = M ′ ∆ g P ∆ g Q(P) = C ψPQ .
{ } ( )
Warning: The Hirvonen covariance function is meant for use with (free-
air) gravity anomalies, i.e., quantities obtained by subtracting nor-
mal gravity from the measured gravity. Nowadays anomalies are
often obtained by subtracting from the observations a high-degree
“normal field”, i.e., a spherical-harmonic expansion. Then one uses
Hirvonen’s formula at one’s own risk!
Alternative functions that are also often used in local applications are the
covariance functions of first and second order autoregressive processes,
or AR(1) or AR(2) -processes:
C ψ = C 0 e−
ψ/ψ
C ψ = C 0 e−(
ψ/ψ )2 .
( ) ( )
0 or 0
3
Reino Antero Hirvonen (1908 – 1989) was a Finnish physical and mathematical
geodesist.
Ó »îá .
Least-squares collocation 191
1
0.9
0.8
0.7
( ) 0.6
C ψ
0.5
0.4
0.3
0.2
0.1 4
2
−4 0
−2 0 −2 y
x 2 4 −4
20
18
16
14
∆ g 12
10
8
6
4
2
35 40
30
y
25
20 25 30 35
15 15 20
10 10 x
5 5
Figure 9.3. An example of least-squares collocation. Here are given two data
points (stars); the surface plotted gives the estimated value ∆
ˆg P for
each point P in the area. Here we use least-squares collocation for
D inter- and extrapolating gravimetric data.
Ó »îá .
192 Statistical methods
ΣPQ ≈ C PQ − C P i C −i j1 C jQ .
ΣPP = C 0 − C P i C −i j1 C jP .
See figure 9.4. Given two points where gravity has been measured and
gravity anomalies calculated: ∆ g1 = 15 mGal, ∆ g2 = 20 mGal. The co-
ordinates in the x and y directions are in kilometres. Assumed is that
Ó »îá .
Least-squares collocation 193
y
1 (15 mGal)
30
P′
2 (20 mGal)
20
10 P
x
10 20 30
1000 mGal2
C 12 = C 21 = = 666.66 . . . mGal2 ,
1 + 200
400
s21P = (30 − 10)2 + (20 − 10)2 km2 = 500 km2 ,
( )
1000 mGal2
C 1P = = 444.44 . . . mGal2 ,
1 + 500
400
s22P = (20 − 10)2 + (30 − 10)2 km2 = 500 km2 ,
( )
1000 mGal2
C 2P = 500
= 444.44 . . . mGal2 .
1+ 400
We also have
mGal2 .
[ ] [ ]
CP i = C P1 C P2 = 444.44 444.44
Ó »îá .
194 Statistical methods
∆ g1
[ ] [ ]
15
∆g j = = mGal,
∆ g2 20
i.e., √
σ∆ g P = ΣPP = ±27.622 mGal.
Summarizing the result:
∆
ˆg P = 9.333 ± 27.622 mGal.
We may observe here, that the gravity anomaly estimate found is much
smaller than its own uncertainty, and thus does not differ significantly
from zero. In fact, not using the observational data at all would leave us
with the a priori estimate
p
∆
ˆg P = 0 ± 1000 mGal = 0 ± 31.623 mGal,
almost as good.
If, instead, we would choose to locate point P ′ in between points 1 and 2, at
1000 mGal2
location (25 km, 25 km) , then C P ′ 1 = C P ′ 2 = 50
= 888.89 mGal2
1 + 400
and ∆ˆg P ′ = 18.667 ± 7.201 mGal, which is clearly better than the a priori
estimate of zero.
And if we had chosen instead the Gauss–Markov covariance function
C = C 0 e− /d ,
s
we would have obtained the results ∆ ˆg P = 7.663 ± 29.272 mGal for the
original point location, and ∆
ˆg P ′ = 16.460 ± 18.426 mGal for the shifted
point location.
Ó »îá .
Least-squares collocation 195
ℓi = g i + ni ⇐⇒ ℓ = g + n.
Ó »îá .
196 Statistical methods
Here are several points j where gravity is given: let us say, N observa-
tions ℓ j = ∆ g j + n j , j = 1, . . . , N . The number of points to be predicted may
be one, P , or also many. The matrices C i j and D i j are square, and the
inverse of their sum exists. C P i is a rectangular matrix. If there is only
one point P , it is a size 1 × N row matrix.
The prediction error is now the difference quantity5 ∆
ˆg P − ∆ g , and its
P
variance (“variance of prediction”) is ennustusvarianssi
{ }
def
ΣPP = Var ∆
ˆg P − ∆ g =
P
{} { } { }
= Var ∆ g P + Var ∆ g P − Cov ∆ g P , ∆ g P − Cov ∆ g P , ∆ g P .
{ }
ˆ ˆ ˆ
and
{ } { }
Cov ∆
ˆg P , ∆ g
P
= Cov ∆ g P , ∆
ˆg P = C P i (C i j + D i j )−1 C jP .
Here, C T T T
iP (or C jP , or C ℓP ) is the transpose of C P i . The matrix
( )−1
Ci j + Di j is symmetric and its own transpose.
The end result is (remember that the signal variance Var ∆ g P = C PP ):
{ }
ΣPP = C PP + C P i (C i j + D i j )−1 C jP −
− C P i (C i j + D i j )−1 C jP − C P i (C i j + D i j )−1 C jP =
= C PP − C P i (C i j + D i j )−1 C jP .
ΣPP ≈ C PP − C P i C −i j1 C jP . (9.11)
Borderline cases:
Note that here, ∆ g P is the true value of the gravity anomaly at point P, which
5
Ó »îá .
Covariance function and degree variances 197
We write this in the following form using the definition of M ′ {·}, equation
9.3:
K ψPQ = M ′ T P TQ(P) =
( ) { }
ˆ 2π ˆ +π/2 ˆ 2π
1
= 2 T P TQ(P) d λP cos φP d φP d αPQ . (9.12)
8π 0 −π/2 0
astevarianssit The coefficients k n are called the degree variances (of the disturbing
6
because { ( )
( ) P nm cos ψ cos mα m ≥ 0,
Ynm α, ψ = ( )
P n|m| cos ψ sin | m| α m < 0,
Ó »îá .
198 Statistical methods
( )
potential). For isotropic covariance functions K ψ the information
content of the degree variances k n , n = 2, 3, . . . is the same as that of the
function itself, and is in fact its spectral representation.
ˆ +π/2ˆ 2π {ˆ πˆ 2π }
2n + 1 ( )
kn = TP TQ(P) d αPQ P n cos ψPQ sin ψPQ d ψPQ d λP cos φP d φP .
16π2 −π/2 0 0 0
Ó »îá .
Propagation of covariances 199
we obtain
∞ ∞ ( ¨ )
∑ ∑ 1
T n2 d σ
( ) ( ) ( )
K ψ = k n P n cos ψ = · P n cos ψ =
4π σ
n=2 n=2
∞
( n ) ∞
( n
)
∑ ∑[ 2
] ∑ ∑
a2nm + b nm a2nm
( ) ( )
= P n cos ψ = P n cos ψ .
n=2 m=0 n=2 m=− n
i.e.,
The degree variances k n of the disturbing potential can be
calculated directly from the spherical-harmonic coefficients.
The literature offers many alternative notations for the degree variances,
like
def def
k n = σ2n = σTT
i .
Ó »îá .
200 Statistical methods
that
∞ ( ) n+1
( ) ∑ R ( )
T φ, λ, r = T n φ, λ .
r
n=2
Symbolically
( ) { ( )}
T φ, λ, r = L T φ, λ, R ,
Here, L is the linear operator
∞ ( ) n+1
∑ R
L{f } = f n,
r
n=2
Symbolically
∞
∑
L{f } = Ln f n,
n=2
where ( )n+1
n R
L =
r
is the spectral representation of the operator L.
( )
We may still write in a certain point P φP , λP , r P in space:
∞
∑
LP { f } = L nP f n ,
n=2
in which ( )n+1
R
L nP = .
rP
( )
Concretely, for the disturbing potential T φP , λP , r P in point P , this
means
( ) { ( )}
T φP , λP , r P = L P T φ, λ, R =
∞
∑
= L nP T n =
n=2
∞ ( )n+1
∑ R
= Tn.
rP
n=2
Ó »îá .
Propagation of covariances 201
Thus we obtain7
∞
∑
L nP L nQ k n P n cos ψPQ =
( ) ( )
K r P , r Q , ψPQ =
n=2
∞ ( )n+1 ( )n+1
∑ R R ( )
= k n P n cos ψPQ =
rP rQ
n=2
∞ ( )n+1
∑ R2 ( )
= k n P n cos ψPQ . (9.15)
r P rQ
n=2
We know (equation 4.8) that there exists the following relationship be-
tween gravity anomalies and the disturbing potential:
∞ ( )n+1
1∑ R
∆g = ( n − 1) T n ,
r r
n=2
7
This works only this cleanly because in this case the operator L n is of multiplier
( )n+1
type, Rr .
Ó »îá .
202 Statistical methods
∆ g φP , λP , r P = L g,P T φ, λ, R
( ) { ( )}
=
∞
∑
= L ng,P T n =
n=2
∞ ( )n+1
∑ n−1 R
= Tn.
rP rP
n=2
n=2
∞ ( )n+1 ( )n+1
∑ n−1 R n−1 R ( )
= k n P n cos ψPQ =
rP rP rQ rQ
n=2
∞ ( )n+2 (
∑ R2 n − 1 )2 ( )
= k n P n cos ψPQ .
r P rQ R
n=2
Often we write
∞ ( )n+2
) def ∑ R2
C ψPQ , r P , r Q = M ∆ g P ∆ g Q =
( { } ( )
c n P n cos ψPQ ,
r P rQ
n=2
n=2
∞ ( )n+1 ( )n+1
∑ R n−1 R ( )
= k n P n cos ψPQ =
rP rQ rQ
n=2
∞ )n+1
R2
(
∑ n−1 ( )
= k n P n cos ψPQ .
rQ r P rQ
n=2
series expansion:
∑
L n1,P L n2,Q M {T n T n } =
{ { } { }}
Cov L 1 T P , L 2 T Q =
n
∑
L n1,P L n2,Q k n P n cos ψPQ ,
( )
=
n
k n = α n−4 .
By writing
( n − 1 )2
cn =
kn,
R
where c n are the degree variances of gravity anomalies, we obtain
α ( n − 1)2 α
cn = ≈ n−2 .
R2 n4 R2
Here, Rα2 is a planet specific constant, value about 1200 mGal2 for the
Earth.
The Kaula rule does not hold very precisely for very high degree numbers.
It applies, by the way, fairly well for the gravity field of Mars, of course
with a different constant (Yuan et al., 2001).
Another well known rule is the Tscherning–Rapp equation (Tscherning
and Rapp, 1974):
A ( n − 1) ( n − 1 )2
cn = = kn.
( n − 2) ( n + B) R
8
William M. Kaula (1926 – 2000) was an American geophysicist and space geode-
sist who studied the determination of the Earth’s gravity field by means of
satellite geodesy.
Ó »îá .
204 Statistical methods
1010
Kaula EGM96
108 EGM2008
Tscherning–Rapp
GOCE, Gatti et al. (2014)
106 EGM96 error variances
EGM2008 error variances
104 GOCE error variances
)
m4 /s4
(
102
Degree variance k n
100
10−2
10−4
10−6
10−8
0 50 100 150 200 250 300 350 400
Degree n
Figure 9.5. Global covariance functions as degree variances. The GOCE model
D cuts off at degree 280.
9
Arne Bjerhammar (1917 – 2011) was an eminent Swedish geodesist.
Ó .
»îá
Collocation and the spectral viewpoint 205
( )
f ψ are desired in the same points. Then equation 9.9 yields
[ ]−1 ( )
f̂ = C f g C g g + D g g g+n (9.16)
with
[ ] ( ( ) ( )) ( )
Cf g ij
= C f g f ψi , g ψ j = C f g ψi , ψ j ,
[ ] ( ( ) ( )) ( )
Cgg ij
= C g g g ψi , g ψ j = C g g ψi , ψ j ,
[ ] ( ( ) ( )) ( )
D gg ij
= D g g g ψi , g ψ j = D g g ψi , ψ j .
If the physics of the whole situation, including the physics of the mea-
surement process, is rotationally symmetric, we must have
N −1
[ ] { ( ) ( )} 1 ∑ ( ) ( )
Cf g ij
= M◦ f ψ i g ψ j(i) = f ψ i g ψ j(i) ,
N
i =0
which, like the geographic average in section 9.4, replaces the statistical
average.
In the same way we obtain
N −1
[ ] { ( ) ( )} 1 ∑ ( ) ( )
Cgg ij
= M◦ g ψ i g ψ j(i) = g ψ i g ψ j(i) .
N
i =0
= C f g ψ i , ψ j = C f g ∆ψk = C f g [ k] ,
[ ] ( ) [ ]
Cf g ij
= C g g ψ i , ψ j = C g g ∆ψk = C g g [ k] ,
[ ] ( ) [ ]
Cgg ij
def
in which ∆ψk = ψ j − ψ i mod 2π and k = ( j − i ) mod N .
( )
Furthermore
= D g g ψ i , ψ j = D g g ∆ψk = D g g [ k] = E n i n j(i) ,
[ ] ( ) [ ] { }
D gg ij
D g g = σ2 I N ,
10
In fact, the unit or identity matrix is also known as the Kronecker delta, and
as a Toeplitz matrix may be interpreted as a discrete version of the Dirac delta
function. Its discrete Fourier transform is
F {I } = 1
(and yes, this extends to infinite frequencies, and contains an infinite amount of
power. . . ).
Ó »îá .
206 Statistical methods
i 2
1
0
N −1
j N −2
∆ψk
ψj
ψi
Without proof we present that the spectral version of equation 9.16 looks
like this:
F Cf g
{} { }
F f̂ }F g +n =
{ }
=
F Cgg + F D gg
{ } {
F Cf g
{ }
F g+n .
{ }
= (9.17)
F Cgg + σ
{ }
2
This is an easy and rapid way to calculate the solution using FFT. In
the limit in which the observations are exact, i.e., σ2 = 0, by equation
9.17 f̂ follows straight from g + n = g. If for a suitable operator L we have
f = L {g}, the equation simplifies as follows:
{ } F { L } F {C } {
gg
F f̂ = F
}
g + n ,
F Cgg + σ
{ }
2
E.g., if the g are gravity anomalies and the f are values of the disturbing
potential, then12
R
F {L} = .
n−1
11
Otto Toeplitz (1881 – 1940) was a German Jewish mathematician who con-
tributed to functional analysis.
12
In real computation it is not so simple. . . the degree number n, which refers to
global spherical geometry, must first be converted to the Fourier wave number
expressed on the computational grid used.
Ó »îá .
Self-test questions 207
The approach is called Fast Collocation, e.g., Bottoni and Barzaghi (1993).
Of course it is used in two dimensions on the Earth’s surface, though our
example is one-dimensional. As always, it requires that the observations
are given on a grid, and in this case also, that the precision of the
material is homogeneous — the same everywhere — over the area. This
requirement is hardly ever precisely fulfilled.
D Self-test questions
Ó »îá .
208 Statistical methods
6. Why, in the study of the Earth’s gravity field, one uses as the
average of quantities the geographical average rather than the
statistical average?
7. Which two different kinds of covariance functions are used for
gravity anomalies on the Earth’s surface? Give the formulas and
name the free parameters.
8. Explain degree variances. What is the difference between degree
variances k n and c n ?
9. Describe Kaula’s rule.
10. What is a Toeplitz circulant matrix?
ΣPP = C PP − C P i (C i j + D i j )−1 C jP ,
in which the observation points are i = 1, . . . , N . Assume there is only one
observation point, point P . Then
ΣPP = C PP − C PP (C PP + D PP )−1 C PP .
ΣPP ≈ D PP .
Ó »îá .
Exercise 9 – 3: Predicting gravity anomalies 209
2 1
σQQ = CQQ − CQP C −
PP C PQ .
Cov ∆ g P , ∆ g1 Cov ∆ g P , ∆ g2
[ { } { } ]
CP i =
σ2PP = C PP − C P i C − 1
i j C jP .
Ó »îá .
210 Statistical methods
1. Compute C i j , C P i , ∆
ˆg P and σ2 .
PP
2. Compare the result with the previous one. Conclusion?
{ } ∞
∑
L nδ g,P L nδ g,Q k n P n cos ψPQ .
( )
Cov δ g P , δ gQ =
n=2
∂2 T
2. Compute the covariance function of the gravity gradient (i.e.,
∂r2
the vertical gradient of the gravity disturbance!).
k n = α n−4 .
Ó »îá .
Exercise 9 – 7: Underground mass points 211
From these one can derive, using propagation of variances, the degree
variances of gravity anomalies
∞ ∞ (
∑ ∑ n −1)
∆g = L ng T n = Tn
R
n=2 n=2
as follows:
)2 ( n − 1 )2 α
c n = L ng n−2 .
(
kn = kn ≈
R R2
By differentiating the above expansion 9.18 for the disturbing potential
∞ ( ) n+1
∑ R
T ( r, ϕ, λ) = T n (ϕ, λ)
r
n=2
the connection between the disturbing potential and the gravity gradient
in the spectral domain.
On the Earth’s surface r = R , or
∞ ∞
∂2 T ⏐⏐
⏐
∑ ( n + 1) ( n + 2) def
∑
= Tn = L ng g T n
∂ r 2 ⏐ r =R R2
n=2 n=2
with
( n + 1) ( n + 2)
L ng g = .
R2
1. Derive an (approximate) equation for the degree variances for the
gravity gradient, which we may call g n , in an analogue fashion as
above for the gravity anomaly degree variances c n :
g n =? · k n ≈? · n? .
2. Conclusion?
Ó »îá .
D Gravimetric measurement devices
10
D 10.1 History
The first device ever built based on a pendulum was a clock. The pendu-
lum equation, √
ℓ
T = 2π ,
g
tells that the swinging time T of a pendulum of a given length is a con-
stant that depends only on the length ℓ and local gravity g, on condition
that the swings are small. The Dutch Christiaan Huygens1 built in 1657
the first useable pendulum clock based on this (http://en.wikipedia.org/
wiki/Pendulum_clock).
When the young French researcher Jean Richer2 visited French Guyana
in 1671 with a pendulum clock, he noticed that the clock ran clearly
slower. The matter was corrected simply by shortening the pendulum.
The cause of the effect could not be the climatic conditions in the tropics,
i.e., the thermal expansion of the pendulum. The right explanation was
that in the tropics, gravity g is weaker than in Europe. After return to
France, Richer had again to make his pendulum longer. The observation
is described in just one paragraph on pages 87 – 88 in his report “Obser-
vations astronomiques et physiques faites en l’isle de Caïenne”, Richer
(1731).
This is how the pendulum gravimeter was invented. Later, much more
precise special devices were built, e.g., Kater3 ’s reversion pendulum,
and the four-pendulum Von Sterneck4 device, which was also used in
1
Christiaan Huygens (1629 – 1695) was a leading Dutch natural scientist and
mathematician. Besides inventing the pendulum clock, he also was the first to
realize (in 1655) that the planet Saturn has a ring.
2
Jean Richer (1630 – 1696) was a French astronomer. He is really only remem-
bered for his pendulum finding.
3
Henry Kater (1777 – 1835) was an English physicist.
4
Robert von Sterneck (1839 – 1910) was an Austrian-Hungarian scientist.
– 213 –
214 Gravimetric measurement devices
Finland in the 1920’s and 1930’s. We must mention also the submarine
measurements, e.g., in the Java Sea by the Dutch F. A. Vening Meinesz in
which it was observed that above the trenches in the ocean floor there is syvänmeren hauta
a notable shortage of gravity, and that they thus are in a state of strong
isostatic disequilibrium (Vening Meinesz, 1928).
For production gravimetric observations, pendulum gravimeters are how-
ever too hard to operate and too slow. For that purpose the spring gravime-
ter has been developed, see section 10.2.
Pendulum gravimeters are in principle absolute measurement devices,
i.e., gravity is obtained directly as an acceleration. There are, however,
systematic effects associated with the suspension of the pendulum that
make that one cannot trust in the absoluteness of measurement after
all. One tried out solution is the very long wire pendulum, e.g., Hytönen
(1972). However, nowadays absolute measurements are made with bal-
listic gravimeters, cf. section 10.3. It has been observed that the older
measurements made with pendulum apparatus in the so-called Potsdam
system are systematically 14 mGal too large. . .
Ó »îá .
The relative or spring gravimeter 215
Figure 10.2. Autograv CG5 spring gravimeter. Image © 2011 David Monniaux,
D Wikimedia Commons, GNU Free Documentation License.
in which ℓ is the mean length of the spring during the oscillation, and
also the equilibrium length in the absence of oscillations.
When the test mass is disturbed, it starts oscillating about its equilibrium
värähtely-yhtälö position. The oscillation equation, obtained by summing the above two
equations, is
d2 ( ) k( )
ℓ − ℓ = − ℓ − ℓ .
dt2 m
The period is √ √
ℓ − ℓ0 δℓ
√
m
T = 2π = 2π = 2π , (10.2)
k g g
Ó »îá .
216 Gravimetric measurement devices
D 10.2.1 Astatization
Gravity pulling at the mass again is mg, and the corresponding torque
τg = mg p cos ϵ.
By differentiation
mpℓ d g + mg p d ℓ − kbc d ℓ = 0
5
https://en.wikipedia.org/wiki/Spring_(device)#Zero-length_springs.
Ó »îá .
The relative or spring gravimeter 217
Reality
Working range
Force
Hooke’s law
Spring, length ℓ
c ( )
k ℓ − ℓ0 sin β
Length
mg
dℓ
= 2.5 · 10−6 m /mGal,
dg
but for a state of disequilibrium. Then, the test mass will be undergoing
an acceleration a, and we have
m ( g − a) pℓ − kbc (ℓ − ℓ0 ) = 0,
6 p
For comparability we should still multiply by b sin β , if we measure the position
of the test mass.
7
Lucien LaCoste (1908 – 1995) was an American physicist and metrologist, who,
as an undergraduate, together with his physics professor Arnold Romberg (1882 –
1974) discovered the principle of the astatized gravimeter and zero-length spring.
Ó »îá .
218 Gravimetric measurement devices
8
Plastic deformation in a metal crystal is mediated by crystal-lattice defects
called dislocations. As dislocations travel through the crystal lattice under load,
the properties of the metal change, which may eventually result in metal fatigue,
a known problem, i.a., in aviation. https://en.wikipedia.org/wiki/Dislocation.
The art of making metals stronger by inhibiting the motion of dislocations,
e.g., by adding carbon to iron to form steel, forms a large part of metallurgy.
https://en.wikipedia.org/wiki/Strengthening_mechanisms_of_materials.
Ó »îá .
The relative or spring gravimeter 219
δ (ε)
)
F (ε
α F (ε) cos (α + δ + ε)
−ε
mg
.
mg
Figure 10.4. The idea of astatization. The elastic force of an ordinary spring
grows steeply with extension (left), whereas the weight of the
test mass is constant. The lever beam and diagonal arrangement
(right) causes the part of the force of the spring in the direction
of motion of the lever (red) to diminish with extension, while
the spring force itself grows similarly with extension. This near-
cancellation boosts sensitivity. The spring used is a zero-length
D spring.
Ó »îá .
220 Gravimetric measurement devices
Cage transporter
Falling prism
g
”Superspring”
Reference prism
Semi-transparent mirror
Laser
Mirror
9
James E. Faller (1934 – ) is an American physicist, metrologist, geodesist and
student of gravitation. He proposed the installation of laser retroreflectors on
the lunar surface in the context of the Apollo project, in order to measure the
distance to the Moon — LLR, lunar laser ranging.
Ó »îá .
The absolute or ballistic gravimeter 221
Figure 10.6. FG5 absolute gravimeter. Figure © National Oceanic and Atmos-
pheric Administration.
D
d2
z = g ( z) ,
dt2
where it is assumed — realistically — that gravity g depends on the
location z within the drop tube. If we nevertheless take g to be constant,
we obtain by integration
d
z = v0 + gt,
dt
1
z = z0 + v0 t + gt2 ,
2
from which we obtain the observation equations of the measurement
Ó »îá .
222 Gravimetric measurement devices
process
⎡ ⎤
z0
1 2
[ ]
zi = 1 ti 2 ti
· v0 ⎦ + n i .
⎣ (10.6)
g
Here, the unknowns10 are zˆ0 , vˆ0 and ĝ. The quantities z i are the in-
terferometrically measured vertical locations of the falling prism, and
n i are the residuals of the measurements. Determining precisely the
corresponding measurement time or epoch t i is of course essential. The
volume of measurements obtained from each drop is large.
We write the observation equations in matric form:
ℓ = A x + n,
where
t21
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
z1 n1 1 t1
⎢
⎢ z2 ⎥
⎥
⎢
⎢ n2 ⎥
⎥
⎢
⎢ 1 t2 t22 ⎥
⎥ ⎡ ⎤
⎢ .. ⎥ ⎢ .. ⎥ ⎢ .. .. .. ⎥ z0
⎢ . ⎥ ⎢ . ⎥ ⎢ . . . ⎥
ℓ=⎢ ⎥, n=⎢ ⎥, A=⎢ ⎥ and x = v0 ⎦ .
⎣
⎢ zi ⎥ ⎢ ni ⎥ ⎢ 1 ti t2i ⎥
⎢
⎢ ..
⎥
⎥
⎢
⎢ ..
⎥
⎥
⎢
⎢ .. .. ..
⎥
⎥ g
⎣ . ⎦ ⎣ . ⎦ ⎣ . . . ⎦
zn nn 1 t n t2n
From this, the solution follows according to the method of least squares
adjustment, from the normal equations
A T A x̂ = A T ℓ
10
It would be easy (exercise!) to add an unknown representing the vertical
gradient of gravity to this.
Ó »îá .
The absolute or ballistic gravimeter 223
Splitting Readout
beam beam
Mirroring
beam
Without gravity
In gravity field
t
11
This is a quantum theoretically erroneous statement. The matter wave of each
individual atom splits into two!
12
In fact, the spinning of the atom’s phase angle acts like a clock, and the speed
at which time elapses depends on the local geopotential (Vermeer, 1983).
Ó »îá .
224 Gravimetric measurement devices
Ó »îá .
The superconducting gravimeter 225
Middle capac-
itor plate
Levitation coil
Feedback
coil
Levitation coil
13
http://en.wikipedia.org/wiki/Mu-metal.
14
Virtanen (2006) reports how the instrument at Metsähovi detected the change
in gravity as workmen cleared snow from its laboratory roof, including a tea
break! “Weighing” visitors to the lab by their gravitational attraction is also
standard fare.
15
Their periods range from under an hour to over twenty hours, and they are of
Ó »îá .
226 Gravimetric measurement devices
g ≈ 9.8 m /s²,
Ó »îá .
Airborne gravimetry and GNSS 227
Ó »îá .
228 Gravimetric measurement devices
g i = Γ i − 〈a i · n i 〉 .
H ∼ ∆ x = v∆ t,
Ó »îá .
Measuring the gravity gradient 229
where v is the flight speed. The separation between adjacent flight lines
is chosen similarly.
The first major airborne gravimetry project was probably the Greenland
Aerogeophysics Project (Brozena, 1992). In this ambitious American-
Danish project in the summers of 1991 and 1992, over 200,000 km was
flown, all the time measuring gravity and the magnetic field, and the
height of the ice surface using a radar altimeter (Ekholm et al., 1995).
After that, also other large uninhabited areas in the Arctic and Antarctic
regions have been mapped, see Brozena et al. (1996), Brozena and Peters
(1994). Already in subsection 8.6.2 on page 176 we made mention of other
large surveys. Activity continues, see Coakley et al. (2013), Kenyon et al.
(2012). The method is well suitable for large, uninhabited areas, but also,
e.g., for sea areas close to the coast where ship gravimeters would have
difficulty navigating long straight tracks. In 1999 an airborne gravimetry
campaign was undertaken over the Baltic, including the Gulf of Finland
(Jussi Kääriäinen, personal comm.).
In addition to the economic viewpoint, an important advantage of air-
borne gravimetry is, that a homogeneous coverage by gravimetric data is
obtained from a large area. The homogeneity of surface gravimetric data
collected over many decades is difficult to guarantee in the same way.
Also the effect of the very local terrain, which for surface measurements
is a hard to remove systematic error source especially in mountainous
terrain, (see section 5.3 on page 99), does not come into play for airborne
gravimetry.
The operating principle of satellite gravimetry, e.g., GOCE (Geopotential
and Steady-state Ocean Circulation Explorer) is similar. An essential
difference is however, that the instrumentation on the satellite is in a
state of weightlessness: Γ = 0 (in a high orbit, or when using an air
ilmanvastus drag compensation mechanism), or Γ is small and is measured using a
sensitive accelerometer (in a low orbit, where air drag is notable).
The greatest challenge in planning a satellite mission is choosing the
flight height. The lowest possible height is some 150 km. At that height,
ajoaine already a tankload of propellant is needed, or the flight will not last long.
erotuskyky However, the resolution of the measurements on the Earth’s surface is
limited, e.g., the smallest details in the Earth’s gravity field “seen” by the
GOCE satellite are 50 − 100 km in diameter.
Ó »îá .
230 Gravimetric measurement devices
We know that gravity increases going down, at least in free air. Going up,
gravity diminishes, about 0.3 mGal for every metre of height.
In topocentric co-ordinates ( x, y, z), where z points to the zenith, this
matrix is approximately
⎡ ⎤
−0.15 0 0
M≈ ⎣ 0 −0.15 0 ⎦ mGal /m,
0 0 0.3
2
where ∂∂ zW2 = ∂∂z g z ≈ 0.3 mGal /m is the standard value for the free-air verti-
cal gravity gradient: Newton’s law gives for a spherical Earth (the minus
sign is because g points downward while the z co-ordinate increases going
up):
GM
gz = − .
(R + z)2
Derivation gives
∂ GM ∂ (R + z) 2gz
gz = = 2 =− ≈ 3 · 10−6 m /s2 /m = 0.3 mGal /m.
∂z (R + z) 3 ∂z (R + z)
2 2
The quantities ∂∂ xW2 and ∂∂ yW2 again describe the curvatures of the equipo-
tential or level surfaces in the x and y directions, equations 3.3, 3.4:
∂2 W g ∂2 W g
=− and =− ,
∂ x2 ρ1 ∂y 2 ρ2
∂2 W ∂2 W
= ≈ −1.5 · 10−6 m /s2 /m = −0.15 mGal /m.
∂ x2 ∂ y2
Ó »îá .
Self-test questions 231
Note that
∂2 W ∂2 W ∂2 W
+ + ≈ 0,
∂ x2 ∂ y2 ∂ z2
the familiar Laplace differential equation. However, the equation is not
exact here: in a co-ordinate system co-rotating with the Earth, the term
for the centrifugal force, 2ω2⊕ , must be added, equation 3.1.
The gravity-gradient field of Sun and Moon is known on the Earth’s
surface as the tidal field, see section 13.1.
D Self-test questions
Ó »îá .
232 Gravimetric measurement devices
When we use a spring gravimeter in the field, we always take care that
The reason for this is that the functioning of a spring gravimeter depends
on the properties of the spring material, which may change as a result of
careless handling or temperature variations.
Furthermore a gravimeter always has a drift, i.e., the relationship be-
tween measured value and true value changes slowly over time. In a non-
factory fresh gravimeter this drift is however very regular and almost
linear.
As a result of the drift, a spring gravimeter cannot be used for absolute
gravity measurement and is therefore called a relative gravimeter.
Question: how is the relative nature of a spring gravimeter and its drift
taken into account
1. How much does a low-pressure zone of 100 hPa (i.e., the pressure is
100 hPa lower than average air pressure 1015 hPa) affect gravity
measured on the Earth surface? (You may assume the low-pressure
zone to be very extended in area.)
2. How much does sea water rise due to the “upside-down barometer
effect” under a low-pressure zone?
3. How much does the effect from point 2 amount to in local gravity
measured on a ship? On the open sea, with a free-air gravity
Ó »îá .
Exercise 10 – 3: Air pressure and gravity 233
gradient of −0.3 mGal /m, density of sea water 1000 kg /m3 . Analyze the
situation carefully18 .
18
And I mean really carefully.
Ó »îá .
D The geoid, mean sea level,
11 sea-surface topography
– 235 –
236 The geoid, mean sea level, sea-surface topography
“the level surface of the Earth’s gravity field that agrees most
closely with mean sea level.”
The practical problem with this definition is, that determining the correct
level of the geoid requires knowledge of mean sea level everywhere on the
world ocean. This is why many “geoid” models in practice don’t coincide
with global mean sea level, but with some locally defined mean sea level
— and often only approximately.
Mean sea level in its turn is also a problematic concept. It is sea level
from which has been computationally removed all periodic effects — but
who can know if a so-called secular effect in reality is perhaps long
period? A sensible compromise is the average sea level over 18 years —
an important periodicity, saros1 , in the orbital motion of the Moon.
The sea-surface topography again is defined as that part of the difference meritopografia
between mean sea level and the geoid, which is permanent. Also here,
the measure of permanency is the time series that are available, as tide
gauges have been widely operating already for about a century, when mareografi
again many satellite time series — TOPEX/Poseidon and its successors —
are just about a quarter of a century long. See figure 12.1.
1
https://en.wikipedia.org/wiki/Saros_(astronomy).
Ó »îá .
Geoids and national height datums 237
Ó »îá .
238 The geoid, mean sea level, sea-surface topography
Absolute land uplift is the motion of the Earth’s crust relative to the
centre of mass of the Earth. This land uplift is measured when
using satellites the orbits of which are determined in a co-ordinate
reference system tied to the Earth’s centre of mass. E.g., GNSS
positioning of tide gauges.
Relative land uplift is the motion of the Earth’s crust relative to mean
sea level. This motion is measured by tide gauges, also called
mareographs.
Geoid rise: as the post-glacial land uplift is the shifting of masses in-
ternal to the Earth from one place to another, it is clear that also
the geoid must change. The geoid rise is however small compared
to the land uplift, only a few percent of it.
d
Equation (the point above a quantity denotes the time derivative dt ):
in which
Ó »îá .
Methods for determining the sea-surface topography 239
The rise in the geoid as a result of land uplift can be simply calculated
with the Stokes equation:
¨ ( )
dN R d
∆ g d σ.
( )
= S ψ
dt 4πγ σ dt
d
Here, dt ∆ g is the change of gravity anomalies over time due to land
uplift. Unfortunately we do not precisely know the mechanism by which
mass flows into the land uplift area in the Earth’s mantle. We may write
d dh
∆g = c ,
dt dt
in which the constant c may range from −0.16 to −0.31 mGal /m.
Up until fairly recently, the most likely value was about −0.2 mGal /m, with
substantial uncertainty. The latest results (Mäkinen et al., 2010) give
−0.16 ± 0.02 mGal /m (one standard deviation), which would seem to settle
the issue. It looks like the Bouguer hypothesis is closer to physical reality.
The flow of mass happens probably within the asthenosphere.
This problem has been studied much in the Nordic countries. The method
has been gravimetric measurement along the 63◦ N parallel (“Blue Road
Geotraverse” project). The measurement stations extend from the Norwe-
gian coast to the Russian border, and have been chosen so, that gravity
along them varies within a narrow range. In this way, the effect of the
scale error of the gravimeters is avoided. Clearly, absolute gravity is of
no interest here, only the change in gravity differences over time between
the stations.
These measurements have been made over many years using high pre-
cision spring or relative gravimeters. In recent years, there has been a
shift to using absolute gravimeters, obviating the need for measurement
lines.
Ó »îá .
240 The geoid, mean sea level, sea-surface topography
Earth’s crust
Asthenosphere
(a)
Bouguer hypothesis. . .
Earth’s crust
Upper mantle
(b)
. . . and free-air hypothesis.
65 5 35
65
10 30
15 20 25
Föllinge
Meldal Stugun
Kopperå Vaasa Joensuu
Vågstranda Kramfors
Äänekoski
60 60
5 35
10 30
15 20 25
Ó »îá .
Global sea-surface topography and heat transport 241
salinity measurements along vertical profiles are used on the open ocean,
and geostrophic levelling if ocean current measurements are used to
determine the Coriolis effect, generally close to the coast.
All methods should give the same results. The Baltic Sea is a textbook
example, where all three methods have been used. A result was that
the whole Baltic Sea surface is tilted: relative to a level surface, the sea
surface goes up from the Danish straits to the bottoms of the Gulf of
Finland and the Bothnian Bay by some 25 − 30 cm.
Oceanographic model calculations show, that this tilt is mainly due to a
salinity gradient: in the Atlantic Ocean, salinity is 30 − 35 o /oo, when in
the Baltic it drops to 5 − 10 o /oo due to the massive production of sweet
water (Ekman, 1992). Of course on top of this come temporal variations,
like oscillations like in a bathtub, the amplitude of which can be over a
metre.
In Ekman (1992) more is said about the sea-surface topography of the
Baltic and its determination.
a = 2 v×−
→ ,
⟨ ⟩
ω ⊕ (11.1)
If a fluid flows on the Earth’s surface, then, in the above equation 11.1,
only the part of − → normal to the surface will have an effect: this part
ω ⊕
has a length of ω→
⟨− ⟩
⊕ · n = ω⊕ sin ϕ, and the vector equation 11.1 may be
replaced by a simpler scalar equation:
a = 2vω⊕ sin ϕ,
def
where a = ∥a − 〈a · n〉∥, i.e., the length of the projection of a onto the
tangent plane to the Earth, and v = ∥v∥ , ω⊕ = −
def def
→
ω⊕ etc. in the familiar
way. The direction of the Coriolis acceleration is always perpendicular to
the flow velocity: when watching along the flow direction, to the right on
the Northern hemisphere, to the left on the Southern hemisphere.
Ó »îá .
242 The geoid, mean sea level, sea-surface topography
As a result of the Coriolis force, the sea surface in the area of an ocean
current is tilted sideways with respect to the current, at an angle
a 2 v ω⊕
= sin ϕ.
γ γ
∂H ω⊕ ∂H ω⊕
= −2v y sin ϕ, = +2v x sin ϕ, (11.2)
∂x γ ∂y γ
2
A popular, though unofficial, unit for ocean current is the sverdrup (http:
//en.wikipedia.org/wiki/Sverdrup), a million cubic metres per second. All the
rivers of the world together make about one sverdrup, while the Gulf Stream is
30 − 150 Sv.
Ó »îá .
Global sea-surface topography and heat transport 243
Ó »îá .
244 The geoid, mean sea level, sea-surface topography
3
6000 years “before present”, 6 ka BP. BP conventionally means: before 1950.
Nowadays is also used b2k, before the year 2000.
Ó »îá .
The sea-level equation 245
Sea level
drops Sea level
Sea level
rises
drops
Greenland
Antarctica
Figure 11.5. The sea-level equation. Sea level reacts in a complicated way when
D continental ice sheets melt.
where
( )
◦ S = S (ω, t) = S ϕ, λ, t describes the variations of sea level as a
( )
function of place ω = ϕ, λ and time t,
◦ I = I (ω, t) is similarly a function of place and time describing the
geometry of ice sheets and glaciers,
◦ S E is the eustatic term, i.e., the variation in ice volume converted to
“equivalent global sea-level variation”, in an equation
m i ( t)
S E ( t) = ,
ρo Ao
4
http://www.atmosp.physics.utoronto.ca/∼peltier/data.php.
5
https://geodynamics.org/cig/software/selen/selen-manual.pdf.
Ó »îá .
246 The geoid, mean sea level, sea-surface topography
function is multiplied with the ice and sea functions and integrated
over the domain in question. These integrals are by the way very
similar to the ones discussed in section 7.1, e.g.:
ˆ t ¨
G s ψ ω, ω′ , t − t ′ S ω′ , t′ d ω′ dt′ ,
{ ( ) ( )} ( )
{G s ⊛o S } (ω, t) =
−∞ ocean
( )
where ψ ω, ω′ is the geocentric angular distance between evalua-
( ) ( )
tion point ω = ϕ, λ and data point ω′ = ϕ′ , λ′ . The measure of the
( ) ( )
surface integral is d ω = N ϕ M ϕ cos ϕ d ϕ d λ, where N, M ≈ R
are the principal radii of curvature of the Earth ellipsoid. As can
be seen, we have here a convolution applied both over the Earth’s
surface ω and over the time axis t.
◦ The overbar designates averaging over the whole ocean surface,
◦ γ0 is an average acceleration of gravity,
◦ G s is the so-called Green’s function of sea level:
G s = G V − γ0 G u ,
G V ψ, t = G rV ψ, t + G eV ψ, t + G vV ψ, t
( ) ( ) ( ) ( )
G u ψ, t = G eu ψ, t + G vu ψ, t
( ) ( ) ( )
The behaviour of sea level can now be computed in this way, that one
first tries to construct an “ice-load history”, i.e., I (ω, t); Then, from this
one tries to calulate iteratively, using the sea-level equation 11.3, S (ω, t).
Note that S describes relative sea-level variation, i.e., changes in the
relative positions of sea level and the Earth’s solid body or Earth’s crust.
It is a function of place: one may not assume that it would be the same
everywhere. In the article Mitrovica et al. (2001) it is shown how, e.g., the
melting water from Greenland flees to the Southern hemisphere, when
the melting water from Antarctica again comes similarly to the North.
This is a consequence from the change in the Earth’s gravity field and
the geoid, when large volumes of ice melt. Another factor is that also the
physical shape of the Earth changes, when the ice load changes: so-called
Glacial Isostatic Adjustment, GIA.
Ó »îá .
The sea-level equation 247
Figure 11.6. Sea-level rise after the last ice age (Wikimedia Commons, © Robert
D A. Rohde, GNU Free Documentation License).
This also complicates the monitoring of global mean sea level from local
measurements: the problem is familiar in Fennoscandia, as the Earth’s
crust, for now, moves up faster than global sea-level rise. . .
Green’s functions in the sea-level equation are functions of both ψ and
time difference ∆ t; this tells us that GIA is a function of both place and
time. On a spherically symmetric Earth, the functions may be written as
expansions, e.g.6
I
∞
( )
R γ0 ∑ ∑
G vV ψ, ∆ t = H (∆ t) k ℓ i e−sℓ i ∆ t
( ) ( )
Pℓ cos ψ ,
M
ℓ=1 i =1
6
We consider here only the plastic or viscous deformation.
Ó »îá .
248 The geoid, mean sea level, sea-surface topography
D Self-test questions
If the velocity of flow of an ocean current is 0.1 m /s and its width 100 km,
compute:
1. How much at latitude 45◦ N is the height difference between its left
and right edges?
2. If the same current was 200 km broad and the velocity of flow
0.05 m /s (i.e., assuming the same depth, also the transport of water
would be the same), compute for that case the height difference
between the left and the right edges.
3. [For fun] if the depth of the current is 1 km, what is the water
transport in sverdrup?
How does the post-glacial land subsidence observed in the United States
and central Europe support a Bouguer type of land-uplift mechanism
Ó »îá .
Exercise 11 – 2: Land subsidence and the mechanism of land uplift 249
Ó »îá .
D Satellite altimetry and satellite gravity
12 missions
– 251 –
252 Satellite altimetry and satellite gravity missions
1
But read this.
Ó »îá .
Satellite altimetry 253
Figure 12.1. Results from the TOPEX/Poseidon and Jason satellites. Left,
global mean sea-level rise; right, correlation between global
mean sea level and ENSO (“El Niño”). © Colorado University at
Boulder’s Sea Level Research Group http://sealevel.colorado.edu/,
D Nerem et al. (2010).
h sat = h 0,sat + ∆ h,
Ó »îá .
254 Satellite altimetry and satellite gravity missions
True orbit
Calculated orbit
h sat ℓ
footprint
Sea-surface topography H
Geoid
Reference ellipsoid
returns to the receiver) will be larger, and the distance travelled by the
radio waves will on average be longer.
Newer satellites use an interferometric technique that differs somewhat
from the description above.
Of all the corrections related to instrumentation, atmosphere, ocean, and
solid Earth, we mention
Ó »îá .
Crossover adjustment 255
ha = N + H + ∆ h + ϵ + n,
In the figure we have three tracks and two crossing points. The observa-
tion equations, which describe the discrepancies in the known crossover
points as functions of the orbit errors, are
ℓ1 = ∆ h 2 − ∆ h 3 + n 1
ℓ2 = ∆ h 1 − ∆ h 3 + n 2
Ó »îá .
256 Satellite altimetry and satellite gravity missions
∆h2 ∆h3
∆h2 − ∆h3
∆h1
∆h1 − ∆h3
Crossover 1
Crossover 2
2
1
or in matric form2
x
ℓ A⎡ ⎤ n
[ ] [ ] ∆h1 [ ]
ℓ1 0 1 −1 ⎣ n1
= ∆h2 ⎦ + , (12.1)
ℓ2 1 0 −1 n2
∆h3
symbolically
ℓ = A x + n.
If one now tries to calculate the solution with ordinary least squares,
)−1 T
x̂ = A T A
(
A ℓ,
one will notice that this doesn’t work. The normal matrix A T A is singular
(check!). This makes sense, as one can move the whole track network up
or down without the observations ℓk changing. No unique solution can
be found for such a system.
Finding a solution requires that something must be fixed. E.g., one
track, or, more democratically, the mean level of all tracks. This fixing is
2
Note the similarity with the observation equations for levelling! Instead of
benchmarks, we have tracks, instead of levelling lines, crossover points.
Ó »îá .
Crossover adjustment 257
T
)−1 T −1 T −1 T −1
[(
T −1 T
]
ℓ = A −1 ℓ.
( ( ) )
A A A ℓ= A A A ℓ= A A A
3
http://maxima.sourceforge.net/.
Ó »îá .
258 Satellite altimetry and satellite gravity missions
and then multiply both sides, and both terms on the right, with the
diagonal matrix
⎡ ⎤
1 0 0
def
D=⎣ 0 1 0 ⎦.
0 0 1c
The result is
Dℓ DA Dn
⎡ ⎤ ⎤ ⎡ ⎤ ⎡ ⎤
ℓ1 ∆h1
⎡
0 1 −1 n1
⎣ ℓ2 ⎦ = ⎣ 1 0 −1 ⎦ ⎣ ∆h2 ⎦ + ⎣ n2 ⎦,
0 1 1 1 ∆h3 0
∆ h = a + b τ,
where the parameter τ is the location along the track reckoned from its
starting point. The dimension of this location can be time (seconds) or
angular distance (degrees). Now the set of observation equations for the
situation described above is
x
⎡ ⎤
a1
ℓ A n
[ ] [
⎢
⎢ b1 ⎥
⎥ [ ]
ℓ1 1 τ21 −1 −τ31 ⎢
]
0 0 a2 n1
⎢ ⎥
= ⎥+ .
⎥
ℓ2 1 τ12 0 0 −1 −τ32 ⎢ b2 n2
⎢
⎥
⎢ ⎥
⎣ a3 ⎦
b3
Ó »îá .
Crossover adjustment 259
Here the design matrix contains, besides the values 1 and −1, also the
expressions τki , in which i is the number of the track, k that of the
crossover point. These are computable when the geometry of the tracks
is known.
Now there are two unknowns for every track, a and b, a constant and
a trend. Of course also this system will prove to be singular. Removing
the singularity can be done by fixing all three parameters b and one
parameter a4 .
The phenomenon that no solution can be found unless something is fixed,
is called a datum defect. Fixing something suitable will define a certain
datum. Between two different datums exists a transformation formula:
in the case of one orbit error parameter per track, this transformation is
a simple parallel shift or translation of all tracks up or down.
The situation is somewhat similar as when defining a height or vertical
reference system, for a country. One needs to fix one point, e.g., Helsinki
harbour. If alternatively one fixes another point, e.g., Turku harbour,
the result is another datum, in which all height values differ from the
corresponding ones in the first datum by a certain fixed amount.
The argument continues to hold if there is a large number of tracks: say,
ten Northgoing and ten Southgoing tracks, crossing in 10 × 10 crossover
points. Here, for two parameters per track, we would have 40 unknowns
and no less than 100 observations. Still, we must constrain the absolute
level and the various trends and possible other deformations of the whole
network of tracks. It gets complicated, but a simple approach is to attach a
priori uncertainties to the unknowns a i , b i to be estimated, e.g., from the
known uncertainties of the orbit prediction available. The least-squares
adjustment equation then becomes
)−1 T
x̂ = A T A + Σ−1
(
A ℓ,
where Σ is the diagonal matrix containing the a priori variances5 σ2a,i ,
σ2b,i of the parameters of each track i . This is referred to as Tikhonov6
regularization.
4
In order to understand this, build, e.g., a three track “wire-frame model” from
pieces of iron wire, tied together by pieces of string at the crossover points.
Crossover conditions don’t in any way fix the values of the trends b, and the
whole absolute level of the frame continues to be unconstrained.
5
We assume the mean error of unit weight to be 1.
6
Andrey Nikolayevich Tikhonov (1906 – 1993) was a Russian Soviet mathemati-
cian and geophysicist.
Ó »îá .
260 Satellite altimetry and satellite gravity missions
Questions:
Answers:
Ó »îá .
Choice of satellite orbit 261
∆ h = a 00 + a 10 x + a 01 y + a 11 x y
with four degrees of freedom. So, fix one bias and three trends, not
all North- or all Southgoing.
6. Any such choice would be arbitrary. Rather use the method de-
scribed above instead, Tikhonov regularization.
where now τ is an angular measure, e.g., the place along the track mea-
sured from the last South-North equator crossing. See Schrama (1989),
where this problem is treated more extensively. In this model, a repre-
sents the size of the orbit, while ( b, c) describe the offset of the centre
of the orbit from the geocentre. This model is three-dimensional: the
orbital arcs with their crossovers form a spherical network surrounding
the Earth. The degrees of freedom left by the crossover conditions are
now the size of this sphere and the offset of its centre from the geocentre:
with ( X , Y , Z ) geocentric co-ordinates, we have
∆h = a0 + a1 X + a2 Y + a3 Z (12.4)
7
One could argue that, in eq. 12.3, the parameter a should be zero, as Kepler’s
third law allows a very precise determination of the orbital size, see section 12.3.
Then, also a 0 = 0 in equation 12.4.
Ó »îá .
262 Satellite altimetry and satellite gravity missions
z
Earth’s
rotation E
Perigee Nodal line
axis
y
Nodal line
ν ω
ν
Satellite i Apogee Perigee
a
Ω b
ω X
Ω̇
θ
Apsidal
Earth
line ea
x rotation
′
x (Greenwich)
(Vernal equinox)
Ascending node
Apogee
◦ If one wishes to study the precise shape of the mean sea surface,
one chooses a long repeat period, in order to get the tracks as close
together as possible on the Earth’s surface.
◦ If one wishes to study the variability of the sea surface, one chooses
an orbit that returns to the same location after a short time interval.
Then, the grid of tracks on the Earth’s surface will be sparser.
8
This is why it is said that Henry Cavendish was the first to “weigh the Earth”. . .
determining GM was already straightforward back then using the orbital motion
of the Moon. The challenge was separating G and the mass of the Earth M,
obtaining the latter in ordinary units of mass.
Ó .
»îá
Choice of satellite orbit 263
Attraction
Solar
caused by Satellite
(apparent)
Earth’s orbital
daily
flattening motion
motion
X
Also parameters describing the figure of the Earth affect satellite mo-
tion, e.g., the quantity J2 , the dynamic flattening, having a value of
J2 = 1082.6267 · 10−6 . It is just one of many so-called spherical-harmonic
coefficients that describe the figure of the Earth and affect satellite orbits.
In the case of J2 the effect is, that the plane of the satellite orbit rotates
at a certain rate (orbital precession), which make the satellite, if she flies
over the same location the next day, do so several minutes earlier. The
equation is, for a circular orbit of radius a:
√
dΩ 3 GM a2E
=− J2 cos i,
dt 2 a3 a2
in which again a E is the equatorial radius of the Earth reference ellipsoid,
and i the inclination of the orbital plane relative to the equator. If we
substitute numerical values into this, we obtain
dΩ cos i
= −1.318,95 · 1018
[ 3.5 −1 ]
m s ,
dt (a E + h)3.5
where h is the mean height of the satellite orbit, conventionally above a
sphere of size equatorial radius a E . If we substitute into this, e.g., the
satellite height h = 800 km (and use a E = 6,378,137 m) we obtain
dΩ ◦
= −1.331,02 · 10−6 cos i rad s−1 = −6day
.589
[ ]
· cos i. (12.6)
dt
For practical reasons — solar panels! — we often choose the satellite
orbit such, that the orbital plane turns along with the annual apparent
◦
0◦. 9856
motion of the Sun, 365.360
25 days = day . See figure 12.6.
Ó »îá .
264 Satellite altimetry and satellite gravity missions
Spring
Summer Winter
Autumn
D 12.3.1 Example
Questions:
Ó »îá .
Choice of satellite orbit 265
i
3 2 1
Answers:
1. The satellite completes 14 orbits per day, i.e., per 1440 minutes:
P = 1440
14 = 102.857 min.
2. The satellite completes 43 orbits in three days, i.e., per 3 × 1440
minutes: P = 3×43
1440
min = 100.465 min.
3. The satellite completes 502 orbits in 35 days, i.e., per 35 × 1440
minutes: P = 35×502
1440
min = 100.398 min.
4. Execute the octave code in table 12.2. The result is 780.604 km.
Ó »îá .
266 Satellite altimetry and satellite gravity missions
%format long
GM=3986005e8;
ae=6378137;
P=100.465*60; % seconds
fac=4*pi*pi; % four pi square
a=(GM*P*P/fac)^0.33333333;
h = a - ae;
printf(’\n\nOrbital height: %8.3f km.\n’, h/1000);
5. The same code, with P=100.398*60, yields 777.421 km. The differ-
ence with the previous is 3.183 km.
6. There are 43 orbits with different ground tracks. That means a
40,000
separation of 360
43 = 8.372 degrees, or at the equator, 43 = 930 km;
less at higher latitudes.
360 40,000
7. 502 = 0.717 degrees, or 502 = 80 km.
8. (a) The 35-day orbit would be excellent for detailed mapping. The
three-day orbit would be able to see, e.g., tides or weather-
related phenomena, but at poor resolution.
(b) The difference in height being only 3 km, and in period, 4 s,
the change in orbit between the two repeat periods should be
easily within reach of even small on-board thrusters. So, yes.
Ó »îá .
Retracking 267
Travel time
Half-
height
rule
Figure 12.9. Analyzing the altimeter return pulse. The classical return pulse
D time measurement uses the “half-height point”.
D 12.5 Retracking
The results of a satellite altimetry mission are published already during
flight in the form of a so-called geophysical data record (GDR) file, con-
taining everything related to the measurement and, e.g., atmospheric
correction terms, tidal corrections, wave parameters, etc.
It is common practice today to process again altimetry measurements al-
ready collected earlier, in order to extract further useful information. The
complete return pulse is analyzed again in a method called retracking9 .
The standard method of analysis is based on that point on leading edge
of the return pulse, which is at half height from the maximum value of
the pulse. This is according to experience a good way to get the travel
time associated with the point in the centre of the footprint, directly
underneath the satellite. In the back part of the pulse are reflections
from the further-away peripheral areas of the footprint.
There are however two situations where this method doesn’t work prop-
erly during flight, and a more careful a posteriori analysis of the pulse is
worthwhile:
9
http://www.altimetry.info/radar-altimetry-tutorial/data-flow/data-processing/
retracking/
Ó »îá .
268 Satellite altimetry and satellite gravity missions
a
http://psc.apl.washington.edu/research/projects/arctic-sea-ice-volume-
D anomaly/.
1. is constant
2. coincides with a level surface, i.e., is the same as the geoid.
In practice however the ocean surface is variable in time and is also not a
level surface. Therefore, other approaches have appeared.
Ó »îá .
Satellite gravity missions 269
Ó »îá .
270 Satellite altimetry and satellite gravity missions
GPS-2 GPS-3
Magnetometer GPS-4
GPS-1
beam
GPS antenna
Acceleration vector
Nom
inal
orbi
t
”CHAMP” Tru
e
orb
Solar it
cells
Earth internal
mass density variations (example)
Figure 12.11. Determining the Earth’s gravity field from GPS orbital tracking
D of a low flying satellite.
Ó »îá .
Satellite gravity missions 271
Sate
te 1 llite
li 2
Satel Distance 220 km
Precise ranging, wavelength 1.5 cm
km
450
Figure 12.12. The principle of the GRACE satellites: measuring the minute
variations in time of the gravity field using SST (satellite-to-
satellite tracking). The changes are due to mass shifts in the
“blue film” – the atmosphere and hydrosphere – and expressed as
D variations in “total sea-floor pressure” (↓).
10
http://commons.wikimedia.org/wiki/File:Global_Gravity_Anomaly_
Animation_over_LAND.gif
Ó »îá .
272 Satellite altimetry and satellite gravity missions
- Measurement of
GOCE satellite acceleration
1 differences
5 -
4 X
3
6
Gradiometer
Accelerometer
Unknown
density variations
Figure 12.13. Determining the Earth’s gravity field with the gravitational gra-
D diometer on the GOCE satellite.
11
Because of this inclination angle, there was a cap of radius 3.3◦ at each pole
within which no measurements were obtained.
12
http://blogs.esa.int/rocketscience/2013/11/11/goce-burning-last-orbital-view/
Ó »îá .
Self-test questions 273
D Self-test questions
Given are two Northgoing satellite tracks and three Southgoing ones.
There are six crossover sites, see figure 12.14.
Ó »îá .
274 Satellite altimetry and satellite gravity missions
Ó »îá .
Exercise 12 – 3: Kepler’s third law 275
Ó »îá .
D Tides, atmosphere and
13 Earth crustal movements
GMR 2 GMR 2 ( 2
)
W= P 2 (cos ζ ) + . . . = 3 cos ζ − 1 +...,
d3 2d 3
where d is the distance to either the Moon or the Sun, R the radius of
the Earth, and ζ the local zenith angle of the Sun or Moon. P2 (cos ζ) is
the Legendre polynomial of degree two. GM is the mass of the Sun or
Moon multiplied by Newton’s gravitational constant. In the case of Sun
and Moon the extra terms (. . .) can be neglected, because these are such
remote bodies: d ≫ R .
According to spherical trigonometry (cosine rule on the sphere)
cos ζ = sin φ sin δ + cos φ cos δ cos h,
Low tide
High tide
High tide
Low tide
D Figure 13.1. Theoretical tide. ζ is the local zenith angle of the Moon (or Sun).
– 277 –
278 Tides, atmosphere and Earth crustal movements
or for n = 2,
( )
P2 (cos ζ) = P2 sin φ P2 (sin δ) +
1 ( )
+ P21 sin φ P21 (sin δ) cos h+
3
1 ( )
+ P22 sin φ P22 (sin δ) cos 2 h.
12
In this, according to table 2.3,
( )
P21 sin φ = 3 sin φ cos φ,
P22 sin φ = 3 cos2 φ,
( )
and we obtain
1(
3 cos2 ζ − 1 = P2 sin φ P2 (sin δ) +
) ( )
P2 (cos ζ) =
2
+ 3 sin φ cos φ sin δ cos δ cos h+
3
+ cos2 φ cos2 δ cos 2 h =
4
1( )1(
3 sin2 φ − 1 3 sin2 δ − 1 +
)
=
2 2
3
+ sin 2φ sin 2δ cos h+
4
3
+ cos2 φ cos2 δ cos 2 h.
4
From this
3 sin2 φ − 1 3 sin2 δ − 1 +
⎡ ( )( ) ⎤
2
GMR ⎣
W= + 3 sin 2φ sin 2δ cos h+ ⎦ .
4d 3
+ 3 cos2 φ cos2 δ cos 2 h
1
The declination is the geocentric latitude of the Moon.
2
The hour angle is the difference in longitude between the Moon and the local
meridian. It vanishes when the Moon is in upper culmination, i.e., due South
when seen from Northern non-tropical latitudes.
3
http://mathworld.wolfram.com/SphericalHarmonicAdditionTheorem.html.
Ó »îá .
Theoretical tide 279
GMR 2 [(
3 sin2 φ − 1 3 sin2 δ − 1 ,
)( )]
W1 = 3
4d
that still depends on δ and therefore is periodic with a 14-day
(half-month) period. Again by using spherical trigonometry:
1 ) (1
2 2 2 2
sin δ = sin ϵ sin ℓ =⇒ sin δ = sin ϵ sin ℓ = sin ϵ − cos 2ℓ ,
2 2
(13.1)
where ℓ is the longitude of the Moon in its orbit, reckoned from the
ascending node (equator crossing), and ϵ is the inclination of the
Moon’s orbit with respect to the equator, on average 23◦ but rather
variable, between 18◦ .3 and 28◦ .6 . Thus we obtain
GMR 2 [( 2
)( 2
(1 1 ) )]
W1 = 3 sin φ − 1 3 sin ϵ − cos 2 ℓ −1 ,
4d 3 2 2
where we have used result 13.1. We split W1 = W1a + W1b into two
parts, a constant4 and a periodic or semi-monthly (“fortnightly”)
part:
GMR 2 [( 2
)(3 2 )]
W1a = 3 sin φ − 1 sin ϵ − 1 , (13.2)
4d 3 2
GMR 2 [( 2
)(3 2 )]
W1b = − 3 sin φ − 1 sin ϵ cos 2 ℓ .
4d 3 2
GMR 2 [ ]
W2 = 3 sin 2 φ sin 2 δ cos h ,
4d 3
GMR 2 [ 2 2
]
W3 = 3 cos φ cos δ cos 2 h .
4d 3
In both, we have in addition to h, still δ as a “slow” variable. These
equations could be written out as sums of various functions of the
longitude of the Moon ℓ.
Use again basic trigonometry, equation 13.1:
1 ] [1
2 2 2 2 2
cos δ = 1 − sin δ = 1 − sin ϵ sin ℓ = 1 − sin ϵ − cos 2ℓ ,
2 2
1
cos 2ℓ cos 2 h = [cos(2ℓ + 2 h) + cos(2ℓ − 2 h)] ,
2
√
sin 2δ = 2 sin δ cos δ = 2 sin δ cos2 δ =
√ [1
2 1 ]
= 2 sin ϵ sin ℓ 1 − sin ϵ − cos 2ℓ ,
2 2
4
Not precisely, because ϵ is (slowly) time dependent.
Ó »îá .
280 Tides, atmosphere and Earth crustal movements
D Table 13.1. The various periods in the theoretical tide. The widely used symbols
were standardized by George Darwin.
a
Lunar fortnightly
b
Solar semi-annual
5
Paul Melchior (1925 – 2004) was an eminent Belgian geophysicist and Earth
tides researcher.
6
Arthur Thomas Doodson (1890 – 1968) was a British oceanographer, a pioneer
of tidal theory, also involved in designing machines for computing the tides. He
was stone deaf.
7
Sir George Howard Darwin (1845 – 1912) was an English astronomer and
mathematician, son of Charles Darwin of Origin of Species fame.
8
Augustus Edward Hough Love (1863 – 1940) was a British mathematician and
student of Earth elasticity.
Ó »îá .
Deformation caused by the tidal force 281
20 −0.1 20
0
0 −0.2 0 X
−20 −0.3 −20 −0.05
−40 −0.4 −40
−60 −0.5 −60 −0.1
−80 −0.6 −80 −0.15
Obliquity ϵ, 0◦ − 90◦ Ecliptic longitude, 0◦ − 360◦
0.2
20 20
0 0
0 0
−20 −0.2 −20 −0.2
−40 −0.4 −40 −0.4
−60 −60 −0.6
−0.6
−80 −80 −0.8
0 5 10 15 20 0 5 10 15 20
Hour angle Hour angle
Figure 13.2. The main components of the theoretical tide. These values must
D still be multiplied by Doodson’s constant D.
Here r is the distance from the geocentre. It is assumed here that the
Love numbers H n , L n depend only on r , i.e., the elastic properties of the
Ó »îá .
282 Tides, atmosphere and Earth crustal movements
Ó »îá .
The permanent part of the tide 283
GMR 2 [( 2
)(3 2 )]
Wperm = 3 sin φ − 1 sin ϵ − 1 ≈
4d 3 2
3GMR 2 ( 2 1)
≈ − 0.7615 sin φ − .
4d 3 3
Ó »îá .
284 Tides, atmosphere and Earth crustal movements
With the combined Doodson’s constant 13.3 for Sun and Moon equal to
3GMSun R 2 3GMMoon R 2
D= 3
+ 3
= (12.3 cm + 26.75 cm) × g = 39.05 cm × g
4 d Sun 4 d Moon
we obtain (1 )
2
Wperm = 29.74 cm × − sin φ × g.
3
We can express this, with the Bruns equation 4.2, into a permanent tidal
geoid effect:
(1 )
2
Nperm = 29.74 cm × − sin φ .
3
From this, Nperm (0◦ ) = 9.91 cm on the equator, and Nperm (90◦ ) =
−29.74 cm on the poles.
This, the geoid effect of the permanent part of the external potential of
Sun and Moon, is also equal to the difference between the mean geoid
and the zero geoid as defined above:
(1 )
def
∆mean
zero N = Nmean − Nzero = 29.74 cm × 2
− sin φ .
3
For heights H above sea level, with H = h − N , we have
(1 )
def
∆mean
zero H = Hmean − Hzero = −29.74 cm × 2
− sin φ ,
3
and for two different latitudes φ1 and φ2 we have for the effect on the
height difference
∆mean mean 2 2
zero H φ2 − ∆zero H φ1 = 29.74 cm × sin φ2 − sin φ1 .
( ) ( ) ( )
Ó »îá .
Loading of the Earth’s crust by sea and atmosphere 285
∆mean mean 2 2
tidefree H φ2 − ∆tidefree H φ1 = 38.66 cm × sin φ2 − sin φ1 .
( ) ( ) ( )
∆zero zero 2 2
tidefree H φ2 − ∆tidefree H φ1 = 8.92 cm × sin φ2 − sin φ1 .
( ) ( ) ( )
9
Hans-Georg Wenzel (1945 – 1999) was a German physical geodesist and geo-
physicist.
Ó »îá .
286 Tides, atmosphere and Earth crustal movements
very good way to study this phenomenon, because many more local, often
poorly known, factors affect local gravity. Measurement by GNSS is
promising but also challenging.
D Self-test questions
D Exercise 13 – 1: Tide
GMR 2 [( 2
)(3 2 )]
W1a = 3 sin ϕ − 1 sin ϵ − 1 ,
4d 3 2
where ϕ is latitude and ϵ is the obliquity of the Earth’s axis of rotation,
currently about 23◦ .
Ó »îá .
D Earth gravity field research
14
D 14.1 Internationally
In the framework of the IAG, the International Association of Geodesy,
research into the Earth’s gravity field is currently the responsibility of
the International Gravity Field Service. The IGFS was created in 2003
at the IUGG General Assembly in Sapporo, Japan, and it operates under
the IAG’s new Commission 2 “Gravity Field”. The United States National
Geospatial-Intelligence Agency (NGA) serves as its technical centre.
An important IAG service of great reputation is the International Gravity
Bureau, BGI, Bureau Gravimétrique International, located in Toulouse,
France (http://bgi.omp.obs-mip.fr/). The bureau works as an international
broker to which countries can submit their gravimetric materials. If some
researcher needs gravimetric material from another country, e.g., for
geoid computation, he can request it from the BGI, who will provide it
with the permission of the country of origin, provided the country of the
researcher has in its turn submitted its own gravimetric materials for
BGI use.
The French state has invested significant funds into this vital interna-
tional activity.
Another important IAG service in this field is the ISG, the International
Service for the Geoid. It has in fact already operated since 1992 under
the name International Geoid Service (IGeS), the executive arm of the
International Geoid Commission (IGeC). The ISG office is located in
Milano (http://www.isgeoid.polimi.it/) also with substantial support by the
Italian state. The task of this service is to support geoid determination
in different countries, for which purpose existing geoid solutions are
collected in a common data base, and international research schools
are organized to develop knowledgability and skills in the art of geoid
computation, especially in developing countries.
Both services, BGI and ISG, are under the auspices of the International
Gravity Field Service IGFS, as two of the many official services of the
IAG. Other IGFS services are the International Center for Earth Tides
– 287 –
288 Earth gravity field research
(ICET), the International Center for Global Earth Models (ICGEM), and
the International Digital Elevation Model Service (IDEMS).
D 14.2 Europe
In Europe operates the EGU, the European Geosciences Union, in the
frame of which much publication and meeting activities relating to gravity
field and geoid are being co-ordinated. The EGU organizes annually
symposia, where always also sessions are included on subjects related to
gravity field and geoid. Also American scientists participate. Conversely
the American Geophysical Union’s (AGU) fall and spring meetings1 are
also favoured by European researchers.
To be mentioned is the Geodetic Institute (“Institut für Erdmessung”) of
Leibniz University in Hannover, Germany, which since 1990 has acted as
the European computing centre of the International Geoid Commission
(IGeC, and produced high quality geoid models for use in Europe (Denker,
1998); (http://www.ife.uni-hannover.de/gravity.html?&L=1).
D 14.4 Finland
In Finland the study of the Earth’s gravity field has mainly been in the
hands of the Finnish Geodetic Institute, founded in 1918, one year af-
ter Finnish independence. The institute has been responsible for the
national fundamental levelling and gravimetric networks and their in-
ternational connections. In 2001 the Finnish Geodetic Institute’s gravity
1
Fall (autumn) meetings in San Francisco, spring meetings somewhere in the
world. The AGU, though American, is a very cosmopolitan player.
Ó »îá .
Textbooks 289
D 14.5 Textbooks
There are many good textbooks on the study of the Earth’s gravity field.
In addition to the already mentioned classic Heiskanen and Moritz (1967),
which is in large part obsolete, we may mention Wolfgang Torge’s book
Torge (1989). Difficult but good is Moritz (1980). Similarly difficult is
Molodensky et al. (1962). Worth reading also from the perspective of
physical geodesy is Vaníček and Krakiwsky (1986). A new book in this
field is Hofmann-Wellenhof and Moritz (2006).
Ó »îá .
D Field theory and vector calculus —
A core knowledge
∆E = 〈F · ∆r〉 ,
the scalar product of force F and path ∆r. Often, also in the sequel, we
leave the angle brackets 〈·〉 off.
Later we shall see that if the points 1 and 2, ∆r = r2 − r1 , are very close
to each other, we may write
dE = 〈F · d r〉 ,
– 291 –
292 Field theory and vector calculus — core knowledge
Let
def
s = 〈a · b〉
be the scalar product of the vectors a and b. Then
⟨ ⟩ ⟨ ⟩
µa · b = a · µb = µ 〈 a · b〉 ,
〈a · b〉 = 〈b · a〉 ,
and often we call
def
√
∥a∥ = 〈 a · a〉
the norm or length of vector a.
The following also applies:
〈a · b〉 = ∥a∥ ∥b∥ cos α,
where α is the angle between vectors a and b.
Let
def
c = 〈a × b〉
be the vectorial product of the two vectors a and b. Then (µ ∈ R):
⟨ ⟩ ⟨ ⟩
µa × b = a × µb = µ 〈 a × b〉 ,
〈a × b〉 = − 〈b × a〉 ,
and thus 〈a × a〉 = 0.
The resulting vector c is always orthogonal to the vectors a and b; the
length of vector c corresponds to the surface area of the parallellogram suunnikas
spanned by the vectors a and b. In a formula:
∥c∥ = ∥a∥ ∥b∥ sin α,
where again α is the angle between vectors a and b. If the angle is zero,
then also the vectorial product is zero (because then, a = µb for some
suitable value of µ).
Ó »îá .
Vector calculus 293
c = 〈 b × a〉
α
∥c∥
. . b
a
If r is the distance of the body (planet) from the centre of motion (the
dr
Sun), and is the distance travelled in a unit of time, then the product
dt
⟨ ⟩
dr
r× (A.1)
dt
is precisely twice the surface area of the triangle or “area” swept over.
Let us take the time derivative of this product, i.e., the expression A.1:
d2r
⟨ ⟩ ⟨ ⟩ ⟨ ⟩
d dr dr dr
r× = × + r× 2 .
dt dt dt dt dt
Ó »îá .
294 Field theory and vector calculus — core knowledge
Angular momentum
〈r × v〉
Planet
Sun
v
1
2 r×v
∥〈 〉∥
Velocity
vector
Radius vector
r
Figure A.2. Kepler’s second law. In the same amount of time, the radius vector
of a planet will “sweep over” a same-sized area — conservation of
D angular momentum.
Ó »îá .
Scalar and vector fields 295
c = 〈a × b〉 =
⏐ ⏐
⏐ i j k ⏐⏐
⏐
= ⏐⏐ a 1 a 2 a 3 ⏐⏐ =
⏐ b b b ⏐
1 2 3
= (a 2 b 3 − a 3 b 2 ) i + (a 3 b 1 − a 1 b 3 ) j + (a 1 b 2 − a 2 b 1 ) k.
I.e.,
c1 = a2 b3 − a3 b2 ,
c2 = a3 b1 − a1 b3 ,
c3 = a1 b2 − a2 b1 .
r = xi + yj + zk,
def ∂ ∂ ∂
∇=i +j +k .
∂x ∂y ∂z
∂F ∂F ∂F
g = grad F = ∇F = i +j +k .
∂x ∂y ∂z
So, the field g (r) = g ( x, y, z) is a vector field in the same space, the gradi-
ent field of F .
Ó »îá .
296 Field theory and vector calculus — core knowledge
D Figure A.3. The gradient. The level curves of the scalar field in blue.
Interpretation: the gradient describes the slope of the scalar field. The
direction of the vector is the direction in which the value of the
scalar field changes fastest, and its length describes the rate of
change with location. Imaging a hilly landscape: the height if the
ground above sea level is the scalar field, and its gradient is pointing
everywhere uphill, away from the valleys toward the hilltops. The g
arrows are the longer, the steeper is the slope of the ground surface.
The gradient operator (like also the divergence and the curl, see
later) is linear:
∂a 1 ∂a 2 ∂a 3
s = div a = 〈∇ · a〉 = + + .
∂x ∂y ∂z
Ó »îá .
Scalar and vector fields 297
Figure A.4. The divergence. Positive divergences (“sources”) and negative ones
D (“sinks”).
Given a vector field a ( x, y, z). We form the vectorial product c of this and
the nabla operator, producing again a vector field:
⏐ ⏐
⏐ i j k
⏐
∂ ∂ ∂
⏐ ⏐
⏐ ⏐
c = rot a = ∇ × a = ⏐⏐ ⏐=
⏐ ∂x ∂y ∂z
⏐
⏐
⏐ a1 a2 a3
⏐
⏐ ∂ ∂ ⏐ ⏐ ∂ ∂ ⏐ ⏐⏐ ∂ ∂
⏐ ⏐ ⏐ ⏐ ⏐ ⏐
⏐
⏐ ⏐ ⏐ ⏐
= ⏐⏐ ∂ x ∂ y ⏐⏐ i − ⏐ ∂x ∂z ⏐ j + ⏐ ∂ y ∂z ⏐k =
⏐ ⏐
⏐
⏐ a1 a2 ⏐ ⏐ a1 a3 ⏐ ⏐ a2 a3 ⏐
∂a 3 ∂a 2 ∂a 1 ∂a 3 ∂a 2 ∂a 1
( ) ( ) ( )
= − i+ − j+ − k,
∂y ∂z ∂z ∂x ∂x ∂y
Imagine a weather map, where low- and high-pressure zones are drawn.
Our vector field is the wind field. The wind circulates (on the Northern
hemisphere) clockwise around the high-pressure zones, and counterclock-
wise around the low-pressure zones. We may say that the curl of the wind
field is positive at the high pressures and negative at the low pressures.
(This is a poor metaphor as it is two-dimensional. In R2 , the curl is
a scalar, not a vector. Just like we need only one angle to describe a
rotational motion, when in R3 we need the three Euler angles.)
Ó »îá .
298 Field theory and vector calculus — core knowledge
Figure A.5. The curl. Positive (clockwise) and negative (counterclockwise) ed-
D dies.
Note that if
a ( x, y, z) = grad F ( x, y, z) ,
then also
a ( x, y, z) = grad (F ( x, y, z) + F0 ) ,
Ó »îá .
Integrals 299
∂ F0 ∂ F0 ∂ F0
grad F0 = i +j +k = 0.
∂x ∂y ∂z
a = grad F = ∇F,
∂ ∂ ∂ ∂ ∂ ∂
div a = ∇a = ∇∇F = F+ F+ F=
∂x ∂x ∂y ∂y ∂z ∂z
∂2 ∂2 ∂2
( )
def
= + + F = ∆F,
∂ x2 ∂ y2 ∂ z2
def ∂2 ∂2 ∂2
∆= + + .
∂x 2 ∂ y2 ∂ z2
For the potential of a “source free” field — e.g., for the gravitational
potential in vacuum, the electrostatic potential in an area of space free of
electric charges — this Delta, or Laplace, operator vanishes.
D A.3 Integrals
D A.3.1 The curve integral
We saw earlier, that work ∆E can be written as the scalar product of force
F and path ∆r:
∆E = 〈F · ∆r〉 .
The differential form of this is
dE = F · d r,
from which one obtains the integral form, the work integral
ˆ B
∆E AB = F · d r.
A
Here is computed the amount of work needed to move a body from point
A to point B by integrating F · d r along the path AB.
Ó »îá .
300 Field theory and vector calculus — core knowledge
Tangent vector t
rot a
Integral
¸
S 〈a · t〉 ds
Closed
path S
dx dy dz
t= i+ j+ k,
ds ds ds
we may also write
ˆ B
∆E AB = 〈F · t〉 ds,
A
the parametrized version of the integral.
Assume given again some vector field a and a surface in space S . Often,
one needs to integrate over the surface S the normal component of a
vector field, the projection of a onto the normal vector of the surface.
Let the normal vector of the surface be n. Then we must integrate
¨
〈a · n〉 dS,
S
symbolically written ¨
a · d S,
S
where the notation d S is called an oriented surface element. It is a vector
pointing in the same direction as the normal vector n.
Like a curve, also a surface can be parametrized. E.g., the Earth’s surface
(assumed a sphere) can be parametrized by latitude φ and longitude λ:
( )
r = r φ, λ . In this case we write as the surface element
dS = R 2 cos φ d φ d λ,
Ó »îá .
Integrals 301
Let S be a surface in space (not necessarily flat) and ∂S its edge curve.
Assume that the surface and its edge are well-behaved enough for all
necessary integations and differentiations to be possible. Then (Stokes):
¨ ˛
rot a · d S = a · d r.
S ∂S
In words: The surface integral of the curl of a vector field over a surface
is the same as the closed path integral of the field around the edge
of the surface.
Special case: If rot a = 0 everywhere (a conservative vector field) then
˛
a · d r = 0,
∂S
i.e., also
ˆ B ˆ B
a · dr = a · d r.
A, path 1 A, path 2
In other words, the work integral from point A to point B does not
depend on the path chosen. And the work done by a body transported
around a closed path is zero.
This explains perhaps better the essence of a conservative force
field. A conservative field can be represented as the gradient of
a potential: a = grad F, where F is the potential of the field. The
Earth’s gravity vector field g ( x, y, z) is the gradient of the Earth’s
gravity potential W ( x, y, z). At mean sea level — more precisely, at
the geoid — the gravity potential is constant; the gravity vector g
stands everywhere orthogonally on the geoid.
Ó »îá .
302 Field theory and vector calculus — core knowledge
n
a
∂V
div a
Figure A.7. The Gauss integral theorem. n is the normal vector to the exte-
rior surface. Note that the Gauss intregral theorem can also be
formulated with the aid of (Michael Faraday’s) field lines: a field
line starts or terminates on an electric charge (i.e., a place where
D div a ̸= 0) or runs to infinity (i.e., through the surface ∂V ).
Ó »îá .
The continuity of matter 303
d
per unit of volume. The second term again, ρ , describes the change
dt
in the amount of matter inside the volume element over time. The two
terms must balance for the “matter accounting” to close.
If the moving fluid is incompressible, then ρ is constant and
d ( )
ρ=0 and div ρ v = ρ div v,
dt
and we obtain
ρ div v = 0 =⇒ div v = 0.
Remember, however, that not necessarily rot v = 0 — i.e., the flow isn’t
necessarily eddy free —, i.e., a potential F for which v = grad F does not
necessarily exist.
Ó »îá .
D Function spaces
B
– 305 –
306 Function spaces
it is easily verified that the above requirements for a scalar product are
met.
One basis in this vector space (a function space) is formed by the so-called
Fourier basis functions,
−
→ 1p
e0 = 2 ( k = 0)
2
−
→
e = cos kx, k = 1, 2, 3, ... (B.4)
k
−
e→
− k = sin kx, k = 1, 2, 3, ...
This is the familiar way in which the coefficients of a Fourier series are
calculated.
D B.2.2 Example
Ó »îá .
Fourier function space 307
K =1
K = 25 K =3 K =5
1.0
a0
0.8
f (x)
0.6
X
0.4 a b
0.2 b5
b3
0
b1
−0.2
0 1 2 3 4 5 6
Figure B.1. Fourier analysis on a step function. Plotted are the Fourier expan-
def p ∑K ∑K
sions f (K ) (x) = 12 a 0 2 + k=1 b k sin kx = 12 − π2 k=1,3,5,... 1k sin kx,
for values of K of 1, 3, 5, and 25. The inset gives the spectrum of
D the function.
Ó »îá .
308 Function spaces
D B.2.3 Convergence
then ˆ 2π {
1 }2
lim f (K) ( x) − f ( x) dx = 0.
K →∞ π 0
Note that this does not mean that, for every x ∈ [0, 2π), f (K) ( x) → f ( x)
when K → ∞. Looking at figure B.1, there will always remain a small
neighbourhood of x = π where the difference f (K) ( x)− f ( x) will be almost 12 ,
even for arbitrarily large values of K . We say that the Fourier expansion
is convergent, but not uniformly convergent.
The Fourier series converges pointwise “almost everywhere” in x ∈ [0, 2π):
in all points except in the two special points x = 0 and x = π.
Also, note the “shoulder” of the expansion even for K = 25. This shoulder
will get narrower for higher K , but not any lower. It is known as the
Gibbs phenomenon.
Lx − λx = 0,
where the problem consists of determining the values λ i for which solu-
tions x i exist.
In a concrete n-dimensional vector space we may write the vector
n
∑
x= xi ei ,
i =1
Ó »îá .
Sturm-Liouville differential equations 309
when also
n
∑ n
∑ [ ]
λx = λ xi ei = λx j e j . (B.6)
i =1 j =1
〈x · Ly〉 = 〈Lx · y〉 .
〈x · A y〉 = 〈 A x · y〉
i.e.,
n n n
[ ] [ n ]
∑ ∑ ∑ ∑
xi ai j yj = a i j x j yi ,
i =1 j =1 i =1 j =1
which is trivially true if, for all i, j values 1, . . . , n,
a i j = a ji , i.e., A = A T .
In other words:
1
An advantage of the numerical representations is of course that one can really
calculate with them.
Ó »îá .
310 Function spaces
A symmetric matrix is a self-adjoint operator.
Lx p = λ p x p ,
i.e.,
( )⟨ ⟩
λ p − λq x p · x q = 0.
⟨ ⟩
If λ p ̸= λ q , we thus must have x p · x q = 0, or x p ⊥ x q . What was to be
proven.
σ2x σ x y
[ ]
Var (xP ) = ΣPP = ,
σ x y σ2y
a symmetric matrix. Here, σ2x and σ2y are the variances, or squares
of the mean errors, of the x and y co-ordinates, whereas σ x y is the
covariance between the co-ordinates.
The eigenvalues of this matrix ΣPP are the solutions of the charac-
teristic equation
σ2x − λ σx y
[ ]
det = 0,
σx y σ2y − λ
2
Actually the eigenvectors may be arbitrarily re-scaled: if x is an eigenvector,
def
then also e = ∥xx∥ is. Thus we obtain an orthonormal basis.
Ó »îá .
Sturm-Liouville differential equations 311
i.e.,
σ2x − λ σ2y − λ − σ2x y = 0.
( )( )
This yields
1( 2 ) 1 [ √]2
σ x + σ2y ±
[ ]
λ1,2 = σ2x + σ2y − 4 σ2x σ2y − σ2x y =
2 2√
1( 2 1 [ 2 ]2
σ x + σ2y ±
)
= σ x − σ2y + 4σ2x y .
2 2
The variance matrix has a variance or error ellipse. The semi-
p p
lengths of its principal axes are λ1 , λ2 and the directions of
the principal axes are the eigenvectors of ΣPP , x1 , x2 , mutually
orthogonal. If the co-ordinate axes are turned into the directions of
x1,2 , then the matrix ΣPP will assume the form
σ2x′ λ1 0
[ ] [ ]
0
Σ′PP = = .
0 σ2y′ 0 λ2
d2
2
x ( t ) − ω 2 x ( t ) = 0. (B.8)
dt
The solution has the general form (α amplitude, φ phase constant)
( )
x ( t) = α sin tω − φ .
⏐ ⏐
d ⏐⏐ d ⏐⏐
On the interval t ∈ [0, T ] we require x (0) = x (T ) , x⏐ = x , i.e.,
dt x=0 dt ⏐ x=T
periodicity. These boundary conditions are an essential part of being
self-adjoint. Then, a solution is found only for certain values of ω —
quantization.
Note that equation B.8 is an eigenvalue problem, form-wise:
Lx − ω2 x = 0,
Ó »îá .
312 Function spaces
As on the right-hand side the first terms vanish and the second terms are
identical, we have
⟨−
→x · L−
→y = L−
⟩ ⟨ → −
x ·→
⟩
y ,
which was to be proven.
Self-adjoint operators have eigenvalues and eigenvectors, in this case
functions, that are mutually orthogonal for different values of ω3 . For the
oscillation equation with the above periodicity conditions they are just
the solution functions
2π k
( )
( )
sin ωk t − φ = sin t−φ , (B.9)
T
3
In fact, for the same value ω there exist two mutually orthogonal periodic
solutions,
2π kt
sin ωk t = sin ,
T
2π kt
cos ωk t = cos .
T
Any linear combination of these is also a valid solution, and is of the general
form B.9.
4
Jacques Charles François Sturm FRS FAS (1803 – 1855) was an eminent French
mathematician, one of the 72 names engraved on the Eiffel Tower.
5
Joseph Liouville FRS FRSE FAS (1809 – 1882) was an eminent French mathe-
matician.
Ó »îá .
Legendre polynomials 313
for values n = 0 . . . ∞, m = − n . . . n.
In this function space, a scalar product is defined:
¨
⟨−
→ ⟩ 1
v ·−
→ ( ) ( )
w = V φ, λ W φ, λ d σ ,
4π σ
Ó »îá .
314 Function spaces
see also Heiskanen and Moritz (1967). Proving this with the help of
equation 2.13 is a long process.
If we now divide the functions Ynm (or, equivalently, R nm , S nm ) by the
square roots of the above factors, we obtain the fully normalized surface
spherical harmonics Y nm .
With those it is again easy to calculate the coefficients f nm of a given gen-
( )
eral function f φ, λ (the overline means that these are fully normalized
coefficients):
¨
⟨ −−→⟩ 1 ( ) ( )
f nm = f · ynm = f φ, λ Y nm φ, λ d σ. (B.10)
4π σ
Ó »îá .
Self-test questions 315
in which {
a nm m ≥ 0,
ynm =
b n| m| m < 0.
D Self-test questions
Ó »îá .
D Why does FFT work?
C
in which
N −1 ( )
( ) 1 ∑ jk
F ωj = f ( xk ) exp 2π i , j = 0, . . . , N − 1. (C.1)
N N
k=0
is
N −1 ( )
∑ ( ) jk
f ( xk ) = F ωj exp −2π i , k = 0, . . . , N − 1. (C.2)
N
j =0
1
For simplicity’s sake. An arbitrary interval will work.
– 317 –
318 Why does FFT work?
FFT is just a very efficient method for computing both these formulas
C.1, C.2. A brute-force calculation of these formulas requires of order N 2
“standard operations”, each of them a single multiplication plus a single
addition or subtraction. If N is even, we may write
⎡N ⎤
2 −1 ( ) N −1 ( )
( ) 1 ∑ jk ∑ jk ⎦
F ωj = ⎣ f ( xk ) exp −2π i + f ( xk ) exp −2π i =
N N N
k=0 k= N
2
⎡N N
⎤
2 −1 ( ) ( 2 −1 (
N)∑ ( ′
)
1 ∑ jk j ) jk ⎦
= ⎣ f ( xk ) exp −2π i + exp −2π i 2 f xk′ + N exp −2π i =
N N N 2 N
k=0 k′ =0
⎡N N
⎤
2 −1 ( ) 2 −1 ( ( )
1 ∑ jk ∑ ) jk ⎦
= ⎣ f ( xk ) exp −2π i + exp (−π i j ) f x k+ N exp −2π i =
N N 2 N
k=0 k=0
N
2 −1 [ )] ( ) ( )
1 ∑ j ( jk + if j even
= f ( x k ) ± f x k+ N exp −2π i , (C.3)
N 2 N − if j odd
k=0
Ó »îá .
D Helmert condensation
D
1
Uniform convergence means that, given r, r ′ , for every ϵ > 0 there is an Nmin
for which
⏐ ⏐
N ( )
r n+1
⏐1 1 ∑ ( )⏐⏐
⏐ − P n cos ψ ⏐ < ϵ
⏐
⏐ℓ r r′ ⏐
n=0
for all N > Nmin , and for all values of ψ. This is a stronger property than mere
convergence.
– 319 –
320 Helmert condensation
¨ ∑
⎡ −(n−2) ⎤
∞ n R −
r ⎣
− (R + H )−(n−2) +
( )
= Gρ ⎦ P n cos ψ d σ.
σ n=0 n−2
n̸=2
+ r 2 ln R +
R
H
H ( n − 2) ( n − 1) H 2
⎡ ⎤
1 − ( n − 2) + −
(R + H )−(n−2) = R −(n−2) ⎣
⎢ R 2 R2 ⎥ ⎦.
3
( n − 2) ( n − 1) n H
− + . . .
2·3 R3
Also the special case n = 2,
H 1 H2 1 H3 1 H4
[ ]
2 R+H 2
r ln = r − + − +... =
R R 2 R2 3 R3 4 R4
H n − 1 H 2 ( n − 1) n H 3
⎡ ⎤
r ⎢ R − 2 R2 + 2 · 3 R4 − ⎥
n
= ⎦,
R n−2 ( n − 1) n( n + 1) H 4
⎣
− +...
2·3·4 R4
is cleanly included into the following expression obtained by substitution:
H n − 1 H2
⎡ ⎤
¨ ∑
∞
r ⎢ R − 2 R2 +
n
Ttint r, φ, λ = G ρ
( ) ⎥ ( )
n−2 ⎣ 3 ⎦P n cos ψ d σ.
σ n=0 R ( n − 1) n H
+ −...
6 R3
(D.1)
Ó »îá .
The exterior potential of the condensation layer 321
H ( n + 3) ( n + 2) H 2
⎡ ⎤
⎢ 1 + (n + 3) R + 2
+
R2 ⎥
(R + H )n+3 = R n+3 ⎣ 3 ⎦.
( n + 3) ( n + 2) ( n + 1) H
+...
2·3 R3
Substitution yields
H n + 2 H2
⎡ ⎤
¨ ∑
∞ ( ) n+1
R ⎢ R + 2 R2 +
Ttext = G ρ R 2
⎥ ( )
3 ⎦ P n cos ψ d σ.
r ( n + 2) ( n + 1) H
⎣
σ n=0
+...
6 R3
(D.2)
This is thus the exterior potential of the topography, or, inside the topo-
graphic masses, the harmonic downward continuation of the exterior
potential, assuming that this is mathematically possible (in the case of
mountainous topography, generally not) and doesn’t diverge.
ext
THelmert = Ttext − Tcext =
Ó »îá .
322 Helmert condensation
⎡ {
n+2 } H2 ⎤
¨ ∑
∞ ( ) n+1 −1 +
R 2 R2
= GρR2
⎥ ( )
⎦ P n cos ψ d σ =
⎢
3
r ( n + 2) ( n + 1) H
⎣
σ n=0
+ +...
6 R3
¨ ∑∞ ( ) n+1 [
n 2 ( n + 2) ( n + 1) H 3
]
R ( )
= Gρ H + + . . . P n cos ψ d σ.
σ r 2 6 R
n=0
(D.4)
Then
∞ ( ) n+1 [
( n + 2) ( n + 1) H n3
]
ext
∑ R n
= 4πG ρ
THelmert H n2 + +... .
r 2 (2 n + 1) 6 (2 n + 1) R
n=0
(D.6)
If the topography is constant, then all terms for which n ̸= 0 vanish; in
the above expression, also the first term vanishes, and the second and
further terms are very small. So, in practice2
ext 4 3 Gρ
THelmert = πH · +... ≈ 0
3 r
as was to be expected. The exterior field remains in this special case
almost unchanged as a result of Helmert condensation.
Let us calculate gravity anomalies from the Helmert potential, but only
using the first term of equation D.6:
∂T 2
∆ g Helmert = + T=
∂r r
∞ [ ]( )n+1
∑ n − ( n + 1) 2 R
= 2πG ρ + H n2 =
2n + 1 r r r
n=0
∞ ( )n+1
1 ∑ n ( n − 1) R
= −2πG ρ H n2 .
r 2n + 1 r
n=0
2
As a curiosity, this result can be interpreted as the potential of a sphere of
crustal matter with radius H (the average height of the topography) located at
the geocentre. Even for a topographic mean height of 10 km the effect on the
geoid would only be 12 mm (exercise!).
Ó »îá .
The total potential of Helmert condensation 323
Ó »îá .
324 Helmert condensation
By substituting into this the equation D.8 for the µQ double mass density
layer, we obtain, by taking the limit r P , r Q ↓ R :
∞ ¨
1 ∑ ( ) ( )
T = n 2πGH ρ HP n cos ψ d σ =
4π σ
n=0
∞ ¨
1 ∑ ( )
= n A B HP n cos ψ d σ.
4π
n=0
Here, we have left the designations P,Q again off as no longer needed for
clarity.
The symbol A B denotes the attraction of a Bouguer plate of thickness H
and density ρ .
Let us develop the quantity [ A B H ] into a spherical-harmonic expan-
sion (Heiskanen and Moritz, 1967 equation 1-71). then, each degree
constituent is (equation 2.16):
¨
2n + 1 ( )
[ A B H ]n = [ A B H ] P n cos ψ d σ,
4π σ
3 1
In fact, a better place for this layer would be level 4 H. This is one of the
approximations made here.
Ó »îá .
The dipole method 325
T 1 A B H πG ρ H 2
δ NHC = ≈ = .
γ 2 γ γ
∞
4πG ρ 1 ∑ n + 1 2 πG ρ H 2
δ NHC ≈ · H ≈ ,
γ 2 2n + 1 n γ
n=0
Ó »îá .
D The Laplace equation in spherical
E co-ordinates
D E.1 Derivation
Consider a small volume element with sizes in co-ordinate directions of
def
∆ r , ∆φ, ∆λ. Look at the difference in flux of vector field a = ∇V between
what comes in and what goes out through opposite faces.
We do the analogue of what was shown in section 1.12.4, using a body
with surfaces aligned along co-ordinate lines, allowing the size of the
body to go to zero in the limit, and exploiting the Gauss integral theorem
1.18. The quantity div a = ∆V is a volume density, and its average value
multiplied by the volume of a body must equal the total flux through the
surface of the body.
Part of the difference in flux f between opposing faces is due to a change
in the normal component of a between the faces, part is due to a difference
in face surface area ω:
f + − f − ≈ ω a + − a − + a ω+ − ω− .
( ) ( )
difference
r − ω r ≈ 2 r ∆ r cos φ∆φ∆λ.
ω+ −
– 327 –
328 The Laplace equation in spherical co-ordinates
r ∆φ
r cos φ ∆r
r r ∆λ cos φ
∆φ
∆λ
λ
Equatorial plane
φ = r cos φ + ∆φ ∆ r ∆λ,
ω+
( )
φ = r cos φ∆ r ∆λ,
ω−
difference
φ − ωφ ≈ − r sin φ∆φ∆ r ∆λ.
ω+ −
∂V
Multiply by a φ = ∂rφ
and divide by body volume r 2 cos φ∆ r ∆φ∆λ,
yielding
tan φ ∂V
∆φ′ V = − .
r 2 ∂φ
This of course in addition to the regular contribution
a+ −
φ − aφ ∂V + ∂V − 1 ∂2 V
( )
1
∆φ V = = − ≈ .
r ∆φ r ∆φ ∂rφ ∂rφ r 2 ∂φ2
Ó »îá .
Solution 329
∆V = ∆r V + ∆λ V + ∆φ V + ∆r′ V + ∆φ′ V =
∂2 V 1 ∂2 V 1 ∂2 V 2 ∂V tan φ ∂V
= + + + − 2 ,
∂r2 r 2 cos2 φ ∂λ2 r 2 ∂φ2 r ∂ r r ∂φ
equivalent to equation 2.7.
D E.2 Solution
D E.2.1 Separating the radial dependency
∂2 R ∂R 1 ∂2 Y ∂2 Y ∂Y
( ) ( )
1 1
r2 2 + 2r =− + − tan φ .
R ∂r ∂r Y cos φ ∂λ
2 2 ∂φ 2 ∂φ
This must again apply for all values r, φ, λ and thus can only be a constant,
p. This yields two equations
2
2∂ ∂R
( )
R
r + 2r − pR = 0,
∂r 2 ∂r
1 ∂2 Y ∂2 Y ∂Y
( )
+ − tan φ + pY = 0.
cos φ ∂λ
2 2 ∂φ 2 ∂φ
For the first equation we try a power law,
R = rq,
yielding
q ( q − 1)) r q + 2 qr q − pr q = 0
with the solution
p = q ( q + 1) .
( )
Solving the second equation for Y φ, λ ,
1 ∂2 Y ∂2 Y ∂Y
( )
+ − tan φ + q ( q + 1) Y = 0, (E.1)
cos φ ∂λ
2 2 ∂φ 2 ∂φ
is trickier. It turns out that q must be an integer. One finds, for n ∈
N0 , that there are positive solutions1 q = n and negative solutions q =
− ( n + 1), with n = 0, 1, 2, . . . With this, the full set of special solutions is
( )
Y n φ, λ
Vn,1 = r n Yn φ, λ , Vn,2 = , n ∈ N0 .
( )
r n+1
i.e., equations 2.8.
1
The above derivation logic doesn’t work for the exponent values q = 1 or q = 0,
but the result holds also for these values.
Ó »îá .
330 The Laplace equation in spherical co-ordinates
Note that both solutions q, the positive and the negative one, produce on
substitution into equation E.1 the same equation for n.
1 ∂2 Y ∂2 Y ∂Y
( )
+ − tan φ + n ( n + 1) Y = 0.
cos φ ∂λ
2 2 ∂φ 2 ∂φ
cos2 φ
Substitution and multiplication by FL yields
cos2 φ ∂2 F ∂F 1 ∂2 L
( )
− tan φ + n ( n + 1) F =− .
F ∂φ2 ∂φ L ∂λ2
Both sides must be again equal to the same constant, which we shall
assume positive and call m:
∂2 F ∂F m2
( )
− tan φ + n ( n + 1) − F = 0,
∂φ2 ∂φ cos2 φ
∂2 L
+ m 2 L = 0.
∂λ2
The first equation is known as Legendre’s equation. Its solutions are the
( )
Legendre functions P nm sin φ , with the integer m = 0, 1, 2, . . . , n. The
second is the classical harmonic oscillator, with solutions2
2
This also explains why m must be an integer.
Ó »îá .
Bibliography
Ole Balthasar Andersen, Per Knudsen, and Philippa Berry. The DNSC08GRA
global marine gravity field from double retracked satellite altimetry. Journal
of Geodesy, 84(3), 2010. DOI: 10.1007/s00190-009-0355-9. 176
Tulu Besha Bedada. Absolute geopotential height system for Ethiopia. PhD
thesis, University of Edinburgh, 2010. URL: https://www.era.lib.ed.ac.uk/
bitstream/1842/4726/3/Bedada2010_small.pdf. 176
Gian Paolo Bottoni and Riccardo Barzaghi. Fast collocation. Bulletin géodésique,
67:119- 126, 1993. 207
– 331 –
332 Bibliography
Martin Ekman. Postglacial Rebound and Sea Level Phenomena with Special
Reference to Fennoscandia and the Baltic Sea. In Juhani Kakkuri, editor,
Geodesy and Geophysics, lecture notes, NKG Autumn School. Finnish Geodetic
Institute Publication 115, 1992. 241, 283
Óîá .
Bibliography 333
René Forsberg and Jānis Kaminskis. Geoid of the Nordic and Baltic region from
gravimetry and satellite altimetry. In J. Segawa, H. Fujimoto, and S. Okubo,
editors, Proceedings, IAG International Symposium on Gravity, Geoid and
Marine Geodesy (GraGeoMar96), Tokyo, Sep. 30 – Oct. 5, 1996, number 117 in
International Association of Geodesy Symposia, pages 540 – 547, Tokyo, 1996.
International Association of Geodesy, Springer Verlag. 288
Markku Heikkinen. Solving the Shape of the Earth by Using Digital Density
Models. Report 81:2, Finnish Geodetic Institute, Helsinki, 1981. 73, 75, 227
R.A. Hirvonen. The continental undulations of the geoid. PhD thesis, Helsinki
University of Technology, 1934. Finnish Geodetic Institute Publication 19.
289
Óîá .
334 Bibliography
Erkki Hytönen. Absolute gravity measurement with long wire pendulum. PhD
thesis, University of Helsinki, 1972. Finnish Geodetic Institute Publication
75. 214
J. Kakkuri, editor. Geodesy and Geophysics, lecture notes, NKG Autumn School
1992, Finnish Geodetic Institute Publication 115, 1981. 165
Juhani Kakkuri. Surveyor of the Globe – Story of the Life of V.A. Heiskanen.
Maanmittauslaitos, 2008. URL: https://readymag.com/u95015526/508134/.
107, 289
P. Melchior. The Tides of the Planet Earth. Pergamon Press, Oxford, 1978. 280
Óîá .
Bibliography 335
Silja Märdla. Regional Geoid Modelling by the Least Squares Modified Hotine
Formula Using Gridded Gravity Disturbances. PhD thesis, Tallinn University
of Technology, 2017. ISBN 978-9949-83-186-9. 96, 151
R.L. Parker. The rapid calculation of potential anomalies. Geophys J R Astr Soc,
31:447 – 455, 1972. http://topex.ucsd.edu/geodynamics/parker.pdf. 179
N.K. Pavlis, S.A. Holmes, S.C. Kenyon, and J.K. Factor. An Earth Gravitational
Model to Degree 2160: EGM2008. In EGU General Assembly 2008, Vienna,
Austria, April 13-18, 2008, 2008. http://earth-info.nga.mil/GandG/wgs84/
gravitymod/egm2008. 55
W.R. Peltier. Closure of the budget of global sea level rise over the GRACE
era: the importance and magnitudes of the required corrections for
global glacial isostatic adjustment. Quaternary Science Reviews, 28(17 –
18):1658 – 1674, 2009. ISSN 0277-3791. doi: http://dx.doi.org/10.1016/j.
quascirev.2009.04.004. URL http://www.sciencedirect.com/science/article/pii/
S0277379109001218. Quaternary Ice Sheet-Ocean Interactions and Land-
scape Responses. 244
John Henry Pratt. The attraction of the Himalaya mountains upon the plumbline
in India. Philosophical Transactions of the Royal Society of London, CXLV:
53 – 100, 1855. 106
Óîá .
336 Bibliography
John Henry Pratt. The deflection of the plumb-line in India and the com-
pensatory effect of a deficiency of matter below the Himalaya mountains.
Philosophical Transactions of the Royal Society of London, CXLIX:745 – 778,
1858. 106
John Henry Pratt. Speculations on the constitution of the Earth’s crust. Pro-
ceedings of the Royal Society of London, XIII:253 – 276, 1864. 106
W. Torge. Gravity and Tectonics. In Geodesy and Geophysics, lecture notes, NKG
Autumn School, Finnish Geodetic Institute Publication 115, 1992. 285
Óîá .
Bibliography 337
Peter Vaníček and Edward Krakiwsky. Geodesy – The Concepts. Elsevier Science
Publishers, Amsterdam, 1986. 289
Martin Vermeer. Geoid studies on Finland and the Baltic. PhD thesis, University
of Helsinki, 1984. 26, 83
Martin Vermeer. Terrain Reduction and Gridding Techniques for Geoid Deter-
mination. In Geodesy and Geophysics, lecture notes, NKG Autumn School,
Finnish Geodetic Institute Publication 115, 1992. 179
Heikki Virtanen and Jussi Kääriäinen. The installation and first results from
the superconducting gravimeter GWR20 at the Metsähovi station, Finland.
Report 95:1, Finnish Geodetic Institute, Helsinki, 1995. 15 p. 225
Phuong Lan Vu, Frédéric Frappart, José Darrozes, Vincent Marieu, Fabien
Blarel, Guillaume Ramillien, Pascal Bonnefond, and Florence Birol. Multi-
Satellite Altimeter Validation along the French Atlantic Coast in the Southern
Bay of Biscay from ERS-2 to SARAL. Remote Sensing, 93(10), 2018. URL:
http://www.mdpi.com/2072-4292/10/1/93/pdf. 266
Lianxing Wen and Don L. Anderson. Layered mantle convection: A model for
geoid and topography. Earth and Planetary Science Letters, 146(3 – 4):367 –
377, 1997. ISSN 0012-821X. doi: http://dx.doi.org/10.1016/S0012-821X(96)
00238-5. http://222.195.83.195/wen/Reprints/WenAnderson97EPSL.pdf. 114
L. Wong and R. Gore. Accuracy of geoid heights from modified Stokes kernels.
Geophys J R Astr Soc, 18:81 – 91, 1969. 155
Óîá .
338 Bibliography
Óîá .
Index
– 339 –
340 Index
Óîá .
Index 341
Óîá .
342 Index
Óîá .
Index 343
Óîá .
344 Index
Óîá .
Index 345
Óîá .
346 Index
Óîá .
Index 347
Óîá .
348 Index
Óîá .
Index 349
Óîá.
350 Index
Óîá .
Index 351
Óîá .
352 Index
Óîá .
Index 353
Óîá .
354 Index
Óîá .
Index 355
Óîá .
356 Index
Óîá .
Index 357
Óîá .