Geophysics & Geodynamics

In 1672 Jean Richer shipped a pendulum clock from Paris to Cayenne, five degrees north of the equator, to observe Mars at opposition and pin down the solar parallax. The clock had been regulated in Paris to beat seconds. In Cayenne it ran slow by two minutes and twenty-eight seconds a day, every day. Richer shortened the pendulum by one and a quarter lignes, about 2.8 millimetres out of 994, and the clock kept time. He recorded the correction, sailed home, and was told he had been careless.

He had not. Those 2.8 millimetres are the first quantitative measurement of the shape of the Earth, and this lesson spends them.

What a shortened pendulum measures

A seconds pendulum swings from one side to the other in one second, so its full period is two seconds. For small amplitude, T = 2 pi sqrt(L/g), and inverting at T = 2 s,

L = g T2/(4 pi2) = g x 0.101321 metres per unit of g in m/s².

Gravity at the Paris Observatory is 9.8094 m/s², asking L = 0.99386 m; at Cayenne it is 9.7807, asking L = 0.99095 m. The pendulum must lose 2.91 millimetres. The Paris ligne is a one hundred and forty-fourth of the pied du roi, 2.256 mm, so that is 1.29 lignes. Richer removed 1.25, within four percent of a difference of three parts in a thousand, in 1672, with instruments shipped into a climate that swelled every joint of the clock case.

His other number disagrees. A clock set to the Paris length beats as sqrt(g/L), so a shortfall dg/g costs (1/2)(dg/g) of its ticks: 0.5 x 2.93 x 10-3 x 86400 = 127 seconds a day. Richer reported 148. Trust the length, a null adjustment checked against local solar time day after day; the rate is only as good as a regulation carried out of Paris months before.

The point: gravity varies across the Earth by half a percent from equator to pole, and a seventeenth-century clock could resolve it.

Newton's number, and Huygens's

Fifteen years later Isaac Newton explained it in Book III, Proposition 19 of the Principia. Imagine two canals of water bored from the surface to the centre, one from the pole and one from the equator, meeting there. The Earth is fluid, so the two columns must balance. The equatorial column is lighter, because rotation removes part of its weight, so it must be longer. The Earth bulges.

The size of the effect is set by one dimensionless number, the ratio of centrifugal acceleration at the equator to gravity there:

m = omega2 a / ge.

Put in the sidereal rotation rate omega = 2 pi / 86164.1 s = 7.2921 x 10-5 rad/s, the equatorial radius a = 6 378 137 m, and ge = 9.7803 m/s²:

omega2 a = (7.2921 x 10-5)2 x 6.378137 x 106 = 0.033916 m/s², so m = 0.033916/9.7803 = 3.468 x 10-3 = 1/288.

Newton used 1/289. For uniform density the canal balance gives a flattening of f = (5/4) m, and 1.25/288 = 1/230.7. Newton wrote 1/230, and predicted from the same calculation that a seconds pendulum at the equator must fall short of the Paris length by roughly what Richer had measured. Richer stopped being an embarrassment and became evidence.

Why 5/4 rather than 1? The bulge is rock, and rock attracts: the new equatorial mass pulls the equator further out, which deforms it more, and the series converges on five quarters. Christiaan Huygens worked the same problem in 1690 with all the mass at the centre, so that the bulge attracted nothing, and got f = m/2 = 1/577. Same rotation rate, answers a factor of 2.5 apart, and the whole difference is an assumption about where the mass sits.

The measured value, from the World Geodetic System ellipsoid of 1984, is f = 1/298.257223563 = 3.3528 x 10-3, between Huygens and Newton and closer to Newton.

The Cassini objection, and why anyone had to sail anywhere

Instead Giovanni Domenico Cassini and his son Jacques Cassini extended the Paris meridian survey and in 1718 reported a degree of latitude longer in the south of France than in the north: the signature of a prolate Earth, a lemon rather than an orange.

The logic inverts most intuitions. On an oblate Earth the surface is flattest near the poles, so you must travel further along it to swing the local vertical through one degree. A degree of latitude gets longer toward the pole. To first order,

M(phi) = M0 (1 + 3 f sin2 phi),

where M0 is a degree at the equator. The quarrel came down to a few hundred metres in a quantity of about 111 kilometres, measured by triangulation over ground that was neither flat nor cooperative.

Torneå and Quito

In 1735 the Academy sent two expeditions to arcs far enough apart that the difference could not hide. Pierre Bouguer, Charles Marie de La Condamine and Louis Godin went to the Viceroyalty of Peru, to the plain south of Quito, spent nine years there, quarrelled bitterly and buried one of their own. Pierre Louis Maupertuis took the Lapland party up the frozen Torne river in 1736 and was home within a year.

One toise is 1.9490 m. Do what the Academy did: take the Lapland and Peru degrees, which bracket the whole range of latitude, and solve M(phi)/M0 = 1 + 3 f sin2 phi for f:

57438/56753 = 1.01207, and sin2(66.33 deg) = 0.8393, so f = 0.01207/(3 x 0.8393) = 4.79 x 10-3 = 1/209.

Oblate, decisively, and close enough to Newton's 1/230 to end the argument. The value is too large because the Lapland degree ran about 440 metres long; Jöns Svanberg remeasured that arc in 1801 to 1803 and shortened it. No error of a few hundred metres could turn 1/209 negative, so the prolate Earth was finished. Notice how good the Peru arc was: 110 614 m against a modern 110 575 m, four parts in ten thousand, with wooden quadrants at 3 000 metres altitude.

Clairaut's theorem: from shape to gravity and back

Alexis Clairaut, who had gone to Lapland with Maupertuis, published in 1743 the result that ties the two measurable quantities together. For a rotating body in hydrostatic equilibrium, whatever its internal density profile,

f + beta = (5/2) m, where gravity varies as g(phi) = ge(1 + beta sin2 phi).

This is Clairaut's theorem: a surveyed quantity plus a pendulum quantity add to a fixed multiple of the rotation ratio, the unknown interior cancelling out. Test it with m = 3.4677 x 10-3 and f = 3.3528 x 10-3:

beta = (5/2)(3.4677 x 10-3) - 3.3528 x 10-3 = 8.6693 x 10-3 - 3.3528 x 10-3 = 5.3165 x 10-3.

The measured value, from the WGS84 polar and equatorial gravity values 9.83218 and 9.78033 m/s², is (9.83218 - 9.78033)/9.78033 = 5.3024 x 10-3. Clairaut's first-order prediction is high by 0.27 percent, which is what you expect from a theorem carried to first order in a quantity of size 3 x 10^-3.

Modern practice replaces the two-term formula with the normal gravity of a reference ellipsoid. The 1980 International Gravity Formula reads

gamma(phi) = 9.780327 (1 + 0.0053024 sin2phi - 0.0000058 sin22phi) m/s².

At Cayenne, latitude 4.94 N, it returns 9.78071 m/s²; at the Paris Observatory, 48.84 N, 9.80967. The difference of 0.02896 m/s², through L = 0.101321 g, is a pendulum 2.93 mm shorter. Richer's 1.25 lignes was 2.82 mm.

Why this matters: weaker equatorial gravity has two causes and you need both. Rotation removes m = 3.47 x 10-3 directly. The other 1.83 x 10-3 is shape: the equator sits 21.4 km further from the centre, which weakens gravity, partly offset by the mass of the bulge beneath you, which strengthens it.

What the flattening says about the inside

A homogeneous fluid Earth would be flattened by 1.25 m; the Earth is flattened by 0.967 m. It resists deformation more than uniform rock would, and there is only one way to do that while staying fluid: put the mass in the middle, where rotation is slow and the equatorial bulge gains nothing.

The Radau-Darwin relation makes this quantitative, converting f and m into the moment of inertia factor C/(M a2):

C/(M a2) = (2/3)[1 - (2/5) sqrt(5m/(2f) - 1)].

Work it. 5m/(2f) = 5(3.4677 x 10-3)/(2 x 3.3528 x 10-3) = 0.0173385/0.0067056 = 2.5857. Subtract one: 1.5857. Its square root is 1.2593. Two fifths of that is 0.50372. So C/(M a2) = (2/3)(1 - 0.50372) = 0.3308.

A uniform sphere gives exactly 0.4; all the mass at the centre approaches 0. The Earth's 0.3308, from a rotation rate and the shape of the surface alone, says the interior is far denser than the outside: mean density 5514 kg/m³ against 2700 to 3000 for surface rocks. The rest of this course says what is down there, and where its boundaries lie.

Common misconceptions

• "Centrifugal force throws the equator outward, and that is the whole story." It supplies only m = 1/288. A uniform fluid Earth reaches 1.25 m, because the bulge attracts the bulge, and that self-gravitation is the quarter Huygens left out.
• "Newton's 1/230 was a mistake." It is the right answer for the problem he posed, a homogeneous fluid. The Earth's 1/298 differs because the Earth is not homogeneous, which makes the gap a measurement of the interior rather than an error.
• "A degree of latitude gets shorter toward the pole, because the meridians converge." Converging meridians shorten a degree of longitude. A degree of latitude follows the curvature of the meridian, which is flattest at the poles, so it lengthens: 110.57 km at the equator, 111.69 km at the pole.
• "The Earth is a sphere to any reasonable precision." The equatorial radius exceeds the polar by 21.4 km, two and a half Everests, which is why everything here is quoted against an ellipsoid.

Looking back

Richer's 2.82 mm of pendulum records three parts in a thousand of gravity and agrees with the modern normal gravity formula to four percent. From m = 1/288 Newton got (5/4)m = 1/230 for a homogeneous fluid, Huygens m/2 = 1/577, and the Earth's 1/298 lies between. Lapland and Peru gave f = 1/209 and killed the prolate Earth. Clairaut predicts beta = 5.3165 x 10-3 against a measured 5.3024 x 10-3, and Radau-Darwin returns a moment of inertia factor of 0.3308, well below the 0.4 of uniform rock: the first evidence here for a dense core.

Next, gravity comes off the ellipsoid and onto the real surface, where the corrections you must apply before a reading means anything are larger than the signal.

Sources

1. Wikipedia contributors. (n.d.). Figure of the Earth. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). French Geodesic Mission to Lapland. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Theoretical gravity. Wikipedia. en.wikipedia.org
4. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 5. Cambridge University Press.
5. Lowrie, W., and Fichtner, A. (2020). Fundamentals of geophysics (3rd ed.), Chapter 2. Cambridge University Press.

A gravity meter sitting on a survey mark west of Denver reads 979 639.1 milligals. That is the whole observation. It is also, on its own, worthless: it tells you that you are at 39.74 degrees north on a continent at 1609 metres, and almost nothing about the rock below. This lesson takes that single number and subtracts from it, one term at a time, everything you already know about, until what remains is a statement about the crust.

A milligal is 10^-5 m/s², one millionth of surface gravity. Modern relative gravimeters resolve about 0.01 mGal, so every correction below has to be carried to a hundredth of a milligal or it will swamp the signal.

Correction one: what the ellipsoid already predicts

The first thing to remove is everything Lesson 1 explained. Normal gravity on the 1980 reference ellipsoid at latitude phi is

gamma = 978 032.7 (1 + 0.0053024 sin2phi - 0.0000058 sin22phi) mGal.

At 39.74 N, sin phi = 0.63944 so sin2phi = 0.40888, and sin22phi = 0.96684. The bracket is 1 + 0.0021681 - 0.0000056 = 1.0021625, so

gamma = 978 032.7 x 1.0021625 = 980 147.7 mGal.

Subtract: 979 639.1 - 980 147.7 = -508.6 mGal. The station is short of the ellipsoid prediction by half a gal, and almost all of that deficit is simply the fact that the meter is a mile above the surface the formula refers to.

Correction two: the free-air term

Gravity falls with distance from the centre. Differentiate g = GM/r2 and you get dg/dr = -2g/r, which at g = 9.81 m/s² and r = 6 371 km is -3.08 x 10-6 s^-2, or -0.3086 mGal per metre once the ellipsoid's latitude dependence is included. The free-air correction adds that back:

FAC = 0.3086 x 1609 = 496.6 mGal.

So the free-air anomaly is -508.6 + 496.6 = -12.0 mGal.

Notice what the correction did and did not do. It moved the station down to the datum through empty space. It said nothing about the 1609 metres of granite and sandstone that are actually between the meter and sea level, and which are pulling upward on the meter the whole time. That rock is the next correction, and it is fifteen times larger than the anomaly you are hunting.

Correction three: the plate, and then the hills

Model the intervening rock as an infinite horizontal slab of thickness h and density rho. The attraction of such a slab is famously independent of how far above it you are:

gplate = 2 pi G rho h.

With G = 6.674 x 10-11, 2 pi G = 4.1934 x 10-10, and the standard crustal density rho = 2670 kg/m³, the coefficient is 1.1196 x 10-6 s^-2 per metre, that is 0.11196 mGal per metre. So

gplate = 0.11196 x 1609 = 180.1 mGal,

and the simple Bouguer anomaly is -12.0 - 180.1 = -192.1 mGal.

The slab is a lie in one specific direction. Real topography has valleys where the slab put rock, and those missing masses were pulling the meter sideways and slightly down; it also has peaks above the station whose mass pulls the meter upward. Both errors have the same sign. Whether the ground rises or falls away from you, the terrain correction is always positive. For a station on the eastern flank of the Front Range a typical value is +6.4 mGal, giving a complete Bouguer anomaly of -185.7 mGal.

What matters here: a gravity anomaly is not an anomaly in gravity. It is the residual after a model has been subtracted, and it is only as meaningful as the model. Change rho from 2670 to 2800 and the Bouguer anomaly at this station moves by 8.8 mGal without anything underground having changed.

Reading the two anomalies against each other

The pair carries more information than either alone. Free-air near zero, Bouguer strongly negative, is the signature of a highland in isostatic balance: the extra topographic mass above the datum is cancelled, at long wavelength, by a deficit of mass below it. Something light is filling space that denser mantle would otherwise occupy.

Run the arithmetic backwards to see how much. A Bouguer deficit of 185.7 mGal, read as a slab of crust replacing mantle with a density contrast of 3300 - 2800 = 500 kg/m³, needs a thickness

t = 185.7 / (0.04193 x 0.500) = 185.7 / 0.02097 = 8.9 km,

using the milligal form 2 pi G rho = 0.04193 rho with rho in g/cm³. So this station sits over roughly nine kilometres of extra crustal root beyond a reference column. That is the number the next lesson derives from first principles and checks against the Himalaya.

The contrast with an uncompensated feature is sharp. A buried salt dome or a young seamount that the lithosphere is holding up by its own strength produces a large free-air anomaly, because nothing offsets it at depth. Free-air anomalies over the world's ocean trenches run to -250 mGal and over the outer rise to +50, and those are not compensated: they are a plate being bent.

The geoid, and why your handheld altitude is wrong

All of the above quietly assumed a datum. Which surface is "sea level" under Denver? The answer is the geoid: the equipotential surface of the Earth's gravity field that coincides with mean sea level over the oceans, continued under the continents. It is not the ellipsoid. Mass excesses pull it up, deficits let it sag, and the departures reach +85 m near New Guinea and -106 m in the Indian Ocean south of Sri Lanka.

The link between the two is Bruns's formula: if T is the disturbing potential, the difference between the real potential and the ellipsoid's at the same point, then the geoid height is N = T/gamma. A disturbing potential of 500 m²/s² lifts the geoid by 500/9.81 = 51 m. Potential differences of a few hundred SI units are tens of metres of surface.

This is not bookkeeping. A GNSS receiver measures height above the ellipsoid, h. Levelling, and water, care about height above the geoid, H. They differ by N:

H = h - N.

At the Denver benchmark a receiver reads h = 1662.3 m and the national geoid model gives N = -17.4 m, so H = 1662.3 - (-17.4) = 1679.7 m. Seventeen metres is the difference between a canal that drains and a canal that does not.

The upshot: there are three surfaces in play, and confusing any two of them costs tens of metres. The ellipsoid is a smooth mathematical fiction, the geoid is a physical equipotential you can measure, and the topography is where you are standing.

Chimborazo: the mountain that did not pull

Return to Peru in 1738. Pierre Bouguer, waiting out the survey, set up on the flank of Chimborazo, a volcano rising 6263 m and standing well clear of its surroundings. His idea was Newton's: a mountain has mass, mass attracts, so a plumb line hung beside it should tilt toward it, and the tilt shows up as a difference between astronomical latitude and the latitude the survey gives.

He estimated the mountain's attraction from its shape and an assumed density and predicted a deflection near 1 arcminute 43 arcseconds. What he measured was about 8 arcseconds, roughly one twelfth of it, on a mountain whose mass he could see. Bouguer, who had no framework for a deficit at depth, concluded instead that the Earth as a whole must be far denser than the rock around him, and derived a mean density of about 4.7 times water. That is not a bad number: the modern value is 5.51, and Nevil Maskelyne's better-controlled repetition at Schiehallion in Scotland in 1774 gave 4.5.

But the shortfall was real and it did not go away. It is the same shortfall as the 185.7 mGal above, seen with an eighteenth-century instrument: mountains sit on roots of light material, and from outside the two nearly cancel. Bouguer had measured isostasy 117 years before anyone named it. He simply had no reason to look downward for the missing pull.

Common misconceptions

• "The free-air correction removes the rock between the station and sea level." It removes only the distance. That is precisely why the Bouguer plate exists as a separate term, and at Denver the plate is 180 mGal against a 12 mGal anomaly.
• "The terrain correction can have either sign." It is always positive. A valley is missing mass that would have pulled the meter down, and a peak above the station pulls it up; both make the observed gravity too small relative to the slab model.
• "The geoid is mean sea level." Close, but the sea surface departs from the geoid by up to about two metres because of currents, winds and temperature, a field called dynamic ocean topography that satellite altimetry maps separately.
• "A negative Bouguer anomaly means there is less mass down there." Less than the model assumed, which was rock of density 2670 continuing to the datum. State the reference density or the number means nothing.

What to carry forward

An observed 979 639.1 mGal became a complete Bouguer anomaly of -185.7 mGal through four steps: normal gravity 980 147.7 removed, a free-air correction of 0.3086 mGal per metre restored over 1609 m, a Bouguer plate of 0.11196 mGal per metre removed, and a terrain correction of 6.4 mGal added back. A free-air anomaly near zero over a strongly negative Bouguer anomaly says the highland is compensated, and the arithmetic converts 185.7 mGal into about 9 km of crustal root. The geoid, an equipotential surface running from -106 m to +85 m against the ellipsoid, is the datum all of this refers to, and H = h - N is the sentence that keeps a GNSS height honest. Bouguer's plumb line at Chimborazo, predicted at 103 arcseconds and observed near 8, is the same result with no instrument more modern than a telescope and a weight on a string.

The next lesson asks what the missing mass actually is, and finds two rival answers with the same surface signature.

Sources

1. Wikipedia contributors. (n.d.). Bouguer anomaly. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Geoid. Wikipedia. en.wikipedia.org
3. National Geodetic Survey. (n.d.). GEOID model homepage. NOAA. geodesy.noaa.gov
4. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 5. Cambridge University Press.
5. Lowrie, W., and Fichtner, A. (2020). Fundamentals of geophysics (3rd ed.), Chapter 2. Cambridge University Press.

Between 1847 and 1855 the Great Trigonometrical Survey of India measured the distance from Kaliana, on the plain near Delhi, to Kalianpur, about 600 kilometres south. Triangulation gave the difference in latitude as 5 degrees 23 minutes 42.29 seconds. Astronomical observation, taken by pointing a telescope at stars and reading the plumb line, gave 5 degrees 23 minutes 37.06 seconds. The two disagreed by 5.236 arcseconds, which at 30.9 metres per arcsecond of latitude is 162 metres of ground.

