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 another debate over a great scientific question -- how old is the earth? At the beginning of the European renaissance in the 1400's, 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, some 5,654 years before the publication. Ussher's contribution was just one among many in an ongoing debate among biblical scholars.

Around the same time, naturalist scholars had began to write about their observations of the geological record found in layers of rock around the world. Inspection of geological formations revealed seashells on mountain tops, thick rock layers that alternated back and forth between aquatic and terrestrial deposition, and canyons carved deep through layers of rock by rivers. At first, the observations were easily reconciled with biblical chronologies. Seashell fossils on mountain tops were evidence of Noah's flood. It was natural then to attempt to reconcil other parts of the geological chronology with other catastrophes like earthquares and volcanoes.

However, the diversity and pattern of the events observed in layers of rock around Europe required many different catastrophes to reconciliation chronological theory and geological observation. Their relationships became so complicated that the whole methodology became unpalatable to some. These "uniformitarians", including James Hutton (1726-1797), John Playfair (1748 - 1819), and Charles Lyell (1797-1875), proposed that the world was and had always been in a state of slow and steady change, continually trading off the uplifting of new mountains against their erosion. From the uniformitarian perspective, the earth was essentially ageless. Lyell's student, Charles Darwin (1809-1882), used uniformitarian theory as a foundation for his own theory of the evolution of life through natural selection.

William Thomson (1824-1907) believed that everything should have recourse to mechanical principles in science. He found the lack of rigor in uniformatarianism... unsatisfactory. The recently developed laws of thermodynamics implied that the earth must age and that this aging must be observable. Thomson thought there must be a physical way to calculate the age of the earth. Thomson considered several ideas and mechanisms including the tidal effects of the moon on the earth, and the temperature of the earth itself. We will focus on the last.

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), widely remembered for his needle problem, took up Newton's 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 published his results in 1778 estimating the earth's age as 74,832 years old.

Fourier's heat equation

In 1807, Joseph Fourier (1768-1830) published his book "The analytic theory of heat". The precocious Thomson read Fourier's book while doing undergraduate studies in 1840 at the age of 15. He soon recognized it's relevance to Buffon's method for estimating the earth's age.

It will be useful to review a little of the modern theory of heat at this point. The history of the theory of heat is interesting and somewhat complicated history, but here are the basics. Heat is a measure of the internal energy of an object. Since heat is a form of energy, it obeys the law of conservation of energy -- in the absence of mechanisms introducing or removing it, total heat must be conserved.

Heat is not measured directly, but in terms of temperature. Temperature is what we sense when we say something is hot or cold, and is a representation of heat's ability to induce change and do work. In modern thermodynamics terminology, temperature is the rate of change in energy per change in entropy. Mathematically, it is more natural to think in terms of reciprocal temperature, which is the change in entropy per change in energy: adding small amounts of energy to cold things at low temperature leads to a large increase in entropy and disorder, and adding a small amount of heat to a hot thing at high temperature causes a small increase in entropy because it is already strongly disordered. The rate of heat added per unit mass divided by the rate of temperature increase is what we call the "specific heat". In practice, we often confound heat and temperature. We can get away with this because for many systems not undergoing a phase-transition like melting, the specific heat is roughly constant, implying a linear relationship between local temperature and local heat.

Now, consider a large metal bar which is heated at some spots and cooled at others. Fourier's law of heat states that the rate of heat flow at each location is proportional to the gradient of the temperature at that location, but the direction of flow is opposite to the direction of the temperature gradient. So, cold flows towards hot and hot flows toward cold. Fourier's law predicts that temperature smooths out over time -- hot spots become cooler, while cold spots warm up. If the specific heat is constant, then, the rate of heat flow is proportional to the gradient of the heat itself, \begin{gather*} J \propto \nabla H \end{gather*}

where \(J\) is the flux and \(H\) is the local heat. It then follows from vector calculus that the rate of change of the local heat will be proportional to the divergence of the flux, or \[\dot{H} \propto \nabla \cdot \nabla H.\] Exploiting the linear relationship between heat \(H\) and temperature \(\theta\) (\(H \propto \theta\)), we arrive at Fourier's heat equation \[\dot{\theta} = \kappa \nabla^2 \theta,\] where \(\kappa\) is the thermal diffusivity of the metal bar (or whatever substance). In three dimensional Cartesian spatial coordinates, this is equivalent to \[\frac{\partial \theta}{\partial t} = \kappa \left( \frac{\partial^2 \theta}{\partial x^2} + \frac{\partial^2 \theta}{\partial y^2} + \frac{\partial^2 \theta}{\partial z^2}\right),\] while in one spatial dimension, \[\frac{\partial \theta}{\partial t} = \kappa \frac{\partial^2 \theta}{\partial x^2}.\]

