Potential Field Methods

+ natural source methods

+ non-invasive
+ inexpensive
+ fast
+ easy data collection, reduction, but...
- non-straightforward interpretation
- low resolution
- ambiguous
- not always applicable

Method

Advantages

Disadvantages

Cost Ratio

Magnetics

Very fast, very cheap

Poor resolution, not always

applicable

Gravity

Fast, cheap

Poor resolution

10

Seismic

Fine detail, good correlation to

geology

$$$

100

From a 1999 Edcon brochure advertising their aerogravity/magnetic surveys:

"The cost of conducting an aerogravity/magnetic survey over a 5,000 square kilometer
concession in South America is in the order of $200,000 to $300,000. The cost of a 3-D
seismic survey over only 250 square kilometers can be ten times that amount."
Pat Millegan, Marathon Oil, on use of G&M in industry:
Pat stresses the importance of diversifying your skills: " seismic does NOT answer all
the questions, all the time...there are MANY seismic failures (e.g., one current Marathon
project). The main reason G&M does not see more use is true "ignorance". My job is 10100 times harder when my "clients" (the exploration groups...I'm in a service group)
know nothing about G&M. Please stress geophysical integration to your students. It is
the smart way to explore, but you don't just throw G&M at everything...don't bother if the
geology isn't conducive to geophysical results."
1972 Costs of Acquisition and Processing of Geophysical Data (Telford et
al.)
x $106

802

89.7

17

1.9

0.7

Petroleum Exploration
seismic
surface grav/mag
airborne mag

Mineral Exploration
airborne mag

19

2.1

ground mag

12

1.5

Other

34

3.8

Total

894

100

Gravity and Magnetics in a Nutshell

Gravity is useful wherever the formations of interest have densities that are
appreciably different from those of surrounding formations.
Some examples:

mapping sedimentary basins, where sedimentary rocks consistently have

lower density than basement rocks
salt bodies: low density of salt
groundwater studies (e.g., Cayman Islands)
burial chambers in pyramids

Magnetics is useful whenever object of investigation has a contrast in magnetic

susceptibility or remanence
Some examples:

mapping structure on basement

mapping sedimentary basins
direct location of ores containing magnetite

History of Gravity Method

Man has always recognized its force: fear of falling; up & down
Galileo, 1590: pendulum period; force on body proportional to weight;
acceleration of g independent of mass; gal = 1 cm/s2
After sun recognized as center of universe, Tycho Brahe (1546-1601)
made extensive measurements of the "peculiar motion" of planets
Johannes Kepler (1571-1630): Kepler's Laws history

1. The planets move in elliptical orbits with the sun at one focus

where a (distance CA below) and b (distance CB below) are the major and
minor semiaxes. The eccentricity, e, is given by c/a, where c is the
distance from the center of the ellipse to one of the foci, and x and y
represent coordinates of points on orbit. (examples: Earth = 0.01673;
Mercury = 0.2056; Pluto = 0.250)

1. A line drawn from the sun to a planet will sweep out equal areas in equal
times (conservation of angular momentum)

2. The square of a planet's period of revolution is proportional to the cube of

the length of the major semiaxis of the orbital ellipse (conservation of
kinetic and potential energy)

Newton, 1687, Philosophiae Naturalis Principia Mathematica: force of

gravity is a property of all matter, Earth included
Jean Richer, 1672: pendulum clock, accurate in Paris, lost a few minutes
per day in Cayenne, French Guiana

Seen as tool to measure variation in geopotential. Newton correctly interpreted as due to

oblateness. French believed otherwise; French Academy of Sciences sent two expeditions, one
to high latitudes of Sweden, other to equatorial Ecuador (included Pierre Bouguer) to compare
length of degree of arc at both sites.

Pierre Simon, Marquis de Laplace: gravity obeys simple differential eq.

(early to mid 1800s)
Lord Cavendish, 1798, determined G, hence mass of Earth (estimate of G was
6.754x10-11)
Torque required to twist quartz fiber:
Torque provided by gravity:
Set equal and solve for G:

(current value 6.6720x10-11 MKS)

Cavendish experiment leads to mass, bulk density of Earth:

When mass causing acceleration, M, is Earth, we use g to represent

acceleration

We know R = 6371 km (how?), g = 9.8 m/s2, G = 6.67x10-11 MKS (what

are the units?), so M = 6.0x1024 kg
(Bold numbers: memorize!)
Bulk density:

But how was the Newton defined? Improvement in accuracy of G (and

hence mass of Earth) over time:

Gravity as Geophysical Tool

Kater, 1818, reversible pendulum: absolute g

Earliest efforts to locate oil-bearing structures involved gravity: just before
1900, Baron Roland von Etvs, Hungary
o torsion, or Etvs balance
o measures distortion in g field from buried bodies
o slow, cumbersome to operate
1915/16 torsion balance survey at 1-well oil field at Egbell, now in
Czechoslovakia; highly successful
1917 Schweider: salt dome in Germany
1922 Shell: Horgada field in Egypt
1922 Spindletop field in Texas - salt structure
Vening Meinesz, 1928, shipborne pendulum
1930s - Gulf Research & Development, 1st gravimeter (direct readings of
g differences; oil boom, LaCoste-Romberg, Worden meters patented

Potential Fields
Fields
A field is a set of functions of space and time.
We are concerned with 2 kinds of fields:
1. Material fields describe some property at a point of the material and at a
given time (intensive quantity)
Examples: density, porosity, magnetic susceptibility, temperature; not a
material property: mass, heat; these are extensive quantities (depend on
extent of material)
2. Force fields describe forces that act at each point of space at a given
time
Examples: gravity, magnetic field, electrostatic field
Fields can be scalar or vector or tensor
A vector field can be described in terms of field lines (or lines of flow, or lines of
force or flux lines). These are lines that are tangent at every point to the vector
field.

Potential Theory

Concept of potential
Example: Consider map of ski area: put arrows everywhere giving magnitude
and direction of slope; It is easier just to give elevation at each point!

In 1-D
In 3-D:

2-D example of relationship between scalar potential and vector field

[Note: is the "del" operator or gradient operator; it is always a vector
quantity; sometimes it is written with an arrow over it, or boldfaced, to
indicate that it is a vector operator] Thus we see that a scalar field (elevation)
can give rise to a vector field (slope)
Another example: temperature field (scalar), heat flow field (vector), where

Conservative Fields
For force fields (vector fields) it can be shown that if the force field is
conservative, it may be (and must be) represented as the gradient of a scalar
field.
1. All force fields derived from scalar field are conservative
2. All conservative fields can be derived from scalar

Let's show that a force field derived from scalar is

conservative:
Conservative:

Stokes Theorem: [Kaplan, Advanced Calculus, 2nd

Ed., p. 344 ff.]

If there are no singularities in F, then U must be continuous and differentiable, so

order of differentiation doesn't matter, and

and therefore F is conservative.

Newton's Universal Law of Gravitation

Newton realized that all of Kepler's laws regarding the motion of the planets
could be explained if a) the planets could be treated as point masses (and he
went off and invented integral calculus to prove this was a good approximation),
and b) if the gravitational force between two objects was proportional to the
product of their masses and inversely proportional to the square of the distance
between them:

In geophysics, we are interested in the force exerted on a "test body" by the

Earth:

Since this force depends on the mass of the body, we can divide both sides by
the mass of the test body (equivalent to using a test mass of 1 unit mass):

For any point mass:

Since any mass distribution can be broken into

infinitesimal ''point masses" and since g is linear in m:

Is gravity conservative?
For point mass, m, observation point P, at a distance r from mass,

Since 1/r has same value at beginning and end of loop, g is conservative. In
fact, by this same reasoning, all central force fields (like f(r), which only
depends r) are conservative.

Potential as Work
Usually define potential as work required to bring unit mass from infinity to
distance r from infinitesimal mass causing the potential:

Defined this way, potential is positive, and tends to zero as r goes to infinity (as
we get an infinite distance from the mass). [Note from Parasnis, 5th Ed., p. 60:
"This is the same definition as that adopted by Kellogg (1953) in his classic
book on potential theory and (implicitly) by Jeffries (1976), among others.
As defined in this manner, V has the units of [m2s-2] and represents the
work done by the field per kilogram of a point mass m0 when m0 moves
from infinity to a distance r from m."]

Computing Gravity from Potential Field

Finding g from U (Cartesian coordinates). With the definition of potential
given above, the acceleration of any point mass towards a mass, m, namely
Gm/r2, is given by:

Finding g from U (spherical coordinates):

Note on signs: defined this way, g will be negative, because it points in the
opposite direction of the unit radial vector. For this reason, you sometimes see g
defined as the positive gradient of potential, so that g (and |g|) will be a positive
number, for convenience.

Integrating over masses to find total field

Because gravity is linear in mass (dm), we could find the gravitational
acceleration due to an extended body by vectorially adding (integrating) the
gravity due to the infinite infinitesimal masses that make up that body, but this
would be complicated. Because potential also depends linearly on mass (dm),
and is scalar, integrating the potential over a body is easier. The potential due to
several (even infinite) dm's is the sum (integral) of the potentials due to individual
dm's. In Cartesian coordinates, for example,

For an arbitrary mass distribution (Cartesian

coordinates)

For an arbitrary mass distribution (spherical coordinates)

Example: what is potential due to sphere of density ?

Deriving Poisson's and LaPlace's equations

First we will derive the Divergence Theorem and Gauss's Theorem

Divergence Theorem and Gauss's Theorem

Consider flux (flow) of material (force lines) through an infinitesimal box:

The flux out of a volume V equals the divergence throughout volume V

Example: For incompressible fluid, flux is zero (no place for fluid to go), so

and since this is true for any arbitrary volume,

Poisson's and Laplace's Equations

Pierre-Simon, marquis de Laplace, born
March 23, 1749, Beaumount-en-Auge,
Normandy, France, died March 5, 1827,
Paris. French mathematician, astronomer,
and physicist who is best known for his investigations
into the stability of the solar system.
Spherical harmonics or Laplace's coefficients
During the years 1784-1787 he produced some
memoirs of exceptional power. Prominent among
these is one read in 1784, and reprinted in the third
volume of the Mchanique cleste, in which he
completely determined the attraction of a spheroid on
a particle outside it. This is memorable for the
introduction into analysis of spherical harmonics or Laplace's coefficients, as also
for the development of the use of the potential - a name first given by Green in
1828.
Simon Denis Poisson, 1781 - 1840. Poisson's most important works were a
series of papers on definite integrals and his advances in Fourier series. This
work was the foundation of later work in this area by Dirichlet and Riemann.
In 1812 discovered that Laplace's equation is valid only outside of a solid.

For gravity,

Consider the net flux out of (or into) a closed volume:

If M is inside the volume, the surface surrounding the mass takes up the entire
"field of view", which is another way of saying that the total solid angle subtended
by the surrounding surface is 4 steradians (a steradian is the 3D equivalent of
a radian; the circumference of a unit circle is 2, hence 2 radians in a circle.
Similarly, the surface area of a unit sphere is 4 steradians).
More on solid angle...

(Laplace's Equation in Integral Form)

From Gauss's theorem,

Since this holds no matter how the volume is chosen,*

If M is outside the volume, total solid angle is 0 (2 ways to look at this: the
surface presents just as much of its front as its back, so they cancel, or notice
that the flux lines which go in one side of the volume bounded by the surface
come out the other side, so the net flux is zero), so

Note that Laplace's equation is just the special case of Poisson's equation where
density is zero.

Applications of Poisson's Equation in Integral Form

1. Gravity due to spherically symmetric body: put imaginary surface
("Gaussian surface") around the sphere

where M is the mass contained within the Gaussian surface.

2. Gravity inside a spherically symmetric hollow shell: put imaginary

surface ("Gaussian surface") anywhere within the hollow region around the
sphere

Since the mass contained within the Gaussian surface is zero,

Optional Assignment: Read Edgar Rice Burrough's At the Earth's Core

Homework problem: find gravity (g(r)) inside and outside a homogeneous
sphere.
3. Gravity due to an infinite slab of thickness h and density : Bouguer's
Formula

Consider "pill box" or cylindrical Gaussian surface

no flux out of sides of cylinder, by symmetry
g through top and bottom must be constant and perpendicular to top and
bottom (again, symmetry), so:

General Solution to LaPlace's Equation in

Spherical Harmonics (Spherical Harmonic
Analysis)

LaPlace's equation is
coordinates,

In spherical coordinates, where r is distance from the origin of the

coordinate system, is the colatitude, and is azimuth or
longitude:

Solutions to LaPlace's equation are called harmonics

In spherical coordinates, the solutions would be spherical
harmonics

Example: show that

, and in rectangular (cartesian)

for point mass (

Solving LaPlace's Equation

Assume variables are separable:

Multiply through by

Last term on LHS depends only on , yet first two do not depend
on , so last term must be constant (and first two must add up to
negative of that constant).

, so

This could be rearranged like this:

This is of the form

where L() has been replaced by x(t)
and the constants "renamed." This is just the ordinary differential
equation for the simple harmonic oscillator problem.

This an ODE, with solution

an integer
Going back to the first two terms, we have

, where m is

Multiply through by

Again, terms must be independent, so both must be constant,

giving this ordinary differential equation:

or
which has the solution
where l is any integer greater than or equal to 0

Finally,

This is another ordinary differential equation, known as

Legendre's Equation, and has solutions of the form
, where
Polynomials,

are the Associated Legendre

are constants

Also, this is kind of neat: The Intel(R) Philanthropic Peer-to-Peer

Program

The general solution to LaPlace's Equation, then, is:

Examples of

Any Legendre polynomial can be found from this generating

function:

Spherical Harmonic Analysis consists of determining values for

(and significance of) constants

o
o

where

for rotating Earth, might neglect dependence, i.e., allow

only m = 0 terms:

are Legendre polynomials

or, for convenience

1. Since a body that is finite in three dimensions (x, y, z) will "look

like" a point mass at infinity, the gravity must tend to GM/r2 as r
goes to infinity, so the potential will go to -GM/r. This eliminates
the C'lm, S'lm terms, because they depend on rl
2. For l = 0, m = 0, the legendre polynomial Plm(cos()) (remember,
this is a function, not a constant times cos(r)) is 1, so C'00 is
identically equal to GM/r, where G is the Univ..., M is the mass of
the body, and r is the distance from it. This term represents the
"sphere" part of the potential.
3. If we set the origin at the center of mass of the body, there will be
as much mass east and west of the center of mass, north and
south of the center of mass, and in front and behind the center of
mass. Therefore, the l = 1, m = 0 term must be zero, because it is
asymmetrical between the northern and souther hemispheres.
So, C10 = 0. This is because P00(cos()) = cos(), which is positive
in the N and negative in the S (or vice versa, since C10, if it
weren't zero, could be negative).

if we pick origin to be center of mass

, n odd, if equator is plane of symmetry, only true on
largest scale...

Other than the coefficients above which can be found from

"common sense" boundary conditions, values for all the other
coefficients are determined by combined satellite and ground
gravity data. The larger "longitudinally symmetrical" terms are:
o
o
o
o

(oblateness)
(pear-shapedness)

o
o

measurements of Earth's gravity field show that the biggest effect

is due to Earth's rotation and bulge

Why Do We Care About Spherical Harmonic Analysis

of Earth's Gravity?

The most complete model for the earths gravitational field,

based on an expansion in a Laplace series, is given by the GEMT2 model. It contains 600 coefficients above degree 36:
More than you'd ever care to know about the GEM-T2 model...
The coefficients define the Earth's gravitational potential at any
point in space (outside the Earth).

One way to visualize the potential field is to look at the shape of

an equipotential surface, usually (and conveniently) the one
corresponding to mean sea level. However...
The l = 0, m = 0 term is the part of the Earth's potential that can
be explained by a perfectly spherically symmetric body (or point
mass). Think of it as a constant term.
The l = 1, m = 0 term expresses the effect of oblateness (as
measured by satellite, but not rotation, since satellites are not
affected by Earth's rotation, although surface gravity
measurements clearly are). Although a thousandth as big as the l
= 0, m = 0 term, it is about a thousand or more times bigger than
the next biggest term.
In order to even see the smaller terms on an equipotential
surface, the oblateness of the Earth, known to have a flattening of
1/298.25, must be subtracted from the "contour map."

Here is what the equipotential surface (geoid) looks like just using some
of the lower degree and order terms:

International Gravity Formula

accounts for variation of gravity with distance from equator

2 effects:
o

rotation of Earth (centripetal acceleration):

where

oblateness of Earth (caused by rotation)

This is the differential equation for the Simple Harmonic Oscillator

(SHO), or a mass on a spring:

where t is time, x is displacement of the spring, m is the mass and k is

the spring constant. The general solution to this equation is:

The undetermined coefficients A and B are determined by initial

conditions (think of them as boundary conditions in time), namely the
position, x, and velocity v, of the mass, when t = 0.

Measuring Gravity
Absolute Measurements

Why?
o
o

ballistics, defense
tectonic studies (e.g., glacial rebound)

mass of Earth
tying relative measurements together
more difficult to achieve high precision than relative
pendulum would work, but can't determine pendulum constant
accurately enough
therefore, use free-fall method
o photograph finely-etched meter stick illuminated by shortperiod, high-intensity flashes at precisely controlled time
intervals (~1 mgal)
o track time of fall of corner-cube reflector with laser
interferometer
o commercial instrument now available (see paper, Carter et al.,
o
o

o
o

EOS)
See also

[Link]
Iodine Stabilized Helium-Neon Laser

Note: You don't have to assume an initial velocity and position of

zero. Given 3 positions and 3 times, you could solve for the
acceleration (how?). In the actual experiments, they gather many
times and positions and then have an overdetermined system in
which they can both improve the accuracy of their acceleration

estimate but also estimate an error on that value.

IGSN71 - network of absolute values

National Geodetic Survey, 1988, Leonard, Oklahoma:

Name

Location

Type

g, microgals

Leonard AA

Seismometer vault

Primary absolute 979,720,911.7

Leonard CA

Old non-magnetic building pier relative

979,720,997.0

Leonard NCMN NCMN

relative

979,721,080.4

Leonard CB

relative

979,738,110.0

Leonard School

Measuring Gravity
Relative instruments

pendulum measurements (still used at

sea)

for point mass on massless string,

where I is moment of inertia about point of suspension; h is distance from
point of suspension to center of mass. But, since K doesn't vary, can
measure change in g:
o Used by Gulf R&D in 1930s; 1-second period (how long),
thermostat, vacuum
o used in late 50s, early 60s by George Woolard at airports, seaports,
large cities, to establish worldwide gravity network

mass on spring

static mass-spring system; k is spring constant:

for 0.1 mgal (10-6 m/s2) accuracy:

interferometer, wavelength of light about 5x10-7 m!
N.B. A 1oC change in T would change length of quartz spring 5.5
microns! (How much of a change in gravity would this appear to be?)

consider system as SHO:

but from above,

, so
thus, increasing the period increases sensitivity (but slows readings)
for ordinary mass-spring system, 20-s period requires 100 m spring/mass
system!

or

LaCoste & Romberg Zero-Length Spring

History of the LaCoste & Romberg meter(s)

moment balance about pivot gives

solving for g:

we want
to be small, and
so a spring with an
unstretched length of zero gives (theoretically) infinite sensitivity.

Worden Gravimeter

Scintrex CG-3 Automated Gravity Meter

"microgal" meter
manual levelling
automatic reading
automatic long-term drift correction

automatic tide correction

stores data

Underwater Gravimeters [From L&R web site]

underwater meters operate on ocean bottom; no averaging as with

shipboard meters
usually shallow; can be modified for use at almost any depth (deep water
operations slow, $$$)
useful in swamps, on muskegs, frozen lakes, ice islands
if tranpsorted by helicopter, can hover over the gravity station while taking
reading
accuracy decreases with depth due to errors in measuring water depth
and position of meter
inherent precision meter about 0.01 mgal
in actual sea operation, base station checks => about 0.1 mgal
water depth usually measured with pressure gauges; accuracy ~1/2% (0.6
meter error in depth => ~0.1 mgal)
overall accuracy of about 0.2 mgal is considered good in a survey in water
160 meters deep

Units

Systeme International (SI; MKS): m/s2

cgs: cm/s2 (gal)
milligal = mgal = 10-3 gal; microgal = gal = 10-6 gal
gravity unit = gu = 0.1 mgal
1 microgal = 0.7 mph/year!

typically desire survey accurate to 0.1 mgal (100 gals);

g ~ 980,000,000 microgals!
to resolve anomalies, must adjust observations for several effects

Our Meters

History of the LaCoste & Romberg Meters

Dr. Arnold Romberg (1882-1974)
[photo courtesy David Dillon, Jr.]

The LaCoste & Romberg Company was begun

in 1939 by a graduate student at the University
of Texas, Lucien LaCoste, and his faculty
advisor, Dr. Arnold Romberg. An interesting
sidenote is that Dr. Romberg's grandson lives
here in Norman. He scanned this photo of his
grandfather for me from a family scrapbook!
For a complete history of the development of
LaCoste & Romberg meters and their company,
see:
[Link]

Airborne Gravity and Gravity Gradiometry

From the Grav-Mag mailing list:
Those of you who access the CSIRO web site should be wary of comparisons with the US Air Force
Gravity Gradiometer Survey System (GGSS), which we developed here in one of our previous incarnations
as Bell Aerospace. This technology is now more than fifteen years old and was an early exposure to using a
gravity gradiometer in an airborne environment. Since those initial flights in a C130 Hercules, much
experience has been gained operating different instruments in a variety of environments, which has resulted
in improved performance. In addition BHP and Bell Geospace have developed or improved upon software
techniques to optimize survey performance for their specific applications. Airborne results from the BHP
Falcon system give a better idea of current capabilities. Although we cannot make any independent
statements about BHP's Falcon gravity gradiometer performance, we believe that the specifications which
can provide a more accurate assessment of current performance can be found on the BHP web site:
[Link]
As described there, the Falcon system has been developed to operate in the higher turbulence experienced
at low terrain clearance (80m), thus improving gradient signal amplitude, which decays as the cube of
distance from source. This substantially changes the CSIRO conclusions about the Eotvos sensitivity
needed to identify targets, since they anticipate 300m as a nominal altitude.
The Falcon instrument measures two curvature gradients. To create a full tensor, five gradient, airborne
instrument, the successful full tensor marine system is now being further developed to optimize it for
airborne applications. Lockheed Martin is currently inviting enquiries from organizations who would be
interested in participating in flight trials of such a system.
Andy Grierson

Data Acquisition, Reduction of Gravity Data

Establishing base stations

tie to IGSN71 or other absolute reading if possible; below are some OK

stations; click to see values:

establish local base

o close to or within survey area
o stable, permanent, easy to find
o good control on x, y, z (benchmark, road intersection, airport)

can establish new base by looping from old base

Determining Elevations

surveying (costs >= gravity survey!)

topo sheet
inertial guidance system
altimeter
GPS, DGPS: GPS Primer
RTK dual-frequency DGPS (2 units with real-time
communication/correction between them) gives cm-level accuracy in
seconds; cost (Feb. 2001) is about $45K
DEMs 7.5-Minute DEM: 30x30-meter data spacing
More on DEMs (10 m vs. 30 m)

Determining Horizontal Position

same as elevation
digitizer handy for getting latitudes

Adjusting Observed Gravity

Tidal correction

secular variations in g are (generally) undesirable

effect of Sun about 50% of Moon
each has 12-hour period (front and back bulge),
but tidal inequality...
2 effects:
o pull of bodies on meter
o distortion of Earth (solid earth tide); adds about 12% to this effect
total magnitude about 0.2, 0.3 mgal (refer to tidal correction lab)
to correct:
o fixed recording gravimeter

located in or near survey area

subtract variations from survey data

probably most accurate correction, but $$
o tide tables (gravity)

read tidal correction for given time and location

many not apply well near water

no longer published
o calculate tidal effect

computer program yields correction for time and location

incorporated into Scintrex meter
o include in drift correction

data from Sandia's Lacoste and Romberg Model G meter, sitting in

Gilbert's lab, connected to chart recorder (light line); dark line from my tide
computer program

Drift Correction

secular variations in g are undesirable

instrument drift (DT, DP, creep), tides
assumptions
o changes are smooth and slow
o changes are independent of location
drift estimated by reoccupation of (base) station
must reoccupy every 2 hours or so:

using same base for entire survey not practical (or necessary):

Station

Time

Dial Divisions Drift Rate

Base

11:20

762.71

GN1

11:42

774.16

GN2

12:14

759.72

GN3

12:37

768.95

GN4

12:59

771.02

Base

13:10

761.18

Elapsed
Time

Correction

Corrected
Reading

shift correction; the next day...

Dial
Divisions Shift

Station

Time

Base

10:20

763.68

GN5

10:42

775.16

GN6

11:14

765.42

GN7

11:37

765.35

GN8

11:59

770.32

Base

12:10

760.28

Drift

Elapsed

Rate

Time

Corrected
Correction Reading

Solution to above; don't look at this until you've tried it yourself!

3-point drift correction
2
o assume g = at + bt + c
o evaluate a, b and c using g known at one station at t1, t2 and t3
o Excel spreadsheet example

Latitude correction

Geodetic Reference System Formulae refer to theoretical estimates of the

Earth's shape
From these GRS formulae we obtains International Gravity Formulae
(IGF)
Several different formulae have been adopted over the years

In these equations,
is geographic latitude and
is commonly referred
to as theoretical gravity or normal gravity
o First internationally accepted IGF was 1930:

o
o

This was found to be in error by about 13 mgals; with advent of

satellite technology, much improved values were obtained.
The Geodetic Reference System1967 provided the 1967 IGF:
Most recently IAG developed Geodetic Reference System 1980,
leading to World Geodetic System 1984 (WGS84); in closed form it
is:

The IGF value is subtracted from observed (absolute) gravity data. This
corrects for the variation of gravity with latitude

International Gravity Formula "Calculator"

Latitude correction: short form

note that over small range, curve is nearly straight-line slope

where is a typical latitude for the (small) field area

miles north, in above formula
subtract this amount for each mile north

Burger Table 6-1

Free-air correction

"flag-pole" correction
accounts for decrease in g(r) due to change (increase) in elevation (r)
does not take in into account mass which may be elevating you

could pick datum (e.g., sea level), compute "exact" difference, but over
small changes in elevation, h, change is nearly linear

note that, once IGF is removed, we take elevations relative to sea level,
not center of Earth!

Taylor Series

this is the free-air effect: g decreases (hence sign) with elevation

alternative "derivation":

how big an effect? take g = 9.83 m/s2; r = 6371 km,

correction: add 0.3086 times elevation in meters

Bouguer correction

"bulldozer" correction
accounts for mass (rock) responsible for elevation change (between
observation point and sea level, usually)
depends on density of material holding you up (Bouguer density)
gravity due to infinite slab, thickness h (m), density (g/cm3):

slab holding you up increases g; must subtract this effect out

requires knowledge of h and

Burger Table 6-2

Selecting Reduction Density

while sea level may be datum, variation in elevation occurs in near

subsurface
methods for selecting Bouguer density:

1. Standard Bouguer density = 2670 kg/m3

nominally average crustal density

ensures continuity between surveys

2. Direct measurement

collect samples, core, drill samples, hand samples

inaccessibility; may not be representative

3. Geologic map to get rock type; get values from tables, graphs, etc.

handbooks (physical properties of rocks and minerals)

sedimentary rock density histograms
rock density ranges
rock density means and ranges
salt density vs. sediment

Rock Type

Density

Ice

880 - 920

sea water

1010 - 1050

Shale

1950 - 2700

limestone, dolomite

2500 - 2850

sandstone

2100 - 2600

soil & alluvium

1650 - 2200

rock salt

1850 - 2150

felsic igneous rocks

2550 - 2750

mafic igneous rocks

2700 - 3000

ultramafic rocks

3000 - 3300

4. Density profile (Nettleton method)

collect closely-spaced g readings over topographic feature

make latitude, free-air correction

make Bouguer correction, with various values of

find which gives least correlation with topography

5. Logs

gamma-gamma density logs

o -ray source (usually Cobalt), geiger counter
o when -rays (photons) collide w/ electrons in rock, Compton
scattering causes them to be reflected
o intensity of reflected beam depends on density of electrons, which
depends on density
neutron density logs
seismic velocity logs (due to rough correlation of density w/ velocity)
borehole gravity meter (see [Link])
o measures g to 0.01 mgal

must overcome small hole size, hostile environment, self-leveling

"width" of investigation ("penetration into sides of borehole") ~ point

separation
best method for getting true formation density

6. Linear regression (least squares) method

assumes no correlation between topography and subsurface density (i.e.,

anomalies are randomly distributed with respect to elevation)
therefore correlation between topography and g will be due to Bouguer
slab
plot gfa vs. elevation, h
fit line through points
slope will approximate 2G; solve for (Bouguer density)

Least Squares Fit

data pairs: x1, y1; x2, y2;...xn, yn, where n is the total number of data points
straight line: y = mx + b, where m is the slope, b the intercept

minimize the sum of the squares of the residuals (SSR):

setting derivatives = 0 gives the normal equations:

these two equations are used to solve for unknowns, m and b. For
homework, show that:

Problems with linear regression method for getting density:

if density varies (decreases in this example) with elevation, will get curve

density may correlate with elevation; e.g., Arbuckles: limestones are hillformers, shales are valley-formers...

Bouguer correction at sea, underground

surface survey

first term replaces water with crust

second term Bouguer corrects to sea level
underwater survey (for accurate value at sea)

underground survey

o
o

Other corrections to gravity at sea

not a stable platform as in land gravity

accelerations can be of order

, >> accuracy desired
corrections also apply to airborne gravity

Acceleration correction

ship can't accelerate up or down very long, so average over t eliminates az

(not so for aircraft)

must know ax, ay (accelerometers, gyroscope to keep accelerometers

level)
natural period of land g meter typically ~10s, but waves also ~10s, so ship
gravimeter must have much longer period

Eotvos correction (Nettleton, p. 116 - 118)

important for shipborne and airborne gravity

component of velocity in E direction increases apparent angular velocity
(Coriolis Force)
biggest single source of error in shipborne gravity comes from error in V,
(although GPS, particularly DGPS, probably helps significantly)
centrifugal acceleration

component in vertical direction

effect in change in angular velocity, ( will be small compared to on

ship or even plane)

if velocity is V, E component is

angular velocity due to this motion is

so, daV is given by

there is also a simple acceleration in vertical direction due to total velocity,

V

Total Eotvos Correction

For V in knots, correction in mgals

Example: 10 knots E at equator => 75 mgals!

Terrain correction

Bouguer correction assumes infinite slab; terrain correction corrects for

this erroneous assumption!
o Bouguer correction works for gently sloping surfaces, like
topography on basement
o error < 3% for slope < 1/5 (see Adams and Hinze, vol. 3, SEG Geotech. &

Terrain correction
o always positive
o requires detailed info on elevation around station, not just at station

Environ. Gphy., p.99)

size of terrain correction depends on relief and its proximity to

station

Burger Table 6-3

Hammer Terrain Correction Chart
Terrain Correction with DEMs

advent of DEMs has made medium to far-field correction much easier

short-field correction may be done using "newly available, reflectorless
laser rangefinders. Such rangefinders permit a detailed digital terrain data
to be acquired in the vicinity of a gravity station within only 2-3 minutes,
permitting the gravity meter operator to acquire the terrain data needed ...
at the same time as the gravity measurements are made."
Geodesy Group at Curtin University - Terrain Effect

Geophysical Software
Sample DEM: Contours, Color Shaded, Shaded Relief, Map

Free-air, Bouguer, Isostatic Anomaly

latitude and free-air correction virtually always made, giving the Free-Air
Anomaly (FAA)

for environmental/exploration work, some Bouguer correction is made

for gravity in mgals and elevation in meters, then

The Standard Bouguer Anomaly uses Bouguer density of 2.67 g/cm3

note that, once is chosen, FAC and BC can be combined into one
correction: the elevation correction

the Complete Bouguer Anomaly also includes the terrain correction

sometimes make a correction for isostasy, giving Isostatic Anomaly

Estimating Survey Error

dependent errors add algebraically

independent errors add vectorially
sources of error:
o meter dial
o meter consistency
o drift
o latitude
o elevation (Free-air and Bouguer partially cancel)
o terrain
o others?

determining error of calculated quantity

Error: A Calculus/Physics Refresher

Suppose you want to find the mass of a homogeneous planet of density and
radius R. You know the volume of a sphere, and that the mass would be the
volume time the density, so you have:

Now, you know that you can never determine neither the density nor the radius
exactly, so you'd like to know how much error there would be the mass
calculation if you make a (small) error in the density of radius.
Let's start with density. What we're really asking is, "how much of a change will
there be in mass if we have a small (technically, infinitessimal) change in density.

This should ring a bell from some math class you took. The quantity
represents the (infinitessimal) change in mass due to a (infinitessimal) change in
density, with R constant. Performing the derivative, we get

This first equation will have the units of mass/density (which happens to be
volume). So if we know density to an accuracy of 100 kg/m3, and the radius of
the planet is 6371 km,

and our error in mass turns out to be that amount. (For reference, the bulk
density of the Earth is about 5500 kg/m3, and the mass is about
The error of 100 kg/m3 is about 1 part in 55, then, and so the mass is off by
about the same ratio.)
Now let's deal with an error in radius. Again, we want to know the expected
change in mass for a small change in radius, so

Notice that the error in M due to an error in R now depends on R! In other words,
if we make an error of 1 km for a big planet, the error in mass will be much
greater than for a small planet. (Physically, you can picture the shell of additional
material 1 km thick; its surface area will vary as R squared.)
Now, you might say, "if we don't know R precisely, how do we know what R to
use in the formula?" And the answer is, you don't, but you can still estimate the
error, even though you won't know the amount of error exactly!

Anomaly Separation and Filtering

applies to gravity and magnetics!

gravity and magnetic fields are each the sum of many anomalies
core, mantle, relief on Moho, deep crustal density anomalies, all produce
long-wavelength anomalies
short-wavelength effect of shallow crustal bodies may overlap
N.B.:
o deep bodies produce only long wavelengths; shallow bodies can
produce long or short wavelengths!
o shortest wavelength that can be produced at given depth is due to
point mass (or spherically symmetric mass)
in exploration surveys, usually interested in shallow features
removal of long-wavelength anomalies (usually due to deep sources)
enhances short-wavelength (i.e., necessarily shallow) features
mustn't "throw out baby with bath water"

regional field: long-wavelength "background" field

residual field: total field - regional field

removal of regional corresponds to high-pass filter (pass high

wavenumber anomalies, reject small)
methods:
o graphical smoothing
o polynomial fitting
o averaging
o upward continuation
o wavelength filtering
an example of the need for regional removal/anomaly separation:
o LA Basin, regional removal

Graphical Smoothing

graphically estimate regional field

subtract regional from total to produce residual
time-consuming, very subjective (good and bad)

Profile

sketch regional field

subtract regional from total to create residual

profile, regional removal

another profile, regional removal

Contour Map

sketch contours representing regional field

estimate regional values at control points
subtract regional values at control points from (total) values at control
points
re-contour residual control point

contour data, smoothing

ring method, contour data
picking ring size

Polynomial Fitting, Trend Surfaces (method of

least squares)

fit a smooth (polynomial) surface to data to represent regional

subtract calculated regional value from observed value at each point to get
residual field
zeroth-order polynomial (constant, g0)
o minimize Sum of Squares of Residuals (SSR)

first-order polynomial in x and y

3 equations
3 unknowns (a, b, c)
second-order polynomial in x and y

north-central Iowa trend surfaces

o
o

Gridding irregularly spaced points

most filtering and averaging schemes require regularly spaced (gridded)

points
can overlay grid on contour map, interpolate values at grid points (groan)
or, use computer program (like Surfer)

Inverse Distance Squared

compute value at grid point based on neighbors

use nearest N points, or points within R of grid location
value at P is (inverse distance squared) weighted average of selected
neighbors

Kriging

fit analytical surface to data points

use surface to compute values at grid points
"Kriging provides a means of interpolating values for points not physically
sampled using knowledge about the underlying spatial relationships in a
data set to do so. Semivariograms provide this knowledge. Kriging is
based on regionalized variable theory and is superior to other means of
interpolation because it provides an optimal interpolation estimate for a
given coordinate location, as well as a variance estimate for the
interpolation value."

Cubic Spline

cubic spline is shape an elastic rod would take if constrained to fit at

control points
bicubic spline is shape an elastic plate would take if constrained to fit at
control points

Minimum curvature

find surface with least curvature that fits points withing a certain tolerance
requires odd number of grid points
Briggs, Ian C., 1974, Machine contouring using minimum curvature.
Geophysics, v. 39, pp. 39-47.
works well with profile data
works well with digitized contour data

Smoothing by averaging

Running or moving average

n-point running average

here weighting factor is 1.0 for all points

Weighted averaging

may want to weight closer points more (or less)

in 2 dimensions (25-point average, for example):

grid method of regional removal

running average, profile
matrix smoothing (running average), contour data

Vertical Derivatives

often take first or second vertical derivative of field

enhances short wavelength anomalies relative to long wavelength

anomalies (?-pass filter)
note, for modelling purposes, that result does not have units of g
tends to delineate edges of anomalous body

4-point average and the 2nd vertical derivative

form residual by subtracting average of 4 closest points from center point

now consider Laplace's equation

in discrete form, we use finite differences

the first backward, first forward, and second central differences are

similarly

so for the second vertical derivative (assume x = y)

which differs from the 4-point grid average only by a constant

thus, second vertical derivative amounts to horizontal "curvature" of field
effect of 2nd derivative
radius of curvature and 2nd derivative
buried river channel [modelled gravity; Burger Fig. 6-31]
salt domes, Texas
LA Basin
Cement field, OK

Upward Continuation

Stokes's Theorem: If gravity values are known everywhere on

Earth's surface, gravity at any higher point can be calculated from
these values
knowing field at one elevation, can compute what field would look like at a
higher elevation (upward continuation) or lower elevation (downward
continuation

this amounts to 1/r3 weighted averaging

downward continuation enhances near-surface bodies more than deeper
bodies, hence lessens effect of regional
consider effect on bodies at different depths, d1 and d2

original photo
cropped
resized
Lanczos filter resampling

Wavelength Filtering - the Fourier Transform

High-Pass Filtering
Unfiltered: Queen: Bohemian Rhapsody
Passes high wavenumber components; attenuates low wavenumber or long
wavelength components
Queen: Bohemian Rhapsody, high-pass at 2000 Hz

Low-Pass Filtering
Passes low wavenumber components; attenuates high wavenumber or short
wavelength components
Queen: Bohemian Rhapsody, low-pass at 1000 Hz

Bandpass Filtering
Passes a band of wavenumbers, i.e., frequencies between u1 and u2 are kept,
rest eliminated (or at least attenuated).
Queen: Bohemian Rhapsody, band-pass between 1000 Hz and 2000 Hz

consider, for simplicity, a 1-D function, g(x) (e.g., a gravity profile), in the
space domain, g(x), and wavenumber domain, G(u):

high-pass filtering consists of passing (keeping) high wavenumber (short

wavelength) terms

low-pass filtering consists of passing (keeping) high wavenumber (short

wavelength) terms

band-pass filtering consists of passing (keeping) only that portion of the

wavenumber spectrum in the band u1<u<u2
the filter function, H(u) will be

the bandpass filter is sometimes called a "boxcar" filter because of its

shape:

performing point-by-point multiplication, then, gives the filtered spectrum:

the inverse Fourier transform is then applied to F(u) to obtain f(x)

2-D Bandpass Filtering

practically, a "ramped" bandpass filter must be employed

this mitigates Gibbs phenomenon - a "ringing" that results from too sharp
a filter

test
test low
test high
sara (key clicks)

eastern Sierra Nevada, complete Bouguer, low-pass>50 km

eastern Sierra Nevada, high-pass>50 km
Garber oil field, Oklahoma
Regional field used
Residual Garber gravity

Other Filter Operations

Filter

H(u,v)

Upward (z<0), downward (z>0) continuation,

where z is continuation height (depth)
1st vertical derivative

2nd vertical derivative

[Note: The 2 terms above may be absent in some texts, depending on whether
the Fourier transform has 1/(2) term]

Directional Derivatives

comments from Dr. Lyatsky on directional derivatives

Wavelength Filtering - Wavelet Processing

Localized analysis of frequency content

used, e.g., for better compression of images (JPEG2, or Lurawave)

notice that analysis of period is different for different time periods

Interpretation of Gravity Data

Uniqueness (ambiguity) Problem

applies to all potential field methods, and, indeed, all geophysical methods
there is an inherent ambiguity in interpretation of gravity data
even if you had gravity at every point on Earth's surface, there are multiple
models that would produce those values because of integral nature of gravity, it
can be proven that any anomaly can be result of an infinite number of density
distributions!

Constraining Interpretations

not hopeless: gravity eliminates even "more infinite" number of density

distributions
combine gravity data with other constraints:
o density of crustal rocks, particularly in local area
o configuration of rocks: well data, regional geology, etc.
o other geophysics: magnetics, seismic, etc.

Interpretation approaches
Direct Interpretation: Inverse method

only possible if many constraints (artificial?) imposed

1. assume general class of model (e.g., buried sphere)

2. analyze anomaly (anomalies) to define specific model

Indirect Interpretation: Forward modelling

1.
2.
3.
4.
5.

assume specific initial subsurface density model

calculate gravity (always do-able, at least numerically)
compare with data
adjust density model as necessary
repeat steps 2 through 4

Role of interpretation in survey planning

should the survey be conducted?!

how should it be conducted?

Simplifying density models

because we are interested in (and indeed only measure) change in g, we

are only interested in changes in density (density contrast)
background density can always be subtracted
furthermore, horizontal slab doesn't contribute to gravity anomalies
in some cases (e.g., sphere, horizontal cylinder), mass
excess/deficiency, is determinable quantity
a datum shift may be made to compare model to data

Another example illustrating the ideas of gravity anomaly and density

contrast:

Vertical g

gravity of Earth >> any anomaly

gravity defines vertical
therefore gravity meters only measure vertical component:

Gravity due to simple bodies

easiest, most versatile approach to survey planning, interpretation

even complex structures produce anomalies similar to simple shapes; for
example, horizontal circular cylinder versus square cylinder:

Infinite Slab

, where d is the thickness

works for gently sloping surfaces

example: topography on basement
error < 3% for slope < 1/5 (see Adams and Hinze, vol. 3, SEG Geotech. & Environ. Gphy., p.99
use estimate magnitude of anomaly for many flat-layer situations:
o relief on density contrast boundary: basement, bedrock, etc.
o dip-slip fault in horizontal strata
o laterally extensive mines
o water removal/recharge in horizontal aquifer

Sphere

applicable to approx. equidimensional bodies (longest dimension <<

depth)
gravity due to sphere

vertical component:

Since,

, we get

Where is g 1/2 of maximum?

"Inversion"

technique
1. find distance (x1/2) from peak of anomaly where anomaly is half maximum
anomaly:

2. depth of body:
3. since

Notes

non-uniqueness:

(can't determine radius AND density contrast)

anomaly size is relative to "baseline", or "background"
in general, half-width to left and right are unequal

Example:
gmax = 35 mgals, 1/2-width at 17.5 mgals = 7.5 m; therefore z = 10 m; assuming
density contrast of 1.0, find radius of sphere.

Infinite Horizontal Cylinder

applicable to bodies much longer in one horizontal direction than in

vertical or other horizontal direction

tunnels, river channels, horst or graben block, etc.

This is a max at x=0, or

where we've dropped the subscript v.

Find the relationship between half-width and depth:

Vertical Cylinder

no simple expression exists for gravity off-axis; but on axis:

special case: infinite slab

special case: semi-infinite cylinder

Note that finite cylinder formula can be derived by superposition of two semiinfinite cylinders

Narrow Vertical Cylinder

while no simple expression exists for gravity off-axis of a thick cylinder, for
thin cylinder an approximate solution exists
good for z > 2a

semi-infinite case

depth criterion:

knowing z, we can now find

finite case

use superposition:

Semi-infinite Horizontal Slab

finite horizontal slab

thin semi-infinite sheet

for a line,

so for a sheet,

Note that

Depth criterion

thin finite sheet

Computing depth:

Compute

at center of anomaly:

Compute

to one side:

On graph with no vertical exaggeration, find depth which yields these

angles:

2D Grids
The vertical gravity component due to a line element of mass per unit length is:

Now consider an arbitrary 2-D body of density :

The area of the shaded element is dz*dx, so

For , z constant, each block contributes the same to g

"Computer" Methods of Interpretation

Talwani, 1973, in Bolt: Computational Methods in Geophysics

3D laminar bodies

Talwani and Ewing, 1960, Geophysics, v. 25, 203-225

Pluoff, 1976, Geophysics, v. 41, 727-739

uses solid angle approach

for thin horizontal sheet,

2-D polygon method

Talwani, Sutton and Worzel, 1959, JGR, 64: 1545 - 1555
Talwani, Worzel and Landisman, 1959, JGR, 64: 49-59

uses line integral approach

This is the method used in most 2D computer modelling programs, like GM-SYS.

3D vertical prisms (method of Cordell and Henderson)

Three-dimensional iterative method of Cordell and Henderson

(Cordell, L., and Henderson, R. G., "Iterative three-dimensional
solution of gravity anomaly data using a digital computer,"
Geophysics 33 (1968), 596-601). Block heights are relative to a
predetermined reference surface; density contrast is also
predetermined. Initial block height might be determined by using
the gravity value above the block and using the Bouguer slab
approximation. Then slab heights are adjusted to give a best fit to
the measured gravity values (or a gridded gravity field derive from
measured gravity data). Figure from Blakely, 1996.

Danes, 1960, Geophysics, 25: 1215-1228

square prisms, infinite bottom depth; get finite prisms by using another set
of prisms
iterative method:
o pick set of prisms
o find g at center of each prism
o adjust heights to match actual field
o close to direct approach (inversion), given assumptions

General 3D Bodies
GRVMAG message from Manik Talwani re: his 3D G&M inversion program
(5/10/2001)

Bhattacharyya, B. K., Navolio, M. E., 1976, A Fast Fourier Transform

method for rapid computation of gravity and magnetic anomalies due to
arbitrary bodies: Geophys. Prosp., 24, 633-649.
Gerard, A., Debeglia, N., 1975, Automatic three-dimensional modeling for
the interpretation of gravity or magnetic anomalies: Geophysics, 40 (6),
1014-1034.
Talwani, M., Ewing, M., 1960, Rapid computation of gravitational attraction
of three-dimensional bodies of arbitrary shape: Geophysics, 25 (1), 203225.
Okabe, M. 1979, Analytic expressions for gravity anomalies due to
homogeneous polyhedral bodies and translations into magnetic
anomalies. Geophysics v44, p730-744.

Interpretation Examples

Infinite Slab: Bedrock depths, Reading, Mass.

Subsurface voids, Medford Caves, Florida and without vertical
exaggeration
Valley geometry, Pine Valley, central Nevada
o location, Bouguer gravity map

o
o

gravity, geologic profiles

o
o

horizontal cylinder model

double-cylinder model

o
o

2-D polygon model