That gap was not error. The survey's own repeat measurements were good to a fraction of an arcsecond, and 5.236 is an order of magnitude larger. Something was tilting the plumb line, and the obvious candidate stood on the northern horizon.

Working the discrepancy the wrong way first

Archdeacon John Henry Pratt did the calculation properly in 1855. He took the mass of the Himalaya and the Tibetan Plateau as their visible volume times ordinary rock density, integrated the horizontal attraction at each station, and predicted that the plumb line at Kaliana should be dragged north relative to Kalianpur by 15.885 arcseconds.

The observed figure was 5.236. The mountains pull about one third as hard as the mountains are there. This is the same failure Bouguer met at Chimborazo in the previous lesson, but now with an instrument good enough that the shortfall could not be blamed on the instrument.

Two responses arrived in the same year, both to the Royal Society, and they are the two models we still use.

Airy: the mountain has a root

Airy's proposal treats the crust as a raft of uniform density floating on a denser fluid mantle. A block that stands higher must also reach deeper, exactly as a thicker iceberg both rises further above the water and sinks further below. Balance the pressure at a depth below the deepest root, where every column must weigh the same, and for a block of height h above the reference surface the root r satisfies

rhoc (h + t + r) = rhoc t + rhom r, so r = h rhoc/(rhom - rhoc),

with t the reference crustal thickness, which cancels. The whole model is that one fraction. Put in a crustal density of 2800 kg/m³ and a mantle density of 3300:

r = h x 2800/500 = 5.6 h.

The Tibetan Plateau stands at a mean elevation of about 4.96 km, so its root is 5.6 x 4.96 = 27.8 km. Adding a 35 km reference crust, the Moho should lie 62.8 km below sea level and 67.8 km below the plateau surface.

Key idea: the root is not extra mass. It is light crust occupying space that dense mantle would otherwise fill, so its gravitational effect is a deficit, and from far away that deficit nearly cancels the visible mountain.

Pratt: the mountain is made of lighter stuff

Pratt's own model keeps the bottom of the crust flat, at a fixed depth of compensation D, and lets density vary sideways. Columns that stand higher are less dense, by exactly the amount that keeps the mass per unit area constant:

rhoh (D + h) = rho0 D, so rhoh = rho0 D/(D + h).

Take D = 100 km and a reference density of 2800 kg/m³. Under the same 4.96 km plateau,

rhoh = 2800 x 100/104.96 = 2668 kg/m³.

A reduction of 132 kg/m³, under five percent, is enough. That is the trouble with the two models: at the surface they are nearly indistinguishable, because both are constructed to hold the total mass in a column fixed, and gravity outside a compensated load mostly sees the total.

The measurement that settled it

Gravity alone will not choose between the two, but a seismic wave will, because a root is a boundary and a density gradient is not. Receiver-function and refraction studies across Tibet find the Moho at 65 to 80 km depth, deepening northward beneath the plateau, against 35 to 40 km on the Indian shield to the south. That is an Airy root, and a slightly deeper one than 5.6 h predicts, because Tibet has also been thickened by tectonic shortening rather than merely floated up.

Pratt is not wrong; it is right somewhere else. Oceanic lithosphere cools as it ages, contracts, and becomes denser, so the ocean floor deepens as the square root of age. That is lateral density variation above a compensation depth, which is Pratt's model with temperature doing the work. A mid-ocean ridge stands 3 km above the abyssal plain for the same reason a Pratt column stands high: it is hot and light, not thick.

Real crust does neither purely. A seamount 60 km across sits on a lithosphere with elastic strength, and the load is carried regionally by flexure of a plate rather than locally by a root. That is Vening Meinesz's amendment, and it is why small loads show large free-air anomalies while wide ones do not: compensation is efficient at wavelengths long compared with the flexural parameter, and absent at short ones.

Fennoscandia: isostasy caught in the act

Everything above is a statement about equilibrium. The best evidence that the mantle really does flow to reach it comes from a place where equilibrium has not yet arrived.

At the last glacial maximum the Fennoscandian ice sheet was roughly 3 km thick over the Gulf of Bothnia. The equilibrium depression under that load is

w = hice rhoice/rhom = 3000 x 917/3300 = 833 m.

The ice was gone by about 10 000 years ago. The land is still coming up. GNSS and levelling give a present maximum rate near 10 mm per year at the head of the Gulf of Bothnia, raised shorelines record several hundred metres of uplift already achieved, and roughly 100 m remain. The free-air anomaly over the region is about -30 mGal, a real mass deficit that has not yet been refilled.

Now extract a number from it. Uplift after an instantaneous unloading decays exponentially, and for a load of horizontal wavelength lambda on a viscous half-space the relaxation time is

tau = 4 pi eta/(rho g lambda).

Fennoscandian shoreline data give tau = 4400 years, that is 1.389 x 1011 s, over a load some 3000 km across. Rearranged,

eta = tau rho g lambda/(4 pi) = (1.389 x 1011)(3300)(9.81)(3.0 x 106)/12.566.

The numerator is 1.349 x 1022, so eta = 1.07 x 1021 Pa s. This is essentially the number Norman Haskell obtained in 1935 from the same data, and it remains the standard reference viscosity for the upper mantle. For scale: water is 10^-3 Pa s, honey about 10, pitch about 10⁸. The mantle is thirteen orders of magnitude stiffer than pitch and still flows fast enough to lift Scandinavia a centimetre a year.

Bottom line: a plumb line that leaned 5.236 arcseconds instead of 15.885 says mountains are compensated; a shoreline rising 10 mm a year says the compensating flow has a measurable viscosity; and one number, 10²¹ Pa s, connects them.

Common misconceptions

• "Isostasy means the Earth is in balance." It means it tends toward balance on a timescale set by mantle viscosity. Fennoscandia is 100 m out of balance right now, and ocean trenches are held tens of kilometres out of it by plate forces.
• "The root under a mountain is dense, which is why it holds the mountain up." The reverse. The root is crust, less dense than the mantle it displaces, and it holds the mountain up by buoyancy exactly as the submerged part of an iceberg does.
• "Airy and Pratt are competing theories and one of them lost." Both mechanisms operate. Continental mountain belts are compensated by thickness, and oceanic lithosphere by thermal density variation, which is Pratt's arithmetic with temperature supplying the density change.
• "Post-glacial rebound is the crust springing back elastically." The elastic part finished as the ice left. What is still happening is viscous flow of mantle material back under the depression, which is why the timescale is thousands of years rather than seconds.

Where this leaves us

A 5.236 arcsecond discrepancy between two Indian survey stations, against Pratt's predicted 15.885, showed that the Himalaya pull about a third as hard as their visible mass requires. Airy explains the shortfall with a root r = h rhoc/(rhom - rhoc), which is 27.8 km under a 4.96 km plateau at densities of 2800 and 3300. Pratt explains it with a density of 2668 instead of 2800 above a 100 km compensation depth. Gravity cannot separate them, because both fix the column mass; seismology can, and the Moho at 65 to 80 km beneath Tibet says Airy. Fennoscandia then turns the idea into a rheology: a 3 km ice sheet depresses the surface by 833 m, the recovery has a relaxation time of 4400 years, and tau = 4 pi eta/(rho g lambda) returns a mantle viscosity of 1.07 x 10²¹ Pa s.

Gravity has now given us a compensated crust and a viscous mantle, but only in outline. To get depths, boundaries and densities we need a probe that goes through the Earth rather than around it, and Module 2 starts with the two kinds of wave that do.

Sources

1. Wikipedia contributors. (n.d.). Isostasy. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Post-glacial rebound. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). John Henry Pratt. Wikipedia. en.wikipedia.org
4. Turcotte, D. L., and Schubert, G. (2014). Geodynamics (3rd ed.), Chapters 2 and 6. Cambridge University Press.
5. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 5. Cambridge University Press.

On 17 April 1889 an earthquake shook Tokyo. Eight thousand kilometres away in Potsdam, Ernst von Rebeur-Paschwitz was using a delicate horizontal pendulum to look for tidal tilting of the ground, and his photographic trace that day carried a disturbance he could not explain. When the Japanese bulletin arrived weeks later he compared times and realised that his instrument had recorded an earthquake on the far side of the planet. Nobody had done that before. The Earth, it turned out, rings.

Everything in this module follows from that: if elastic waves cross the whole planet, their arrival times are a survey of the inside. But there are two kinds of body wave, they carry different information, and one of them stops dead at a depth of 2891 km. This lesson is built around the comparison, because the differences are where the physics lives.

Two waves from two moduli

An elastic solid resists two distinct kinds of deformation. Squeeze it and it pushes back: that stiffness is the bulk modulus K, the ratio of pressure change to fractional volume change. Shear it, sliding one face past another without changing volume, and it also pushes back: that is the shear modulus mu. Each restoring force supports its own wave:

VP = sqrt((K + 4mu/3)/rho) and VS = sqrt(mu/rho).

The P wave is a compression, particles moving along the direction of travel, and it needs both moduli because compressing a block in one direction without letting it spread sideways involves shear too. The S wave is a shear, particles moving across the direction of travel, and it needs only mu.

Work three materials. Granite: K = 56 GPa, mu = 28 GPa, rho = 2700 kg/m³. Then K + 4mu/3 = 56 + 37.3 = 93.3 GPa, and

VP = sqrt(9.33 x 1010/2700) = sqrt(3.456 x 107) = 5.88 km/s, VS = sqrt(2.8 x 1010/2700) = 3.22 km/s.

Upper-mantle peridotite: K = 130, mu = 70 GPa, rho = 3300. Then K + 4mu/3 = 223.3 GPa and VP = 8.23 km/s, VS = 4.61 km/s. Those are the numbers that will identify the mantle in the next lesson.

A liquid is the third case, and it is special: a fluid by definition cannot support a static shear stress, so mu = 0 exactly. Then VS = 0 and VP = sqrt(K/rho). There is no shear wave, not a slow one, not an attenuated one. The restoring force that a shear wave exists to oscillate against is absent.

In short: P waves measure K + 4mu/3, S waves measure mu alone, and anywhere mu vanishes the S wave simply does not exist.

Why P always wins the race

Divide the two expressions and the density cancels:

VP/VS = sqrt(K/mu + 4/3).

Since K and mu are both positive, the ratio can never fall below sqrt(4/3) = 1.155, so the P wave is always first, everywhere, in every material. For most crustal rock K/mu is near 2 and the ratio is about 1.73, the square root of three. Granite above gives 5.88/3.22 = 1.826.

The same ratio is usually quoted as Poisson's ratio:

nu = (VP2/VS2 - 2)/(2(VP2/VS2 - 1)).

For granite, VP2/VS2 = 3.334, so nu = 1.334/4.668 = 0.286. As VS goes to zero the ratio goes to infinity and nu goes to 0.5, the fluid limit. Poisson's ratio is therefore a fluid detector: values above about 0.30 in a reservoir often mean pore fluid, and 0.5 means no solid framework at all.

The boundary at 2891 km, read from four numbers

Here is the heart of the lesson. The Preliminary Reference Earth Model tabulates properties on both sides of the core-mantle boundary:

Check the last row, because it is the whole argument. Base of mantle: rho VP2 = 5566 x (13716)2 = 1.047 x 1012 Pa, and 4mu/3 = 392 GPa, so K = 1047 - 392 = 656 GPa. Top of core: mu = 0, so K = rho VP2 = 9903 x (8065)2 = 644 GPa.

The bulk modulus is the same within two percent across the boundary. The density almost doubles. And yet the P velocity falls by 41 percent. Run it the other way to be sure: if the core had the mantle's rigidity as well as its own density, VP = sqrt((644 + 392) x 109/9903) = 10.2 km/s rather than 8.06. Every bit of the missing speed is missing rigidity.

The core of it: the core-mantle boundary is not a place where the Earth becomes less stiff to compression. It is a place where it stops resisting shear, and the S wave simply has nothing left to propagate on.

Four waves, four jobs

Body waves go through; surface waves run along the boundary and decay with depth. Each tells you something the others cannot.

Surface waves are dispersive: long periods sample deeper, faster material and therefore arrive first, spreading a single pulse into a train that can last many minutes. That dispersion is not a nuisance. Measuring group velocity as a function of period and inverting gives a shear velocity profile with depth, which is how most of what we know about the upper 300 km of the mantle under the oceans was obtained.

Using the gap between P and S

Because the two waves leave together and travel at different speeds, their separation grows with distance. For near-surface paths with VP = 8.0 and VS = 4.6 km/s,

1/VS - 1/VP = 0.2174 - 0.1250 = 0.0924 s per km,

so distance = (tS - tP)/0.0924 = 10.8 (tS - tP) km. A station recording S 60 seconds after P is about 650 km from the epicentre. Three such circles intersect at one point, which is how every earthquake was located before computers, and how the arithmetic inside a modern locator still begins.

The method also fails in an instructive way. It assumes a single constant velocity, and real velocity climbs with depth, so rays bottom out and arrive early. Correcting that is exactly the travel-time curve of the next lesson.

Common misconceptions

• "S waves cannot cross a liquid because liquids are less dense." Density is irrelevant, and the outer core is in fact nearly twice as dense as the mantle above it. The reason is mu = 0: a liquid has no shear restoring force, so a shear wave has nothing to oscillate against.
• "The P velocity drops at the core because the core is softer." Its bulk modulus is 644 GPa against the mantle's 656, a two percent difference. The drop from 13.72 to 8.06 km/s is the loss of mu and the near doubling of rho.
• "P stands for powerful." It stands for primary, meaning first. On most seismograms S and the surface waves are considerably larger in amplitude, which is why damage correlates with them and not with the P arrival.
• "Seismic rays travel in straight lines." Velocity increases with depth almost everywhere, so rays curve continuously and turn back toward the surface. A ray that reaches 25 degrees of distance bottoms near 600 km depth rather than travelling a chord.

Putting it together

Two elastic moduli give two body waves: VP = sqrt((K + 4mu/3)/rho) and VS = sqrt(mu/rho). Granite at K = 56 and mu = 28 GPa gives 5.88 and 3.22 km/s; peridotite at 130 and 70 GPa gives 8.23 and 4.61. The ratio sqrt(K/mu + 4/3) can never fall below 1.155, so P is always first. At the core-mantle boundary PREM gives densities of 5566 and 9903, P speeds of 13.72 and 8.06, and bulk moduli of 656 and 644 GPa: the stiffness to compression barely changes, the density nearly doubles, and the entire velocity drop is the disappearance of a 294 GPa shear modulus. That is what the absence of shear waves in the outer core proves, and it is the strongest single statement in geophysics about the state of matter 3000 km beneath your feet. Finally, 10.8 (tS - tP) turns an arrival interval into kilometres, approximately.

Next, that approximation gets repaired, and the repair is what located the Moho, the core and the inner core.

Sources

1. Wikipedia contributors. (n.d.). S wave. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Preliminary reference Earth model. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Core-mantle boundary. Wikipedia. en.wikipedia.org
4. Stein, S., and Wysession, M. (2003). An introduction to seismology, earthquakes and Earth structure, Chapters 2 and 3. Blackwell.
5. Dziewonski, A. M., and Anderson, D. L. (1981). Preliminary reference Earth model. Physics of the Earth and Planetary Interiors, 25(4), 297-356.

At 09:03 local time on 8 October 1909 a magnitude 6 earthquake struck the Kupa Valley, about forty kilometres southeast of Zagreb. Andrija Mohorovicic, who ran the Zagreb observatory, wrote to every station in Europe that might have felt it and assembled arrivals from twenty-nine of them. Plotting arrival time against distance, he found something the single-velocity picture forbids: beyond a few hundred kilometres there were two P arrivals, and the later one at short range became the earlier one at long range.

That crossing is the subject of this lesson, and the same trick, applied at planetary scale, produces the core.

Two paths, two straight lines

Suppose a crust of thickness h and speed V₁ lies on a mantle of speed V₂, with V₂ greater. A wave leaving a shallow source can go directly through the crust, or it can strike the interface at the critical angle, run along the top of the faster layer, and leak energy back up. The two arrival times are

tdirect = x/V1 and thead = x/V2 + 2h cos(ic)/V1, with sin(ic) = V1/V2.

The first is a line through the origin with slope 1/V1. The second is a shallower line with slope 1/V2 and a positive intercept. Two lines of different slope must cross. Setting them equal and solving,

xc = 2h sqrt((V2 + V1)/(V2 - V1)).

Take Mohorovicic's velocities, V1 = 5.6 and V2 = 7.9 km/s, and his crustal thickness of 54 km. Then (7.9 + 5.6)/(7.9 - 5.6) = 13.5/2.3 = 5.870, whose square root is 2.423, so

xc = 2 x 54 x 2.423 = 262 km.

Now run the inversion the way a seismologist actually does it, from the intercept rather than the crossover, because the intercept is measured on a fitted line rather than at a single crossing point. With sin(ic) = 5.6/7.9 = 0.7089, cos(ic) = 0.7054, and the intercept

ti = 2h cos(ic)/V1 = 2 x 54 x 0.7054/5.6 = 13.60 s.

So a measured intercept of 13.60 s returns h = ti V1/(2 cos ic) = 13.60 x 5.6/1.411 = 54.0 km. That layer boundary is the Mohorovicic discontinuity, universally shortened to the Moho, and it is the base of the crust everywhere on Earth.

Two honest caveats. Modern refraction and receiver-function work under the Dinarides puts the Moho nearer 40 km than 54, because Mohorovicic had to assume constant velocities in each layer and the real crust speeds up with depth. And the method is blind to a low-velocity layer: a slow zone under a fast one generates no head wave at all, so it leaves no trace on the curve. Refraction sees only boundaries where speed jumps upward.

Remember: the observable is a slope and an intercept. Slopes give velocities, the intercept gives a depth, and everything else in this lesson is that idea applied at a larger radius.

From slope to turning point

Over the whole Earth the travel-time curve is not straight, because velocity rises with depth and rays bend continuously. The governing quantity is the ray parameter

p = r sin(i)/v(r),

which is constant along a ray in a spherically symmetric Earth. At the deepest point of a ray, the turning point, the ray is horizontal, so sin(i) = 1 and

p = rturn/v(rturn).

The reason this is useful is that p is also the slope of the travel-time curve, p = dT/dDelta, which you can measure with a ruler. Read the slope of the P curve near 100 degrees and you get about 4.43 seconds per degree. Convert to radians: 4.43 x 57.296 = 253.8 s.

That ray bottoms just above the core, where the mantle P velocity is 13.72 km/s. So

rturn = p v = 253.8 x 13.72 = 3482 km,

and the depth is 6371 - 3482 = 2889 km. The modern value for the core-mantle boundary is 2891 km. A ruler on a graph, one velocity, and one line of algebra put the largest boundary inside the Earth within two kilometres.

The dark ring

Beyond about 103 degrees the direct P wave stops arriving. Between 103 and 143 degrees there is a shadow zone, and past 143 degrees P returns, late and strong, as the phase PKP that has passed through the core.

The cause is the velocity drop worked out in the previous lesson: 13.72 km/s falls to 8.06 at the boundary. Snell's law at a drop in velocity bends a ray toward the normal, so rays entering the core are deflected sharply downward and emerge far around the far side, skipping a band of the surface entirely.

