Earth's gravity field

For the irrigation district's long leveling lines, Ruth Albrecht wanted gravity values at the benchmarks so that the heights could be computed rigorously. Theo Martel spent two days with a borrowed relative gravimeter, starting and ending each day on a station of the national gravity network maintained by NGS, and reading the meter on every benchmark and on pillar CR1. The readings looked almost identical, about 9.81 m/s² everywhere, yet the differences of a few hundredths of a percent between them carry the information that turns leveled differences into heights. This lesson describes the Earth's gravity field, the normal gravity model that serves as its smooth reference, and the quantities that geodesists compute from gravity measurements.

Gravitation, rotation and gravity

Newton's law of gravitation says that every mass attracts every other mass with a force proportional to the product of the masses and inversely proportional to the square of the distance. For the Earth as a whole, the attraction at distance \(r\) from a spherically symmetric Earth would be \(GM/r^2\), where \(GM = 3.986\,005 \times 10^{14}\ \text{m}^3/\text{s}^2\) is the geocentric gravitational constant. Geodesists use \(GM\) rather than \(G\) and \(M\) separately, because the product is known from satellite orbits to about one part in a billion, while \(G\) alone is known only to about one part in ten thousand.

A gravimeter on the ground rotates with the Earth, so it also feels a centrifugal acceleration directed away from the rotation axis, of magnitude \(\omega^2 p\), where \(\omega\) is the angular velocity and \(p\) the distance from the axis. Gravity is the vector sum of gravitation and centrifugal acceleration. At the equator the centrifugal term reaches about 0.034 m/s², reducing gravity by roughly 0.3 percent; at the poles it vanishes. Combined with the flattening, which puts polar points closer to the center, this makes gravity at the poles about 0.53 percent larger than at the equator.

Gravity is measured in meters per second squared, but geodesy and geophysics use the gal, named after Galileo: 1 Gal = 1 cm/s² = 0.01 m/s², and 1 milligal (mGal) = \(10^{-5}\) m/s². Gravity at the Earth's surface is about 980,000 mGal. Modern absolute gravimeters measure it to a few microgals, a precision of a few parts per billion.

Potential and level surfaces

Working with a vector field is cumbersome, so geodesy uses the gravity potential \(W\), a scalar whose gradient is the gravity vector. It is the sum of the gravitational potential \(V\) and the centrifugal potential:

The sign convention in geodesy makes \(W\) positive and increasing toward the Earth. A surface on which \(W\) is constant is an equipotential or level surface. Water at rest settles on such a surface, and a spirit bubble centers when the instrument is tangent to one. Gravity is everywhere perpendicular to the level surfaces, so a plumb line is a curve that crosses each level surface at right angles. Moving along a level surface requires no work against gravity; moving from one level surface to another requires work equal to the difference in potential. That is the physical meaning of height, and lesson 7 builds on it.

Level surfaces are not parallel. Where gravity is stronger, level surfaces are closer together, because the same potential difference corresponds to a shorter distance: \(dW = -g\,dH\). Since gravity increases toward the poles, level surfaces converge toward the poles. This simple fact is why leveled height differences depend on the route taken, and why rigorous heights need gravity.

The normal gravity field

The actual gravity field is irregular, so geodesists compare it with a smooth reference, the normal gravity field. It is generated by a level ellipsoid: an ellipsoid of revolution that is itself a level surface of its own gravity field, with the same mass, rotation rate and flattening as the Earth. GRS80 defines such an ellipsoid. Its four defining constants, \(a\), \(GM\), \(J_2\) and \(\omega\), fix the geometry, the potential on the surface (\(U_0 = 62\,636\,860.85\ \text{m}^2/\text{s}^2\)) and gravity everywhere outside it.

On the surface of the level ellipsoid, normal gravity depends only on latitude and is given in closed form by Somigliana's formula:

with \(\gamma_e = 9.780\,326\,7715\ \text{m/s}^2\) at the equator and \(\gamma_p = 9.832\,186\,3685\ \text{m/s}^2\) at the poles for GRS80, so that \(k = 0.001\,931\,851\,35\). An older series form, the International Gravity Formula 1980, agrees with Somigliana's formula to better than 0.1 mGal and is still found in textbooks and software.

