Further reading

- William Thomson’s 1863 paper on the age of the earth.
- John Lienhard’s podcast on the controversy.
- A historical overview
- An excellently detailed historical account by Hall and Wilkie (2015) of Thomson’s calculation and its context, exerpted from the book
*Kelvin, Thermodynamics, and the Natural World*, edited by Collins et al., 2016. - The magazine
*Scientific American*played an important role at the time, with an article in 1877, a note from Kelvin in 1895, and a special historical review issue of the history of attempts to determine how hold the earth is.

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.

William Thomson (1824-1907), eventually to become known as Lord Kelvin, believed that everything should have recourse to mechanical 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 to consider – the effects of the moon on the earth, and 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.

Hot and cold are not substances, but opposite extremes of the measure of heat. The only theoretical difference between heating and cooling is a sign change.

Unless something is introducing or removing it, heat is conserved.

Heat is not measured directly, but in terms of temperature. Temperature is a representation of heat’s ability to induce change. In modern thermodynamics, temperature is the rate of change in energy per change in entropy. Mathematically, it is more natural to think in terms of inverse temperature, which is the change in entropy per change in energy: Adding small amounts of energy to cold things leads to a large increase in entropy and disorder.

Adding a small amount of energy to a hot thing causes a small increase in entropy because it is already strongly disordered.Fourier’s law states that the rate of heat flow is proportional to the gradient of the temperature, 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.

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\]

where

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

The earth is a very large sphere, with a radius \(r = 6.37 \times 10^6\) meters.

Earth was initially molten rock of uniform temperature

The earth has remained homogeneous (meaning that it’s substance is the same everywhere) and isotropic (meaning that heat moves at the same speed in all directions) for it’s whole lifetime.

The temperature at the earth’s surface is fixed at approximately 0 degrees Celsius

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.

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

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}(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\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)
```

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.

Daily variation in temperature – but these are small, and average out over the year. Similarly, annual fluctuations in temperate are rather small compared to high temperatures we expect to encounter at the earth’s core, and are also self-cancelling.

We can now try to assess the validity of our flat-earth assumption. We find that only the outer shell of the earth shows significant temperature change under the flat earth assumption, which suggests that the solution would be very similar for a spherical model.

uncertainty in thermal diffusivity – different rock types, soil vs rock We now believe there are rapid changes in temperature at the boundary layers between crust and mantle, and mantle and core.

uncertain initial temperature: Thomson’s initial temperature estimate may be on the high end, but if the earth started off cooler, then the earth would be even younger than estimated.

There might be additional sources of heat. George Darwin, a successful student of Thomson’s, thought that the heat in the earth’s mantle might be due to tidal effects from the moon.

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.

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

- Schrodenger equation in quantum mechanics
- Navier-Stokes equation in fluids
- Black–Sholes equation in finance
- Fisher-KPP equation in genetics

William Thomson’s flat-earth heat equation solution \(\theta(x,t) = \theta_0 \operatorname{erf}\left( \frac{x}{2 \sqrt{\kappa t}}\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=200,000\) from \(x=0\) to \(x=7 \times 10^6\) when \(\kappa = 2.8 \times 10^{3}\).

Is the flat-earth approximation reliable?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)

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?

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

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\).

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?

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

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

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.

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.

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

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?

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.

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\).

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

Minimum size of a marine reserve.

Optimal fishery management always contains a reserve.