Lecture 34: On the age of the earth (2 lectures)

(previous, next)


References for further reading


Catastrophists vs. uniformitarians:

Early natural philosophy at the beginning of the European renaissance derived from Aristotilian and biblical sources.

James Ussher's biblical chronology (1581-1656) placed creation at the 23rd of October, 4004 BC (\(6 \times 10^3\) years)

Such theories had a hard time explaining observations from the natural world like seashells found on the tops of mountains and monstrous creatures that no longer seemed to exist.

The first philosophies attempting to explain these irregularities used recourse to catastrophic events. In particular, the biblical flood described by Noah was often invoked. However, the diversity and pattern of the events observed in layers of rock around Europe required many diverse catastrophies to reconciliation theory and observation.

Uniformitarians proposed that the world was and had always been in a state of slow and steady change, with the creation and erosion of mountains. From the uniformitarian perspective, the earth was essentially ageless.

A law for heat

Joseph Fourier (1768-1830) published an important work on the theory of heat in 1807. This became the basis of our modern understanding.

Hypothesize that the rate of heat-flow is proportional to heat gradient. Then the local change in heat is proportional to the heat gain minus heat loss.

\[\frac{\partial H}{\partial t} = D \frac{\partial^2 H}{\partial x^2}\] Green function \[H(x,t) = \frac{1}{\sqrt{4 \pi Dt }} e^{-x^2/4Dt}\] \[\sin(\omega x) e^{-a t}, \quad D \omega^2=a\]

Applications:

Kelvin

William Thomas believed that everything should have recourse to physical principles in science. In 1862, building on an idea of Buffon, he published his ideas on using the heat equation of Fourier to attempt to calculated the age of the earth. He derived an estimate of around 100 million years. While this seems to split the difference among existing theories of the time, it was much shorter than Darwinians and uniformitarians believed at the time.

A toy problem

Let's consider a toy-problem first, to get a feel for things. Suppose we have a long steel bar that was heated in a furnace and then removed and insulated so that the heat energy only escapes though the un-insulated ends. The bar has length \(L\) and thermal diffusivity \(\kappa\). We know the furnace had a temperature of a \(1,000\) degrees above room temperature, but the steel bar has been cooling down for a while now. We'd like to know how long, based on just the temperature gradiant at the end of the bar.

Let \(\theta(t,x)\) represent the excess temperature above room-temperature of location \(x\) allong the steel bar at time \(t\). The bar is \(1\) meter long, and has thermal diffusivity of \(\kappa = 4.2 \times 10^{-6}\) square meters per second. The current temperature gradient at the end of the bar is \(20\) degrees Celcius per meter.

\[\theta(0,x) = \theta _ 0, \quad \theta(t,0) = \theta(t,L) = 0, \quad \frac{d \theta}{d t} = \kappa \frac{d^2 \theta}{dx^2}\]

\[\theta(t,x) = \sum _ {n=1}^{\infty} a_n \sin \left( \frac{n\pi}{L} x \right) e^{ -\lambda_n t}\]

\[\theta(t,x) = \frac{4 \theta_0 }{\pi} \sum_{n=0}^{\infty} \frac{e^{- \frac{\kappa t}{L^{2}} \pi^{2} \left(2 n + 1\right)^{2}}}{\left(2 n + 1\right)} \sin{\left (\frac{\pi x}{L} \left(2 n + 1\right) \right )}\]

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Here is some simple code to check that our
# Fourier series terms actually solve our heat equation.
from sympy import *
var('x,n,L,kappa,t')
# on the interval [0,L]
f = 4/(2*n+1)/pi*sin((2*n+1)*pi/L*x)*exp(-kappa*pi**2*(2*n+1)**2/L**2*t)

# Let's check quickly that our Fourier terms solve the ODE and BCs
f.diff(t) - kappa * f.diff(x,x), f.subs(x,0), f.subs(x,L)
theta = Sum(f, (n,0,oo))

# for illustration that the series matches our expected initial
# conditions, here are the first 20 terms when L = 1 and kappa = 1.
# In a plot, you'll see strong Gibbs phenomena because of the initial
# discontinuity at the boundaries.
print Sum(f, (n,0,20)).subs(L,1).subs(kappa,1).doit().evalf()

The leading, slowest-decaying term in this solution is \[\theta(x,t) = \frac{4 \theta _ 0}{\pi} e^{- \frac{\kappa t}{L^{2}} \pi^{2}} \sin{\left (\frac{\pi x}{L} \right )}\]

Taking a derivative and solving for \(x\) and evaluating at the rod's end (\(x=0\)), we can determine \(t\), how long the bar has been cooling.

\[t = \frac{L^{2}}{\kappa \pi^{2}} \log{\left (\frac{4 \theta_0}{L \frac{\partial \theta(x=0,t)}{\partial x} } \right )}\]

So, in principle, at least, we can use the temperature gradient to approximate how long the bar has been out of the oven.

Calculation from spherical cooling

But more realistically, the earth is a sphere in 3 dimensions, so we have to work with a more complicated heat equation.

\[\frac{d \theta}{dt} = \kappa \left( \frac{d^2 \theta}{dx^2} + \frac{d^2 \theta}{dy^2} + \frac{d^2 \theta}{dz^2} \right)\]

In spherical coordinates \((r,\theta,\phi)\), the heat equation becomes (using subscripts to represent derivatives, for convenience)

\[\theta _ t = \kappa \frac{1}{r^2} \left( r^2 \theta_r \right)_r = \kappa \left( \frac{2}{r} \theta_r + \theta_{rr} \right)\]