Gravity changes with height

Above the ellipsoid, normal gravity decreases. Near the surface the rate is the free-air gradient, about −0.3086 mGal per meter. A more precise second-order formula is

The decrease of gravity with height is substantial: climbing 1000 m reduces gravity by about 309 mGal, far more than the variations caused by geology. That is why raw gravity values cannot be compared directly; they must be referred to a common surface or compared with normal gravity at the right height.

Disturbances and anomalies

The difference between actual and normal gravity carries the information about the Earth's irregular mass distribution. Two definitions are in use, and confusing them is a common error. The gravity disturbance compares gravity at a point with normal gravity at the same point: \(\delta g = g_P - \gamma_P\). The gravity anomaly compares gravity at a point with normal gravity at a corresponding point on the reference that is displaced by the geoid height: \(\Delta g = g_P - \gamma_Q\). In practice, the free-air anomaly uses normal gravity at a height equal to the orthometric (or normal) height above the ellipsoid, rather than the ellipsoidal height. Because GNSS now gives ellipsoidal heights directly, the disturbance has become easier to compute, but the anomaly remains the quantity used in classical geoid computation by Stokes's formula.

Geophysicists go one step further. The simple Bouguer anomaly subtracts the attraction of the rock between the station and the geoid, modeled as an infinite flat slab of standard density 2670 kg/m³, which attracts \(2\pi G\rho H = 0.1119\) mGal per meter of height. The complete Bouguer anomaly also corrects for the surrounding terrain. Bouguer anomalies over the continents are strongly negative, because mountains and plateaus are supported by low-density roots, the principle of isostasy.

# Normal gravity on GRS80: Somigliana's formula, height correction and simple anomalies
import math

a, f = 6378137.0, 1 / 298.257222101
GM, omega = 3.986005e14, 7.292115e-5
ge, gp = 9.7803267715, 9.8321863685          # normal gravity at equator and pole, m/s^2
b = a * (1 - f); e2 = f * (2 - f)
k = b * gp / (a * ge) - 1                    # Somigliana constant
m = omega ** 2 * a * a * b / GM
print("k =", format(k, ".12f"), "  m =", format(m, ".12f"))

def somigliana(lat):
    s2 = math.sin(math.radians(lat)) ** 2
    return ge * (1 + k * s2) / math.sqrt(1 - e2 * s2)

def series(lat):                             # the 1980 International Gravity Formula series
    s = math.sin(math.radians(lat))
    return 9.780327 * (1 + 0.0053024 * s * s - 0.0000058 * math.sin(math.radians(2 * lat)) ** 2)

def gamma_h(lat, h):                         # normal gravity at ellipsoidal height h (second order)
    s2 = math.sin(math.radians(lat)) ** 2
    return somigliana(lat) * (1 - 2 / a * (1 + f + m - 2 * f * s2) * h + 3 * h * h / (a * a))

mGal = 1e-5                                  # 1 mGal = 1e-5 m/s^2
for lat in (0, 30, 45, 47.506, 60, 90):
    g0 = somigliana(lat)
    print("lat", str(lat).rjust(7), " gamma0 =", format(g0, ".7f"), "m/s^2   series - closed =",
          format((series(lat) - g0) / mGal, "+.3f"), "mGal")

lat = 47.506
w2p = omega ** 2 * (a / math.sqrt(1 - e2 * math.sin(math.radians(lat)) ** 2)) * math.cos(math.radians(lat))
print()
print("centrifugal acceleration at the site =", round(w2p / mGal), "mGal (", round(w2p, 5), "m/s^2 )")

h = 1004.8
gh = gamma_h(lat, h)
print("normal gravity at h =", h, "m:", format(gh, ".7f"), "m/s^2")
print("vertical gradient near the ellipsoid =", round((gamma_h(lat, 1.0) - somigliana(lat)) / mGal, 4), "mGal per m")