Try to get the core radius from the shadow edge with straight rays and watch it fail, because the failure is instructive. A straight chord from the source tangent to a sphere of radius r_c emerges at distance Delta with cos(Delta/2) = rc/R. At Delta = 103 degrees, cos(51.5 deg) = 0.6225, giving rc = 0.6225 x 6371 = 3966 km. That is 14 percent too large. Real rays are concave upward because velocity increases with depth, so they reach a given distance from a shallower turning point than a straight line would; assuming straightness therefore puts the obstacle too far out. The ray-parameter method above has no such flaw, which is why it lands on 3482 km.

S waves make the same point more bluntly. There is no S beyond 103 degrees at all, anywhere, ever. Not a delayed arrival, not a weak one. Richard Dixon Oldham published the argument in 1906 from records of the 1897 Assam earthquake, and Beno Gutenberg nailed the depth to 2900 km in 1913 by timing reflections and diffractions around the boundary.

Lehmann's P prime

The shadow should be dark. It was not quite.

Inge Lehmann, running the Danish seismic network, kept finding faint P arrivals inside the shadow zone, well before PKP could possibly appear. Her clearest data came from the Buller earthquake in New Zealand on 16 June 1929, recorded across Europe. In 1936 she published a three-page paper with the shortest title in the literature, P prime, proposing that the core has an inner core with a higher velocity, whose boundary reflects and refracts energy into the forbidden band.

Her model put that boundary near 1400 km radius. Normal-mode and body-wave work later refined it to 1221 km, and PREM gives the inner core a P velocity of 11.03 km/s against 10.36 at the base of the outer core: a jump upward, exactly the opposite of the drop at the core-mantle boundary, which is why the inner core focuses energy rather than deflecting it away.

The inner core is solid. The evidence is the phase PKJKP, an S wave crossing the inner core, which is extremely weak and was only convincingly reported decades later; more robustly, the Earth's free oscillations after very large earthquakes have periods that require a finite shear modulus at the centre. PREM assigns the top of the inner core a shear velocity of 3.50 km/s, which corresponds to mu = 12 764 x 35002 = 156 GPa.

So what?: every one of these is the same procedure. Measure arrival time against distance, find where the curve does something a uniform Earth cannot do, and turn that feature into a radius.

Common misconceptions

• "The shadow zone is dark because the core absorbs the energy." It is a geometric effect. Snell's law at a velocity drop from 13.72 to 8.06 km/s bends rays steeply downward, so they emerge beyond 143 degrees instead of within the band.
• "The S shadow and the P shadow have the same cause." They do not. P is merely refracted away from 103 to 143 degrees and comes back as PKP. S is absent beyond 103 degrees permanently, because it cannot enter a liquid at all.
• "Refraction surveys find every layer." They find only boundaries where velocity increases downward. A low-velocity zone generates no head wave and is invisible, which is why Mohorovicic's method could not have found the asthenosphere.
• "Mohorovicic measured the crust to be 54 km thick, so it is." That number depends on assuming two uniform layers. Modern work under the same region gives about 40 km. The discontinuity is real; the depth was model-dependent, and saying so is part of quoting it.

The short version

A refracted arrival that overtakes the direct one gives a crossover xc = 2h sqrt((V2 + V1)/(V2 - V1)), which for 5.6 and 7.9 km/s over a 54 km crust falls at 262 km, and an intercept of 13.60 s returns the same 54 km. For the deep Earth the ray parameter p = rturn/v(rturn) equals the measured slope dT/dDelta: 4.43 s per degree is 253.8 s per radian, and multiplying by the 13.72 km/s at the base of the mantle gives a turning radius of 3482 km, or a boundary 2889 km down. Straight rays would have said 3966 km, which is 14 percent wrong and shows why curvature cannot be waved away. S vanishes past 103 degrees because the outer core is liquid; P reappears past 143 as PKP; and the faint arrivals Lehmann found in between require an inner core, now placed at 1221 km radius with a shear modulus of about 156 GPa.

Velocities are not the goal, though. Density is, and the next lesson converts one into the other.

Sources

1. Wikipedia contributors. (n.d.). Andrija Mohorovicic. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Shadow zone. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Inge Lehmann. Wikipedia. en.wikipedia.org
4. Stein, S., and Wysession, M. (2003). An introduction to seismology, earthquakes and Earth structure, Chapter 3. Blackwell.
5. Lowrie, W., and Fichtner, A. (2020). Fundamentals of geophysics (3rd ed.), Chapter 8. Cambridge University Press.

In 1923 Leason Adams and Erskine Williamson, working at the Carnegie Institution's Geophysical Laboratory in Washington, published a way to turn seismic velocities into densities. Their assumption was simple enough to state in one sentence: suppose a region of the Earth is chemically uniform, at the same temperature gradient a rising parcel would follow, and compressed only by the weight of what lies above. Then density and velocity are locked together, and a curve of velocity against depth becomes a curve of density against depth.

The assumption is false in about a fifth of the mantle. Finding out exactly where it fails turned out to be more useful than the places where it works, and this lesson follows that trail.

Building the equation

Start with hydrostatic equilibrium. In a shell at radius r, the pressure gradient balances weight:

dP/dr = -rho g(r), with g(r) = G m(r)/r2 and m(r) the mass inside r.

Now bring in the elasticity. The bulk modulus is defined by K = rho dP/drho, so drho/dP = rho/K. And from the previous two lessons, K/rho = VP2 - (4/3)VS2, a combination worth its own name, the seismic parameter:

Phi = VP2 - (4/3)VS2 = K/rho.

Chain the two together, drho/dr = (drho/dP)(dP/dr), and the Adams-Williamson equation falls out:

drho/dr = -rho g/Phi.

Everything on the right is measurable. Seismology gives V_P and V_S at every depth, so Phi is known; g follows from the mass already accumulated; and you integrate downward from a starting density at the top of the region. What you must assume is that no other cause of density change is at work: no change in chemistry, no change in crystal structure, no departure from an adiabatic temperature profile.

The lower mantle, where it works

Take the region from 771 km depth, the top of the lower mantle, to the core at 2891 km. Reference values there: at the top, rho = 4443 kg/m³, VP = 11.07, VS = 6.24 km/s, so

Phi = 122.5 - (4/3)(38.94) = 70.6 km²/s².

At the base, VP = 13.72 and VS = 7.265 give Phi = 188.2 - 70.4 = 117.9. Use mid-shell values: Phi = 94.2 km²/s² = 9.42 x 107 m²/s², rho = 5005 kg/m³, and g = 10.2 m/s², since gravity rises from about 9.9 to 10.65 across the lower mantle. Then

drho/dr = -(5005)(10.2)/(9.42 x 107) = -5.42 x 10-4 kg/m³ per metre, that is 0.542 per kilometre of descent.

Over the 2120 km of the lower mantle that predicts a density increase of 0.542 x 2120 = 1149 kg/m³. The observed increase, 4443 to 5566, is 1123. The equation is high by 2.3 percent, which for a one-line model with mid-shell averages is a triumph.

Worth holding on to: the lower mantle behaves like one substance being squeezed. Two thousand kilometres of rock, an eighth of the Earth's radius, and self-compression explains almost all of the density change.

The transition zone, where it collapses

Now do the same arithmetic from 400 km to 700 km depth. At 400 km, rho = 3540, VP = 9.13, VS = 4.93, so Phi = 83.4 - 32.4 = 51.0. At 700 km, VP = 10.75, VS = 5.95, so Phi = 115.6 - 47.2 = 68.4. Mid values: Phi = 5.97 x 107, rho = 3960, g = 9.97.

drho/dr = -(3960)(9.97)/(5.97 x 107) = -6.61 x 10-4, or 0.661 kg/m³ per km.

Over 300 km that is a predicted increase of 198 kg/m³. The observed increase is 3540 to 4380, which is 840.

Self-compression accounts for less than a quarter of it. This is the discrepancy Francis Birch pressed hard in 1952, and it admits only two explanations: either the composition changes at these depths, or the same material rearranges itself into a denser crystal structure. Birch argued for the second, and he was right.

What the something else is

Mantle rock is dominated by olivine, magnesium iron silicate. Under pressure it does not simply compress; it reorganises.

• At about 13.5 GPa, near 410 km, olivine converts to wadsleyite, roughly 6 percent denser in the same chemistry.
• Near 520 km, wadsleyite gives way to ringwoodite, a spinel structure, in a weaker and broader step.
• At about 23.5 GPa, near 660 km, ringwoodite breaks down entirely into bridgmanite plus ferropericlase, a change of about 8 percent in density and the largest single step inside the mantle.

These are not guesses. Each was reproduced in a diamond anvil cell at the right pressure and temperature, and the depths match the seismic discontinuities to within the uncertainty of the mantle geotherm.

The two main boundaries behave differently under temperature, and the difference matters for the rest of this course. Their Clausius-Clapeyron slopes have opposite signs: about +3 MPa/K at 410 and about -2 MPa/K at 660. Work the consequence for a cold slab, 500 K below ambient. The pressure gradient in the mantle is rho g, about 34.7 MPa per km near 410 km.

At 410, a positive slope means colder material transforms at lower pressure: DeltaP = 3 x 500 = 1500 MPa, so the boundary rises by 1500/34.7 = 43 km. At 660, a negative slope pushes the transformation to higher pressure in cold material: DeltaP = 2 x 500 = 1000 MPa, and with rho g = 39.9 MPa/km the boundary is depressed by 1000/39.9 = 25 km. Seismic imaging finds exactly that pattern, an elevated 410 and a depressed 660, over subducting slabs.

The sign at 660 has a dynamical consequence too. A slab arriving there is denser than its surroundings, but the phase boundary it must push down is buoyant relative to it, so the transition resists penetration. Some slabs punch straight through, some flatten and lie horizontally for a thousand kilometres before sinking. That argument returns in Lesson 9.

PREM, and the layer with the odd name

All the numbers above come from the Preliminary Reference Earth Model, published by Adam Dziewonski and Don Anderson in 1981. PREM is a one-dimensional model: density, the two velocities and attenuation as functions of radius alone, fitted simultaneously to body-wave travel times, surface-wave dispersion, free-oscillation periods, and the two integral constraints from Lesson 1, the total mass 5.972 x 10²⁴ kg and the moment of inertia factor 0.3308. It has been the reference against which everything is compared for over forty years.

Near the bottom the model runs out of road. The lowermost 200 to 300 km of the mantle, called D double prime, has velocity gradients unlike the rest of the lower mantle, strong lateral variation, and in places thin patches where V_S drops by 10 to 30 percent, the ultra-low velocity zones. The name is an accident of bookkeeping: Keith Bullen labelled the Earth's shells A through G in the 1940s, region D was the lower mantle, and when it had to be split the pieces became D prime and D double prime.

Part of the explanation arrived in 2004, when Murakami and colleagues found that bridgmanite transforms again, near 125 GPa and 2500 K, into post-perovskite. That depth is right for the top of D double prime, and the transition explains the velocity jump and some of the anisotropy. The rest of the layer's character is thermal: it is the boundary layer above a core more than 1000 K hotter than the mantle, and Module 3 treats it as such.

The point: Adams-Williamson is not a failed model. It is a null hypothesis, and each place it fails is a discovery: 410, 520, 660 and the base of the mantle are all named because self-compression could not account for what the waves reported.

Common misconceptions

• "The 410 and 660 km discontinuities are changes in composition." Both are phase changes in essentially the same material. The mantle above and below 660 is chemically similar; ringwoodite has merely broken down into bridgmanite and ferropericlase.
• "Density is measured directly by seismology." It is not. Body waves constrain velocities, and density enters only through models like Adams-Williamson plus the integral constraints of total mass and moment of inertia. That is why PREM's density is much less well resolved than its velocities.
• "D double prime is the outer core." It is the bottom of the mantle, solid silicate, sitting above the boundary at 2891 km. The odd name comes from Bullen's alphabetical labels, not from any property of the layer.
• "If the equation fails, the seismic data must be wrong." The reverse. The velocities are among the best-determined quantities in the Earth sciences; it was the assumption of a chemically uniform, adiabatic shell that broke, and its breaking is how the transition zone was found.

What you now know

Hydrostatic balance plus the definition of the bulk modulus gives drho/dr = -rho g/Phi with Phi = VP2 - (4/3)VS2. Through the lower mantle, mid-shell values of Phi = 94.2 km²/s², rho = 5005 and g = 10.2 predict a density rise of 1149 kg/m³ against an observed 1123, so 2120 km of rock is explained by compression alone. Across the transition zone the same calculation predicts 198 and the Earth delivers 840, and the shortfall is made up by olivine converting to wadsleyite near 410 km, to ringwoodite near 520, and breaking down to bridgmanite and ferropericlase near 660. The Clapeyron slopes of plus 3 and minus 2 MPa/K raise the 410 by about 43 km and depress the 660 by about 25 km in a slab 500 K cold. PREM ties all of it to a total mass of 5.972 x 10²⁴ kg and a moment of inertia factor of 0.3308, and stops being adequate only in D double prime, where post-perovskite and a hot core-mantle boundary take over.

That boundary layer is a thermal object, so the next module stops asking what the Earth is made of and starts asking how hot it is and what the heat is doing.

Sources

1. Wikipedia contributors. (n.d.). Adams-Williamson equation. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Transition zone (Earth). Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Post-perovskite. Wikipedia. en.wikipedia.org
4. Dziewonski, A. M., and Anderson, D. L. (1981). Preliminary reference Earth model. Physics of the Earth and Planetary Interiors, 25(4), 297-356.
5. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 8. Cambridge University Press.

In 1862 William Thomson, not yet Lord Kelvin, published a paper called "On the secular cooling of the Earth". He had a differential equation, three numbers, and a conclusion: the Earth had been a solid body for something like 100 million years, and no geologist was entitled to more. For the next thirty-five years that number, refined downward to between 20 and 40 million years, was the hardest constraint in the Earth sciences, and it was wrong.

What makes the episode worth a lesson is that the calculation is correct. Every step of it survives. The error is in an assumption so natural that it took an engineer with a grudge to notice it.

The calculation, done properly

Thomson's model is a half-space that starts at a uniform temperature T₀ and has its surface suddenly held at zero. Heat leaves by conduction alone. The solution to the diffusion equation is

T(z,t) = T0 erf(z/(2 sqrt(kappa t))),

where kappa is thermal diffusivity. Differentiate at the surface, where the error function is linear, and the gradient is

(dT/dz)0 = T0/sqrt(pi kappa t).

That gradient is measurable in a mine. Invert for the age:

t = T02/(pi kappa (dT/dz)02).

Put in what Thomson had. The melting temperature of rock, then thought to be about 3870 K above the surface; the gradient in British mines, about 38 K per km, so 0.038 K/m; and a diffusivity of 1.25 x 10-6 m²/s measured on rock samples in his own laboratory. Then

T02 = 1.4977 x 107, pi kappa = 3.927 x 10-6, (dT/dz)2 = 1.444 x 10-3,

so the denominator is 5.671 x 10-9 and

t = 1.4977 x 107/5.671 x 10-9 = 2.641 x 1015 s = 84 million years.

Notice how the answer depends on the inputs. The age goes as the square of the assumed melting temperature and inversely as the square of the gradient, so a 20 percent error in either moves the answer by 40 percent or more. That sensitivity is why Thomson's published range, 20 to 400 million years, was so wide, and why he kept tightening it toward the low end as estimates of rock melting points fell.

How deep the calculation can see

Here is the question nobody asked for thirty years. The cooling wave has penetrated to a depth of roughly sqrt(kappa t):

sqrt(1.25 x 10-6 x 2.641 x 1015) = sqrt(3.30 x 109) = 57 km.

Fifty-seven kilometres, on a planet with a radius of 6371. Thomson's solution says nothing whatever about the other 99 percent of the Earth. It is a statement about a thin skin, and it converts a measured surface gradient into an age only if you assume that the temperature at the bottom of that skin has been falling in exactly the way pure conduction from an initially uniform body requires.

Why this matters: the result is not a measurement of the Earth's age. It is a measurement of how long the outer 57 km has been losing heat on the assumption that nothing resupplies it from below.

The geologists, and then the engineer

Geologists had been saying for decades that the number was too small, and they had evidence, though of a kind Thomson did not weigh heavily. Charles Darwin, estimating the time to strip the chalk from the Weald in the first edition of the Origin in 1859, arrived at about 300 million years for that one erosional episode. Stratigraphers counting sediment thicknesses against modern deposition rates kept producing hundreds of millions of years. Thomson's reply was consistent: rates of erosion and deposition are guesses, whereas conduction is physics.

The objection that finally worked came from John Perry, professor of engineering at Finsbury Technical College and once Thomson's own assistant. In a series of letters in 1894 and a paper in Nature in 1895, Perry made a precise point. Suppose the deep interior is not a rigid conductor but able to move, or simply has a much larger effective diffusivity. Then heat is delivered to the base of the conducting lid faster than diffusion alone would deliver it, the lid's gradient does not decay as the inverse square root of time, and the inferred age is no longer bounded by 100 million years.

The arithmetic is unforgiving. In the half-space solution, t is proportional to 1/kappa for a fixed gradient, but if the interior's effective diffusivity is larger by a factor f while the lid's stays the same, the surface gradient is maintained for a time longer by roughly f. Perry's own numbers: an interior ten times more conductive allows 840 million years; a hundred times allows 8.4 billion. He concluded that ages of two to three billion years were entirely compatible with the same measured gradient.

Thomson conceded that the mathematics was right and denied that the physics applied, on the grounds that the interior is rigid because it transmits shear waves and supports the tides. He was answering the wrong objection. A material can be elastic on a timescale of seconds and flow on a timescale of millions of years, which is exactly the contrast between the shear modulus of 294 GPa and the viscosity of 10²¹ Pa s already met in this course. Perry had no name for that, but he had the consequence right.

The term neither of them had

Radioactivity was discovered in 1896. By 1903 Pierre Curie and Albert Laborde had measured the heat output of radium, and in 1906 John Strutt, later the fourth Baron Rayleigh, measured radium concentrations in common rocks and reached a conclusion that ended the argument in one step: if the whole Earth contained radium at the concentration found in surface rocks, it would produce far more heat than flows out of it. A cooling Earth was not the only source of the gradient, and might not even be the main one.

Here is the modern accounting. Total heat loss, from roughly 38 000 borehole and marine measurements, is 47 terawatts, give or take 2. Spread over 5.10 x 10¹⁴ m² of surface that is a mean flux of 92 mW/m², with continents averaging about 65 and the ocean floor about 96, because oceanic lithosphere is young and still cooling.

The ratio of radiogenic production to total loss, about 20 over 47 or 0.43, is the Urey ratio. It says the Earth is still cooling: it loses rather more than it makes, at a rate that corresponds to the mantle dropping roughly 100 K per billion years.

Work the crustal term yourself, because it is the one you can check. Take continental crust 40 km thick, with heat production of 2.0 microwatts per cubic metre in the upper 10 km, where granites concentrate uranium, thorium and potassium, and 0.4 in the lower 30 km. The contribution to surface heat flow is

(2.0 x 10-6)(10 000) + (0.4 x 10-6)(30 000) = 0.020 + 0.012 = 0.032 W/m² = 32 mW/m².

Against a continental mean of 65, that leaves about 33 mW/m² arriving from the mantle below. Roughly half the heat under your feet on a continent was made in the rock immediately beneath you, and the other half came from 2900 km down. Strutt's result in one line.

Since 2005 the radiogenic term has been measured a second way, independently. Uranium and thorium decay chains emit antineutrinos above the detection threshold, and KamLAND in Japan and Borexino in Italy have counted them. The inferred geoneutrino signal gives a radiogenic power of about 20 TW with an uncertainty still of order 8, consistent with the table and obtained without any assumption about the Earth's composition at all.

In short: Thomson's equation was right, Perry's objection was right, and the decisive term was one that did not exist as a concept when either of them started.