Thomson's model

In 1862, Thomson finally published his ideas on applying the Fourier's heat equation of Fourier to the problem of the earth's age. Thomson set up his model as follows. Imagine that the earth is was a molten homogeneous sphere of uniform temperature similar to the sun at the birth of our solar system. Since that time, it has cooled from the outside in, with the surface temperature quickly approaching what we see today at the earth's surface (approximately 0 degrees Celsius), while the inside has cooled more slowly, and remains molten today.

Assuming the temperature profile is spherically symmetric and we use spherical coordinates with \(r\) being the radius from the center of the earth, the heat equation with constant conduction rate \(\kappa\) takes the variable-coefficient form \begin{gather} \frac{\partial \theta}{\partial t} = \kappa \left( \frac{2}{r} \frac{\partial \theta }{\partial r} + \frac{\partial^2 \theta}{\partial r^2} \right) \end{gather} with boundary conditions \begin{gather} \theta(r _ {\oplus},t) = 0, \\ \theta(0 < r < r _ {\oplus},0) = \theta_0, \\ \theta'(0,t) = 0, \end{gather}

where \(r _ {\oplus}\) is the earth's radius, and \(\theta _ 0\) is the earth's initial temperature.

Analyzing the dimensions of our symbols, we find there are three dimensionless groups, \begin{gather*} \left( \frac{\theta}{\theta_0}, \frac{r}{r _ {\oplus}}, \frac{t \kappa}{r _ {\oplus}^2}\right), \end{gather*} and so, all solutions can be expressed as \(\theta(r,t) = \theta_0 \psi(r/r _ {\oplus}, \kappa t / r _ {\oplus}^2)\), where \begin{gather} \frac{\partial \psi}{\partial t} = \frac{2}{r} \frac{\partial \psi }{\partial r} + \frac{\partial^2 \psi}{\partial r^2} \end{gather} with boundary conditions \begin{gather} \psi(1,t) = 0, \quad \psi(0<r<1,0) = 1. \quad \psi'(0,t) = 0. \end{gather}

Flat earth again

We will make one additional modelling approximation -- a flat-earth approximation. If the earth is young and its cooling has been modest over it's lifetime and not penetrated very deeply, only the surface temperature should be effected, and the curvature of the earth should have only a small effect on the cooling process over this spatial scale.

To make this idea precise, we introduce a new positive dimensionless scaling parameter \(\epsilon\), and new independent coordinates \(x := (1 - r)/\epsilon\) and \(s = t/\epsilon^2\). Then under this change of coordinates, \begin{gather} \frac{\partial \psi}{\partial s} = \frac{2 \epsilon}{\epsilon x-1} \frac{\partial \psi }{\partial x} + \frac{\partial^2 \psi}{\partial x^2}. \end{gather} If only the surface of the earth has cooled, then we should be able to perform our calculations over just a thin slice of the earth, in which case we can assume \(\epsilon \ll 1\). When \(\epsilon\) is small, it follows that \begin{gather} \frac{\partial \psi}{\partial s} \approx \frac{\partial^2 \psi}{\partial x^2} , \quad \psi(x=0,t) = 0, \quad \psi(x>0,0) = 1. . \end{gather}

We've fudged things a little here with our initial condition. Strictly speaking, we should only consider a finite domain where \(x \leq 1/\epsilon\), rather than a semi-infinite domain. But if our hypothesis of a young earth is correct and the flat-earth approximation is accurate, the temperature of the core of the earth should not have departed measurably from it's initial temperature of 1, and our fudged semi-infinite domain should give the same answer as the exact finite-domain version. But this all depends on the actual age of the earth, which we have not yet calculated. So the validity of our flat-earth approximation is something we'll have to revisit after we've made our estimates of age.

Let's take a moment to consider what we've done.

Method of images

