Geophysics 210 September 2008

B1 – Size and shape of the Earth

B1.1 Is the Earth flat?

Speculations …….

• Once believed that the Earth is flat and that ships could sail over the edge.
• View persisted into the Middle ages and was an issue in recruitment for Columbus.
• The Flat Earth Society still alive and well http://theflatearthsociety.org/ and describe
the Earth as “being a disk with a circumference of about 78225 miles and a diameter of 24900
miles. The sun and moon are both disks about 32 miles in diameter …and are about 3000 miles
above the Earth, and the stars about 100 miles above the sun and moon. The Flat Earth Society
also maintains that the Earth is accelerating upward at a rate of 9.8 m/s², thereby simulating
gravity. This upward momentum is caused by the "Universal Accelerator", a vague term used by
the Society to describe a force that originated at the Big Bang and caused the Earth to speed
upwards. Gravity cannot exist on a flat Earth since the disc shape would eventually collapse on
itself. However, other planetary bodies such as the moon and the sun have gravitational pulls,
causing the gravitational force on an object to decrease as it increases in altitude. This also
allows spacecraft to orbit.” http://en.wikipedia.org/wiki/Flat_Earth_Society

B1.2 Historical perspective

• Early Greek view was that the world was surrounded by Oceanus, origin of all rivers.
• Anaximander (600 B.C.) – cylindrical earth surrounded by celestial sphere

• Pythogoras (582-507 B.C.) believed

the Earth was a sphere, which was
considered the most harmonious
geometric shape.

• Aristotle (384-322 B.C.) described observations that supported the theory that the
Earth was a sphere. These included (1) the fact that the shadow of the moon is

1
 Geophysics 210 September 2008

circular in lunar eclipses (2) constellations were higher in the sky as one traveled
south.
• More historical details of ancient Indian, Armenian and Islamic studies can be found
in the article : http://en.wikipedia.org/wiki/Spherical_Earth

Observations that suggest the Earth is a sphere

• Mountain peaks lit by the Sun after

sunset.
• Ships disappear below the horizon as
they sail across ocean.
• The moon looks like a disk. Is the
Earth the same shape?
• The Earth casts a circular shadow
during lunar eclipses.

B1.3 Quantitative approach to computing the radius of the Earth

Eratosthenes (275-195 B.C.)

Estimated size of Earth from

observations that the elevation of the
sun varied with position on the
Earth’s surface in Egypt.

Measurement 1: Angular distance Aswan - Alexandria

On the summer solstice, the sun was overhead in Aswan, since it illuminated a
deep well. On the same day, the sun was at an angle of 7° 14’ to the vertical in
Alexandria.

Measurement 2 : Distance on ground Aswan – Alexandria

The distance Aswan-Alexandria was estimated at 5000 stadia = 925 km

Combining these measurements, can show that the circumference of the Earth = 360
*925/7.23 = 46058 km. Correct value = 40030 km. Error = 15%

2
 Geophysics 210 September 2008

Jean Picard (1620-1682) : Measured the length of ~1 degree of latitude in France in