Who was actually wrong about what

It is tempting to file this as physics overruling geology and then being corrected. The record is less tidy. Thomson's method was sound, and the same half-space solution, applied to the ocean floor rather than the whole planet, is still the standard model of how oceanic lithosphere cools: seafloor depth really does increase as the square root of age, and heat flow really does fall as its inverse square root, for 80 million years. Perry was right about convection before anyone could demonstrate it, but his suggested ages of 2 to 3 billion years were a plausibility argument, not a measurement. And radiogenic heat, which settled the dispute, turned out to supply less than half the outgoing flux, so the Earth is cooling, much as Thomson said. The answer, 4.55 billion years, came from a completely different direction, which is Lesson 16.

Common misconceptions

• "Kelvin made an arithmetic mistake." He did not. Put his inputs into his formula today and you get 84 million years. The defect is the assumption that conduction is the only transport mechanism below the lid.
• "Radioactivity means the Earth is not cooling." Radiogenic heat is about 20 TW against 47 TW of loss. The deficit is real, and the mantle is cooling at something like 100 K per billion years.
• "Ocean floor is hotter because it is nearer the mantle." It is hotter because it is younger. Heat flow falls with the inverse square root of crustal age, from several hundred mW/m² at a ridge axis to about 50 at 100 million years.
• "A solid cannot convect." Solidity is a statement about behaviour on a timescale. The mantle transmits shear waves with a modulus of 294 GPa at periods of seconds and flows like a fluid of viscosity 10²¹ Pa s over millions of years, and both are true at once.

The takeaway

Half-space cooling gives t = T02/(pi kappa (dT/dz)2), and Thomson's inputs of 3870 K, 38 K/km and 1.25 x 10^-6 m²/s return 84 million years. The same numbers give a thermal penetration depth of only 57 km, which is why the result constrains a skin and not a planet. Perry's 1895 objection, that a mobile or more conductive interior keeps resupplying the lid, scales the permitted age with the enhancement factor: ten times gives 840 million years, a hundred times gives 8.4 billion. Strutt's 1906 measurement of radium in rocks supplied the missing source, and the modern budget is 47 TW out against about 20 TW of radiogenic production, a Urey ratio of 0.43, with continental crust alone supplying 32 of the 65 mW/m² measured on land. Geoneutrino counting at KamLAND and Borexino now confirms the radiogenic term independently.

Both men assumed a temperature profile without checking whether it was stable. The next lesson checks.

Sources

1. Wikipedia contributors. (n.d.). Age of the Earth. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Earth's internal heat budget. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Geoneutrino. Wikipedia. en.wikipedia.org
4. Turcotte, D. L., and Schubert, G. (2014). Geodynamics (3rd ed.), Chapter 4. Cambridge University Press.
5. Davies, J. H., and Davies, D. R. (2010). Earth's surface heat flux. Solid Earth, 1(1), 5-24.

A borehole in the Canadian Shield gives a temperature gradient of about 25 K per kilometre. Take it seriously and extrapolate. At the base of the mantle, 2891 km down,

T = 25 x 2891 = 72 275 K.

The surface of the Sun is 5800 K. Iron boils at 3130 K. A mantle at 72 000 K would not be rock, or liquid rock, or even a liquid: it would be a plasma, and the Earth would have no solid interior to transmit the shear waves of Lesson 4. So the gradient cannot continue. Somewhere below the boreholes the temperature profile changes character completely, and this lesson finds where and why.

Where conduction has to stop

Conduction moves heat down a gradient at a rate q = -k dT/dz. With a thermal conductivity of 4 W/(m K), the Canadian gradient carries 4 x 0.025 = 0.100 W/m², which is 100 mW/m², in the right range for measured continental values once crustal radiogenic production is included.

Now ask the reverse question. If the whole mantle were conductive, what gradient would it need to deliver the heat that actually comes out? The mantle supplies about 40 of the 47 TW, the remaining 7 being radioactivity in the continental crust, so over 5.10 x 10¹⁴ m² that is 78 mW/m². A conductive mantle delivering that flux over 2891 km needs a total temperature difference

DeltaT = q d/k = 0.078 x 2.891 x 106/4 = 56 400 K.

The same absurdity from the other direction. Conduction is simply not capable of moving the Earth's heat through 2900 km of rock at any temperature the rock could survive. Something else carries it, and the only candidate is motion of the rock itself.

Key idea: the impossibility of the conductive geotherm is not a detail. It is the argument that the mantle must convect, and it needs no observation beyond a borehole gradient and a conductivity.

The gradient a moving mantle actually has

If material circulates, most of the mantle is neither heating nor cooling in place: a parcel carries its heat with it and changes temperature only through compression and expansion. That is the adiabatic gradient,

(dT/dz)ad = alpha g T/cp,

with alpha the volume coefficient of thermal expansion and c_p the specific heat. Evaluate it in the upper mantle: alpha = 3 x 10-5 per K, g = 9.9 m/s², T = 1600 K, cp = 1250 J/(kg K):

(dT/dz)ad = (3 x 10-5)(9.9)(1600)/1250 = 0.475/1250 = 3.8 x 10-4 K/m = 0.38 K per km.

Sixty-five times gentler than the borehole gradient. Carried over the 2600 km of convecting mantle between the boundary layers, and allowing for alpha falling with pressure toward 1 x 10^-5, the adiabat contributes only about 800 to 1000 K of the total rise. The geotherm therefore has three parts, and almost all of the temperature change happens in two thin ones.

That is the shape to remember: steep, flat, steep. Nearly all the resistance to heat flow lives in two boundary layers a few percent of the mantle's thickness, and the enormous middle is nearly isothermal once the adiabat is removed.

Does it convect? Compute the number

Convection begins when buoyancy overcomes the two things that oppose it, viscous drag and thermal diffusion. The ratio is the Rayleigh number:

Ra = alpha rho g DeltaT d3/(kappa eta).

Every term is now in hand. Use whole-mantle averages: alpha = 2 x 10-5 per K, since expansivity falls with pressure; rho = 4500 kg/m³; g = 9.8 m/s²; the superadiabatic temperature difference DeltaT = 1500 K, which is the part of the total that is available to drive flow once the adiabat is subtracted; d = 2.9 x 106 m; kappa = 1 x 10-6 m²/s; and the viscosity from Fennoscandia, eta = 1 x 1021 Pa s.

Build it up so you can see where the size comes from. d3 = 2.4389 x 1019 m³. Then

alpha rho = 0.09, times g gives 0.882, times DeltaT gives 1323, times d3 gives 3.227 x 1022.

The denominator is kappa eta = 10-6 x 1021 = 1015. So

Ra = 3.227 x 1022/1015 = 3.2 x 107.

Compare that with the critical value. Linear stability analysis of a layer heated from below gives Ra_c = 657.5 for free-slip boundaries and 1708 for rigid ones. The mantle exceeds the harder of those by a factor of

3.2 x 107/1708 = 18 700.

This is not a marginal case. Note also that the d3 dependence is what makes the answer so large: cut the layer to 700 km, the upper mantle alone, and Ra falls by a factor of 71, to 4.5 x 10⁵, still supercritical by more than two hundred times. Convection in the mantle is not a hypothesis that needs delicate conditions. It is hard to avoid.

How much good the convection does

The Nusselt number is the ratio of the heat actually transported to what conduction alone would carry across the same layer with the same temperature difference. Conduction across 2900 km with DeltaT = 1500 K and k = 4 W/(m K) gives

qcond = 4 x 1500/2.9 x 106 = 2.07 x 10-3 W/m² = 2.07 mW/m².

The observed mantle flux is 78 mW/m², so

Nu = 78/2.07 = 38.

Convection carries thirty-eight times what conduction could. And because all the resistance sits in the boundary layers, the top one must be thin enough to conduct the whole flux:

delta = d/Nu = 2 900/38 = 76 km.

The observed thickness of oceanic lithosphere, from the flattening of seafloor depth with age and from the depth of the seismic low-velocity zone, is 90 to 120 km. A scaling argument with no free parameters has predicted the thickness of a tectonic plate to within 30 percent. That is the single best check that the whole picture hangs together.

The top boundary layer also reveals itself directly. A plate cooling as a half-space thickens as 2.32 sqrt(kappa t), which at 80 million years is 2.32 x sqrt(10-6 x 2.52 x 1015) = 2.32 x 50 200 = 116 km, and seafloor depth follows 2500 + 350 sqrt(t) metres with t in millions of years: 2500 m at a ridge crest and 6000 m at 100 million years. Both of those are Lesson 7's half-space solution, now recognised as a boundary layer rather than a whole planet.

Two numbers that say what kind of fluid this is

Mantle convection is unlike any flow in ordinary experience, and two dimensionless numbers show how far.

The Prandtl number compares momentum diffusion with heat diffusion:

Pr = eta/(rho kappa) = 1021/(4500 x 10-6) = 2.2 x 1023.

The Reynolds number compares inertia with viscosity. Take a plate speed of 5 cm per year, which is 1.58 x 10-9 m/s:

Re = rho u d/eta = (4500)(1.58 x 10-9)(2.9 x 106)/1021 = 2 x 10-20.

Inertia is irrelevant by twenty orders of magnitude. There is no turbulence, no eddies, no vortex shedding. Stop driving the mantle and it stops moving essentially instantly in geological terms. Flow patterns are set entirely by the balance of buoyancy against viscous resistance, which is why numerical mantle convection solves the Stokes equations with the inertial terms deleted rather than the full Navier-Stokes equations.

One more timescale. A parcel moving at 5 cm a year crosses the mantle's 2900 km in 2.9 x 106/0.05 = 5.8 x 107 years, so a full circuit takes a few hundred million years. That is why the deep mantle retains chemical heterogeneity: the stirring time is comparable to the age of the ocean basins, not short compared with the age of the Earth.

The upshot: the mantle is a solid at infinite Prandtl number and zero Reynolds number, supercritical by four orders of magnitude, moving a centimetre or two a year, and transporting thirty-eight times the heat conduction could manage.

Common misconceptions

• "The mantle convects, so it must be molten." It is solid, with a shear modulus of 294 GPa at its base. Flow happens by the slow migration of crystal defects, which gives an effective viscosity of 10²¹ Pa s, and solidity at seismic periods is perfectly compatible with that.
• "The geotherm is the borehole gradient continued downward." Extrapolating 25 K/km reaches 72 000 K at the core. The real profile is steep through two thin boundary layers and nearly adiabatic, at 0.3 to 0.4 K/km, through the 2600 km between them.
• "The lithosphere is a chemically distinct layer." It is primarily the conductive thermal boundary layer at the top of the convecting system, which is why its thickness grows as the square root of age, reaching about 116 km at 80 million years.
• "A larger Rayleigh number means faster but otherwise similar convection." It also means thinner boundary layers and more vigorous small-scale structure. Since delta scales as d/Nu and Nu grows with Ra, increasing Ra makes the plates thinner, not merely quicker.

Summing up

A conductive geotherm is impossible: 25 K/km extrapolated reaches 72 275 K at the core-mantle boundary, and a conductive mantle carrying the observed 78 mW/m² would need a temperature difference of 56 400 K. What the mantle has instead is an adiabat, alpha g T/cp = 0.38 K/km in the upper mantle, sandwiched between a 100 km conductive lithosphere and a 200 to 300 km conductive D double prime. The Rayleigh number alpha rho g DeltaT d3/(kappa eta) evaluates to 3.2 x 10⁷, which is 18 700 times the rigid-boundary critical value of 1708, and even the upper mantle alone exceeds it two hundredfold. The Nusselt number is 38, so the top boundary layer must be about 76 km thick to conduct the whole flux, against an observed 90 to 120 km. A Prandtl number of 2.2 x 10²³ and a Reynolds number of 2 x 10^-20 say that inertia plays no part at all, and a transit time of 58 million years across the mantle says the stirring is slow enough to leave heterogeneity behind.

Which is a prediction. If the mantle fails to mix completely, the leftovers should be visible, and the next lesson goes looking for them.

Sources

1. Wikipedia contributors. (n.d.). Mantle convection. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Rayleigh number. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Geothermal gradient. Wikipedia. en.wikipedia.org
4. Turcotte, D. L., and Schubert, G. (2014). Geodynamics (3rd ed.), Chapter 6. Cambridge University Press.
5. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 7. Cambridge University Press.

A P wave leaving a Tonga earthquake and recorded in Sweden should arrive, according to the one-dimensional model of Lesson 6, at a time you can predict to within a second. Sometimes it arrives four seconds early. Sometimes it is three seconds late. Those residuals are not noise. Collect a few million of them, from thousands of earthquakes to thousands of stations, and the pattern in them is a map of where the mantle is fast and where it is slow.

This lesson is about that map, and about the argument it has failed to settle.

What a residual is worth

Travel time along a ray is t = integral of ds/v. Perturb the velocity slightly and the change in arrival time is

dt = -integral of (dv/v2) ds.

Take a concrete case: a ray crossing 1000 km of mantle where the velocity is 10 km/s and is 1 percent fast. Then

dt = -(0.01)(1000/10) = -1.0 second.

One second per thousand kilometres per percent. Since arrival times can be picked to a few tenths of a second on a good record, tomography can in principle see anomalies of a few tenths of a percent, which for mantle rock corresponds to a temperature difference of roughly 100 K or a small change in composition.

The difficulty is not sensitivity but geometry. Earthquakes happen at plate boundaries and seismometers stand on land, so rays sample the mantle unevenly: beneath western North America, Japan and Europe the coverage is dense; beneath the southern Indian Ocean it is thin. An inversion asked to fit residuals with more unknowns than well-constrained directions will invent structure, so every tomographic model is damped, and the damping smears real anomalies and shrinks their amplitudes.

The standard honesty check is the checkerboard test: put a synthetic pattern of alternating fast and slow blocks into the model, compute the arrival times the real ray set would have measured, invert those, and see how much of the checkerboard comes back. In well-sampled regions blocks of 200 to 300 km are recovered; under the mid-ocean gaps, features smaller than 1000 km simply do not return. Read any tomographic image with that in mind, because amplitude and sharpness are model choices as much as observations.

What matters here: a tomographic image is the output of an inversion with a regularisation parameter, not a photograph. The right question about a feature is never only whether it is there, but whether the ray coverage could have resolved it.

The thing everybody agrees about: slabs

The clearest result in seismic tomography is that subducted lithosphere is visible as fast material, and that some of it reaches the bottom of the mantle.

The best-imaged example is the Farallon Plate, which spent 150 million years sliding eastward beneath North America and was almost entirely consumed by about 30 million years ago. It has not gone anywhere. Tomography finds a coherent sheet of fast mantle beneath eastern North America and the western Atlantic, mostly between 800 and 2000 km depth, in the right place and with roughly the right area to be that plate.

Turn it into a rate. If a slab that left the surface 150 million years ago now lies 2000 km down, its mean sinking speed is

2.0 x 106 m/1.5 x 108 yr = 1.3 cm per year,

comparable to plate speeds at the surface, and implying that a slab needs about 2.9 x 106/0.013 = 220 million years to reach the core. That single number kills strict layered convection: if slabs cross 660 km and keep going, the upper and lower mantle exchange material, and the two cannot be separate chemical reservoirs.

Not every slab does it. Beneath the Izu-Bonin arc the fast anomaly flattens and lies horizontally in the transition zone for over a thousand kilometres before descending, which is the negative Clapeyron slope of Lesson 6 doing exactly what it was predicted to do. Slab behaviour at 660 is a spectrum, not a rule.

The thing at the bottom

Sitting on the core-mantle boundary are two enormous regions where shear velocity is 2 to 3 percent low, one beneath Africa and one beneath the central Pacific, roughly antipodal, each around 1000 km tall and together covering about 20 percent of the boundary. They are the large low-shear-velocity provinces, and they appear in every global model built since the 1980s, whatever the data or the damping.

Are they simply hot? Three observations say probably not, or at least not only.

• Their edges are sharp. Waveform modelling of S waves grazing the African province requires a boundary a few tens of kilometres wide, and a purely thermal anomaly in a convecting fluid would diffuse into a gradual one.
• Bulk sound speed and shear speed are anti-correlated inside them: shear velocity drops while bulk sound speed rises slightly. Heating lowers both. A change in composition, for instance more iron, can lower one and raise the other.
• The temperature excess needed to explain a 3 percent shear drop thermally is of order 1500 K, which would put the material well above its solidus and produce a great deal more melt than is observed.

The current consensus, and it is a soft one, is that they are thermochemical piles: hot, denser than average because of composition, swept into two heaps by the pattern of slab arrival at the core-mantle boundary, and long-lived because the density excess keeps them from rising.

The dispute: what rises from the bottom of the mantle

In 1971 W. Jason Morgan proposed that narrow columns of hot material rise from the deep mantle, stay fixed while plates move over them, and produce volcanic chains such as Hawaii. Fifty-five years later the proposal is still contested, and the disagreement is a good example of two credible reconstructions of the same evidence.

The case that plumes are real and deep. Ocean island basalts carry helium with ³He/⁴He ratios up to 35 times atmospheric at Hawaii and near 50 at Iceland, against about 8 for mid-ocean ridge basalt. Helium-3 is primordial and is not made in the Earth, so a high ratio means a reservoir that has never been thoroughly degassed, which is easiest to place deep. Full-waveform tomography published by French and Romanowicz in 2015 resolved broad low-velocity conduits, several hundred kilometres wide rather than the narrow pipes of the original model, beneath about a dozen hotspots, extending to the base of the mantle. And the surface locations of most hotspots sit above the margins of the two provinces rather than at random, which is what you would expect if the province edges are where upwellings are launched.

The case that they are not. Don Anderson and Gillian Foulger argued for decades that the plume hypothesis explains too much and predicts too little. Of the several dozen features called hotspots, only a handful show a clean age progression along a chain. Iceland's low-velocity anomaly is imaged confidently in the upper mantle and much less so below 660 km. Melting anomalies can be produced without any deep source: a fertile patch of recycled crust in the shallow mantle melts more readily than ordinary peridotite at the same temperature, and lithospheric extension can crack the plate and let melt through.

Then there is the measurement that made the fixity argument much harder to sustain. Ocean Drilling Program Leg 197, in 2001, cored Detroit Seamount at the northern end of the Emperor chain, 81 million years old, and measured the magnetic inclination locked into its lavas. The paleolatitude that comes out is about 36 degrees north. Hawaii sits at 19 degrees north today. If the hotspot had been fixed, the two would agree. Instead the source appears to have moved south by roughly 15 degrees between 81 and 47 million years ago, at something like 40 to 50 mm per year, comparable to a plate speed.

What would settle it is not rhetoric but resolution. Full-waveform inversion using entire seismograms rather than picked arrival times, fed by ocean-bottom seismometers that fill the sampling gaps, is steadily improving the images. A conduit continuously resolved from the province margin to the surface beneath a chain, in models with different damping and different data, would close the argument for that chain. So far, the strongest cases are Hawaii, Iceland, Reunion and Samoa, and the weakest are the many features that were called hotspots because they were volcanoes in an unexpected place.

Bottom line: slabs reaching the core-mantle boundary and two thermochemical piles sitting on it are established. Narrow plumes rising from those piles are the best available explanation for a small number of chains and an over-used explanation for many others.

Common misconceptions

