Kelvin’s age of the Earth (1862)

In 1862 William Thomson, later Lord Kelvin, treated the Earth as a semi-infinite solid that began at a uniform molten temperature and has cooled from its surface ever since. The temperature below the surface is then \(T_s + (T_i - T_s)\operatorname{erf}(x/2\sqrt{\alpha t})\) (semi_infinite_temperature()), and the surface gradient \((T_i - T_s)/\sqrt{\pi\alpha t}\) falls with time. From the gradient measured in mines, about 1 °F per 50 ft, he found an age of about 98 million years. This example redoes his arithmetic in SI units and checks the gradient against the profile itself.

import matplotlib.pyplot as plt
import numpy as np

from engineeringkit.heat_transfer import semi_infinite_temperature

Kelvin’s numbers

An initial excess of 7000 °F, a gradient of 1 °F per 50 ft and a diffusivity of 400 ft² per year, converted to SI.

year = 365.25 * 86400
dT = 7000 * 5 / 9  # K
gradient = (5 / 9) / (50 * 0.3048)  # K/m
alpha = 400 * 0.3048**2 / year  # m^2/s
age = dT**2 / (np.pi * alpha * gradient**2)
print(f"diffusivity {alpha:.3e} m²/s, gradient {gradient * 1e3:.1f} mK/m")
print(f"age of the Earth: {age / year / 1e6:.0f} million years")
assert 95e6 < age / year < 100e6

diffusivity 1.178e-06 m²/s, gradient 36.5 mK/m
age of the Earth: 97 million years

The gradient from the profile

Differentiate the error-function profile numerically at the surface and compare with the closed form.

Ts, Ti = 273.0, 273.0 + dT
dx = 1.0
numeric = (semi_infinite_temperature(Ti, Ts, dx, alpha, age) - semi_infinite_temperature(Ti, Ts, 0.0, alpha, age)) / dx
print(f"surface gradient from the profile {numeric * 1e3:.3f} mK/m")
assert np.isclose(numeric, gradient, rtol=1e-4)

depth = np.linspace(0, 400e3, 400)
fig, ax = plt.subplots(figsize=(6, 4))
for myr in [1, 10, 98, 1000]:
    ax.plot(semi_infinite_temperature(Ti, Ts, depth, alpha, myr * 1e6 * year) - 273.0, depth / 1e3, label=f"{myr} Myr")
ax.invert_yaxis()
ax.set(xlabel="temperature [°C]", ylabel="depth [km]", title="Kelvin's cooling Earth")
ax.legend()
fig.tight_layout()

surface gradient from the profile 36.454 mK/m

Why he was wrong

Radioactive heating (discovered in 1896) and, above all, convection in the mantle, argued by John Perry in 1895, keep the interior hot. The radiometric age, 4.54 billion years, would need a surface gradient about \(\sqrt{4540/98}\) times smaller than the one measured.

print(f"the gradient for 4.54 Gyr would be {gradient * np.sqrt(age / (4.54e9 * year)) * 1e3:.1f} mK/m")

plt.show()

the gradient for 4.54 Gyr would be 5.3 mK/m