To solve his flat-earth model, Thomson used a clever technique called the "method of images". The function \[G(x-a,s) = \frac{1}{\sqrt{4 \pi s }} e^{-x^2/4s}\] solves the one-dimensional heat equation on an infinite domain with initial condition \[G(x-a,0) = \delta(x-a).\] Such a solution is called a Green function, in honor of George Green -- an ordinary miller whose self-taught ideas revolutionized the fields of applied mathematics to which he applied himself. The Green function can reduce the problem of solving the one dimensional heat equation initial-value problem on an infinite domain to quadratures. For any initial condition \(\psi _ 0(x)\) on an infinite domain, the solution of the heat equation is \[ \psi(x,s) = \int _ {-\infty}^{\infty} G(x-a,s) \psi_0(a) da.\] We can confirm this by differentiation and substitution. For some details on Green functions and Delta functions, see this appendix

However, Thomson was faced with a slightly different problem -- he wanted to solve the heat equation on a semi-infinite domain with boundary condition \(\psi(0,t) = 0\). This is where the method of images enters our conversation. Suppose that for each Green function to the right of zero, we place another one at the same distance to the left of zero with equal magnitude but opposite sign. The linearity of the heat equation implies that these two will cancel each other out exactly at the midpoint between the -- the origin. Suppose, then, we consider the function \[\Psi(x,s) = \int _ 0 ^{\infty} G(x-a,s) - G(x+a,s) da.\] This function satisfies the boundary condition of being zero at the earths surface because when \(x=0\), then \(\Psi(x,s) = G(-a,s) - G(a,s) = 0\) for all shifts \(a\). The function \(\Psi(x,s)\) also satisfies the initial condition \(\Psi(x,0) = 1\) for \(x \in (0,\infty)\). Thus, \(\Psi(x,s)\) is the solution of Thomson's model. Applying some calculus, this solution can be simplified to show that \[\Psi(x,s) = \operatorname{erf}\left( \frac{x}{2 \sqrt{s}}\right),\] where \[\operatorname{erf}(z) := \frac{2}{\sqrt{\pi}} \int _ {0}^{z} 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\psi(0,t)}{dx} = \frac{1}{\sqrt{\pi t}}.\] Solving for time, Thomson found the surprisingly simple formula \[s = \frac{1}{\pi \psi'(0)^2},\] or in re-dimensionalized form, \[t = \frac{\theta_0^2}{\pi \kappa \theta'(0)^2}.\]

Thanks mining work of the time, it was well known that the earth's initial temperature increases \(\theta'(0) \approx 1/28\)th of a degree per meter. But Thomson acknowledged there was uncertainty as to which values to use for the initial temperature \(\theta_0\) and the thermal conductance \(\kappa\). Still, reasonable values like an initial temperature of the earth \(\theta_0 \approx 4,000^{\circ}\) Celsius, and a thermal diffusivity \(\kappa \approx 2.8 \times 10^{3}\) implied an age of \(1.4\) million years. Thomson himself concluded in 1869 that the earth was not more than 1 billion years old, and most likely 100 million years old. While far far larger than biblical chronologies, Thomson's estimate was far short of the modern estimate of 4 billion years.

Estimating the importance of assumptions

Our flat-earth approximation relied on the earth's lose of heat being confined to a thin boundary-layer at the earth's surface. Based on an estimate of the earth's radius \(r = 6.34 \times 10^6\) meters and our other parameter values, Thomson's model predicts the cooling will not yet have affected the earth any deeper than 4% of the earth's radius, implying \(\epsilon = 0.04\) in our flat-earth approximation. This seems small enough that the flat-earth approximation we made was reasonable.

The disagreement between Thomson's 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 -- nobody had a 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.


  1. William Thomson's flat-earth heat equation solution \(\psi(x,s) = \theta_0 \operatorname{erf}\left( \frac{x}{2 \sqrt{\kappa s}}\right),\) 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=1.4 \times 10^{6}\) from \(x=0\) to \(x=6 \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} + r e^{\theta},\] based on the Frank-Kamenetskii approximation of the Arrhenius law. For the following, assume \(r=1\).

    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. In light of your results for our spaceless fire model earlier in an earlier chapter, 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 \[\frac{dn}{dt} = D \frac{d^2n}{dx^2} + 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 )