• "Tomography photographs the mantle." It solves an underdetermined inverse problem with damping. Anomaly amplitudes are systematically reduced by that damping, and features smaller than the local resolution length can be artefacts. Checkerboard tests are how you find out which is which.
• "Blue means cold and red means hot." The colours mean fast and slow. Velocity responds to temperature, composition, melt fraction and grain size. The anti-correlation of shear and bulk sound speed in the low-velocity provinces is exactly a case where the thermal reading fails.
• "Slabs stop at the 660." Some do stall and lie flat for a thousand kilometres, but the Farallon slab is imaged down to 2000 km and other fast anomalies reach the core-mantle boundary. Whole-mantle exchange is not in doubt.
• "Hotspots are fixed reference points for plate motion." Leg 197 measured a paleolatitude of 36 degrees north for an 81 million year old Emperor seamount against Hawaii's present 19 degrees. The source moved, at a rate comparable to plate speeds.

Pulling it together

A travel-time residual of one second buys you a one percent anomaly over a thousand kilometres at 10 km/s, which is sensitive enough to see 100 K of temperature difference; what limits tomography is ray coverage, and the checkerboard test is how the limit is quantified. The Farallon slab, imaged between 800 and 2000 km beneath eastern North America, gives a sinking rate near 1.3 cm per year and a mantle transit time of about 220 million years, which rules out strictly layered convection. On the core-mantle boundary sit two provinces, beneath Africa and the Pacific, 2 to 3 percent slow in shear, about 1000 km tall, covering a fifth of the boundary, with sharp edges and anti-correlated bulk sound speed that together argue for a compositional as well as thermal origin. Whether narrow plumes rise from their margins is genuinely open: helium-3 ratios of 35 and 50 times atmospheric and full-waveform images of broad conduits support it, while the absence of age progressions at most hotspots and the 15 degrees of southward motion measured at Detroit Seamount undermine the fixity the original model assumed.

Heat and flow are now in place. The next module turns to the field the flowing core generates, which supplies both a clock and a record of where every plate has been.

Sources

1. Wikipedia contributors. (n.d.). Seismic tomography. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Large low-shear-velocity provinces. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Mantle plume. Wikipedia. en.wikipedia.org
4. French, S. W., and Romanowicz, B. (2015). Broad plumes rooted at the base of the Earth's mantle beneath major hotspots. Nature, 525(7567), 95-99.
5. Tarduno, J. A., et al. (2003). The Emperor Seamounts: Southward motion of the Hawaiian hotspot plume in Earth's mantle. Science, 301(5636), 1064-1069.

In 1839 Carl Friedrich Gauss published a general theory of terrestrial magnetism, and with it a method that has never been superseded. He wrote the magnetic potential as a sum of spherical harmonics, split that sum into terms that could only come from inside the Earth and terms that could only come from outside, and fitted both to the magnetic observations then available. The external coefficients came out small. The field you measure with a compass is made below your feet, not above your head, and Gauss proved it with arithmetic rather than assertion.

Everything since has been about where below, and how.

The shape of the field, in numbers

To first approximation the geomagnetic field is a dipole at the centre of the Earth, tilted from the rotation axis. For a dipole, at magnetic colatitude theta, the two components at the surface are

Br = -2 B0 cos(theta) and Btheta = -B0 sin(theta),

so the magnitude is B = B0 sqrt(1 + 3 cos2theta). With B0 = 30 microtesla, the field is 30 microtesla at the magnetic equator and 60 at the poles: a factor of two from end to end, which is why a compass needle that balances level in Ecuador dips hard in Norway.

The relation that matters most for the rest of this course is the inclination. Divide the two components:

tan(I) = Br/Btheta = 2 cot(theta) = 2 tan(lambda),

where lambda is magnetic latitude. At 45 degrees, tan I = 2, so I = 63.4 degrees. Measure the inclination frozen into a rock and you can solve for the latitude at which it formed. Lesson 11 lives on that equation.

The dipole moment follows from B0 = mu0 m/(4 pi a3). Inverting with B0 = 29.8 microtesla and a = 6371 km gives

m = 4 pi a3 B0/mu0 = 7.7 x 1022 A m².

Now the real field. The International Geomagnetic Reference Model gives the first three Gauss coefficients for 2020 as g10 = -29 405, g11 = -1451 and h11 = 4653 nT. The dipole magnitude is the root sum of squares,

sqrt(29 4052 + 14512 + 46532) = sqrt(8.704 x 108) = 29 502 nT,

and the tilt from the rotation axis is

tan(tilt) = sqrt(14512 + 46532)/29 405 = 4874/29 405 = 0.1657, so tilt = 9.4 degrees.

That places the geomagnetic pole near 80.6 degrees north. About 90 percent of the field's energy at the surface sits in these three coefficients; the rest, the non-dipole field, is what makes local declination vary from place to place in ways no globe can capture.

Remember: the geomagnetic pole, from the tilted dipole fit, and the magnetic dip pole, where a needle stands vertical, are different points hundreds of kilometres apart. Neither is the geographic pole, and only the first is what paleomagnetism reconstructs.

The field will not hold still

London's declination was 11 degrees east of north in 1580. By 1820 it was 24 degrees west. It is about 1 degree east today. A compass corrected with a four-hundred-year-old chart would put a ship on the rocks.

This is secular variation, and its statistics are the strongest constraint on where the field is made. Four facts:

• The dipole moment has fallen from about 8.5 x 10²² A m² in 1840 to 7.7 x 10²² now, roughly 9 percent in 180 years.
• Non-dipole features drift westward at about 0.2 degrees of longitude per year, which at the core-mantle boundary is a speed near 15 km per year.
• The north magnetic dip pole crept across the Canadian Arctic at some 10 km per year through most of the twentieth century, then accelerated to 50 or 60 km per year after 1990 and crossed the date line.
• Over the South Atlantic there is a region where the surface field is 30 percent below the dipole value, and it is growing. Satellites in low orbit take measurable radiation damage crossing it.

A permanently magnetised rock cannot do any of this. Magnetite loses its magnetisation above its Curie temperature of 580 degrees Celsius, and hematite above 680, and those temperatures are reached at 20 to 30 km depth on a continental geotherm. So the only permanently magnetisable part of the Earth is a thin shell, and a thin rigid shell cannot rearrange its field in decades. The source must be somewhere that moves, is electrically conducting, and is large. That is the outer core, and nothing else qualifies.

Two numbers that define a dynamo

Consider the field simply sitting in the core and being left alone. Currents decay by ohmic dissipation, and the magnetic field diffuses like heat, with a magnetic diffusivity

etam = 1/(mu0 sigma).

Liquid iron alloy at core conditions has an electrical conductivity of about 5 x 105 S/m, so

etam = 1/((1.2566 x 10-6)(5 x 105)) = 1.59 m²/s.

The free decay time for the largest mode in a sphere of radius R is tau = R2/(pi2 etam). With R = 3.48 x 106 m,

tau = 1.211 x 1013/(9.87 x 1.59) = 7.71 x 1011 s = 24 000 years.

Yet magnetised rocks record a field of comparable strength 3.5 billion years ago. Twenty-four thousand years against three and a half billion is a factor of 145 000. The field is not a leftover. It is being made continuously, right now, and would be gone in a geological instant if the machine stopped.

What keeps it going is fluid motion stretching field lines faster than diffusion can smooth them out. The ratio of the two effects is the magnetic Reynolds number:

Rm = u L/etam.

Take the flow speed inferred from westward drift, u = 5 x 10-4 m/s, that is half a millimetre per second or 15 km per year, and the outer core thickness L = 2.26 x 106 m:

Rm = (5 x 10-4)(2.26 x 106)/1.59 = 1130/1.59 = 710.

Numerical and laboratory dynamos start to self-sustain somewhere between Rm of 10 and 100, depending on the flow geometry. At 710 the core is far above threshold. Compare this with the ordinary Reynolds number of the mantle, 2 x 10^-20, computed in Lesson 8: the two numbers look alike and describe opposite worlds. In the mantle, viscosity wins totally. In the core, advection of the magnetic field wins comfortably.

The core of it: a decay time of 24 000 years says the field must be regenerated; a magnetic Reynolds number of 710 says the core is capable of regenerating it.

How the machine is powered

A dynamo converts mechanical work into magnetic energy, so something must be doing work on the fluid. Two sources are available, and they are unequal.

Thermal convection carries heat out of the core into the mantle, roughly the 9 TW of Lesson 7. Only a fraction of that can be converted, because a heat engine between 4000 K at the inner core boundary and 3700 K at the top is limited by its Carnot factor to about 7 percent, and the ohmic dissipation actually required is estimated at 0.2 to 0.5 TW.

Compositional convection is the more efficient half. As the inner core freezes onto itself, it rejects light elements, probably sulphur, silicon and oxygen, into the liquid just above it. That buoyant fluid rises without needing any heat to be transported at all, so it is not limited by a Carnot factor. Most current models make it the dominant power source, which has a striking implication: the geodynamo may have got substantially stronger when the inner core began to freeze, somewhere between one and 1.5 billion years ago.

The rotation of the Earth organises the result. Coriolis forces are enormous compared with viscous ones in the core, and they align convective motion into columns parallel to the rotation axis. That alignment is why a field generated by turbulent fluid comes out looking like a dipole aligned with the spin axis to within 9.4 degrees, rather than pointing in a random direction.

Common misconceptions

• "The Earth contains a giant magnet." The core is far above the Curie temperature of iron, 770 degrees Celsius, so it cannot be permanently magnetised at all. The field is made by electric currents in moving conducting fluid and would decay in 24 000 years without them.
• "The magnetic pole is where the dipole axis meets the surface." Two different points. The geomagnetic pole comes from fitting the dipole and sits near 80.6 degrees north; the dip pole, where a needle stands vertical, is a local feature of the full field and moves at tens of kilometres a year.
• "The field is weakening, so a reversal is starting." The dipole has fallen 9 percent since 1840, but it has been higher than average for the last two thousand years and the present rate is not exceptional in the paleomagnetic record. A decline of this size is ordinary secular variation.
• "Rocks in the crust generate the field." Crustal magnetisation contributes short-wavelength anomalies of tens to hundreds of nanotesla, useful for mapping geology, but the main field is 50 000 nT and changes on decadal timescales that a rigid magnetised shell cannot produce.

What to remember

Gauss's 1839 harmonic separation put the source inside the Earth. A dipole gives B = B0 sqrt(1 + 3 cos2theta), so 30 microtesla at the equator and 60 at the poles, and tan I = 2 tan(lambda), which is the whole basis of paleolatitude. The 2020 coefficients of -29 405, -1451 and 4653 nT give a dipole of 29 502 nT tilted 9.4 degrees, a moment of 7.7 x 10²² A m², and about 90 percent of the surface field energy. Secular variation, London's declination swinging 35 degrees since 1580 and the dip pole accelerating to 50 km per year, rules out a crustal source, as does a Curie depth of 20 to 30 km. In the core, a conductivity of 5 x 10⁵ S/m gives a magnetic diffusivity of 1.59 m²/s, a free decay time of 24 000 years, and, at a flow speed of half a millimetre per second, a magnetic Reynolds number of 710, comfortably above the threshold for self-sustaining dynamo action. The power comes mostly from light elements released as the inner core freezes, and rotation is what makes the output a dipole rather than a mess.

Sometimes the machine reverses. The record of those reversals is written on the ocean floor, and reading it is the next lesson.

Sources

1. Wikipedia contributors. (n.d.). Earth's magnetic field. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Dynamo theory. Wikipedia. en.wikipedia.org
3. National Centers for Environmental Information. (n.d.). Geomagnetism. NOAA. ncei.noaa.gov
4. Lowrie, W., and Fichtner, A. (2020). Fundamentals of geophysics (3rd ed.), Chapter 5. Cambridge University Press.
5. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 3. Cambridge University Press.

In 1956 Keith Runcorn published a set of paleomagnetic directions from British and European rocks of many ages and drew the conclusion that seemed forced: the magnetic pole had moved. Not a little. From somewhere near the present equator in the Precambrian, north through the Pacific, arriving at its present position only recently. The path had a name, apparent polar wander, and for a few years the obvious reading of it was the literal one.

That reading is wrong, and this lesson follows it until it breaks, because the exact point of failure is one of the most elegant arguments in the Earth sciences.

What a rock actually records

Begin with the measurement, since everything later depends on its reliability. When basalt erupts at 1200 degrees Celsius it has no magnetisation; iron oxide grains are above their Curie temperature and thermally disordered. As the flow cools through the blocking temperature, a few tens of degrees below the Curie point of 580 for magnetite, the grains lock their moments into alignment with the ambient field. That is thermoremanent magnetisation, and in a stable grain it survives for billions of years, because the energy barrier to flipping it is large compared with thermal energy at surface temperature.

Sediments do it differently. Magnetic grains settling through water rotate toward the field before they are buried, giving a detrital remanence that is weaker, slightly shallowed in inclination, and averaged over the time the sediment took to accumulate. That averaging is a feature: a lava flow records an instant, including whatever non-dipole wobble was present, while a sediment core records a mean.

Either way the measurement gives two angles: declination D, the azimuth of the horizontal component, and inclination I, the dip below horizontal. From Lesson 10, tan I = 2 tan(lambda), so inclination alone gives paleolatitude. Longitude is unrecoverable from a single site, because a dipole field is axially symmetric and has no way to label meridians.

Key idea: a rock remembers its latitude and its orientation, never its longitude. Every paleomagnetic reconstruction of the continents is therefore free to slide in longitude, and that freedom is the single largest source of ambiguity in the field.

Turning two angles into a pole

The standard product of a paleomagnetic study is a virtual geomagnetic pole: the position the geomagnetic pole would have had, assuming a centred dipole, to produce the direction you measured where you measured it. Work one.

A site at 40 degrees north, on the prime meridian, yields D = 30 degrees and I = 60 degrees. First the paleolatitude:

tan(lambda) = tan(60)/2 = 1.7321/2 = 0.8660, so lambda = 40.9 degrees,

and the paleomagnetic colatitude, the angular distance from site to pole, is p = 90 - 40.9 = 49.1 degrees. Now walk that distance from the site along the azimuth D, using the spherical cosine rule:

sin(latpole) = sin(latsite) cos(p) + cos(latsite) sin(p) cos(D)

= (0.6428)(0.6543) + (0.7660)(0.7562)(0.8660) = 0.4207 + 0.5016 = 0.9223,

so the pole sits at 67.3 degrees north, at a longitude obtained from a companion formula. Repeat this for rocks of many ages at one continent and the poles trace a path: an apparent polar wander curve.

Where the obvious reading breaks

Suppose the curve means what it says, and the pole really migrated. Then there is exactly one pole at each moment in the past, so every continent must give the same curve. That is a hard, checkable prediction, and it is the one that fails.

Runcorn compared the European path with the North American one. The two have the same shape, the same sequence of bends, the same ages at each bend. They are not in the same place. The North American curve lies about 30 degrees of longitude west of the European one for the whole Paleozoic and Mesozoic, and the separation closes toward the present.

Now trace the failure precisely. There are only three ways out.

• The measurements are wrong. But the paths agree in shape to a precision far better than their separation, and independent laboratories reproduce them. An error that preserved shape while translating position by 30 degrees is not a plausible error.
• There were two poles. The field is 90 percent dipole, and a dipole has one axis. A quadrupole strong enough to fake a second pole would have left other signatures, and would have to have persisted for 400 million years.
• There was one pole, and the continents were not where they are now. Rotate North America eastward about 30 degrees, closing the Atlantic, and the two curves lie on top of each other.

The third option is the only survivor, and it arrives with a number attached: the Atlantic is about 30 degrees wide and was not there in the Paleozoic. Note what has happened to the original conclusion. The pole did not wander; the observer did. This is why the curves are called apparent polar wander, and the word was added after the argument, not before.

So what?: the failure of a prediction is what converted continental drift from a map-shaped intuition into a measurement with an error bar. Wegener had the idea in 1912 and no instrument; Runcorn had an instrument and no need for the idea until his own data forced it.

Reversals, and a clock made of them

There is a second thing the rocks record, and at first it looked like a problem. In 1906 Bernard Brunhes, sampling lavas in the Massif Central, found flows whose magnetisation pointed south and up rather than north and down. In 1929 Motonori Matuyama showed in Japanese and Manchurian basalts that the anomaly was not random: reversed flows were systematically older than normal ones.

Once potassium-argon dating could put ages on lavas, in the early 1960s, the reversal sequence became a timescale.

Reversals are irregular. Over the last ten million years they have come at four or five per million years; during the Cretaceous superchron the field held one polarity for 38 million years. A transition itself takes one to ten thousand years, during which the dipole drops to perhaps a tenth of its usual strength and the direction swings erratically. Remember that duration: it becomes important in the next section.

The prediction that tested everything at once

In 1963 Fred Vine and Drummond Matthews, and independently Lawrence Morley whose paper was rejected by two journals, put two ideas together. Suppose Harry Hess was right that new ocean floor forms at a ridge and spreads away. Suppose the field reverses. Then the basalt cooling through its blocking temperature at the axis is magnetised in whatever polarity is current, and as it spreads it carries that polarity outward. The ocean floor is a tape recorder.

The hypothesis made three predictions that could all fail.

• Magnetic anomalies over ocean floor should form stripes parallel to the ridge axis.
• The pattern should be symmetric about the axis, because crust is added to both sides.
• The sequence of widths should match the sequence of chron durations from dated lavas on land, with one scale factor, the spreading rate.

All three held. The Eltanin-19 profile across the Pacific-Antarctic ridge, published in 1966, is symmetric about the axis to a degree that is uncomfortable to look at if you do not believe in spreading, and the Reykjanes Ridge survey showed the same.

Now extract the number the whole of plate tectonics needs. Suppose a survey across a slow ridge finds reversal boundaries at these distances from the axis:

Check the first: 14.0 km in 0.78 Myr is 14.0/0.78 = 17.9 km per million years, and one km per Myr is one mm per year. The three agree to better than one percent, so the full spreading rate is 2 x 17.8 = 35.6 mm per year, typical of the South Atlantic. Run the same arithmetic on the East Pacific Rise, where the Brunhes-Matuyama boundary lies 55 km out, and you get a half-rate of 55/0.78 = 70.5 mm per year and a full rate of 141, four times faster.

Consistency across three widely separated boundaries is the real result. It means the rate has been steady for six million years, and it means the chron ages derived from potassium-argon dating of continental lavas correctly predict widths on an ocean floor thousands of kilometres away. Two independent clocks, one radiometric and one geometric, agreeing.

And now the detail flagged earlier. A reversal takes one to ten thousand years. At a half-rate of 17.8 mm per year, 10 000 years builds 178 metres of crust. Stripe widths here are 14 000 metres and up. The transition zones are narrow compared with the stripes by a factor of a hundred, which is why the pattern is sharp rather than smeared, and why this method works at all at slow ridges.

Common misconceptions

• "Apparent polar wander shows that the magnetic pole has moved thousands of kilometres." It shows relative motion between the continent and the pole. Because two continents give paths of identical shape in different places, the motion must be the continents', and closing the Atlantic by 30 degrees merges the curves.
• "Magnetic stripes are alternating bands of different rock." The rock is the same basalt throughout. What alternates is the direction of its remanent magnetisation, locked in as it cooled through the blocking temperature.
• "A positive magnetic anomaly lies directly over normally magnetised crust." Only roughly. The measured anomaly is the vector sum of the ambient field and the crustal field, so at low magnetic latitudes it is skewed and the peaks are offset from the block edges; reduction to the pole is applied before widths are read.
• "Paleomagnetism gives the full former position of a continent." It gives paleolatitude and orientation only. An axially symmetric dipole field cannot encode longitude, so reconstructions need an additional constraint, such as hotspot tracks or the fit of continental margins.

Recap

