Early natural philosophy at the beginning of the European renaissance derived from Aristotilian and biblical sources.
James Ussher's biblical chronology (15811656) 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.
Joseph Fourier (17681830) published an important work on the theory of heat in 1807. This became the basis of our modern understanding.
Hypothesize that the rate of heatflow 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\]
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.
Let's consider a toyproblem 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 uninsulated 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 roomtemperature 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 )}\]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

The leading, slowestdecaying 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.
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 SturmLouville 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 zeroderivative 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}\]
1 2 3 4 5 6 7 8 9 10 11 12 13 

Again, dropping all but the slowestdecaying 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)\]
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). 
The sun's outer temperature is about 5800 degrees Kelvin. All silicates melt by 1500 degrees Kelvin (1200 Celcius).
Silicon's thermal diffusivity is about \(0.8{\rm cm}^2/{\rm s}=2.523\times 10^3{\rm m}^2/{\rm yr}\). although it is a function of temperature. Earth radius is \(6.378 \times 10^6\) meters.
The geothermal gradient in stable parts of the crust is an increase of 25 to 30 degrees celcius per kilometer of depth.
1 degree Farenheit per 50 feet is about \(3.64 \times 10^{2}\) degrees Celcius per meter.
So, now we can do the calculations, ... but the answers come out a bit odd.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 

And if we think hard about these formulas, we discover that they are a little fishier than we would like.
If, instead, we approximate the earth as a semiinfinite 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}}\]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

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