Heat and the Age of the earth

( previous, home, next )

Further reading

History of the history of the earth

In the early 1800’s, scientists across Europe were embroiled in debate over a great scientific question – how old is the earth and how has it been shaped.

At the beginning of the European renaissance, natural philosophy derived from Aristotle’s writings and biblical sources. The most common practice for determining the age of the earth was to count generations backward through the bible to Adam and Eve. In 1650, James Ussher’s history of the world was published. Ussher’s history was based on biblical chronology, and placed the beginning of the world at the 22rd of October, 4004 BC, only 5,654 years before the publication. Ussher’s contribution was one in an ongoing debate among Christian scholars over the age of the earth. However, as scientists began to closely observe the geological record found in layers of rock around the world, they began to refine their stories of world history to match their observations. Sea shells were found on the tops of mountains. Thick rock layers alternated back and forth between aquatic and terrestrial deposition. Canyons carved deep through layers of rock by rivers.

The first philosophies attempting to explain these irregularities used recourse to catastrophic events, hence the name Catastrophists. In particular, the biblical flood described by Noah was often invoked, in addition to well-remembered earthquakes and volcanic eruptions. However, the diversity and pattern of the events observed in layers of rock around Europe required many diverse catastrophes to reconciliation theory and observation. In the process, their complexities became unpalatable to some.

Uniformitarians proposed that the world was and had always been in a state of slow and steady change, with the uplifting and subsequent erosion of mountains. From the uniformitarian perspective, the earth was essentially ageless. Uniformatarianism was started by James Hutton (1726-1797), based on observations of erosion, refined by John Playfair (1748 - 1819) and popularized by Charles Lyell (1797-1875). Lyell’s student, Chuck Darwin (1809-1882) engaged in the debate: geological, uniformitarian estimates would eventually be consistent with the evolutionary time scales that he would propose. When it came time to fight Kelvin, Thomas Henry Huxley (1825 - 1895) took up Darwin’s cause.

Enter William Thomson

William Thomson (1824-1907), eventually to be come Lord Kelvin, believed that everything should have recourse to physical principles in science. He found the lack of rigor in uniformatarianism… unsatisfactory. Thomson thought there must be a physical way to calculate the age of the earth. He found two mechanisms worth considering.

  1. The effects of the moon on the earth.
  2. The temperature of the earth.

You can read more about the first, if you are interested. We will focus on the second.

Thomson was not the first person to attempt to use the cooling of the earth to estimate its age. One of the first attempts to provide a mechanistic alternative to biblical chronology came from the religiously devote Isaac Newton (1642-1726). Newton had developed a theory of heat based on the idea that the rate of cooling of a hot object was proportional to it’s temperature difference from its background. He proposed that the earth was once molten rock and had cooled to its spheroid shape. Using his law of cooling and crude observations of a one-inch sphere of iron, he estimated that the earth was 50,000 years old, about 10 times older than Ussher had proposed. Georges-Louis Leclerc, Comte de Buffon (1707-1788) took up his idea and ran experiments on 10 spheres, from 0.5" to 5" in 0.5" increments, finding the cooling time from white hot to handle-able. Extrapolating from data with some corrections for the sun’s heating, Buffon estimated the earth’s age as 74,832 years old.

In 1862, building on an idea of Buffon, Thomson published his ideas on applying the heat equation of Fourier (which was NOT yet a well known tool) to the problem of the earth’s age. 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 \theta}{\partial t} = \kappa \frac{\partial^2 \theta}{\partial x^2}\]

In three dimensions, the curvature in each direction can contribute to the rate of change, so

\[\frac{\partial \theta}{\partial t} = \kappa \nabla^2 \theta\]


\[\nabla^2 \theta = \frac{\partial^2 \theta}{\partial x^2} + \frac{\partial^2 \theta}{\partial y^2} + \frac{\partial^2 \theta}{\partial z^2}.\]

Thomson’s model

In spherical coordinates for a spherically symmetric solution, the heat equation with constant conduction rate \(\kappa\) takes the variable-coefficient form

\[\frac{\partial \theta}{\partial t} = \frac{\kappa}{r} \frac{\partial \theta }{\partial r} + \kappa \frac{\partial^2 \theta}{\partial r^2}\]