Basalt locks in thermoremanent magnetisation as it cools through a blocking temperature just below 580 degrees Celsius, and tan I = 2 tan(lambda) converts the recorded inclination into paleolatitude: a site at 40 degrees north with D = 30 and I = 60 degrees gives a paleolatitude of 40.9 degrees and a virtual geomagnetic pole at 67.3 degrees north. Runcorn's European and North American polar wander paths have the same shape and differ by about 30 degrees of longitude, and since one dipole cannot have two axes, the continents must have moved; closing the Atlantic by that amount superposes the curves. The reversal sequence, from Brunhes at 0 to 0.78 Ma back through Matuyama, Gauss and Gilbert, and including the 38 million year Cretaceous superchron, provides a dated polarity timescale. Vine, Matthews and Morley predicted that spreading would print that timescale on the seafloor in parallel, symmetric stripes, and the widths deliver spreading rates: 14.0, 46.0 and 105.0 km to three boundaries give half-rates of 17.9, 17.8 and 17.8 mm per year, while the East Pacific Rise gives 70.5. Because a reversal takes only 178 metres of crust to record at slow rates, the stripes stay sharp.

Those rates are velocities at one place. Turning a set of them into the motion of a rigid plate over a sphere is the next lesson.

Sources

1. Wikipedia contributors. (n.d.). Paleomagnetism. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Vine-Matthews-Morley hypothesis. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Geomagnetic reversal. Wikipedia. en.wikipedia.org
4. Vine, F. J., and Matthews, D. H. (1963). Magnetic anomalies over oceanic ridges. Nature, 199(4897), 947-949.
5. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 3. Cambridge University Press.

Leonhard Euler proved in 1775 that any displacement of a rigid body on a sphere, however complicated it looks on a map, is a single rotation about an axis through the centre. Plate tectonics inherited that theorem intact. A plate does not slide across the Earth; it turns about an axis, and the axis pierces the surface at a point called the Euler pole. Two numbers for the pole and one for the rate specify the entire motion of a plate forever.

This lesson is a procedure: take those three numbers, produce a velocity you could measure with a survey mark, and then check it against three instruments that did not exist when the theory was written.

Three numbers into one velocity

For a rotation at angular rate omega about a pole, the linear speed at a point an angular distance Delta from that pole is

v = omega R sin(Delta).

The sine is the whole geometry. At the Euler pole itself the plate spins but does not translate, so v = 0. At 90 degrees from it, on the equator of the rotation, the speed is maximum. Transform faults on that plate must lie along small circles about the pole, which is how Euler poles were first located: fit circles to the transforms.

Work the Pacific relative to North America. The pole sits near 48.7 degrees north, 78.2 degrees west, in Quebec, with an angular rate of 0.78 degrees per million years. San Francisco is at 37.8 north, 122.4 west. Angular distance first, from the spherical cosine rule:

cos(Delta) = sin(48.7)sin(37.8) + cos(48.7)cos(37.8)cos(44.2)

= (0.7513)(0.6129) + (0.6600)(0.7902)(0.7167) = 0.4605 + 0.3738 = 0.8342,

so Delta = 33.5 degrees and sin(Delta) = 0.5513. Convert the rate to radians: 0.78 x pi/180 = 0.013614 rad per million years. Then

v = (0.013614)(6371)(0.5513) = 47.8 km per million years = 48 mm per year.

Geodesy measures about 50 mm per year of Pacific relative to North American motion across California, of which the San Andreas fault itself carries roughly 35 and the rest is distributed across the Eastern California Shear Zone and the Basin and Range. A rigid-plate model built from magnetic anomalies and transform azimuths, containing no Californian data at all, predicts the total to within a few percent.

The point: plate motion is a rotation, not a translation, so the same plate moves at different speeds in different places, and the number you quote is meaningless without a location and a reference plate.

Where the three numbers come from

Global models such as NUVEL-1A and MORVEL are built from three kinds of observation, and it is worth seeing that they are independent.

• Spreading rates from magnetic anomaly widths, exactly as computed in Lesson 11, give the magnitude of relative motion across every ridge.
• Transform azimuths give the direction, because a transform must be parallel to the local small circle about the pole.
• Earthquake slip vectors at subduction zones give direction where there is no ridge to measure.

Then comes the closure condition. Around any circuit of plates the rotation vectors must sum to zero: omegaAB + omegaBC + omegaCA = 0. That is a strong test, because it is over-determined. If the Pacific-Antarctic, Antarctic-Nazca and Nazca-Pacific poles do not close, one of them is wrong or a plate is not rigid. In practice the misfits are small, which is the quantitative statement that plates really do behave as rigid bodies to within a few millimetres per year over their interiors.

A chain of volcanoes as a velocity record

The Hawaiian-Emperor chain runs 6000 km from Kilauea, erupting now, northwest to Midway at 27.7 million years, and then turns sharply north and continues to Detroit Seamount at 81 million years, near the Aleutian trench.

Take the first segment as a speedometer. Midway lies 2432 km from Kilauea and its lavas are 27.7 million years old:

2432 km/27.7 Myr = 87.8 km per million years = 88 mm per year.

That is the Pacific plate's speed over whatever is making the volcanoes, and it is among the fastest motions on Earth.

The bend is the interesting part. At about 47 million years the chain changes azimuth by roughly 60 degrees, from a north-northwesterly Emperor trend to a west-northwesterly Hawaiian one. The textbook reading for thirty years was that the Pacific plate changed direction, and the date was used to time a global tectonic reorganisation.

Lesson 9 supplied the complication: Ocean Drilling Program Leg 197 measured a paleolatitude near 36 degrees north for 81 million year old lavas at Detroit Seamount, against Hawaii's 19 degrees today. If the source drifted south at 40 to 50 mm per year while the plate moved north, then some of the bend is the source moving and only some is the plate turning. Current reconstructions split it, with hotspot motion dominating before 47 Ma and plate motion change after. The chain is still a velocity record; it is a record of relative velocity between plate and source, which is not the same thing as absolute plate motion, and treating the two as identical was the error.

Measuring it directly, three ways

Everything above is inferred from geology millions of years old. Since about 1990 the motions have been watched happening.

GPS. A geodetic receiver does not use the navigation code; it tracks the phase of the carrier wave, whose wavelength is 19 cm, and after days of observation and corrections for satellite orbits, clocks, the ionosphere and the wet troposphere, a station position is known to a few millimetres horizontally. Repeat for five years and the slope of the position time series is a velocity with an uncertainty under 1 mm per year. Continuous stations on plate interiors reproduce the geological rates: the Pacific plate near Hawaii moves about 70 mm per year west-northwest in a global reference frame, stable Eurasia about 25, and the Nazca plate over 60 toward South America.

InSAR. Two radar images of the same ground from the same orbit at different times can be differenced in phase. Each cycle of phase difference, one fringe, corresponds to half a radar wavelength of displacement along the line of sight, because the signal travels out and back. For Sentinel-1, wavelength 5.55 cm, one fringe is 2.77 cm. An interferogram showing twelve fringes across a valley therefore records

12 x 2.77 = 33 cm of line-of-sight motion,

mapped continuously over the whole scene rather than at scattered survey marks. That is the trade: GPS gives three components at a point, InSAR gives one component everywhere.

Seafloor geodesy. Radio does not penetrate seawater, so a ship or buoy fixes itself by GPS and then ranges acoustically to transponders on the seabed. The precision is centimetres rather than millimetres, and each campaign is expensive, but it reaches the part of a subduction zone that matters most: the offshore fault.

The 11 March 2011 displacement field

Japan's GEONET network of about 1200 continuous GPS stations was running when the magnitude 9.0 Tohoku earthquake broke. The coseismic field it recorded is the best-documented in history.

The pattern itself carries the physics. Land subsided while the seafloor rose, which is the signature of slip on a shallowly dipping thrust: the hanging wall moves up and seaward near the trench and the region behind it drops. And the fact that displacement grows from 5.3 m on land to 24 m offshore says the slip was concentrated near the trench, in the shallowest part of the interface, which was precisely the part most models had assumed was too weak to store elastic strain. That assumption is why the tsunami was underestimated.

One arithmetic check ties it to the plate model. The Pacific plate converges with northeast Japan at about 83 mm per year. Fifty metres of slip represents

50 m/0.083 m per year = 600 years of accumulated convergence.

The last comparable event on that segment is the Jogan earthquake of 869 CE, 1142 years earlier, whose tsunami deposits had been mapped inland on the Sendai plain in the years before 2011. The geodesy, the geology and the plate rate agree to within a factor that matters far less than the fact that all three pointed the same way.

Worth holding on to: an Euler pole predicts a velocity, GPS measures it, InSAR maps its spatial pattern, and the difference between what a fault has stored and what it has released is the earthquake budget of Module 5.

Common misconceptions

• "A plate moves at one speed." Speed is omega R sin(Delta), so it varies from zero at the Euler pole to a maximum 90 degrees away. Quoting a plate velocity without a location is like quoting a wind speed without a place.
• "The Euler pole is where the plate is moving fastest." The exact opposite. At the pole the plate rotates in place and translates not at all, which is why transform faults trace small circles centred on it.
• "The Hawaii-Emperor bend dates a change in Pacific plate motion." Only partly. The paleolatitude of Detroit Seamount requires the source itself to have moved south, so the bend records relative motion between plate and hotspot, and separating the two takes independent evidence.
• "One InSAR fringe equals one radar wavelength of ground motion." Half a wavelength, because the pulse makes a round trip. For Sentinel-1's 5.55 cm that is 2.77 cm per fringe, and using the wrong factor doubles every displacement you report.

What to carry forward

Euler's theorem reduces the motion of a rigid plate to a pole and a rate, and v = omega R sin(Delta) turns those into a velocity anywhere: a Pacific-North America pole at 48.7 N, 78.2 W rotating at 0.78 degrees per million years gives 48 mm per year at San Francisco, against about 50 measured. Global models are built from spreading rates, transform azimuths and slip vectors, and tested by the requirement that rotation vectors close around every plate circuit. The Hawaiian chain gives 88 mm per year from Kilauea to Midway, and its 60 degree bend at 47 Ma records plate motion and hotspot motion together rather than plate motion alone. Modern geodesy checks all of it: GPS to a few millimetres, InSAR at 2.77 cm per fringe for Sentinel-1, and GPS-acoustic on the seabed. In 2011 those instruments recorded 5.3 m of eastward motion and 1.2 m of subsidence on the Oshika Peninsula, 24 m east and 3 m up on the seafloor, and about 50 m of slip near the trench, which is 600 years of convergence at 83 mm per year released in three minutes.

Which raises the question the next module opens with: what pushes the plates in the first place?

Sources

1. Wikipedia contributors. (n.d.). Euler pole. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Hawaiian-Emperor seamount chain. Wikipedia. en.wikipedia.org
3. Jet Propulsion Laboratory. (n.d.). GPS time series. NASA. sideshow.jpl.nasa.gov
4. Wikipedia contributors. (n.d.). 2011 Tohoku earthquake and tsunami. Wikipedia. en.wikipedia.org
5. Fowler, C. M. R. (2005). The solid Earth: An introduction to global geophysics (2nd ed.), Chapter 2. Cambridge University Press.

Lay the twelve major plates out in a table with two columns: how big each one is, and how fast it moves. There is no relationship. Africa covers 78 million square kilometres and moves 21 mm per year; Eurasia covers 68 million and moves 7; the Cocos plate covers 2.9 million, less than four percent of Africa's area, and moves 86.

Now replace area with the fraction of each plate's perimeter that is a subduction zone. The relationship appears immediately, and it is not subtle.

Every plate with a substantial subducting edge moves at 60 mm per year or more. Every plate without one moves at 27 or less. Donald Forsyth and Seiya Uyeda made that observation quantitative in 1975, and the argument about what it means has been running ever since.

The candidate forces, sized

Three forces are usually proposed, and they can be estimated from quantities already established in this course.

Slab pull. A descending slab is colder than the mantle around it, so it is denser, and the excess weight pulls. Take a density excess of 80 kg/m³, from thermal contraction plus the elevated 410 transition inside the cold slab, a slab 100 km thick, and a descending length of 600 km. The force per metre of trench is

Fsp = (80)(9.8)(1 x 105)(6 x 105) = 4.7 x 1013 N/m.

Ridge push. This is badly named. There is no engine at the axis shoving plates apart. What there is, is a horizontal pressure gradient inside a plate that thickens and deepens as it cools, so the integrated weight of the elevated young lithosphere pushes the older, deeper lithosphere sideways. For a plate of age t the result is approximately g alpha rhom DeltaT kappa t. At 100 million years, that is 3.156 x 1015 s, and with alpha = 3 x 10-5, rhom = 3300, DeltaT = 1300 K and kappa = 10-6:

Frp = (9.8)(3 x 10-5)(3300)(1300)(10-6)(3.156 x 1015) = 4.0 x 1012 N/m.

So slab pull exceeds ridge push by a factor of about 12. Both are real, and one is an order of magnitude larger.

Basal drag. Viscous shear between plate and asthenosphere. With an asthenospheric viscosity of 10¹⁹ Pa s, a plate speed of 5 cm per year, that is 1.58 x 10-9 m/s, and a 100 km channel,

tau = eta u/h = (1019)(1.58 x 10-9)/105 = 1.6 x 105 Pa, that is 0.16 MPa.

Integrated over the Pacific plate's 1.03 x 1014 m², the total is 1.6 x 1019 N. Compare slab pull along the Pacific's roughly 20 000 km of trench: 4.7 x 1013 x 2 x 107 = 9.4 x 1020 N, some 58 times larger. And note the sign: if the plate moves faster than the mantle beneath, drag opposes the motion. Forsyth and Uyeda's regression found exactly that, which is also why plate area does not predict plate speed. A bigger plate has more drag as well as more of whatever drives it.

In short: slab pull of 4.7 x 1013 N/m, ridge push of 4.0 x 1012 N/m, and a basal drag that mostly resists. That is the arithmetic both sides of the dispute accept.

The slab-pull account, and the number that embarrasses it

The straightforward reading is that slabs pull the plates behind them, the way a chain hanging over the edge of a table pulls the rest of the chain. Forsyth and Uyeda inverted twelve plate velocities against six candidate force terms by least squares and found the coefficient on the trench-length term dominant and the others small. Clint Conrad and Carolina Lithgow-Bertelloni reached a similar conclusion in 2002 by a completely different route, driving a mantle flow model with the density field inferred from slabs and finding that it reproduces most of the observed plate motion.

Here is the difficulty. Suppose the full 4.7 x 1013 N/m were transmitted as tension through oceanic lithosphere 100 km thick. The stress would be

sigma = 4.7 x 1013/105 = 4.7 x 108 Pa = 470 MPa.

Oceanic lithosphere cannot carry that. Its integrated strength, from laboratory friction and creep laws, corresponds to a few tens of megapascals of deviatoric stress averaged over the plate, and the World Stress Map, assembled from boreholes and earthquake focal mechanisms, finds plate interiors at 20 to 30 MPa, often in compression rather than tension. A plate under 470 MPa of tension would simply tear.

So most of the slab's negative buoyancy is not delivered to the surface plate as tension. It is balanced locally: by viscous resistance as the slab pushes mantle aside, by bending resistance at the trench, by the negative Clapeyron slope at 660 km from Lesson 6. Estimates of the fraction reaching the plate run from 10 to 30 percent, which brings the transmitted stress down to 50 to 140 MPa, still high but no longer impossible once the plate's thickness is taken as its full mechanical thickness rather than 100 km.

The mantle-flow account

The alternative is not that slabs are unimportant but that they drive the plates indirectly. On this reading, laid out most persistently by Don Anderson and developed in many convection models since, plates are the cold upper boundary layer of mantle convection, and asking what drives them is like asking what drives the top of a pot of boiling water. The density anomaly of a slab drives flow in the whole mantle, and that flow exerts tractions on the base of every plate, including plates with no subducting edge at all.

Africa is the test case this account points to. It has essentially no subducting boundary, so the slab-pull term for it is zero, and yet it moves at 21 mm per year while rifting internally along the East African Rift. Something is moving it, and the candidates are basal traction from mantle flow, in particular from the upwelling above the African low-velocity province of Lesson 9, and the integrated ridge push from the long spreading boundaries that surround it.

The observation that bears most directly on the question is seismic anisotropy. Olivine crystals align with the direction of shear, so the splitting of SKS phases beneath a plate records the direction in which the asthenosphere is being sheared. Beneath most fast-moving oceanic plates the fast direction is parallel to absolute plate motion, which is what you expect if the plate is dragging the mantle and drag is resistive. Beneath parts of the western United States and beneath Africa the fast directions do not follow plate motion, which is what you expect if mantle flow is organised by something other than the plate above it.

Why this matters: the two accounts are not mirror images, and they make different predictions. The slab-pull picture predicts that plate speed should be almost fully determined by trench length, and it nearly is. The mantle-flow picture predicts residuals that correlate with deep structure, and Africa is one. A complete answer needs both, which is why the modern literature no longer asks which force wins but how much of the slab's energy reaches the surface as plate-parallel tension and how much as basal traction.

What would settle it

Three lines of work are narrowing the question. Stress measurements in plate interiors test how much tension a plate actually carries, and the World Stress Map already excludes the naive 470 MPa. Global flow models that predict plate velocities, geoid, dynamic topography and anisotropy simultaneously from one density field are a strong test, because a model tuned to fit plate motions alone can be wrong about everything else. And seafloor seismometers are beginning to give anisotropy under the oceans, where coverage has always been worst and where the drag argument is cleanest.

Common misconceptions

• "Ridge push is the ridge pushing." There is no engine at the axis. The force is the horizontal gradient of pressure inside a plate that thickens and subsides as it cools, and it grows with plate age: about 4.0 x 10¹² N/m at 100 million years.
• "Convection currents drag the plates along like conveyor belts." Basal drag computes to 0.16 MPa and, for fast plates, mostly resists motion. That is why plate speed is uncorrelated with plate area, since area scales the drag as well as any driving traction.
• "Slab pull of 4.7 x 10¹³ N/m acts as tension in the plate." That would be 470 MPa in a 100 km plate, far beyond its strength, and plate interiors measure 20 to 30 MPa. Most of the slab's weight is balanced by resistance to its own descent.
• "Africa has no subduction, so plate tectonics cannot explain its motion." It explains it with the other terms: long surrounding ridges and basal traction from large-scale mantle flow. Africa is not a counterexample, it is the case that requires the second mechanism.

The takeaway

Plate speed is uncorrelated with plate area and strongly correlated with the fraction of a plate's boundary that subducts: Cocos at 0.42 and 86 mm per year against Eurasia at 0.00 and 7. Slab pull, from an 80 kg/m³ density excess over a 100 km by 600 km slab, is 4.7 x 10¹³ N per metre of trench; ridge push at 100 million years is 4.0 x 10¹², twelve times smaller; basal drag is 0.16 MPa and, integrated over the Pacific, 58 times smaller than its slab pull, with a sign that resists. Forsyth and Uyeda's 1975 regression on twelve plates found the trench term dominant, and flow models driven by slab density reproduce most plate motion. But the full slab pull would require 470 MPa of tension in the plate, against 20 to 30 MPa measured, so only 10 to 30 percent of it can be reaching the surface as tension and the rest is balanced by resistance to descent and delivered, if at all, through the mantle. Africa, with no subducting edge and 21 mm per year, is the case that needs basal traction, and seismic anisotropy is the measurement that distinguishes the two.

Whatever the force, plate boundaries do not slide smoothly. They stick, and then they do not, and the next lesson quantifies the difference.

Sources