with boundary conditions \(\theta(r _ {earth},t) = 0\) and \(\theta _ r(0,t) = 0\) by radial symmetry. Decompose into a spatial and a temporal part by separation of variables \(\theta(r,t) = X(r) Y(t)\), and we get \[\begin{gather} Y_t + \kappa \lambda^2 Y = 0, \\ X_{rr} + \frac{2}{r} X_r + \lambda^2 X = 0. \end{gather}\] The first of these has the usual solution \(Y(t) \propto e^{- \kappa \lambda^2 t}\). The second is trickier to deal with. It's a homogeneous variable coefficient second order ordinary differential equation (odd how those adjectives work better backwards). In general, these equations require special functions or series solutions. In our case where \(\lambda\) is unknown, we have an irregular Sturm--Louville problem. But we are actually lucky. It turns out the general homogeneous solution of this equation is \[\alpha \frac{\sin(\lambda r)}{r} + \beta \frac{\cos(\lambda r)}{r}\] where \(\alpha\) and \(\beta\) are free constants. Our central boundary condition of zero-derivative at the origin implies \(\beta=0\) Then \(X(r_{earth}) = 0\) means \(\lambda_n = n \pi / r_{earth}\) where \(n\) is a positive integer. Thus, our eigenfunction expansion of the solution is \[\begin{gather} \theta(r,t) = \sum_{n=1}^{\infty} \alpha_n \frac{\sin(\lambda_n r)}{r} e^{- \kappa \lambda_n^2 t}, \quad \lambda_n = \frac{ n \pi }{ r_{earth}}. \\ a_n = (-1)^{n+1} \frac{\theta_0 r_{earth}}{n \pi} \end{gather}\]

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
from sympy import *
from sympy.abc import R, kappa, r, t, n
L = lambda n : n*pi/R
a = lambda n : (-1)**(n+1)*R/n/pi
theta = Sum(a(n)*sin(L(n)*r)/r*exp(-kappa*L(n) ** 2*t),(n,1,oo))
pprint(theta)

# For illustration that the series matches our expected initial
# conditions, here are the first 20 terms when L = 1 and kappa = 1.
# In a plot, you'll see strong Gibbs phenomena because of the initial
# discontinuity at the boundaries.
theta = Sum(a(n)*sin(L(n)*r)/r*exp(-kappa*L(n) ** 2*t),(n,1,20))
print theta.subs(R,1).subs(kappa,1).doit().evalf()

Again, dropping all but the slowest-decaying term, \[\begin{gather} \theta(r,t) \approx \frac{\theta_0 r_{earth}}{\pi r} \sin\left(\frac{\pi r}{r_{earth}} \right) e^{- \kappa \pi^2 t / r_{earth}^2}, \end{gather}\] Differentiating and evaluating at the surface of the sphere, \[\theta_r(r_{earth},t) = - \frac{\theta_0}{r_{earth}} e^{- \kappa \pi^2 t / r_{earth}^2}.\] Solving, \[t \kappa = \left( \frac{r_{earth}}{\pi}\right)^2 \ln\left( \frac{\theta_0}{-\theta _ r(r _ {earth},t) \, r _ {earth}} \right)\]

Data and calculation

Symbol Value Description
\(\theta _ 0\) \(1,200\) Uniform initial temperature in Kelvin above surface
\(\theta_r(0,t)\) \(3.6 \times 10^{-2}\) Temperature gradient in the earth (degrees per meter)
\(\kappa\) \(2.8 \times 10^{3}\) | Thermal diffusivity (square meters per year)
\(r_{earth}\) \(6.378 \times 10^6\) Radius of the earth (meters).

So, now we can do the calculations, ... but the answers come out a bit odd.

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
from math import log, pi

# Here are our parameter estimates ...
grad = 2.5e-2
T0 = 1.2e3
re = 6.4e6
k = 2.8e3

# 1-dimension calculation
L = 2*re
log(4*T0/(abs(grad)*L))*L**2/(pi**2*k)

# Spherical calculation
log(T0/(abs(grad)*re))*re**2/(pi**2*k)

And if we think hard about these formulas, we discover that they are a little fishier than we would like.

Kelvin's Semi-infinite approximation

If, instead, we approximate the earth as a semi-infinite slab that starts off at constant temperature and then cools from the surface, we get an easier problem (and the one Kelvin acctually used -- based on Green functions rather than Fourier series). The solution

\[\theta(x,t) = \theta_0 \operatorname{erf}\left( \frac{x}{2 \sqrt{\kappa t}}\right)\]

which uses a special function called the error function because of it's role in probability theory.

\[\operatorname{erf}(x) := \frac{2}{\sqrt{\pi}} \int _ {0}^{x} e^{-u^2} du\]

After evaluating, we find the magic formula that William Thomas found ...

\[\frac{d\theta(0,t)}{dx} = \frac{\theta _ 0}{\sqrt{\pi \kappa t}}\]

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sympy import *
from sympy.abc import *
theta_0 = Symbol('theta_0')

theta = lambda x, t : theta_0*erf(x/2/sqrt(t*kappa))

# check that this is really a solution of the heat equation
assert 0 == theta(x,t).diff(t) - kappa * theta(x,t).diff(x,x)
# check the boundary condition
assert 0 == theta(0,t)

# So, if we calculate t based on our parameters above, ...
slopeFormula = theta(x,t).diff(x).subs(x,0)
tFormula = Eq(2.5e-2, slopeFormula).subs([(theta_0, 1.2e3),(kappa, 2.8e3)])
print solve(tFormula, t)

With these estimates, we get a young-earth prediction of only \(200,000\) years old.