Note that if the radius of the earth is large and the conduction rate is small, then over short times, we might be able to make a flat-earth approximation,

\[\frac{\partial \theta}{\partial t} = \kappa \frac{\partial^2 \theta}{\partial r^2}.\]

We must check the flat-earth approximation’s validity at the end of our calculation. If the temperature at the earth’s core has significantly decreased over the time in question, then we can improve our calculation accuracy significantly by returning to the spherical form of our heat equation.

Green function solutions

For some details on Green functions and Delta functions, see this appendix

Thomson’s solution

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

Green function \[G(x,t) = \frac{1}{\sqrt{4 \pi \kappa t }} e^{-x^2/4\kappa t}\] solves the one-dimensional heat equation.

\[\Psi(x,t) = \theta _ 0 \int _ 0 ^{\infty} G(x-a,t) - G(x+a,t) da\]

The function \(\Psi(x,t)\) satisfies the boundary condition of being zero at the earths surface because when \(x=0\), \(G(-a,t) - G(a,t) = 0\) for shifts \(a\). The function \(\Psi(x,t)\) also satisfies the initial condition \(\Psi(x,0) = \theta _ 0\) for \(x \in (0,\infty)\). Thus, \(\Psi(x,t)\) is the solution of Thomson’s model, in integral form. Applying some calculus, this solution can be simplified to show that \[\theta(x,t) = \Psi(x,t) = \theta _ 0 \operatorname{erf}\left( \frac{x}{2 \sqrt{\kappa t}}\right),\] where \[\operatorname{erf}(x) := \frac{2}{\sqrt{\pi}} \int _ {0}^{x} e^{-u^2} du.\] The special function \(\operatorname{erf}\) is called the “error function” because of it’s role in probability theory.
After using Liebnitz’s formula to differentiate, \[\frac{d\theta(0,t)}{dx} = \frac{\theta _ 0}{\sqrt{\pi \kappa t}}.\] Solving for time, we find the seemingly magic formula that William Thomson found … \[t = \frac{\theta _ 0^2}{\pi \kappa \theta'(x=0)^2}.\]

[Show code]
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)

Data and calculation

Symbol Value Description
\(\theta _ 0\) \(1,500\) Uniform initial temp. in Celsius above surface temp.
\(\theta'(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).

Thomson argued that the initial temperature of the earth should have been the same as that of the sun, and did experiments to estimate the conduction of heat by stone. With these estimates, William Thomson get a young-earth prediction of only \(200,000\) years old. This split the difference between bible-based calculations of thousands of years and geological estimates of many millions of years, but is well short of today’s estimate of 4 billion years old.

Estimating the importance of assumptions

If Thomson had really been thinking well, he might have realized the disagreement between his estimate and geological estimates pointed to a mysterious, yet unknown phenomena that was keeping the center of the earth warm. But at the time, the general feeling was that physics was mostly a solved subject, and that there was little new left to discover. (They had no clue about atomic theory, general relativity, and quantum phenomena yet). Thomson’s model did not account for radioactivity, which was only discovered towards the end of Thomson’s life. The natural fission of heavy atoms like uranium generates heat at the earth’s core. When Lord Rayleigh adjusted Thomson’s calculations to include radioactivity, he estimated the age of the earth about 2 Billion years.

However, even Lord Rayleigh was wrong. His theory requires way more radioactivity than it turns out we think there actually is within the earth! The real source of error in Thomson’s model was his assumption of homogeneity and isotropy – this fails to allow for convection processes. The earth’s mantel is solid, but moves slowly – on the order of centimeters a year. It takes 10’s of millions of years for a convection cycle to complete. But even this slow movement is enough to change the earth’s interior heat flow. Mantel convection effectively shrinks the crust and allows a steep temperature gradient to persist over much longer time scales than would be predicted otherwise.

Further applications to diffusion

By the early 1900’s, scholars had begun applying Fourier’s heat equation to the study of other problems involving diffusion and random walks.


  1. William Thomson’s flat-earth heat equation solution [(x,t) = _0 ( ),] can only be a good approximation for the earth if the changes in temperature are confined to the upper surface of the earth. Using scipy.special.erf, plot this solution at time \(t=200,000\) from \(x=0\) to \(x=7 \times 10^6\) when \(\kappa = 2.8 \times 10^{3}\).
    Is the flat-earth approximation reliable?

  2. A legend has been passed down through the years that (as a good Frenchman), Jean-Baptiste Fourier was motivated to carry out his research into the natural philosophy of heat primarily because of a love of food and wine. Specifically, he was having trouble keeping his bottles of wine at the right temperature. In the summer, the bottles would get too hot and go sour. In the winter, the bottles would get too cold, freeze, and crack. The clear solution was to store the bottles in a cellar where the temperature changes would be less, but how deep should that cellar be? The deeper the cellar, the less temperature flucturation we’d expect over the year, but the deeper you go, the more expensive it is to build and the more trouble to get the wine up and down. And if possible, you’d like a cellar that gave you cooler wine in the summer and warmer wine in the winter. The thermal diffusivity of ground around the wine-cellar is about 0.01 centimeters square per second. Solve the 1-dimensional heat equation using separation of variables to answer the questions below. (Hint. Use the complex exponential \(e^{i\omega t}\) and take the real part at the end to simplify your algebra.) (Tung, 2007, Lin & Segel, 1974)

    1. If you consider daily temperature variations of plus or minus 8 degrees, how deep does the cellar have to be for the daily variation in temperature to be less than 2 degrees?

    2. If the annual mean daily temperatures oscillate between 0 and 30 degrees Celsius, determine the optimal depth for Fourier’s wine cellar.

schematic of cross-scection of Fourier's wine cellar

  1. As nuclear physics predicted and history has shown, standard nuclear power plants have an inherent vulnerability. Under stable operation, the heat a reactor generates is dissipated by it’s enveloping cooling system as quickly as it is generated. However, as the reactor temperature increases, the fission rate increases. There is a critical tipping point at which heat is generated faster than the cooling system can compensate, fission becomes a run-away chain reaction, and we get a meltdown. A simple model of this phenomena can be made by adding a nonlinear growth term to the equation for temperature. Let \(\theta(x,t)\) represent the reactor’s temperature over a one-dimensional cross-section from \(-L/2\) to \(L/2\). At the boundaries, the temperature \(B\) is set by the cooling system. Within the reactor, the temperature \(\theta\) obeys the non-dimensionalized partial differential equation \[\frac{d\theta}{dt} = \frac{d^2\theta}{dx^2} + e^{\theta}.\]

    1. If the react is operating at steady-state and the temperature distribution is not changing, what is the differential equation the temperature distribution will solve?

    2. Where, between \(-L/2\) and \(L/2\), does the reactor maintain its maximum temperature?

    3. Rewrite your steady-state equation as a first-order system.

    4. Numerically solve your first-order steady-state system for \(x \in [0,2.5]\) for maximum temperatures \(m \in \{-1,0,1,2,3,4\}\), and plot your solutions.

    5. If the cooling system keeps the boundary temperature \(B=1\), what is the largest reactor size \(L^*\) for which a steady-state temperature profile exists.

    6. If \(L = 3\), what is the cooling system’s maximum allowed temperature \(B^*(L)\) for which a steady-state temperature distribution can exist?

    7. A simplification of this reactor model is to remove space entirely and to treat heat loss as a linear effect, so \[\frac{d\theta}{dt} = -a \theta + e^{\theta}.\] Use the basic graphical methods you learned for ordinary differential equations to characterize to possible reactor behaviors when \(a=5\) and when \(a = 1\).

    8. In light of your results for our spaceless model, how do you interpret the fact that for \(B < B^*(L)\), there are two different steady-state solutions satisfying the boundary conditions?

  2. The linear Fisher-KPP equation [ = D + r n ] describes the invasion of an exponentially growing population into new territory. The speed of these invasions can be approximately determined by finding the slowest travelling wave solution.

    1. Assume \(n(x,t) = exp(w (x-ct) )\), where \(w\) is the wave shape and \(c\) is the speed. Substitute and determine a relationship between \(c\) and \(w\) in terms of \(D\) and \(r\).

    2. Find the value of \(w\) that minimizes \(c\), and determine the corresponding minimum speed.

  3. Minimum size of a marine reserve.

  4. Optimal fishery management always contains a reserve.

( previous, home, next )