1. Wikipedia contributors. (n.d.). Slab pull. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Ridge push. Wikipedia. en.wikipedia.org
3. United States Geological Survey. (n.d.). This dynamic Earth: The story of plate tectonics. pubs.usgs.gov
4. Forsyth, D., and Uyeda, S. (1975). On the relative importance of the driving forces of plate motion. Geophysical Journal International, 43(1), 163-200.
5. Turcotte, D. L., and Schubert, G. (2014). Geodynamics (3rd ed.), Chapter 6. Cambridge University Press.

The State Earthquake Investigation Commission that reported on the 1906 San Francisco earthquake had something no previous commission had ever had: three geodetic surveys of the same ground, from 1851 to 1865, from 1874 to 1892, and again immediately after the shaking. Harry Fielding Reid laid the three sets of triangulation points on top of each other and found something nobody expected. Marks far from the fault had been moving steadily for fifty years. Marks near the fault had barely moved at all until 18 April 1906, when they jumped several metres and landed back on the smooth trend the distant points had been following all along.

The ground had been storing the motion, and then it gave it back. Reid called it elastic rebound, and this lesson turns his picture into numbers.

What is actually stored, and where

Elastic rebound has three stages, and keeping them separate prevents most of the confusion in the subject.

• Loading. The plates move at their steady rate. The fault surface is locked by friction, so the rock on either side deforms elastically instead of sliding. Strain energy accumulates in a volume tens of kilometres across, not in the fault plane.
• Rupture. Shear stress somewhere on the locked patch exceeds the frictional strength. That patch slips, which raises the stress on its neighbours, which slip in turn. The failure front runs along the fault at 2 to 3 km per second.
• Relaxation. The two sides have returned toward their unstrained shape, having released most but not all of the accumulated strain, and loading begins again.

Note where the energy lives. The fault does not store anything; it is a surface. The elastic strain is in the surrounding rock, which is why geodetic networks tens of kilometres wide can watch a fault load, and why the width of the deforming zone tells you the locking depth.

Remember: an earthquake is not the fault breaking. The fault is already there. It is the fault ceasing to hold, and the surrounding rock springing back.

The one equation that matters

The size of an earthquake, measured physically rather than by how much a needle wiggled, is its seismic moment:

M0 = mu A D,

where mu is the shear modulus of the rock, A is the ruptured area, and D is the average slip. Each factor is a physical quantity you could in principle measure with a tape and a laboratory sample. Take the 1906 rupture.

Area first: A = (4.77 x 105)(1.5 x 104) = 7.155 x 109 m². Then

M0 = (3.0 x 1010)(7.155 x 109)(4.5) = 9.66 x 1020 N m.

Now convert. Moment magnitude, defined by Hiroo Kanamori in 1977 and put in this form by Thomas Hanks and Kanamori in 1979, is

Mw = (2/3)(log10 M0 - 9.1), with M₀ in newton metres.

log10(9.66 x 1020) = 20.985, so Mw = (2/3)(11.885) = 7.92.

The catalogue value for 1906 is 7.9. Four numbers, one logarithm, and you have reproduced it.

Test the sensitivity, because that is what tells you how much to trust a magnitude. Use W = 12 km and D = 4.0 m instead, both defensible readings of the same field data: A = 5.724 x 109, M0 = 6.87 x 1020, and Mw = (2/3)(20.837 - 9.1) = 7.82. The moment changed by 40 percent and the magnitude moved by 0.1. Magnitude is a logarithm of a logarithm's worth of physical detail, which is why arguments about the second decimal place are pointless and why a difference of 0.3 is a factor of 2.8 in energy.

How hard the rock was squeezed: stress drop

Moment tells you how much slip happened over how much area. It does not tell you how much stress was released, and that is a separate and more revealing quantity. For a rupture much longer than it is wide, which describes 1906 exactly, the stress drop is approximately

Delta sigma = (2/pi) mu D/W.

With the numbers above,

Delta sigma = (0.6366)(3.0 x 1010)(4.5)/(1.5 x 104) = 5.7 x 106 Pa = 5.7 MPa.

Fifty-seven bar. About the pressure in a scuba tank at half charge. Now here is the fact that organises all of earthquake physics: stress drops measured for earthquakes from magnitude 2 to magnitude 9 almost all fall between 1 and 10 MPa. Nine orders of magnitude in moment, one order of magnitude in stress drop.

That near-constancy is called self-similarity, and it has a consequence you can derive. If Delta sigma is fixed, then D scales with the rupture dimension, and since A scales with the dimension squared, M0 = mu A D scales with the dimension cubed. Double the length of a rupture and you multiply its moment by eight, which is 0.6 of a magnitude unit. Every scaling relation between magnitude and rupture length in engineering seismology descends from this.

The stress drop also raises a puzzle worth stating. If earthquakes release only 1 to 10 MPa, but the frictional strength of rock at 10 km depth under lithostatic load should be of order 100 MPa, then faults are either much weaker than laboratory friction predicts, or they release only a small fraction of the stress they carry. The absence of a measurable heat flow anomaly along the San Andreas, which a strong fault sliding at 100 MPa would produce, favours the first, and the mechanism is still argued about.

Why the old magnitudes stop counting

Charles Richter's 1935 scale was a practical device: read the largest amplitude on a Wood-Anderson torsion seismometer, correct for distance, take the logarithm. It works because most small earthquakes radiate their energy near the instrument's natural period of 0.8 seconds. The surface-wave magnitude M_s does the same at 20 seconds.

Both fail for large events, and the reason is a period. A rupture of length L propagating at v_r radiates most strongly at frequencies below the corner frequency

fc = vr/L.

For 1906, fc = 3.0/477 = 0.0063 Hz, a period of 159 seconds. The 20 second waves that M_s measures are far above the corner, on the part of the spectrum that no longer grows as the earthquake gets bigger. So M_s saturates near 8.2 and M_L near 6.5 to 7, and the 1960 Chilean earthquake and the 1964 Alaskan earthquake both came out as M_s 8.3 to 8.5 despite differing in moment by a factor of three.

The upshot: M_w does not saturate because it is not an amplitude at a chosen period. It is the low-frequency limit of the spectrum, which is M0 itself, and that limit keeps rising however long the rupture takes.

Closing the loop with the plate rate

One last piece of arithmetic ties this lesson to Lesson 12. If 1906 released 4.5 m of slip and the fault is loaded at the San Andreas slip rate of about 34 mm per year, the time to reload is

4.5/0.034 = 132 years.

Paleoseismic trenching at Vedanta Marsh and elsewhere on the northern San Andreas gives recurrence intervals of roughly 200 to 250 years for events of this size. The model under-predicts, and the gap is informative rather than embarrassing: not every rupture releases the same slip, some plate motion is taken up on the Hayward and Calaveras faults and in distributed deformation, and part of the loading is released aseismically by creep. A simple recurrence estimate is a lower bound on the interval and should be quoted as one.

Common misconceptions

• "The fault stores the energy." The strain energy is elastic, in the rock on both sides, over a volume tens of kilometres wide. That is why geodetic networks can watch a fault load, and why the width of the deforming zone measures the locking depth.
• "Each magnitude unit is ten times bigger." Ten times the amplitude, but 101.5 = 31.6 times the energy, since moment goes as the three halves power. Two units is a factor of 1000 in energy.
• "Richter magnitude is what news reports mean." Almost never, above magnitude 4. Richter's scale saturates near 6.5 and has not been used for large earthquakes since the 1970s; reported values are moment magnitudes.
• "A bigger stress drop means a bigger earthquake." Stress drop is nearly independent of size, between 1 and 10 MPa from magnitude 2 to magnitude 9. What makes an earthquake big is the area that slips, not how hard each square metre was squeezed.

Putting it together

Reid's three triangulation surveys showed distant marks moving steadily while near-fault marks lagged and then jumped, which is elastic rebound: strain stored in the rock around a locked fault, released when friction fails. Seismic moment M0 = mu A D for 1906, with L = 477 km, W = 15 km, D = 4.5 m and mu = 3.0 x 10¹⁰ Pa, gives 9.66 x 10²⁰ N m, and Mw = (2/3)(log10M0 - 9.1) = 7.92 against a catalogue 7.9; changing W to 12 km and D to 4.0 m moves the moment by 40 percent and the magnitude by 0.1. The stress drop (2/pi) mu D/W is 5.7 MPa, within the 1 to 10 MPa band that holds from magnitude 2 to magnitude 9, and that self-similarity is why moment scales as the cube of rupture dimension. The corner period of 159 seconds explains why M_s at 20 seconds saturates near 8.2 and M_w does not. Reloading 4.5 m at 34 mm per year takes 132 years, against a paleoseismic 200 to 250, and the difference is slip taken up elsewhere.

All of which describes one earthquake after it has happened. Whether it could have been forecast beforehand is the next lesson, and the answer is not encouraging.

Sources

1. Wikipedia contributors. (n.d.). Moment magnitude scale. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Elastic-rebound theory. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). 1906 San Francisco earthquake. Wikipedia. en.wikipedia.org
4. Hanks, T. C., and Kanamori, H. (1979). A moment magnitude scale. Journal of Geophysical Research, 84(B5), 2348-2350.
5. Stein, S., and Wysession, M. (2003). An introduction to seismology, earthquakes and Earth structure, Chapter 4. Blackwell.

Parkfield, California, population 18, sits on a 25 km stretch of the San Andreas that produced moderate earthquakes in 1857, 1881, 1901, 1922, 1934 and 1966. The intervals are 24, 20, 21, 12 and 32 years, averaging 21.8. In 1985 William Bakun and Allan Lindh published, in Science, a formal prediction: a magnitude 5.5 to 6 earthquake on that segment between 1988 and 1993, at 95 percent confidence. The US Geological Survey instrumented the ground more densely than any other place on the planet and waited.

The earthquake arrived on 28 September 2004, magnitude 6.0, eleven years outside the window. It also ruptured from southeast to northwest, the opposite direction to 1966.

This lesson starts from that failure and traces it back to the exact assumption that broke, because two of the three statistical laws involved are excellent and only the third was the problem.

The law that works: how many, of what size

Count earthquakes in any region over any long enough interval and the distribution of sizes obeys the Gutenberg-Richter relation:

log10 N = a - b M,

where N is the number of events of magnitude M or larger, a measures how seismically productive the region is, and b measures the relative proportion of large to small. Globally b is close to 1.0, which is the most robust number in the whole of earthquake statistics.

Check it on the global catalogue. Roughly 1300 earthquakes of magnitude 5 or more occur per year, about 130 of magnitude 6 or more, and about 15 of magnitude 7 or more. Each unit of magnitude divides the count by about ten, so b = 1.0, and from log10(1300) = a - 5 we get a = 8.11 for the whole Earth. Use it to predict: N(M >= 8) = 108.11 - 8 = 1.3 per year, and the observed long-run average is about one. Use it for a region with a = 5.0 and b = 1.0 and you get 10 events per year at magnitude 4 and up, one every ten years at magnitude 6 and up, one per century at magnitude 7.

Now combine it with Lesson 14's energy scaling, because the result is counterintuitive and important. Energy goes as 101.5M while frequency goes as 10-M, so energy released per magnitude bin goes as 100.5M: each bin releases 100.5 = 3.16 times as much as the bin below. Sum the geometric series downward and every earthquake smaller than the largest one, added together, releases only about 46 percent of what that largest one released. Seismic energy budgets are dominated by their biggest event, and no amount of small-earthquake activity relieves the strain that a great earthquake will take.

One practical matter that Lesson 17 will need. A real catalogue does not follow the straight line down to zero magnitude, because small events are missed by sparse networks. The magnitude at which the plot departs from linearity is the magnitude of completeness, M_c, and fitting b below it gives a wrong answer every time.

The core of it: Gutenberg-Richter tells you the rate of earthquakes of each size over the long run. It says nothing whatever about when.

The second law that works: after the fact

Fusakichi Omori, studying aftershocks of the 1891 Mino-Owari earthquake, found in 1894 that their rate decays as the reciprocal of time. In the form Tokuji Utsu gave it in 1961, Omori's law is

n(t) = K/(t + c)p,

with p usually between 0.9 and 1.4 and c a small constant covering the first hours, when the catalogue is incomplete because arrivals overlap.

Work it. Suppose 100 aftershocks above the completeness threshold occur on the first day, and take p = 1, c negligible. Then K = 100 per day at t = 1 day, so on day 10 the rate is 10 per day and on day 100 it is 1 per day. The cumulative number between day 1 and day 100 is

N = K ln(100/1) = 100 x 4.605 = 461.

A logarithm, which means the sequence has no natural end: aftershocks of the 1811 to 1812 New Madrid earthquakes are arguably still being recorded. Add Bath's law, that the largest aftershock averages 1.2 magnitude units below the main shock, and you can tell a community after a magnitude 7 that their most likely worst aftershock is around magnitude 5.8 and that the rate will halve roughly each time the elapsed time doubles. That is a genuine, useful, quantitative forecast.

The model that broke

The Parkfield prediction rested on neither of those. It rested on the characteristic earthquake model plus a time-predictable assumption, which together say: this fault segment repeatedly produces the same earthquake, it fails when the accumulated stress reaches a fixed threshold, and it reloads at a constant rate. Given those, the next event is due one average interval after the last.

Take the three assumptions separately and test each against what happened.

• The same earthquake each time. The 1966 rupture nucleated at the northwest end of the segment and ran southeast. The 2004 rupture nucleated at the southeast end and ran northwest. Same segment, same magnitude, opposite rupture. Whatever is characteristic about this piece of fault, the rupture process is not.
• A fixed failure threshold. Faults are not uniform. Strength varies along strike with the distribution of asperities, fluid pressure and gouge, and it changes after every rupture. There is no physical reason for a heterogeneous surface to fail at the same integrated stress twice.
• Constant loading. This is the assumption that fails most concretely. In 1983 the magnitude 6.4 Coalinga earthquake and the magnitude 6.0 Nunez event broke 25 km northeast of Parkfield and changed the stress on the Parkfield segment by a few tenths of a megapascal. Against a stress drop of a few megapascals, that is a perturbation of order ten percent, capable of moving a 22 year cycle by years. Neighbouring faults talk to each other, and a single-segment model cannot hear them.

Look again at the intervals with those three failures in mind: 24, 20, 21, 12, 32. The coefficient of variation is about 0.35, which is not a clock. The 12 year gap from 1922 to 1934 and the 32 year gap from 1934 to 1966 differ by a factor of 2.7, and the prediction window was built on their mean as though the scatter were measurement noise rather than physics.

What the experiment did produce

Calling Parkfield a failure without qualification would be a second mistake. The instruments were in place when the earthquake finally came, and the resulting dataset is the most complete record of a moderate earthquake ever assembled: borehole strainmeters within a few kilometres of the rupture, creepmeters across the surface trace, dense seismic arrays, and the San Andreas Fault Observatory at Depth, which between 2004 and 2007 drilled through the fault at 3.2 km and brought up core from an actively creeping strand.

The headline result was negative, and it is worth stating plainly. No short-term precursor was observed. Not anomalous strain, not foreshock patterns distinguishable from ordinary seismicity, not radon, not electrical or magnetic signals. Strainmeters within 10 km of the nucleation point recorded nothing above their noise in the hours before. If a precursor of the kind decades of proposals had assumed exists, it did not appear here, on the best-monitored fault segment in the world, at the moment it was watched for.

Bottom line: Parkfield did not show that prediction is hard. It showed that the specific model on which deterministic prediction rested, a fault that repeats itself on a schedule and announces itself beforehand, is not how this fault behaves.

What replaced prediction

Three things took its place, and none of them is prediction in the original sense.

• Probabilistic seismic hazard analysis. Combine Gutenberg-Richter rates for every mapped fault and background zone with ground-motion models, and produce the probability of exceeding a given shaking level in 50 years. That is what building codes use, and it needs no knowledge of when.
• Operational earthquake forecasting. Omori plus Gutenberg-Richter, implemented as epidemic-type aftershock sequence models, gives a running probability that updates with every event. After a magnitude 6 the chance of a larger shock in the next week is a few percent rather than the background one in ten thousand, and saying so honestly is useful.
• Earthquake early warning. Not prediction at all: detection. A network reads the P wave, estimates location and magnitude in a few seconds, and alerts places the damaging S wave has not yet reached. At 100 km epicentral distance the S wave takes 100/3.5 = 28.6 seconds, and with about 7 seconds for detection and processing that leaves roughly 21 seconds of warning: enough to stop a train, close a valve, or move away from a window. Directly beneath the epicentre it leaves nothing, which is the system's permanent limitation.

Common misconceptions

• "Small earthquakes relieve stress and prevent large ones." With b = 1, each magnitude bin releases 3.16 times the energy of the one below, so everything smaller than the largest event contributes only about 46 percent of what that event released. A thousand magnitude 4s do not substitute for one magnitude 7.
• "Aftershocks stop after a few weeks." Omori's law decays as 1/t, whose integral is a logarithm, so the sequence has no natural end. What ends is the point at which the aftershock rate drops below the background rate for the region.
• "Parkfield showed that earthquakes are unpredictable in principle." It showed that one specific model failed on one specific segment. The stronger claim, that no precursor of any kind exists, is not established, though thirty years of monitoring have made large reliable precursors much less likely.
• "Early warning predicts earthquakes." It measures one that has already started, and races the S wave using the difference between 6 km/s and 3.5 km/s plus the speed of electronics. There is no warning at all directly above the epicentre.

Where this leaves us

Gutenberg-Richter, log10N = a - bM with b near 1.0, reproduces the global rates of 1300, 130 and 15 events per year above magnitudes 5, 6 and 7, and combined with the 101.5M energy scaling shows that every event below the largest contributes only 46 percent of its energy. Omori's law, n(t) = K/(t + c)p with p near 1, turns 100 aftershocks on day one into 461 by day 100 and never quite reaches zero. Both are excellent. The Parkfield prediction used neither: it used a characteristic earthquake recurring on a schedule, and all three of its assumptions failed, most clearly when Coalinga in 1983 perturbed the segment's stress by a tenth of a megapascal. The earthquake came in 2004, eleven years late and running the other way, with no detectable precursor on the best-instrumented fault on Earth. What survives is probabilistic hazard, aftershock forecasting, and early warning that buys about 21 seconds at 100 km.

The last module leaves the present entirely, for a clock that runs on lead isotopes, and then returns to put a real catalogue in your hands.

Sources

1. Wikipedia contributors. (n.d.). Gutenberg-Richter law. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Omori's law. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Earthquake prediction. Wikipedia. en.wikipedia.org
4. Bakun, W. H., and Lindh, A. G. (1985). The Parkfield, California, earthquake prediction experiment. Science, 229(4714), 619-624.
5. Stein, S., and Wysession, M. (2003). An introduction to seismology, earthquakes and Earth structure, Chapter 4. Blackwell.

Clair Patterson arrived at the University of Chicago in 1948 as a graduate student and was handed what sounded like a contained problem: measure the lead in a single crystal of zircon and, from the uranium that had decayed into it, calculate the crystal's age. He spent the next five years failing. Every measurement gave a different answer, and all of them were too high. The lead he was detecting was not coming out of his sample.

That failure is the subject of the second half of this lesson, because it turned into one of the most consequential public health results of the twentieth century. First, the clock itself.

Two uranium isotopes, two clocks, one mineral

Uranium-lead dating is the most powerful geochronometer available, for a reason that is easy to state and hard to overvalue: uranium has two long-lived isotopes with very different half-lives, and both end up as lead.

²⁰⁴Pb has no radioactive parent, so it serves as the reference against which radiogenic additions are measured. The present-day ratio ²³⁸U/²³⁵U is 137.88 in every terrestrial and meteoritic sample, which is what makes the two clocks comparable.