1669-70 and obtained a value of 6329 km for the radius of the Earth. Polar radius
known to be 6357 km which is an error of 0.4%
(http://en.wikipedia.org/wiki/Jean_Picard)

Another approach - the double sunset

• A seated observer watches the sunset over the horizon (A) at the Equator. As soon
as the sun sets, he jumps to his feet and enjoys a few extra seconds of sunlight
before the sun sets again (B).

(A) (B)

The observer was 2 m tall and the two sunsets were separated by 7.5 s. Through
what angle (θ) does the Earth rotate in this time?

θ = ______ degrees = ______ radians

h
The radius of the Earth can be shown to be R = where the angle θ is measured
θ2
in radians and the observer is h m tall.

R = _______ m

B1.4 Non-spherical Earth

It was soon realized that the Earth was not a perfect sphere. Jean Richer (1620-
1682) was a French astronomer who made observations in 1671 from the island of
Cayenne in French Guyana. In addition to measuring the distance of planet Mars from
Earth, he also noted that a pendulum swung slower on Cayenne Island than in Paris.
He deduced that this was because gravity was weaker on Cayenne, which implied that
Cayenne was further from the centre of the Earth than Paris.

Isaac Newton (1642-1727) suggested that the Earth was an oblate spheroid because it
rotates - somewhat flattened at the poles.

3
 Geophysics 210 September 2008

A debate followed about whether this was true, because measurements were not very
accurate. Surveys in Peru (equator, 1735-43) and (equator, 1735-43) showed that
flattening occurred as suggested.

The deformation is now described by

the International Reference Ellipsoid.
This is the shape the Earth would
have if it’s composition was uniform.
A number of ellipsoids have been
defined as data quality has improved.

Ellipsoid Semi-major Semi-minor Inverse flattening

reference axis a axis b (1/f)

GRS 80 6,378,137.0 m ≈ 6,356,752.314 140 m 298.257 222 101

WGS 84 6,378,137.0 m ≈ 6,356,752.314 245 m 298.257 223 563

WGS = World Geodetic System

More Details at

http://en.wikipedia.org/wiki/Figure_of_the_Earth

http://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS84

Density variations mean that there are discrepancies between the surface and the
reference ellipsoid. A surface called the geoid takes this into account and will be
discussed in B1.3.

Implications of a non-spherical Earth : Chimborazo in Ecuador is higher than

Mount Everest, if measured from the centre of the Earth.

4
 Geophysics 210 September 2008

B2. BASICS OF GRAVITY EXPLORATION

B2.1 Newtonian gravitation

Consider two point masses that are a distance r apart. Newton’s theory of gravitation
predicts that they will attract each other with a force F that is given by:

Gm1 m2
F=
r2

The quantity G is called the gravitational constant (or “big G”) but is actually a small
number. Newton deduced this equation from observing the motion of planets and moons
in the solar system. The units are as follows:

F Newton (N)
m kg
r metres (m)
G = 6.67 x 10-11 N m2 kg-2

It can be shown that if one of the masses is finite in size (e.g. a planet or the sun), then F due to
will be the same as if all the mass were concentrated at the centre.

If m1 = m2 = 1,000,000 kg (1000 tonnes) then the variation of F as r increases will show the
“inverse square law”.
 Geophysics 210 September 2008

Consider the mass m2. Newton’s Third Law of Motion predicts that this mass will accelerate
with an acceleration a, where:

F = m2 a

Rearranging this equation gives

F
a=
m2

Now we know the value of F from the first equation so

Gm1 m2 1 Gm
a=g= 2
= 21
r m2 r

where g is called the gravitational acceleration.

Let us now consider that m1 is the Earth and

m2 is a small object that we are going to
drop. This equation tells us that the
acceleration does not depend on the mass
of the object being dropped.

This was proved by Galileo who allegedly

dropped masses from the leaning tower of
Pisa in Italy. This result says that a small
mass and a large mass will fall with the
same acceleration.
 Geophysics 210 September 2008

B2.2 Density of rocks and minerals

The Greek letter rho (ρ) is used to represent density. Rocks and minerals found on Earth have
densities that range from 1000-7000 kg m-3. Often densities are quoted in g cm-3.

To convert, remember that 1 g cm-3 = 1000 kg m-3.

Pure minerals can exhibit a high density since the atoms are closely packed together.

Magnetite ρ = 4.90-5.20 g cm-3

Pyrite ρ = 4.90-5.20 g cm-3
Galena ρ = 7.40-7.60 g cm-3

Sedimentary rocks generally have lower densities since the atoms are not as closely packed
together and pore space is filled with lower density materials.

Water ρ = 1.00-1.05 g cm-3

Alluvium ρ = 1.96-2.00 g cm-3
Shale ρ = 2.00-2.70 g cm-3
Limestone ρ = 2.60-2.80 g cm-3
Dolomite ρ = 2.28-2.90 g cm-3

The range of density values reflects the degree of weathering and the porosity. Why?

The composition is secondary, but explains why dolomite is more dense than limestone.

Most sedimentary rocks increase in density with depth, owing to increased compaction.
This reduces the pore space available for low density materials such as air, water and
hydrocarbons.

Igneous rocks are generally more dense owing to minimal porosity.

Granite ρ = 2.50-2.70 g cm-3

Basalt ρ = 2.70-3.20 g cm-3

The density depends primarily on the rock composition. Mafic rocks are generally
more dense than felsic rocks owing to increased proportion of heavier elements
such as Fe and Mg.

B2.3 Units for gravity measurements (milligals)

The gravitational acceleration at the Earths surface is 9.8 ms-2. Subsurface variations in rock
density produce very small changes in this value, so it is more convenient to use a smaller unit.

9.8 ms-2 = 980 cm s-2

= 980 gal (after Galileo)
= 980,000 milligals

1 milligal = 10-5 ms-2

 Geophysics 210 September 2008

B2.4 Approximate calculation to estimate the size of gravity anomalies

Consider a spherical ore body with density ρ and radius a that is buried at a depth z below the
surface. The extra pull of gravity will be greatest at a point P directly above the ore body.

Remember that the pull of gravity (g) for a sphere is the same as if all the mass were
concentrated at the centre.

4 3
Total mass of ore body = volume x density = πa ρ
3

Excess mass , mE = mass of ore body – mass of rock that was already there
4 3 4
= πa ρ - πa 3 ρ 0
3 3
4 3
= πa ( ρ − ρ 0 )
3

Above ore body, the change in gravity due to the ore body (gravity anomaly) is given by

Gm E 4Gπa 3 ( ρ − ρ 0 )
Δg = =
z2 3z 2

Now guess some values

a = 30 m
z = 40 m
ρ = 4000 kg m-3
ρ0 = 2000 kg m-3

Careful use of a calculator gives Δg = 0.94 mgal

Perspective

This can be compared to value of g = 980,000 mgals

Fractional change = 9.6 10-7

e.g. An 80 kg person would feel 0.08 g lighter!!!!!

MJU 2008
 Geophysics 210 September 2008

B3 Variation of gravity with latitude and elevation

By measuring the subtle changes in the acceleration of gravity from one place to another, it is
possible to learn about changes in subsurface density.

However, other factors can cause gravity to vary with position on the Earth. These effects must
be removed from measurements in order to use gravity data to study the interior of the Earth.

B3.1 Variation of gravity with latitude

It is observed that at the Equator, g E = 978,033 mgal while at the poles g P = 983,219 mgal

This difference is 5186 mgal, which is a lot larger than changes in gravity because of
subsurface density.

Can this observation be explained by the fact that the Earth is a rotating ellipsoid?

(A)The Earth is distorted by rotation

The Earth is an
oblate spheroid.

R E = 6378 km

R P = 6357 km.

Qualitative answer

Since a point on the Equator is further from the centre of the Earth than the poles, gravity
will be weaker at the Equator and g E < g P

Quantitative answer

GM E
For a sphere g (r) = 2
where the mass of the Earth, ME = 5.957 1024 kg.
r
At the North Pole, RP = 6357 km and g P = 983,219 mgal.

If we move up 21 km to the equator, the decrease in gravity will be 6467 mgal

Thus g E = g P - 6467 mgal, which is too much to explain the observed difference
between the Equator and the Poles.

1
 Geophysics 210 September 2008

(B) - Centrifugal forces vary with latitude

The rotation of the Earth also causes gravity to vary with latitude.

Qualitative answer

Imagine you are standing at the North Pole. The rotation of the Earth will not change g, all
that will happen is that you rotate once a day.

Now imagine you are at the Equator. If we could increase the rotation rate of the Earth
enough, you would be ultimately be thrown into space (i.e. become weightless). Thus rotation
makes gravity weaker at the equator.

Quantitative answer

An observer is at a point with latitude θ.

This observer travels around the rotation
axis in circle with radius r = R cos θ.

The rotation rate is ω (radians per sec).

This corresponds to a radial acceleration

a = r ω 2 oriented towards the rotation
axis.

Assume RE = 6378 km

In a “vertical” direction (defined as pointing towards the centre of the Earth) this has a
component :

ar = a cos θ = Rω2 cos2 θ

Now the Earth rotates once per day so ω = 2π / (24 x 60 x 60) = 7.27 10-5 rad s-1

At the North Pole θ = 90º ar = 0

At the Equator, θ = 90º ar = 0.03370 m s-2 = 3370 mgal

g P = g E +3370 mgal

gE < gP

Question : You are standing on the Equator. How fast would the Earth need to rotate to throw
you into space?

2
 Geophysics 210 September 2008

(C) Mass distribution of the Earth

These two factors both make g E < g P so to get the observed difference we need to find a
factor that has the opposite effect. The change in shape from a sphere to an ellipsoid
redistributes the Earth’s mass. Thus results in more mass between points on the Equator and
the centre of the Earth, than between the poles and the centre of the Earth.

Qualitative answer

Consider the case shown above where the flattening is extreme. The observer at the pole
experiences the pull of gravity in all directions, and this almost cancels out. An observer on
the equator only experiences the pull of gravity due to mass located to the left.

Thus this effect will make g E > g P

Quantitative answer

Integration shows that g E ~ g P + 4800 mgal

3
 Geophysics 210 September 2008

Overall variation of g with latitude

Combining these three effects (A,B and C) gives

g P = g E + 6467+ 3370 - 4800 mgal = g E +5037 mgal ( approximately as observed)

These factors are represented in the following equation, that defines the variation of g with
latitude θ

g (θ) = 9.78031846 (1+ 0.0053024 sin ² θ – 0.0000058 sin² 2θ)

This equation is called the Geodetic Reference System for 1967.

More recent revisions are essentially the same, but with ever more decimal places …..

Calculation 1 - What value does this equation predict for Edmonton?

In Edmonton θ = 53˚ 30’ 25” N and the GRS67 equation gives

g = 9.78031846 ( 1+0.003417902-0.000005395) m s −2

= 9.81369388 m s −2

Calculation 2 – How rapidly does gravity vary in a north-south direction?

The variation of g with latitude is important when a survey extends over a significant north-
south distance. Differentiating the GRS67 equation with respect to θ yields

4
 Geophysics 210 September 2008

dg
= 9.78031846 (0.0053024 x 2 sin θ cos θ - 0.0000058 x 4sin 2θ cos 2 θ)
dθ
= 0.049526 m s −2 per radian

= 0.0008655 m s −2 per degree

= 86.550 mgal per degree

= 0.7868 mgal km −1 ( 1 degree latitude = 111 km)

All the these equations define the expected value of theoretical gravity (or normal gravity)
at latitude θ. Differences between this value and what is actually measured are anomalies that
we will analyse for information about subsurface density structure.

Calculation 3 – How much lighter would you feel after driving from Edmonton to
Calgary?

Assume the scales read 80 kg in Edmonton

Change in gravity = 300 x 0.7868 mgals

Fractional change = 300 x 0.7868 / 981369.388 = 2.34 10-4

“
Mass” in Edmonton = 80 kg > “Mass” in Calgary = 80 (1-2.34 10-4) = 79.98 kg

Change = 19 grams!!!!

5
 Geophysics 210 September 2008

210 B3.2 Variation of gravity with elevation

(i) The Free air correction

Newton’s Theory of Gravitation states that at a distance, r, from the centre of the Earth

GM
g(r) =
r2

This means that as you move away from the centre of the Earth, the acceleration of gravity (g)
decreases. In Edmonton, g = 9.81 ms −2 and if you move up a distance, Δh, the acceleration of
gravity will decrease by

Δg = 3.086 Δh x 10 −6 m s −2
= 0.3086 Δh mgal

Consider the exciting topography of a flat topped mountain:

Gravity measurements are made at points A and B. The difference in elevation means that gB will
be less than gB by an amount

Δg = 0.3086 x 100 = 30.86 mgal

When collecting gravity data, our real interest is to determine the density of the rocks below
ground. The change in elevation from ‘A’ to ‘B’ will thus contaminate the data. The Free Air
correction is a mathematical way of undoing the effect of elevation. It allows us to correct the
data collected at ‘B’ in order to make it equivalent to data collected at the same elevation as ‘A’.

In gravity surveys, we always define a reference level for the survey. Free Air corrections are
made relative to this level. In general, any reference level could be chosen, but sea level is
commonly chosen in coastal areas. In Alberta, the average level of the prairies would be a good
choice.

If a gravity measurement was made Δh above the reference level, we must add

CFA = 0.3086 Δh mgal

CFA is called the Free Air correction for a given gravity measurement.

Similarly, if a gravity measurement was made Δh below the reference level, we must subtract

CFA = 0.3086 Δh mgal

Question : to keep data accurate to 0.1 mgal, how accurately must we know the elevation?

6
 Geophysics 210 September 2008

(ii) The Bouguer correction

Unfortunately, this is not the end of story! Compare the gravity measurements at ‘A’ and ‘B’. At
point A, the gravity measurement is solely due to structure below the reference level (blue). At
‘B’ the gravity measurement is due to structure below the reference level, plus the gravitational
pull of the 100 metres of mountain (red). The net result is that gB > gA

From section B2.3 the magnitude of this extra gravitational attraction is approximately

gB - gA = 2πG ρ Δh

where ρ is the density of the mountain.

Thus to remove this effect we need to subtract CB = 2πG ρ Δh from the observed gravity
measurement at ‘B’. This is called the Bouguer correction and

CB = 0.00004193 ρ Δh mgal

Note that to apply the Bouguer correction we need to estimate ρ, the density that lies between ‘B’
and the reference level. Using the value ρ = 2670 kg m −3 this gives

CB = -0.1119 Δh mgal

This value represents an average density for crustal rocks. Other information (e.g borehole gravity
data or Nettleton’s method) may be used to give a better estimate of the density.

Pierre Bouguer

Summary

Measurement above reference level Add Free Air Subtract Bouguer

correction correction

Measurement below reference level Subtract Free Air Add Bouguer

correction correction

7
 Geophysics 210 September 2008

B4 – Gravity anomalies

B4.1 Gravity anomalies of some simple structures

To understand how geological structures can cause gravity anomalies, let us consider
some simple shapes. Obviously these models are too simple to explain real geology, but
they will illustrate some important concepts.

B4.1.1 Buried sphere

Gravity measurements are made on a surface profile across a buried sphere. The sphere
has an excess mass MS and the centre is at a depth z.

To calculate the pull of gravity, we can use the fact that a sphere has the same
gravitational pull as a point mass located at it’s centre. Simple mathematics can be used
to show that at Point P, the vertical component of g is given by the equation below.

GM S z
gz = 3
(x2 + z 2 ) 2

Suppose:
Radius, a = 50 m Depth, z = 100 m
Density contrast , Δρ = 2000 kg m-3 Excess mass, MS = 109 kg

The variation in gz can be plotted on a profile and map

1
 Geophysics 210 September 2008

Note that:

● gz has it’s maximum value directly above the sphere at x = 0 m (Point A)

GM
● The maximum acceleration is g zmax = 2 S
z
● At Point B, gz has fallen to half the peak value. The distance A-B is called the
half-width of the curve (x½).

Can show that x½ = 0.766 z

This is a very useful equation because it means that if we measure x½ , we can

calculate the depth (z) since
z = 1.3 x½
This allows quantities measured at the surface to be used to find out about
subsurface structure. Note that z is the depth to the centre of the sphere.

● Note that a shallow sphere produces a gravity anomaly that is quite narrow
(short spatial wavelength), while a deeper sphere produces a gravity anomaly
that is wide (long spatial wavelength).

●A gravity survey would measure g zmax and x½ values. Once z is computed, the
g zmax z 2
excess mass can be computed as M S =
G
● Far away from the sphere, gz becomes very small

B4.1.2 Buried cylinder

When gravity measurements are made across a buried cylinder, it can be shown that the
variation in gz will be :

2Gπa 2 zΔρ
gz =
(x2 + z 2 )

This curve is drawn below for a cylinder with

radius, a = 50 m depth of axis, z = 100 m

density contrast , Δρ = 2000 kg m-3 horizontal location, x = 0 m

2
 Geophysics 210 September 2008

Note that :

● the maximum value of gz is located directly above the axis of the cylinder (A)
max 2Gπa 2 Δρ
gz =
z
● gzmax for a cylinder is larger than gzmax for a sphere of the same radius? Why?

● For a cylinder can show that the half-width x½ = z

● Cannot distinguish a buried sphere from a cylinder with just a single profile. Need to
collect gravity on a grid and make a map.

B4.1.3 Uniform layer of rock

A layer of rock has an infinite extent, thickness ∆z and a density ρ. The gravitational
attraction of this slab at the point P is:

g z = 2πGρΔz

Note that gz does not depend on the distance from the layer to the point P. Why?

Consider the two density models shown below.

3
 Geophysics 210 September 2008

What can we say about the gravitational acceleration (gz) of the two models?

This is an example of non-uniqueness in geophysics, and occurs when more than one
Earth model can explain the same set of geophysical data.

B4.1.4 Sedimentary basin

Computations for more complicated shapes cannot be done with analytical formula. A
numerical method must be implemented on computer. When ever using a new piece of
software, always be suspicious about the results …. especially if the software was
expensive! See Geophysics 224 notes for details of some tests of this computer program.

Consider now a simple model that represents a sedimentary basin. Note that the density
of the sedimentary rocks in the basin is lower than that of the surrounding (crystalline)
rock.

The horizontal dashed line represents the acceleration of gravity (gB) due to an infinite
layer, with density contrast of 1000 kg m-3 and thickness 1 km. Use the results of 4.1.3 to
verify the result.

The lower panels show the first and second horizontal derivatives (gradients) of gz across
the basin. Note how the gradients define the edges better than the gravity anomaly.

An example of this in real data is presented later in this section (Alberta basement and
Chicxulub impact crater).

4
 Geophysics 210 September 2008

B4.2 Measuring gravity anomalies

B4.2.1 Absolute gravity measurements

An object is dropped and accelerates at a rate g. After time t it will have fallen a distance
gt 2 2x
x where x = . Absolute value of gravity computed from g = 2
2 t

Absolute gravimeters are generally more expensive than a mass-on-a-spring gravimeter

and can be slower to operate. Typical can measure gravity anomalies down to microgal
level (μgal). Micro-g LAcoste FG-5 gravimeter

5
 Geophysics 210 September 2008

For details of instrumentation http://www.microgsolutions.com/

Application to tectonics in BC http://gsc.nrcan.gc.ca/geodyn/wcda/absgrav_e.php

B4.2.2 Relative gravity measurements

Since we have seen that it is the differences in gravity that is generally more important
than absolute values, we do not need absolute measurements of gravity at every survey
location. Often relative gravity measurements can be made over a survey area, and then
tied to an absolute value by using the relative gravimeter at a location that was
previously surveyed with an absolute gravimeter.

Portable pendulum
L
The period of oscillation (T) a pendulum, length (L) is given by T = 2π
g
Note that as gravity gets stronger, the pendulum swings more quickly. Accuracy
around 0.25 mgal when popular in the 1930’s.

Mass-on-a-spring gravimeters

The mass experiences a force of F = mg and the spring stretches an amount s.

Hooke’s Law states that

F = ks = mg

where k is a measure of the stiffness of the spring (the spring constant). If the
gravimeter is then taken to a location where the acceleration of gravity is stronger by
an amount δg, then the spring will stretch a little bit more, δs.

k(s+δs) = m(g+δg)

6
 Geophysics 210 September 2008

Subtracting these two equations gives

kδs
δg =
m

δs can be very small, so various

engineering features are used to amplify
the movement. Typical accuracy is 0.01
mgal.

To put this in perspective, what change in elevation produces 0.01 mgal?

Two widely used relative gravimeters are:

●LaCoste-Romberg gravimeters (http://www.lacosteromberg.com/relativemeters.htm)

This uses a zero length spring and temperature control to measure to 0.01
milligal. The instrument applies a (known) force to keep the spring the same
length. This avoids non-linear elastic effects (i.e. departures from Hookes
Law).

●Worden gravimeter (http://www.gravityservices.com/worden_gravity_meter.html)

7
 Geophysics 210 September 2008

B4.2.3 Gravity survey procedures

● Collect gravity data on a 2-D grid and repeat measurements at cross-over points
where lines intersect. This will give a good idea of the repeatability of the
measurements.
Gravity data must be corrected for a number of factors. This includes (a) a slow
stretching of the spring in the gravimeter (instrument drift) and (b) variations in
gravity due to the tides. These effects can be removed by setting up a series of base
stations that are visited several times a day. See Geophysics 224 notes for more
details.
●We can obtain absolute gravity measurements from a relative gravimeter by
making measurements at pre-surveyed stations where the absolute values is
already known.

A network of these stations has

been established across Canada
and is called the Canadian
Standardized Gravity Network
(CGSN).
http://www.geod.nrcan.gc.ca/aboutus/
gravity/grvstds_e.php

●Often several survey crews needed for each gravimeter crew. It is vital to know the
elevation of each measurement location. Differential GPS is sometimes good
enough, but leveling may be needed.
● Seafloor gravity surveys use a
gravimeter that is lowered to the
seafloor on a cable.

●Marine and airborne gravity

surveys. Ground is covered much
more quickly than with land-based
methods, but in measurements
made further away from targets.
Need to carefully remove the
effects of acceleration caused by
waves (sea) and turbulence (Air).
LaCoste Romberg AirSea meter
shown on right.

8
 Geophysics 210 September 2008

B4.2.4 Gravity measurements from space

Many types of remote sensing surveys can be carried out with satellites. However,
variations in gravity cannot be made from an orbiting satellite. Why

Elevation of sea-surface

However, the sea surface responds to

the rock structure below. If there is
excess mass (high density) then the
sea water will bunch up, since it is
attracted by gravity. Similarly there
will be a dip in the sea surface if the
density is lower.
Radar is used to precisely measure
the elevation of the sea-surface.
Subsurface density structure can be
inferred from variations in sea level.

Note the patterns associated with plate boundaries such as mid-ocean ridges and
subduction zones. These will be discussed in detail later in this course.
Note that the sea-surface also represents the geoid (discuss in B5)

9
 Geophysics 210 September 2008

GRACE (Gravity Recovery and Climate Experiment)

• The GRACE mission (has given new information

about the gravity field of the Earth. Two satellites
are in a low orbit and the distance between them
is accurately measured. Changes in this
separation are caused by increases and decreases
in gravity.

• A map of gravity anomalies (milliGals) highlight short wavelength features better

than a map of the geoid.
• (http://www.csr.utexas.edu/grace/gravity/gravity_definition.html)

• Applications in oceanography (Tapley et al., 2004).

• Temporal change in gravity was observed after the 2004 Sumatra earthquake
(Han et al., 2006). This was due to mass re-distribution of the plates and was
around 15 μgal.
• More details at http://www.csr.utexas.edu/grace/gallery/gravity/).

GOCE (Gravity field and steady-state Ocean Circulation Explorer)

• European Space Agency mission. (http://www.esa.int/esaLP/LPgoce.html)
• Low orbit satellite, streamlined to cut through uppermost atmosphere
• Launch scheduled September 10, 2008
• Will also measure gravity gradients
• Objectives (a) to determine gravity-field anomalies with an accuracy of
1 mGal (where 1 mGal = 10–5 ms–2) (b) to determine the geoid with an
accuracy of 1-2 cm (c) to achieve the above at a spatial resolution better than
100 km.

10
 Geophysics 210 September 2008

B4.2.5 Sample gravity anomaly calculation

A gravity anomaly is the quantity left over after the effects of latitude and elevation
have been accounted for. This typically requires a set of corrections to be made.

Consider a gravity measurement that was made on campus.

milligals
(a) Measurement on gravimeter
(Difference between campus and the CGSN station at the Airport)
49.012

(b) Previously surveyed value at CGSN station at Airport 981117.890

(c) Value on University of Alberta campus is (b)+(a) 981166.902

(d) Value on campus predicted by GRS67 equation (θ = 53.506944˚) 981369.388

(e) Gravity anomaly is (c)-(d) -202.486

(f) Free air correction for 600 m elevation = 300 x 0.3086 = 185.160

(g) Free air anomaly is (e)+(f) -17.326

(h) Bouguer correction for 600 m elevation = 0.1119 x 600 = 67.140

(i) Bouguer anomaly is (g) –(h) -84.466

Compare this value with the map in section B4.3.4

11
 Geophysics 210 September 2008

B4.3 Examples of gravity anomalies

B4.3.1 Caves and cavities

●Cave location in karst terrain. Taken

from Burger Figure 6-37. The caves
produce a decrease in the gravity
anomaly (note that the Bouguer
anomaly is the gravity measurement
after the Free Air and Bouguer
corrections have been made)

● Sand and clay have a lower density

than limestone. Note that the variable
thickness of sand and clay can mask
the effect of the voids.

● Microgravity were recently used by United Nations weapons inspectors to look for
underground bunkers in Iraq prior to the 2003 invasion.

“Microgravity meters -- also called gravimeters -- measure minute differences in gravitational

pull at one site versus another. Large underground voids, such as tunnels or weapons production
facilities, slightly lower Earth's gravitational pull at the surface right above the voids.
Gravimeters can detect these differences, indicating where such facilities might exist. According
to a source familiar with the inspections, gravimeters operate too slowly to efficiently scan large
areas. However, they work well within a single structure, such as a palace or a bunker, where
single and/or multiple basements are suspected” (Geotimes, 2002)

12
 Geophysics 210 September 2008

B4.3.2 Mineral deposits

Ore bodies are often higher density than the host rock and can produce positive
gravity Bouguer anomalies. Gravity can be used to estimate the excess mass of an
ore deposit, using Gauss’s theorem. While non-uniqueness prevents the spatial
distribution being uniquely determined, the total excess mass can be estimated
reliably.

The example below comes from Voisey’s Bay, Labrador. This massive sulphide
deposit has a pronounced positive gravity anomaly. Note that gravity inversion is
an automated procedure that determines a density model that fits the measured
gravity data. An inversion is a solution of the inverse problem, and non-uniqueness
must be taken into account.

Simple half-width calculation can be used to estimate the depth of the body.

High density sulphides again produce a

positive gravity anomaly. Note that the
shallow depth of the ore body gives a
short wavelength anomaly. The deeper
geological structure produces longer
wavelength anomalies that are
sometimes called the regional trend.

Sometimes the regional trend is

subtracted from the gravity to emphasize
the short wavelength features that are
due to shallow structures. The quantity
remaining is called the residual gravity
anomaly.

13
 Geophysics 210 September 2008

B4.3.3 Impact craters

Gravity anomaly over the Chicxulub

impact crater in Mexico. Note the
circular pattern of high and low values.
These anomalies were discovered in the
1970s but the significance was not
understood until the late 1980s. Now
known that a major comet-asteroid
impact occurred here 65 million years
ago. However, role of the impact in the
K-T extinction is still debated.

Horizontal gravity gradient data at

Chicxulub. Note that the gradient shows
the edge of the crater, which is also
defined by a line of sinkholes (cenotes).
The crater filled with sedimentary rocks
after the impact so there is almost no
surface expression.

Now known that a major comet-asteroid

impact occurred here 65 million years
ago. However, role of the impact in the
K-T extinction is still debated. Crater
has a radius close to 100 km and
impacting object believed to be 10-15
km across. Rebound of the crust formed
an uplift in the centre of the crater.

Mass deficit is 2 10e+12 tonnes

Reference:

O. Campos et al, Gauss’s theorem, mass

deficiency at Chicxulub crater and the
extinction of the dinosaurs, Geophysics,
63, 1585-1594, 1998)

14
 Geophysics 210 September 2008

B4.3.4 Bouguer anomaly map of Canada

● A map of Bouguer anomaly should reflect variations in sub-surface density. Variations

due to latitude and elevation should have been removed. Low (negative) values of
Bouguer anomaly indicate lower density beneath the measurement point. High
(positive) values of Bouguer anomaly indicate higher density beneath the
measurement point.

● on the Canadian Shield gravity anomalies are small (10 to -50 mgals). Variations due to
changes in composition and thickness of crust.

● Negative values coincide with mountain ranges (Canadian Cordillera in west and
Laurentides in east). In these areas the crust is thicker than normal. Since crustal rocks
have a lower density than mantle rocks, this gives the upper 100 km a lower average
density. Thus Bouguer anomaly is weaker. See discussion later in class of isostacy.

● Large positive Bouguer anomalies offshore where crust thins.

15
 Geophysics 210 September 2008

B4.3.4 Bouguer anomaly map of Alberta

● What are the dominant features in the Bouguer anomaly map of Alberta? Sketch
a profile from Fort McMurray to the Rockies and try and account for the main
features observed above.

● White lines denote the boundaries of basement blocks. Note correlation between
gravity gradients and these boundaries.

● Data and figure from : M. Pilkington et al, Can. J. Earth Sciences, 37, 1453-
1471, 2000.

16
 Geophysics 210 September 2008

B5 : Gravitational potential and the geoid

B5.1 Gravitational potential energy

We have considered the Earth’s gravity field in terms of the acceleration that a test mass
would experience (g).

An alternative way to visualize the gravity field is by using the gravitational potential
energy (U). This has some advantages because g is a vector while U is a scalar.

Example 1 : Point mass at z = 0

• Each of the circles shown above represents an equipotential surface. That means
no work is required to move a mass along the surface.
• However, moving away from the mass requires that work is done and converted
into gravitational potential energy.
• If you move towards the mass, then some gravitational potential energy is
converted into kinetic energy (falling).
• The acceleration (pull) of gravity (g) is at right angles to the equipotential surface.
• Rate of change of potential with distance is proportional to gravitational
acceleration (g = −∇U )

1
 Geophysics 210 September 2008

Example 2 : g constant as elevation increases. g = −∇U

Earth surface at z=0

Example 3 : Point mass located just below the Earth’s surface

• Examples 1 and 2 can be combined to represent a mass buried in the Earth

• Note that this has the effect of bending the equipotential surfaces upwards.
• g always orthogonal to equipotential surfaces

Equipotentials and the sea surface

• Water will flow downhill until it has reached the lowest possible level, since it is
not rigid. Thus on a calm day, when ocean currents and weather are ignored, the
sea surface is an equipotential surface. Large mounds of water do not persist in
the open ocean!
• Suppose that there is some excess mass (high density) buried under the seafloor.
This will distort the equipotential surface so that it bulges upwards and results in
a small rise in sea level above the mass excess.

2
 Geophysics 210 September 2008

• Can also think of this as the mass excess pulling water towards itself.

Similarly, if there is a mass deficit (low density) then a slight depression will form in the
sea surface.

Question : Would this be a good place for water skiing?

Hint - How is the pull of gravity related to the sea surface?

5.2 The Geoid

• The geoid is defined as the equipotential surface that coincides with mean sea
level. If the Earth had uniform density and no topography, then the geoid and
International Reference Ellipsoid (IRE) would be identical.
• However, non uniform density and topography result in the geoid being relatively
bumpy.
• In the oceans, the geoid is defined by the sea-level and can be measured with
radar altimetry.
• On land the geoid corresponds to the level that a hypothetical ocean would have.
Location must be computed from gravity measurements (land based and satellite
data). See details of GRACE described below.

• Differences between geoid and

IRE are called geoid
undulations.

An interactive program that calculates the height of the Geoid:

http://www.ngs.noaa.gov/cgi-bin/GEOID_STUFF/geoid99_prompt1.prl

Geoid elevation in Edmonton = ______ m

3
 Geophysics 210 September 2008

The geoid was increasingly well defined by satellite data collected from 1957 onwards.
Long wavelength features are shown below.

• Undulations are 100 m above and

below the IRE. Major hole in the
Indian Ocean and a bulge over
Indonesia and Australia.
• Geoid highs are caused by excess
mass.
• Geoid lows caused by mass
deficit

A more detailed view of the geoid has come from the GRACE project (Gravity Recovery
and Climate Experiment). This has been underway since 2002 and uses microwave
measurements of the distance between two satellites to measure the gravity field with
very high precision. (http://www.csr.utexas.edu/grace/)

Geoid from GRACE data (m)

4
 Geophysics 210 September 2008

5.3 Why is the geoid important?

5.3.1 Information about density structure in the interior of the Earth

• The geoid contains a lot of information about the internal structure of the Earth.
• Figure below shows the geoid superimposed on a topography map. Red: -7.7 m
Magenta = -52.8 m.

• The long spatial wavelengths of geoid anomalies indicate that they originate in
deep density contrasts that are associated with mantle convection.
• Observed geoid variations can be explained if density contrasts from subduction
zones persist to mid-mantle depths (Hager, 1984).

5.3.2 Surveying

• The geoid is always orthogonal to the local acceleration of gravity (g).

• Thus elevations measured with levelling (surveying) will be relative to the geoid
and are expressed as heights above sea level.
• In contrast surveys that use GPS satellites will give elevations relative to the IRE.
• To reconcile the two measurements of elevation, the difference between the IRE
and geoid (N) must be known.

H = height above geoid (orthometric height)

h = height above IRE (ellipsoidal height)

N = geoid height (undulation)

To illustrate this difference, consider a ship sailing across the Indian Ocean. GPS
measurements will indicate that it drops into a 70 m deep hole. Conventional surveying
tells us that it stays at sea level the whole time.

More reading : http://en.wikipedia.org/wiki/Geoid

5
 Geophysics 210 September 2008

References
Hager, B., Subducted slabs and the Geoid: constraints on mantle rheology, J. Geophys.
Research, 89, 6003-6015, 1984.
Han, S.C., C.K. Shum, M. Bevis, C.Ji, C.Y. Kuo, Crustal dilation observed by GRACE
after the 2004 Sumatra-Andaman earthquake, Science, 313, 658-662, 2006.
Richards, M. A. and B. Hager, Geoid Anomalies in a dynamic Earth, J. Geophys.
Research, 89, 5987-6002, 1984.
Tapley, B.D., D.P. Chambers, S. Bettadpur and J.C. Ries, Large scale ocean circulation
from the Grace GGMo1 Geoid, Geophys. Res. Lett., 2004.

MJU 2008

6
 Geophysics 210 September 2008

B6 Isostacy

B6.1 Airy and Pratt hypotheses

Himalayan peaks on the Tibet-Bhutan border

• In the 19th century surveyors used plumblines and theodolites to map India. A
plumb line was used when measuring the elevation angle of a particular star, as
described in B1 when Jean Picard measured the radius of the Earth.

• North-South distances in India were measured with two methods (a) changes in
elevation of stars and (b) direct measurements on the ground using triangulation.

• When a plumbline was used, corrections were made for the anticipated attraction of
the Himalaya

• The results of (a) and (b) gave answers that differed by an angle of 5” (0.0014°).
This discrepancy was shown to be because the Himalayan peaks did not deflect
the plumbline as much as predicted (deflection was 1/3 of the expected value).

B6.1.1 Pratt’s Hypothesis

An explanation of this phenomenon was given by John Pratt who was Archdeacon of
Calcutta. Pratt’s hypothesis of isostacy proposed that topography is produced by
crustal blocks with varying density, that terminate at a uniform depth.

1
 Geophysics 210 September 2008

At the compensation depth, pressure is equal at all points (it behaves as a liquid). Thus at
the compensation depth, the pressure below the mountain (B) must equal the pressure
below the Indian plains (A).

ρ c t = ρ1 ( h + t )

Rearranging gives

ρct
ρ1 = which simplifies to ρ1 = ρc t / (h + t)
(h + t )

What crustal density (ρ1) is needed to explain the 5 km high Tibetan Plateau?
Assume ρc = 2800 kg m-3 and ρm = 3100 kg m-3 and t = 30 km.

ANSWER : ρ1 = 2400 kg m-3. What could cause this type of density change?

B6.1.2 Airy’s hypothesis

George Airy accounted for these observations with a different idea. In Airy’s hypothesis
of isotacy, the mountain range can be thought of as a block of lithosphere (crust) floating
in the asthenosphere (lava). Mountains have roots, while ocean basins have anti-roots.

In his 1855 paper, George Airy who was the Astronomer Royal wrote:

“It appears to me that the state of the earth’s crust lying upon the lava may be compared with
perfect correctness to the state of a raft of timber floating upon water; in which, if we remark one
log whose upper surface floats much higher than the upper surfaces of the others, we are certain
that its lower surface lies deeper in the water than the lower surfaces of the others”

2
 Geophysics 210 September 2008

If the system is stable (no external forces) it is said to be in isostatic equilibrium. At the
compensation depth, the pressure due to material above is constant at all locations
(below this depth the Earth behaves as a liquid).

A plateau of height h is supported by a crustal root of depth r. The normal crustal

thickness is t. In this region, the acceleration of gravity is g.

Pressure at a depth h in a medium of density ρ is given by P = ρgh

Thus equating the pressure at the compensation depth at ‘A’ and ‘B’ we can write

t ρc + r ρm = (h + t + r) ρc which simplifies to r ρ m = h ρ c + r ρc

Note that we have assumed that g has the same value at each location. This may seem to
contradict the last few weeks of classes, but is valid as a first order approximation.

Re-arranging this equation gives the thickness of the crustal root,

r = h ρc / (ρm – ρc)

3
 Geophysics 210 September 2008

Question

How deep a root is needed to support the 5 km high Tibetan plateau?

Assume ρc = 2800 kg m-3 and ρm = 3100 kg m-3

Answer r = 47 km and the crustal thickness = h + t + r = 82 km.

Was George Airy correct about crustal thickness in Tibet? Modern seismic exploration
has shown that the crustal thickness in Southern Tibet is in the range of 75-85 km.

Note : George Airy explained the small deflection of the pendulum as

“It will be remarked that that the disturbance (gravity anomaly) depends on two
actions; the positive attraction produced by the elevated table land; and the
diminution of attraction, or negative attraction, produced by the substitution of
a certain volume of light crust for heavy lava”

Question : What Free Air and Bouguer anomalies would be measured?

To consider this, compute the attraction of mass above the compensation depth at points
‘A’ and ‘B’.

gA = 2πG (ρct + ρmr) and gB = 2πG (h+t+r) ρc

Consider the condition we previously derived for buoyancy?

What does this tell us about gA and gB?

Are the values of gA and gB Bouguer or Free Air anomalies?

4
 Geophysics 210 September 2008

B6.2 Under compensation and over compensation

Are topographic features on the Earth’s surface in isostatic equilibrium? This can be
investigated by calculating the isostatic gravity anomaly as follows:

(1) Topography is used to predict the expected depth of the crustal root needed to
support the topography in isostatic equilibrium.

(2) Using this depth, the predicted Bouguer anomaly ΔgR is computed.

(3)The Isostatic Gravity Anomaly is defined as

ΔgI = ΔgB−ΔgR ,

where ΔgB is the measured Bouguer anomaly and ΔgR is the predicted Bouguer
anomaly of the root.

Three situations can arise:

(a) Complete Compensation: Topography and roots in equilibrium.

ΔgB = ΔgR and ΔgI = 0

5
 Geophysics 210 September 2008

(b) Over-compensation: If surface material is removed (e.g. by erosion) then this

results in a crustal root that is too large. To restore equilibrium, upward motion will
occur.

ΔgB < ΔgR and ΔgI < 0 (negative)

(c) Under-Compensation: In this case, the root will be too small to support the
feature the feature. Tectonic forces may give additional support (e.g. plate flexure,
dynamic topography) or else subsidence will occur.

ΔgB > ΔgR and ΔgI > 0 (positive)

6
 Geophysics 210 September 2008

Bouguer gravity anomaly Isostatic gravity anomaly

Details at : http://gdcinfo.agg.nrcan.gc.ca/products/grids_e.html

B6.3 Gravity anomalies and isostatic equilibrium of plate boundaries

B6.3.1 Are mountain belts in isostatic equilibrium?

Canadian Cordillera

• Figure above is taken from Fluch et al., (2003) and shows gravity anomalies
across the Southern Canadian Cordillera on profile AB.
• Note the negative Bouguer anomaly.
• The isostatic gravity anomaly is quite close to zero, indicating that the Canadian
Cordillera is close to isostatic equilibrium.
• This may be the result of a combination of Pratt and Airy hypotheses.
• The crust is thick under the Continental Ranges (Airy) but thinner and with a hot
upper mantle under the Omineca belt (Pratt).

7
 Geophysics 210 September 2008

Clowes et al., CJES, (1995).

Tibet and the Himalaya

• Fowler Figure 10.19 on page 540 shows the Bouguer anomaly measured across
the Himalaya from Cattin et al., (2001). It also shows the predicted Bouguer
anomaly assuming complete isostatic compensation.

• In Ganges Basin ΔgB < ΔgR which represents over-compensation. This occurs
because the Indian Plate is deflected (pulled) downwards by loading in Tibet.
This makes the crustal root thicker than needed to support the observed
topography.

• The Himalaya and southern part of the Tibetan Plateau is under-compensated

with ΔgB > ΔgR which can be thought of not having enough crustal root to
support the observed topography. Thus the Himalaya are partially supported by
flexure of the Indian Plate.

• To the North of the Indus-Tsangbo suture, complete compensation occurs and

the observed topography can be supported by buoyancy forces. George Airy was
correct!

8
 Geophysics 210 September 2008

B6.3.2. Mid-ocean ridges

Figure 9.11 from Fowler (2005) shows gravity anomalies across the mid-Atlantic Ridge
at 45° N. Mid-ocean ridges are large submarine mountain ranges located where plates move
apart to create ocean basins.

• Free air anomaly small, but not

exactly zero. This shows that the
ridge is not quite in isostatic
equilibrium. Topography of ridge
supported by low density upper
mantle (hot and partially molten).

• Figure 9.11 shows a range of

models that all fit the observed
data, giving an illustration of
non-uniqueness.

B6.3.3 Subduction zones

• Figure 9.59 from Fowler (2005) shows the Free Air gravity anomaly across the
Chile Trench and Andes at 23° S, taken from Grow and Bowin, (1975).
• Like many subduction zones, this model shows a characteristic pair of low-high
gravity anomalies.
• Note that the gravity modelling includes the phase transition of basalt to eclogite
in the subducting slab. This corresponds to an increase in density of 400 kg m-3
within the slab.
• Negative gravity anomaly due to deep trench that is filled with low density water
and sediments.
• Positive gravity anomaly on ocean side of the volcanic arc.
• Flexure in the subducted plate can cause under compensation on the overriding
plate. This is because flexural forces give partial support to the topography. In this
situation, Buoyancy forces do not completely support the observed topography.

B6.4 Isostatic rebound

• We can determine if topographic features are in isostatic equilibrium by studying

gravity anomalies.
• An additional perspective on isostacy can be obtained by looking at time
variations that result from rapid changes in the size of ice sheets or large lakes.
• These natural experiments allow the viscosity of the mantle to be determined.
• The basic concept shown in Figure 5.18 of Fowler (2005).

9
 Geophysics 210 September 2008

Canadian Shield
• Strong evidence for post-glacial isostatic rebound comes from the raised
beaches on the coast of Hudson’s Bay.

http://www.earth.ox.ac.uk/~tony/watts/LANDSCAPE/LANDSCAPE.HTM

http://gsc.nrcan.gc.ca/geodyn/gchange_e.php

• Surveying has revealed a pattern of uplift beneath the former location of the
Laurentide Ice Sheet with a maximum rate of around 1 cm per year.
• This is surrounded by a ring of subsidence that is caused by flexure of the
lithosphere. Note that horizontal motions also occur.

www.uwgb.edu/dutchs/STRUCTGE/EarthMvts.HTM http://www.homepage.montana.edu/~geol445/hyperglac/isostasy1/

10
 Geophysics 210 September 2008

• Total uplift since the ice sheet melted is in excess of 100 m in the centre of
Hudson’s Bay.
• Modern uplift rates are much slower than immediately after the ice sheets melted.
• These values have been confirmed by gravity measurements made by the GRACE
satellite (Tamisiea et al., 2007). Gravity changes have a peak change of 1 μGal per
year from 2002-2007.
• Geoid motion was also determined by GRACE and is ~1 mm per year.
• Mantle convection can change gravity over timescales 10 times longer than the
timescales for isostatic rebound. The shorter time scale of post-glacial rebound
allows the phenomena of mantle convection and rebound to be separated.

Fennoscandia

• Very similar observations of isostatic rebound can be seen in Fennoscandia

• Maximum uplifts of 1 cm per year are observed in GPS data.
• Islands are emerging from the Baltic Sea and the border between Sweden and
Finland have been revised on several occasions.
• Kvarken Archipelago in Finland has unique geography

• Note that horizontal motions also occur, and also the yearly variation that is due to
climatic effects on the ground and atmosphere.

11
 Geophysics 210 September 2008

http://www.oso.chalmers.se/~hgs/docent/docans.html

Figure 5.20 in Fowler (2005) shows that an asthenosphere viscosity of 1021 Pa s is

consistent with the uplift history in the Baltic.

For reference, the following values are for the viscosity of some other liquids:

Viscosity (Pa s)
Water 10-2
Honey 2-10
Peanut butter 250
Pitch 108

B6.5 References
Airy, G.B., On the computations of the effect of the attraction of the mountain masses as
disturbing the apparent astronomical latitude of stations in geodetic surveys, Phil. Trans. R.
Soc. London, 145, 101-104, 1855.
Chen, W. P. and S. Ozalaybey, Correlation between seismic anisotropy and Bouguer gravity
anomalies in Tibet and its implications for lithospheric structure, GJI, 135, 93-101, 1998.
Clowes, C. Zelt, J. Amor and R. M. Ellis, Lithospheric structure in the Southern Canadian
Cordillera from a network of seismic refraction lines, CJES, 32, 1485-1513, 1995.
Flück, P., R. D. Hyndman, and C. Lowe, Effective elastic thickness Te of the lithosphere in
western Canada, J. Geophys. Res., 108(B9), 2430, doi:10.1029/2002JB002201, 2003
Grow, J. A. and Bowin, C. O., 1975, Evidence for High-Density Crust and Mantle Beneath the
Chile Trench Due to the Descending Lithosphere, J. Geophys. Res. 80, 1449–1458.
Pratt, J.H., On the attraction of the Himalaya Mountains, and of the elevated regions beyond
them, upon the plumb line in India, , Phil. Trans. R. Soc. London, 145, 53-100, 1855.
Tamisiea, M., J. Mitrivica, J.L. Davis, GRACE gravity data constrain ancient ice geometries and
continental dynamics over Laurentia, Science, 316, 881-883, 2007.

MJU 2008

12
 Geophysics 210 September 2008

B7: Tides, rotation and orbit

B7.1 Tides
• Tides on the Earth arise from the
gravitational influence of the
Moon and Sun. This results in
motion of seawater and also in
changes in gravity that can
influence gravity surveys.
• Effect of sun and moon are
similar, so they will initially be
considered separately.
• Tidal range can exceed 18 m in
coastal locations such as the Bay
of Fundy in Nova Scotia.

B7.1.1 Lunar tides in the ocean

• Earth and moon revolve around their common centre of mass which is located a
distance c from the centre of the Earth.

• If ME = 5.97 10 24 kg, MM = 7.36 10 22 kg and R = 384400000 m

Show that c = 4678 km

• To simplify the computation, must choose an appropriate reference frame.

Answer should be the same in any reference frame, so we can make life simple by
choosing a reference frame that does not rotate and which is centered at the
centre-of-mass of the Earth-Moon system.

• In this frame of reference, the Earth translates and does not rotate. Thus all points
on Earth move in a small circle with the same radius. See animation in Figure 7
at http://www.vialattea.net/maree/eng/index.htm for a clear demonstration.

• Centrifugal motion is thus the same for all points on surface of Earth.

1
 Geophysics 210 September 2008

• However, pull of gravity varies across the surface of the Earth. This can be
illustrated by comparing pull of gravity at ‘A’ and ‘B’ with the value at the centre
of the Earth. At point ‘A’ can show that the net acceleration due to the moon is :
GM M GM M
AA = − = 0.113 milligals
(R − r ) 2
R2
GM M GM M
Similarly at point B, AB = − = -0.107 milligals
(R + r ) 2
R2

• Note the sign difference. This corresponds to a 6% reduction in the acceleration

of gravity due to the moon, from A to B.

• When computed at all points on the surface of the Earth, the following pattern of
tidal forces is computed.

This deforms the sea surface to form two

tidal bulges. As the Earth rotates within
the liquid, two high and two low tides
per day are observed.

This is the principal lunar semidiurnal

constituent of the tide (M2). Period is 12
hours 24 minutes, owing to relative
rotation of moon around the Earth and
rotation of Earth.

• Note that tidal force is due to the spatial gradient (change with distance) of the
gravitational pull of the moon.

• Theoretical amplitude of M2 tide should be 54 cm of vertical motion and is

referred to as a body tide. However seawater moves due to changes in g and
amplifies this effect (1 m).

• The Earth rotates faster than the moon orbits the Earth.

• Friction occurs between seawater and seafloor as the Earth moves beneath the
tidal bulges. This drags the bulges forward and tides lead the moon by 10°.

2
 Geophysics 210 September 2008

• Energy loss slows the Earths rotation by 0.0023 seconds per century (i.e. days get
longer).

• Another effect is that the moon is

accelerated and moves into a
larger orbit. This lunar recession
increases Earth-moon distance by
3.8 cm per year. This has been
measured by laser ranging
experiments from Apollo and
Lunakhod missions to the moon.

B7.1.2 Solar tides in the ocean

• Similar physics to that considered above for the moon

• Using mass of sun, Ms = 1.988 1030 kg and an average Earth-Sun distance

Rs = 1.5 1011 m can show that for the S2 tide:
GM S GM S
AA = − = 0.05012 milligals
( RS − r ) 2
Rs2
and
GM S GM S
AB = − = -0.05011 milligals
( RS + r ) 2
Rs2

• This corresponds to a 0.017 % reduction in the gravitational pull of the sun from
point A to point B.

• Note that the gravitational attraction of the sun is 180 times stronger than that of
the moon. However, the attraction of the moon varies more rapidly with distance
(higher spatial gradient), giving a larger tidal force.

Relative magnitude of M2 and S2 tides

The tides are due to the gradient of gravity forces (i.e. how much they vary per
metre). Illustrated below where the acceleration of gravity (g) produced by the moon
and sun is shown as a function of distance from the Earth Note that:

3
 Geophysics 210 September 2008

(1) g sun > g moon at the Earth.

dg sun dg moon
(2) < . Look carefully and you can see that the slope of gmoon is twice
dr dr
that of gsun at the location of the Earth.

4
 Geophysics 210 September 2008

B7.1.3 Actual tides in the ocean

• Tidal bulges due to the moon and sun can be added, according to the relative
orientation of the moon and sun.

• Neap tides occur when moon and sun are at right angles (first quarter and last
quarter)

• Spring tides occur when sun and moon are collinear (new moon and full moon).

In addition, the simple tidal model listed above is complicated by a number of factors:

5
 Geophysics 210 September 2008

• axis of rotation is inclined at 23°

• variable ocean depth, coastlines
• Earth-moon and Earth-Sun distances vary, since both orbits are elliptical
• Earth is an ellipsoid, not a sphere
• Factors such as rotation of Earth, variable water depth and coastlines can amplify
this even more.

• Coastlines result in a complicated pattern of tides and resonance effects (e.g.

Rossby waves) can develop.

• Note small tidal amplitude in Mediterranean Sea and Gulf of Mexico. Why?

• Tidal phase (high and low) occur at the same time on the cotidal lines. These are
shown in white on the figure above. M2 component dominates other contributions
with 12 and 24 hour period.

• Amphidromic points occur where cotidal lines meet and there is no vertical
motion.

• About 1 terrawatt (1012 W) of frictional energy is dissipated in the ocean. See

Egbert and Ray (2000) for details.

6
 Geophysics 210 September 2008

B7.1.4 Earth Tides

• The solid Earth deforms in response to tidal forces in exactly the same way as the
ocean. However, the Earth does not flow horizontally in the same way as the
ocean, so tides are smaller (50-100 cm vertical displacement).

• Earth tides any point on Earth can be computed at

http://www.taygeta.com/etides/index.html

• Relative motions of Earth and ocean can reduce absolute amplitude of ocean tides

• Figure below shows the Earth tide sequence for (left) Edmonton and (right) the
North Pole in September 2007.

• Why is variation smaller at the North Pole? Is this what you would expect to see?

7
 Geophysics 210 September 2008

B7.1.5 Consequences of Earth tides

Triggering earthquakes

• Cochran et al., (2004) showed a

correlation between occurrence
of shallow thrust earthquakes and
strongest tides.

• Stress due to tide is small,

however if an earthquake is
about to occur, this can be the
straw that breaks the camels
back.

• Factor of three variation is

observed in the number of
earthquakes was associated with
tidal stress.

Eruptions of volcanoes and geysers

Castle Geyser, Yellowstone NP

At high Earth tides, the rock is dilated and eruptions are more frequent. At low tides the
dilation is less and fractures are closed.

8
 Geophysics 210 September 2008

Need to correct gravity survey data for tidal effects

• Comparison of predicted and observed tides in Montreal. Telford Figure 2.6.

• Some modern gravimeters make these corrections automatically.

References

General reference about tides: http://en.wikipedia.org/wiki/Tide

http://www.vialattea.net/maree/eng/index.htm
Good explanation about origin of tides and centrifugal forces. Considers the physics of
tides in a number of reference frames and shows that the same answer can be obtained,
provided that fictitious forces are correctly handled.

Cochran, E.S., J.E. Vidale, S. Tanaka, Earth Tides can trigger shallow thrust fault
earthquakes, Science, 306, 1164-1166, 2004.

Egbert G.D., and R. D. Ray, Significant dissipation of tidal energy in the deep ocean
inferred from satellite altimeter data, Nature, 405, 775-778, 2000.

9
 Geophysics 210 September 2008

B7.2 Changes in length of day (LOD)

The Earth rotates because of angular momentum that was acquired by the solar system as
it condensed from a rotating gas cloud into a protoplanetary disk and then individual
planets and moons. This momentum has been conserved over the life of the solar system.

B7.2.1 Present situation

Variation with latitude and Julian Day

http://en.wikipedia.org/wiki/Day_length

Geographical extent of day and night as seasons vary:

http://www.fourmilab.ch/cgi-bin/uncgi/Earth/action?opt=-p&img=learth.evif

B7.2.2 Variations over geological timescales

Archean Early Earth may have rotated every 8 hours

PreCambrian Tidal rhythmites can record number of days per year and
per month. Rocks from Australia consistent with LOD =
21.9 hours at 620 Ma (Williams, 1997 and 2000). Record
spanned 60 years.

Phanerozoic

Figure above shows the predicted number of days per year, based on extrapolation
of present rates. Wells (1963) showed that the number of daily growth rings in
coral are consistent with this rate of change:

Middle-Devonian New York 400 Ma 400 days per year

Pennsylvanian Texas 385 Ma 385 days per year

Rate of change in LOD = 2 x 10-5 s per year = 2 ms per century

10
 Geophysics 210 September 2008

B7.2.3 variations over historical record

• Historical increase in the length of day (LOD) is around 2.4 ms per

century.
• This rate of change seems small on a timescales of millennia, but the
actual effect is significant because the change is integrated as shown
below.
• Stephenson (2003) shows that this factor is needed to account for the
astronomical records made by many ancient civilizations (Babylonians,
Greeks, Chinese, Arabs).
• The figure below is from Morrison and Stephenson (2004) and shows
constraints on the cumulative time shift (ΔT) based on eclipses. This
figure just uses information on whether an eclipse could have been
observed or not. No record of timing was used.

• The next figure uses records that have times associated with the
observation of an eclipse.

11
 Geophysics 210 September 2008

• Assume LOD changes at rate, r = 2 x 10-5 s per year. Consider a year at a

time T in the past. The time delay in that year = 365 T r seconds.
• Thus total time delay (seconds) from present day back to a time T is

t =T t =T
⎡t 2 ⎤
ΔT = ∫ rt 365dt = 365r ⎢ ⎥ = 182rT 2
t =o ⎣ 2 ⎦ t =0

• Using a value of T = 2000 years gives a delay of 14000 seconds ( 4 hours),

in agreement with observations above.

Explanation for change in length of day (long term trends)

• Angular momentum of the Earth-Moon system is conserved.

• The tidal bulges are not collinear with the Moon (B7.1) and this results in
a torque which slows down the Earth's rotation, i.e. the days are getting
longer.
• While kinetic energy (orbital motion) is dissipated and converted to heat,
the angular momentum of the Earth-Moon system must be conserved. The
loss of angular momentum from the Earth must be balanced by an increase
in momentum for the moon. This results in the moon moving into a higher
orbit at a rate of 3.7 cm per year. This is measured by satellite laser
ranging and lunar laser ranging that measure the Earth-Moon distance to
a precision of 1 part in 10 billion.
• This rate could not have persisted for the entire age of the Earth. Present
day configuration of oceans and continents may be causing a high rate of
tidal dissipation at present.
• The recent rate of change in LOD may be due to mass redistribution
caused by isostatic rebound and sea level changes (see Peltier, 1988).

B7.2.4 Recent records

• Observations with telescopes (since 1620 AD), from satellite laser-ranging

(SLR), lunar laser ranging (LLR) and radio astronomy (very long baseline
interferometry - VLBI) show that the length of day (LOD) changes at a
rate of 1.4 ms/century.
• Significant fluctuations observed from year-to-year and seasonally.
Figure from Hopfner (1998) below.

12
 Geophysics 210 September 2008

Explanation for short term changes in LOD

• Climatic effects - exchanges of angular momentum between the atmosphere and

the solid earth. Average wind speed up - rotation rate down and vice versa.
• Effect of fluid flow in core? Angular momentum exchanged between core and
mantle. Depends on coupling that is caused by topographic irregularities at the
core-mantle boundary? Electromagnetic coupling?
• Change in LOD after the 2004 Sumatra earthquake (Chao and Gross, 2005)

13
 Geophysics 210 September 2008

B7.3 Orbital changes

B7.3.1 Basic orbital parameters http://www.phy.hk/wiki/englishhtm/Motion.htm

B7.3.2 Temporal changes in orbital parameters

• The Earth’s orbit changes over time, mainly because of gravitational interactions
between the Earth and the other planets.
• These orbital changes cause changes in the amount of solar energy received by
the Earth.
• This insolation at high latitudes is believed to control the growth and retreat of
polar ice sheets.

(a) Eccentricity
• This defines the difference between the orbit and a circle.
• A circle has an eccentricity of zero.

Present eccentricity = 0.017

413 kyr cycle 0.005 to 0.058
100 kyr cycle 0.03 to +0.02 (combination of 95 and 125 kyr cycles)

Time variations in eccentricity are due to the gravitational effect of other planets.

(b) Obliquity

• This is angle between the Earth’s rotation axis

and the ecliptic. Present value = 23.44°
• 41 kyr cycle causes obliquity to vary between
22.1° and 24.5°
• Lower values of obliquity are believed to favour
the start of ice ages, as insolation is reduced at
high latitudes in summer. This results in less of
the winter’s snow melting.

14
 Geophysics 210 September 2008

(c) Precession:

• This is the wobble of the axis of rotation relative to the fixed stars.
• Gyroscopic action that is a consequence of the Earth not being a uniform sphere.
• Rotation axis moves in a circle with a period of 21-26 kyr.
• This variation only has an effect if the orbit is elliptical.
• This can produce a temperature difference between the North and South
Hemispheres.

• At present, the Earth is closest to the Sun (perihelion) during the southern
summer. This makes seasonal variations more extreme in the Southern
Hemisphere (6.8% difference in summer radiation).

Milankovitch cycles
• In1920 Milutin Milankovitch computed the combined climatic effects of these
variations and suggested that these changes could have triggered the ice ages.
• However at the time, no detailed datasets were available to prove, or disprove,
this hypothesis.
• http://en.wikipedia.org/wiki/Milankovitch_cycles

15
 Geophysics 210 September 2008

B7.3.3 Evidence in support of Milankovitch’s theory

Climate record in deep ocean sediments

• Deep ocean cores provide a reliable record of past climate through changes in the
concentrations of oxygen isotopes (δ18O). When cooler, ocean water becomes
enriched in the heavier isotope and the value of δ18O will be higher.
• Hays et al.,(1976) showed that these variations could be correlated with variations
in orbital obliquity and precession.
• Recent studies from the Ocean Drilling Program also show a more complicated
picture. Rial (1999) studied δ18O variations of deep sea cores from ODP Site 806
covering the last 2.1 Ma. (B) is a moving spectrogram in which the periodicity in
the oxygen signal is measured over limited time periods of 350 ka. (C) is called
the power spectrum for the entire set of oxygen records and it shows the different
Milankovich cycles but with a resonance phenomena.

Climate records in loess deposits

• Loess is a sedimentary rock deposited by the wind. It can be deposited over long-
periods of time and can record climate variations.

• Corroborating evidence for periods in the geological record are shown on next
page.

• Left column shows the loess-soil pattern. Next column is the geomagnetic
polarity. GSR = grain size distribution. Right column is the corresponding power
spectrum which shows peaks near 100,000 years, 41,000 years, and 23,000 years.

• These periods all correspond to Milankovitch time periods. Plots courtesy of

Professor N. Rutter in EAS.

16
 Geophysics 210 September 2008

Climate records in ice cores

Climate record from the Vostok ice core in Antarctica

(http://en.wikipedia.org/wiki/Vostok_Station)

17
 Geophysics 210 September 2008

Does orbital inclination have an effect?

The plane of Earth’s orbit shifts with a 70,000 year period. This factor was not
considered by Milankovitch. If defined relative to Jupiter, the period is 100, 000 years.
Variations in dust content have been suggested as a mechanism to produce insolation
variations on this timescale. Note that in the last million years the ice ages are dominated
by a periodicity of 100,000 years.

However Kawamura et al., (2007) validated Milankovitch’s hypothesis and showed that
local insolation variations can explain the 100,000 year cycle.

B7.3.4 How did the Egyptians align the pyramids?

• Spence (2000) shows that

precession could have caused a
systematic variation in the
orientation of the pyramids.

• This gives the most precise dates

available for the construction of
each pyramid.

18
 Geophysics 210 September 2008

B7.4 Wobbles
The Earth also undergoes a series of orbital wobbles. These are called free nutations as
they do not require an external driving force.

B7.4.1 Chandler wobble

Animation at “The Hutton Commentaries”

http://huttoncommentaries.com/subs/PSResearch/UandM_PS2001/Undrstnd_Montr_PS2001.htm

• Pole moves by 3-6 m over a

period of 14 months. This was
discovered by Seth Chandler in
1890’s and is due to precession
of an ellipsoidal Earth.

• Motion is damped so it must be

continually excited. Now
believed that the combined
effect of atmosphere and ocean
keeps the Chandler wobble going
(Gross, 2000).

Figure from Gross (2000) is shown below. This shows that a good correlation between
observed wobble and combination of forces from atmosphere and ocean (2 black
curves). Correlation with atmosphere (red curve) is much lower.

Gross (2000)

19
 Geophysics 210 September 2008

B7.4.2 Annual wobble

Gyration of 3 m due to imbalance of density caused by weather system over Siberia.

Destructive interactions occur every few years between Chandler wobble and annual
wobble. With the high precision of modern GPS data, these interactions can be tracked.
e.g. winter 2005-2006.

Lambert et al., 2006 show that:

“centimetre level polar motion

displacements during the 2005–2006
winter season are almost fully explained
by major pressure events on the
continents and on the ocean, especially a
depression over Northern Europe in-
phase with similar events over North
America.”

B7.4.3 Markowitz wobble

Another source of possible polar wandering with 24 year period. This may be due to
coupling of the inner core and mantle. See Battesby (2006) for review of research by Dr.
Mathieu Dumberry (new faculty member in Geophysics at the University of Alberta).

20
 Geophysics 210 September 2008

B7.5 References
Battesby, S., Watch that wobble, New Scientist, 18 February 2006.
Chao, B. F., and R. S. Gross (2005), Did the 26 December 2004 Sumatra, Indonesia, Earthquake
Disrupt the Earth's Rotation as the Mass Media Have Said?, Eos Trans. AGU, 86(1), 1.
Egbert, G.D., R.D. Ray, Significant dissipation of tidal energy in the deep ocean inferred from
satellite altimeter data, Nature, 405, 775-778, 2000.
Gross, R.S., The excitation of the Chandler wobble, Geophysical Research Letters, 27, 2329-
2332, 2000.
Hays, J.D., J. Imbrie, N.J. Shackleton, Variations in the Earth’s Orbit: Pacemaker of the Ice Ages,
Science, 194, 1121-1132, 1976.
Höpfner, J., Seasonal variations in length of day and atmospheric angular momentum
Geophysical Journal International, 135, 407–437, 1998.
Kawamura et al, Northern Hemisphere forcing of climatic cycles in Antarctica over the past
360,000 years, Nature, 448, 912-917, 2007.
Stephenson, F.R., Historical eclipse’s and Earth’s rotation, Geophysics and Astronomy, 2003
Lambert, S.B., C. Bizouard, V. Dehant, Rapid variations in polar motion during the 2005-2006
winter season, Geophysical Research Letters, 33, L13303, doi:10.1029/2006GL026422, 2006
Kawamura et al, Northern Hemisphere forcing of climatic cycles in Antarctica over the past
360,000 years, Nature, 448, 912-917, 2007.
Morrison, L.V., and F.R. Stephenson, Historical values of the Earth’s clock error ΔT and the
calculation of eclipses, Journal of Historical Astronomy, XXXV, 327-336, 2004.
Peltier, W.R., Global Sea Level and Earth Rotation, Science, 240, 895-901, 1988.
Rial, J., Pacemaking the Ice Ages by Frequency Modulation of Earth's Orbital Eccentricity,
Science, 285, 564-568, 23, 1999.
Spence, K., Ancient Egyptian chronology and the astronomical orientation of the pyramids,
Nature, 408, 320-324, 2000.
Wells, J.W., Coral growth and geochronometry, Nature, 197, 948-950, 1963.
Williams,G.E., Precambrian Length of Day and the Validity of Tidal rhythmite paleotidal values,
Geophysical Research Letters 24, 421-424, 1997.
Williams, G.E., Geological constraints on the Precambrian history of Earth's rotation and the
Moon's orbit, Reviews of Geophysics 38, 37-60, 2000.

Satellite Laser Ranging (SLR) http://ilrs.gsfc.nasa.gov/science_analysis/tides.html

Lunar Laser Ranging (LLR) http://sunearth.gsfc.nasa.gov/eclipse/SEhelp/ApolloLaser.html

VLBI http://cddisa.gsfc.nasa.gov/vlbi_summary.html
http://ivs.nict.go.jp/mirror/about/vlbi/whatis.html

Earth Orientation http://maia.usno.navy.mil/eop.html

Recession of the Moon http://www.talkorigins.org/faqs/moonrec.html

Earth Rotation http://en.wikipedia.org/wiki/Earth's_rotation

Tidal Acceleration http://en.wikipedia.org/wiki/Tidal_acceleration
Delta_T http://en.wikipedia.org/wiki/Delta_T

21
 Geophysics 210 September 2008

B8 : Radial density structure of Earth

• Density, ρ(r), varies with radial distance from the centre of the Earth
• Note that where composition changes, there is a jump in density (e.g mantle to outer core
at depth of 300 km)
• In regions where composition stays constant, increasing pressure causes an increase in
density through compaction.
• This model is the PREM (Preferred Reference Earth Model) that was first defined by
Dziewonski and Anderson (1981).

However, this density model cannot be derived from gravity measurements on the Earth’s surface
because of Gauss’s theorem. This states that

∫ g.dS = 4πGM
S
where M = mass enclosed by the surface ‘S’.

In the case of a radially symmetric Earth the surface ‘S’ is the surface of the Earth and g has a
constant value. Thus

4πRE2 g = 4πGM
 Geophysics 210 September 2008

where RE is the radius of the Earth. Note that M is mass enclosed by ‘S’ and the location is not
important. Thus all mass distributions give same value of g at the surface.

• Could map the variation of density with depth if we could drill to the centre of the Earth
and measure g. Figure above shows g (r) for the PREM density model and an Earth with
uniform density (straight line).

• Lowest value of g is at the centre of the Earth. Why?

• The variation of density with depth, ρ(r) was determined from seismic studies,
constrained by knowledge of the total mass of the Earth.

Dziewonski, A.M, D.L. Anderson, Preliminary reference Earth model, PEPI, 25, 297-356, 1981.