H = 1019.9                                   # orthometric height of the pillar, m
g_obs = 9.8055300                            # gravity measured on the pillar, m/s^2
dist = (g_obs - gh) / mGal                   # gravity disturbance: normal gravity at the same point (h)
anom = (g_obs - gamma_h(lat, H)) / mGal      # free-air anomaly: normal gravity at height H above the ellipsoid
boug = 2 * math.pi * 6.674e-11 * 2670 * H / mGal   # Bouguer plate of standard density 2670 kg/m^3
print("gravity disturbance =", round(dist, 1), "mGal;  free-air anomaly =", round(anom, 1), "mGal")
print("Bouguer plate =", round(boug, 1), "mGal;  simple Bouguer anomaly =", round(anom - boug, 1), "mGal")

k = 0.001931851355   m = 0.003449786003
lat       0  gamma0 = 9.7803268 m/s^2   series - closed = +0.023 mGal
lat      30  gamma0 = 9.7932487 m/s^2   series - closed = +0.055 mGal
lat      45  gamma0 = 9.8061992 m/s^2   series - closed = +0.067 mGal
lat  47.506  gamma0 = 9.8084650 m/s^2   series - closed = +0.067 mGal
lat      60  gamma0 = 9.8191784 m/s^2   series - closed = +0.047 mGal
lat      90  gamma0 = 9.8321864 m/s^2   series - closed = -0.016 mGal

centrifugal acceleration at the site = 2295 mGal ( 0.02295 m/s^2 )
normal gravity at h = 1004.8 m: 9.8053655 m/s^2
vertical gradient near the ellipsoid = -0.3085 mGal per m
gravity disturbance = 16.4 mGal;  free-air anomaly = 21.1 mGal
Bouguer plate = 114.2 mGal;  simple Bouguer anomaly = -93.1 mGal

The program first confirms that the Somigliana constant computed from the defining values matches the published one, then tabulates normal gravity. The series formula agrees with the closed form to less than 0.07 mGal everywhere. At the pillar, the centrifugal acceleration is about 2300 mGal, a reminder of how much rotation contributes. Theo's measured gravity, 9.805 530 m/s², exceeds normal gravity at the same point by 16.4 mGal, the gravity disturbance. Using the orthometric height of 1019.9 m, which is 15.1 m greater than the ellipsoidal height because the geoid lies below the ellipsoid here, gives a free-air anomaly of 21.1 mGal; the 4.7 mGal difference is the free-air gradient times the geoid height. Removing the Bouguer slab gives a simple Bouguer anomaly of about −93 mGal, typical of the northern Great Plains near the Rocky Mountain front, where the crust is thick.

Measuring gravity

Absolute gravimeters drop a mirrored test mass in a vacuum chamber and track its fall with a laser interferometer referenced to an atomic clock, measuring \(g\) directly to a few microgals. They are expensive and slow, so national networks use them only at a framework of stations. Relative gravimeters, based on a sensitive spring or, in observatories, a superconducting sphere levitated in a magnetic field, measure differences in gravity between stations. Spring meters drift, so surveys are run in loops that start and end on known stations, like leveling loops, and readings are corrected for Earth tides, which change gravity by up to about 0.3 mGal over a day, and for the height of the instrument.

Large areas are covered from the air and from space. The US GRAV-D project measured airborne gravity over the whole country to support the new geoid, and similar airborne surveys have filled gaps in Canada's coverage. The GRACE and GRACE Follow-On missions measure the changing distance between two satellites in the same orbit to map gravity at long wavelengths and its changes over time, while GOCE carried a gradiometer that measured the gravity gradients directly. Combined global models such as EGM2008, and the newer models that support NAPGD2022, express the field as spherical harmonic coefficients up to very high degree, resolving features about 10 km across.

Gravity also changes with time. Removal of groundwater, melting of glaciers and the slow rebound of the crust after the last ice age all change both heights and gravity. Repeated absolute gravity measurements, together with GNSS, separate vertical motion of the crust from changes in mass beneath it. A station that rises by 1 cm without any change in mass beneath it sees gravity drop by about 3 microgals.