For each system the accumulated daughter relative to surviving parent is

206Pb*/238U = e(lambda238 t) - 1 and 207Pb*/235U = e(lambda235 t) - 1,

where the asterisk means radiogenic. Two equations, one unknown, which is one equation more than you need. That redundancy is the whole method.

The concordia curve

Plot the first ratio against the second and sweep t from zero upward. The locus traced out is concordia, and any mineral that has kept every atom of lead it ever made must lie exactly on it.

Check the last row by hand. lambda238 t = (1.55125 x 10-10)(4.55 x 109) = 0.70582, and e0.70582 = 2.0256, so the ratio is 1.0256. For the other, lambda235 t = (9.8485 x 10-10)(4.55 x 109) = 4.4811, and e4.4811 = 88.35, giving 87.35. Notice how much faster the ²³⁵U clock runs: in 4.55 billion years the ²³⁸U system roughly doubles its parent-normalised daughter while the ²³⁵U system multiplies it by eighty-eight.

Now run the inverse problem, which is what an analysis actually gives you. A zircon returns 206Pb*/238U = 0.5500 and 207Pb*/235U = 14.50. Two ages:

t206 = ln(1.5500)/1.55125 x 10-10 = 0.43825/1.55125 x 10-10 = 2826 Ma,

t207 = ln(15.500)/9.8485 x 10-10 = 2.74084/9.8485 x 10-10 = 2783 Ma.

They differ by 43 million years, 1.5 percent. The point lies slightly below concordia, and the crystal is discordant.

What matters here: discordance is information, not failure. A single clock that gives one number cannot tell you whether to believe it. Two clocks that disagree tell you the sample has lost lead, and by how much.

George Wetherill showed in 1956 what to do with it. Zircons from one rock that all crystallised at the same time and all lost lead during the same later event fall on a straight line in the concordia diagram, a discordia. That line cuts concordia twice: the upper intersection gives the original crystallisation age, the lower gives the age of the disturbance. A set of damaged samples therefore yields two dates instead of none, which is why U-Pb geochronology survives in rocks that have been metamorphosed repeatedly.

Dating a planet you cannot sample

There is a catch that took thirty years to see. The Earth's crust has been melted, mixed and recycled continuously, so no terrestrial rock records the moment of accretion. The oldest known terrestrial mineral, a zircon from the Jack Hills of Western Australia, is 4.404 billion years old; the oldest intact rock, the Acasta gneiss, is 4.03. Both are younger than the planet.

The way around it is to date something that never got recycled, and then show that the Earth belongs to the same system. That is what the lead-lead isochron does.

Plot ²⁰⁷Pb/²⁰⁴Pb against ²⁰⁶Pb/²⁰⁴Pb for a set of objects that formed at the same time from the same starting lead but with different uranium contents. All of them begin at the same point, the primordial lead composition, and each moves away from it at a rate set by its own U/Pb ratio. After time t they lie on a straight line whose slope is

slope = (e(lambda235 t) - 1)/[137.88 (e(lambda238 t) - 1)].

At t = 4.55 Ga that is 87.35/(137.88 x 1.0256) = 87.35/141.41 = 0.6177. Measure a slope, invert this expression, and you have an age without knowing any sample's uranium concentration at all, which is why the method survives even when uranium has been mobilised.

Patterson published the result in 1956 using five meteorites: two iron and three stony, chosen so their U/Pb ratios spanned a wide range. The troilite, an iron sulphide, in the Canyon Diablo iron meteorite contains almost no uranium at all, so its lead has barely evolved since it formed; it defines the primordial composition at ²⁰⁶Pb/²⁰⁴Pb = 9.307 and ²⁰⁷Pb/²⁰⁴Pb = 10.294 and anchors the low end of the line. The age came out at 4.55 billion years, with an uncertainty of about 70 million.

Then the second half of the argument, which is the part usually left out. Patterson also analysed modern ocean sediment, which averages the lead of the whole continental crust, and found that it plots on the same isochron. Earth and meteorites share one starting composition and one age. Without that step the 4.55 would date meteorites and nothing else.

Modern refinements have moved the number very little. Calcium-aluminium inclusions in the Allende meteorite, the oldest solids in the solar system, date to 4.567 Ga; Earth's accretion finished around 4.50 Ga; the Moon-forming impact is placed near 4.51. Patterson's figure, from a mass spectrometer built partly by hand, sits within a percent of all of them.

The lead that would not go away

Return to the five years of failure. Patterson eventually established that the excess lead in his samples was atmospheric and laboratory contamination, at levels hundreds of times higher than the lead he was trying to measure. His response was to rebuild the laboratory: acid-leached everything, distilled his own reagents, filtered the air, and created what amounted to the first ultraclean laboratory. Blank levels fell by orders of magnitude, and the age measurement became possible.

Then he asked why there was so much lead in the first place. Natural weathering could not supply it. He measured lead in ocean water and found the surface enriched relative to deep water, which is the signature of a source from above rather than from rivers. He measured Greenland ice, layer by layer, with Mitsunobu Murozumi and Tsaihwa Chow, and published in 1969 a profile in which lead concentration rises from roughly 0.001 micrograms per kilogram in pre-industrial ice to about 0.2 by the mid-1960s, a rise of more than a hundredfold, with the steepest part beginning in the 1920s.

The 1920s are when tetraethyl lead was introduced as a petrol additive.

In 1965 Patterson published Contaminated and natural lead environments of man, arguing that people then alive carried roughly a hundred times the natural body burden and that the accepted safe levels had been set by comparing contaminated people with other contaminated people. The industry response was sustained: he lost contracts, was kept off a National Research Council panel on atmospheric lead in 1971, and was attacked at length in the literature. The measurements held. Leaded petrol was phased out in the United States from 1975, and mean blood lead in the population fell from about 12.8 micrograms per decilitre in the late 1970s to 2.8 by 1991, a drop of nearly 80 percent that tracks the removal almost exactly.

So what?: the contamination that made the age of the Earth hard to measure was itself a measurement. Patterson could tell it was anthropogenic because he knew what the natural level had to be, and he knew that because he had spent five years failing to reach it.

Common misconceptions

• "Radiometric dating assumes no daughter was present initially." The isochron method assumes nothing of the sort. It solves for the initial composition as the intercept, which is exactly why Canyon Diablo troilite, with almost no uranium, is so valuable: it pins that intercept directly.
• "A discordant zircon is a failed analysis." A discordia line through several discordant grains gives two ages, crystallisation at the upper intersection with concordia and disturbance at the lower. The disagreement between the two clocks is the signal.
• "4.55 billion years is the age of the oldest rock on Earth." The oldest terrestrial mineral is a 4.404 Ga Jack Hills zircon and the oldest intact rock is the 4.03 Ga Acasta gneiss. The 4.55 comes from meteorites, and applies to the Earth only because terrestrial lead falls on the same isochron.
• "Patterson's environmental work was a separate career." It came directly out of the dating problem. Only someone who had established the natural background to a part in a thousand could demonstrate that the modern burden was a hundred times it.

What you now know

Two uranium isotopes decaying to two lead isotopes give two independent clocks in one crystal, with lambda238 = 1.55125 x 10-10 and lambda235 = 9.8485 x 10-10 per year. Concordia is the locus where both agree, running from the origin to 87.35 and 1.0256 at 4.55 Ga. A zircon at 0.5500 and 14.50 gives 2826 and 2783 Ma, discordant by 1.5 percent, and a discordia line through several such grains recovers both the crystallisation age and the age of the later disturbance. Because no terrestrial rock survives from accretion, the age of the Earth comes from a lead-lead isochron whose slope at 4.55 Ga is 0.6177, anchored at the primordial end by Canyon Diablo troilite at ²⁰⁶Pb/²⁰⁴Pb = 9.307, and applied to the Earth because ocean sediment plots on the same line. Patterson's 4.55 Ga still stands against a solar system age of 4.567 Ga from Allende inclusions. And the contamination that delayed the measurement for five years turned into the Greenland ice profile, a hundredfold rise in atmospheric lead beginning in the 1920s, and a fall in mean blood lead from 12.8 to 2.8 micrograms per decilitre once the petrol additive was removed.

One lesson remains, and it hands you the data instead of the conclusion.

Sources

1. Wikipedia contributors. (n.d.). Uranium-lead dating. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Clair Cameron Patterson. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Age of the Earth. Wikipedia. en.wikipedia.org
4. Patterson, C. (1956). Age of meteorites and the Earth. Geochimica et Cosmochimica Acta, 10(4), 230-237.
5. Dalrymple, G. B. (1991). The age of the Earth. Stanford University Press.

A query to the USGS earthquake catalogue search comes back as a comma-separated file whose header begins time,latitude,longitude,depth,mag,magType and runs to twenty-two fields in all. None of those fields is a b-value. Getting one out of the file takes five decisions, and four of them can be made badly without anything on your screen looking wrong.

So this lesson walks the path twice: once on a teaching catalogue of 860 events, binned in steps of 0.2 and built so the right answer is known in advance, and again with one decision changed. That change moves the estimated recurrence of a magnitude 6 earthquake from about once a century to about once a decade. Everything apart from the events themselves is the real machinery: the column names, the estimators, the windows and the arithmetic.

Decision one: which rows are earthquakes

The catalogue of the Advanced National Seismic System does not hold only earthquakes. Its type column also carries quarry blasts and chemical explosions, which in a mining district can outnumber natural events at small magnitudes and are not drawn from the same size distribution. Filter on type before you read a single magnitude, then de-duplicate on the event identifier rather than the origin time, since two networks reporting one shock place it seconds apart.

Decision two: which magnitude scale you are reading

The magType column says which scale produced the number beside it: ml for local magnitude, md for duration, mb for the short-period body wave, mw for moment. These are not interchangeable, and the difference is not a constant offset.

The Gutenberg-Richter relation log10 N = a - b M makes b a slope in magnitude units. Suppose your scale is compressed relative to moment magnitude, so the reported value is M' = k M with k below 1. Substituting gives log10 N = a - (b/k) M': the slope you fit is not b but b/k. A true b of 1.01 read on a scale with k = 0.9 comes back as 1.12, an 11 percent error before any statistics are done. A uniform offset is harmless by the same algebra, since adding a constant moves a and leaves b alone. Stretching hurts; shifting does not.

Worth holding on to: b belongs to a scale as much as to a region, so two b-values from different magnitude types are not comparable, and a catalogue that switches type partway through has a break in it.

Decision three: where the catalogue stops being complete

Bin the surviving magnitudes. Here are the 860 events in bins 0.2 wide, with the count in each bin and the cumulative count at or above it, which is the N in the Gutenberg-Richter relation.

Read the middle column, not the cumulative one. A distribution obeying Gutenberg-Richter rises without limit as magnitude falls, since there are always more small events than large. This one rises to 200 in the bin centred on 2.5 and then turns over: 155, 95, 48, 22. Nothing in the Earth does that. What falls away below 2.5 is the network, not the seismicity, and the turnover is therefore the magnitude of completeness, M_c, that Lesson 15 promised you would need.

Taking that peak is the maximum curvature method, the quickest criterion that works. It has a known bias: where the roll-off is gradual the peak sits slightly below the magnitude at which the catalogue is genuinely complete, so adding 0.2 is common practice. The alternative, from Wiemer and Wyss, fits a line above each trial M_c in turn and takes the lowest whose synthetic counts reproduce 90 or 95 percent of the observed cumulative numbers. It is slower, and on a catalogue this small it returns the same 2.5.

Decision four: which estimator

The obvious move is least squares on the logarithm of the cumulative column. Do not. Each cumulative count contains every count above it, so the points are not independent, and the bins near M_c carry hundreds of events while those near the top carry one or two, which least squares weights alike. Fitted from 2.5 to 4.3 it gives 1.07.

Use maximum likelihood instead. Magnitudes above M_c are exponentially distributed, and Aki showed in 1965 that the likelihood is maximised by

b = log10(e) / (Mmean - (Mc - dM/2)),

where M_mean is the mean magnitude of events at or above M_c and dM is the bin width. That correction is not decoration: binned magnitudes sit at bin centres, so the sample's true lower edge is half a bin below the centre you called M_c.

Work it. The 540 events at or above 2.5 contribute 200 x 2.5 = 500.0, then 126 x 2.7 = 340.2, 80 x 2.9 = 232.0, 50 x 3.1 = 155.0, 32 x 3.3 = 105.6, 20 x 3.5 = 70.0, 13 x 3.7 = 48.1, 8 x 3.9 = 31.2, 5 x 4.1 = 20.5, 3 x 4.3 = 12.9, 2 x 4.5 = 9.0 and 1 x 4.7 = 4.7. The sum is 1529.2, so

Mmean = 1529.2/540 = 2.8319.

The corrected lower edge is 2.5 - 0.1 = 2.4, the denominator is 2.8319 - 2.4 = 0.4319, and with log₁₀(e) = 0.434294,

b = 0.434294/0.4319 = 1.01.

The catalogue was built with b = 1.0, so the procedure recovered the planted value. Now the uncertainty, because a b-value without one is not a measurement. Aki's 1965 paper gives the standard error as b/sqrt(n), and sqrt(540) = 23.24, so

sigmab = 1.01/23.24 = 0.043.

That sets what you may claim: a b of 1.01 plus or minus 0.043 cannot be told apart from the global average of 1.0, nor from a neighbour that came out at 0.95. Reaching sigma = 0.02 would need (1.01/0.02)2 = 2500 events above completeness, which for this network is a century of recording.

Drop the half-bin correction and you divide by 2.8319 - 2.5 = 0.3319 instead, giving 1.31: 30 percent high and far outside the error bar. A term that looks like bookkeeping separates a right answer from a confident wrong one.

The upshot: maximum likelihood, the half-bin correction and the standard error are three lines of arithmetic that decide whether your b-value means anything.

Decision five: declustering, and what it is for

Gutenberg-Richter describes independent events, and a catalogue is not that: every aftershock sequence is a burst of dependent events clustered in space and time. The commonest remedy is the window method of Gardner and Knopoff, who in 1974 asked whether southern California seismicity with aftershocks removed is Poissonian. Around each event they delete every smaller event inside windows that grow with mainshock magnitude:

d = 10(0.1238M + 0.983) km, and t = 10(0.5409M - 0.547) days for M below 6.5.

A magnitude 3.0 clears 10^1.354 = 22.6 km for 10^1.076 = 11.9 days; a magnitude 5.0 clears 40.0 km for 143.7 days; our largest event, at 4.7, clears 36.7 km for 98.9 days. The windows are generous on purpose, which is the method's main criticism: it deletes independent events with the dependent ones.

Suppose declustering removes 188 of the 540 events above completeness. The remaining 352 over 20 years give 17.6 per year rather than 27, so a falls by log10(540/352) = 0.19 while the slope barely moves, since aftershocks follow much the same size distribution as the shocks that triggered them. Declustering is mostly a statement about rate.

From two fitted numbers to a sentence someone can use

The full catalogue holds 540 events above 2.5 in 20 years, or 27 per year, and the intercept follows from that rate at the completeness magnitude:

a = log10(27) + 1.0055 x 2.5 = 1.4314 + 2.5138 = 3.945.

Now predict. For magnitude 5 and above, N = 10(3.945 - 5.028) = 0.083 per year, a return period of 12 years; for magnitude 6, N = 10(3.945 - 6.033) = 0.0082 per year, or 123 years. Convert that to a probability with a Poisson model, which is what declustering was for: the chance of at least one in 50 years is 1 - e-0.41 = 0.33. Rate times interval gives 0.41 and overstates it by a quarter.

One caveat belongs in the same breath. The largest event here is 4.7, and the fit says a magnitude 7 every 1240 years: an extrapolation 2.3 magnitude units beyond anything observed, which no row of data tests. Real hazard models bound the top of the distribution with fault dimensions and palaeoseismology instead.

Change one decision and watch it move

Suppose you skipped the incremental plot and took M_c = 2.1 from the network's advertised detection threshold, the commonest way this analysis goes wrong. The arithmetic runs identically: 790 events at or above 2.1, mean magnitude 2.6395, corrected edge 2.0, and

b = 0.434294/0.6395 = 0.68.

The rate is 790/20 = 39.5 per year, so a = 1.5966 + 0.6791 x 2.1 = 3.023, and for magnitude 6 N = 10(3.023 - 4.075) = 0.089 per year: a return period of 11 years against the 123 you got by doing it properly. That factor of eleven is not noise. Including the incomplete bins flattens the distribution, a flatter distribution means a smaller b, and a smaller b puts far more weight on the large events at the extrapolated end.

In short: the magnitude of completeness is not a tidying step before the analysis. It is the analysis. Every quantity downstream inherits it.

Common misconceptions

• "More data is always better, so fit from the smallest magnitude available." The opposite, below M_c. Adding the incomplete bins took b from 1.01 to 0.68 and the magnitude 6 return period from 123 years to 11. Events you never recorded are not zeros; they are missing rows, and a fit cannot tell the difference.
• "A straight-looking cumulative plot means the catalogue is complete." Cumulative counts are smooth by construction, since each point contains those above it, and the roll-off at the bottom looks like curvature you can fit through. Incompleteness shows in the incremental distribution.
• "b = 1.0 and b = 0.95 are different b-values." Only an error bar can say. On 540 events the standard error is 0.043, so those two are one number; separating them at one sigma needs several thousand events above completeness.
• "Declustering shrinks the catalogue, so it shrinks the hazard." It lowers the rate, which is why a fell by 0.19 here, but hazard models use declustered rates precisely because the Poisson arithmetic behind a 50 year probability requires independence. Skipping it does not make a region safer; it makes the probability wrong.

Looking back

Five decisions stand between a downloaded file and a defensible number, and the file records none of them. Filter the non-tectonic types, de-duplicate on the event identifier, and check that one magnitude scale runs the whole period. Find M_c at the peak of the incremental distribution, 2.5 here, where the counts turn over from 200 to 155 to 95. Estimate the slope by maximum likelihood with the half-bin correction, b = 0.434294/(2.8319 - 2.4) = 1.01, and attach sigmab = b/sqrt(540) = 0.043. Decluster, and expect the rate to fall while the slope holds. Only then convert: a = 3.945, one magnitude 6 every 123 years, a 33 percent chance of one in 50 years. Take M_c = 2.1 instead and that becomes 11 years.

Which is where the course ends, in the place it began. A pendulum clock at Cayenne ran slow, and the interesting question was never the clock but what the number meant. Everything since has been that same move: a gravity reading corrected four times, a travel-time curve crossing at 262 km, a density profile that failed at 410 km, magnetic stripes counted outward from a ridge, a slope of 0.6177 that dates the planet at 4.55 billion years. The Earth does not hand over its interior. It hands over numbers, and what you may say about the inside depends entirely on how carefully you read them.

Sources

1. U.S. Geological Survey. (n.d.). ComCat documentation: Comprehensive catalog of earthquake events and products. USGS Earthquake Hazards Program. earthquake.usgs.gov
2. U.S. Geological Survey. (n.d.). FDSN event web service documentation. USGS Earthquake Hazards Program. earthquake.usgs.gov
3. Aki, K. (1965). Maximum likelihood estimate of b in the formula log N = a - bM and its confidence limits. Bulletin of the Earthquake Research Institute, 43, 237-239.
4. Wiemer, S., and Wyss, M. (2000). Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the western United States, and Japan. Bulletin of the Seismological Society of America, 90(4), 859-869.
5. Gardner, J. K., and Knopoff, L. (1974). Is the sequence of earthquakes in Southern California, with aftershocks removed, Poissonian? Bulletin of the Seismological Society of America, 64(5), 1363-1367.