Chapter 1
The heat of the matter

The Earth is a heat engine, and as with all such devices it is doomed eventually to die! Indeed, it is this very heat loss that is responsible for lithospheric growth and ultimately provides the energy source driving the motions of the lithosphere. Heat loss from the Earth's interior is accompanied both by conduction and advection, which can be simply codified in terms of the requirement for the conservation of energy - and which therefore can serve to illustrate one of the great conservation principles of physics.

1  Conduction and advection

Thermal energy may be transported in two fundamentally different modes; namely by conduction and by advection. In the former heat diffuses through matter while in the latter it is carried by the physical transport of matter.


Everyday experience witnesses the fact that thermal energy is transferred from hot bodies to cool bodies, with the rate of thermal energy transfer, or heat flow, dependent on the nature of the material. Fourier postulated that the heat flow, q, is proportional to the negative of the temperature gradient, because heat flows down temperature gradients, with the constant of proportionality given by a material dependent parameter known as the conductivity, k, a postulate consistent with most subsequent experimental investigation:

q = -k grad  T
Written in this form grad T signifies the gradient of T (often written in the equivalent form T) and as such is a vector operator. That is, grad defines the magnitude and direction of the gradient of a scalar quantity such as temperature.


The physical motion of material results in advective transport of heat. Using a Eulerian description it is easy to see that the rate of heat advected through a fixed point is dependent on the velocity of material transport and the temperature gradient. Thus, as either the tempertaure gradient or the velocity tends to zero, the rate of heat advection must also tend to zero.

2  The thermal energy balance

In order to understand the thermal evolution of material subject to conductive heat transfer it is necessary to derive the thermal energy balance that expresses the fundamental physical constraint provided by the requirement of conservation of energy. Ignoring the advection of heat through the movement of the medium then the thermal energy balance in a control volume can be expressed simply as:

rate of gain of thermal energy
rate of thermal energy flowing into the control volume
rate of thermal energy flowing out of control volume
rate of heat production in control volume
Since the thermal energy of a volume of unit dimensions is given by:
Cp  r  T
where Cp is the heat capacity, the first term in Eqn 2, the rate of gain of thermal energy, can be expressed simply as
( Cp r T )
which for temperature independent r and Cp is equivalent to
Cp  r  T
The difference between the second and third terms in Eqn 2 gives the negative of the divergence of the heat flow across the boundary of the control volume:
- div q
Note that the negative sign arises because divergence actually measures the rate of loss whereas in Eqn 2 we are really interested in the rate of gain (or convergence). Assuming that k is independent of T then the difference between the second and third terms in Eqn 2 can be expressed in terms of T in the following way by substituting the expression for q given in Eqn 1:
k div (grad  T )
= k  ·(T )
= k

+ 2T
+ 2T

The final term in Eqn 2 is given by:
r H
where H is the heat production per unit mass.

Thus the thermal energy balance given by Eqn 2 can be given in the following ways:

= k(div  ( grad T)) + H
= k(·T) + H
= k

+ 2T
+ 2T

+ H
where k is the thermal diffusivity:
k = k
r Cp
The various forms of the thermal energy balance in Eqn 3 form the diffusion equation or, alternatively, the heat equation.

The diffusion-advection equation appropriate to whole rock advection, is obtained by adding an additional term to Eqn 3 which expresses the rate at which heat is advected through a point. This advective contribution is given by the product of the velocity and the temperature gradient, giving:

= k2T - v . T + H
where the second term on the right is the advective term.

For many purposes, it is useful to recast Eqn 4 in a dimensionless form:

= 2T- PeT  T
where the dimensionless variable T and t are given by
T = T/T0
z = z/l
t = t k/ l2
where l is the appropriate length-scale (for example, the thickness of the lithosphere), and T0 is the characteristic temperature difference at the appropriate length-scale. In Eqn 5 PeT is the thermal Peclet number which is a dimensionless variable
PeT = v  l
. The Peclet number expresses the ratio of the advective to diffusive terms. For PeT > 10 the advective term dominates the thermal evolution and diffusion can be largely ignored, while for PeT < 1 diffusion dominates and the advective term can be largely ignored. For 1 < PeT < 10 both diffusion and advection contribute to the thermal evolution.

3  Thermal time constants

It is useful to define the characteristic diffusive response time to thermal perturbations when heat is transferred by conduction. For example consider a bunsen burner applied to the base of a conducting plate (geologically, this is analagous to the impingement of a mantle plume on the base of the lithosphere). We are interested in how long does it takes for the top surface to feel the effect of the perturbed lower boundary condition. Clearly the response time depends on the length-scale (i.e. the thickness of the plate) and (inversely) on the thermal conductivity (or, equivalently, diffusivity). In this problem we can ignore the effects if internal heat production (which does not change the response time to a change in boundary conditions). From Eqn 3 the term on the right, [(T)/( t)], gives the order the temperature divided by the timescale for thermal response, t. Similarly the frist term on the right of Eqn 3, k( div  ( grad T)), is of order kT/ l2, where l is the appropriate length scale. Equating the two quantities gives :

t l2
The constant of proportionality normally used is 1/p2 giving:
t = l2
p2  k
For the lithosphere typical parameter values are l = 125 km and k = 10-6 m 2 s-1 giving t = 50 Ma.

An important feature of the thermal time constant is its dependence on l2. For the lithosphere this implies that short wavelength fluctuations are damped on timescales which are very short compared to lithospheric-scale perturbations (and, incidently, explains why a 3 m deep cellar is insensitive to seasonal atmospheric temperature fluctuations).

4  Continental geotherms

In the steady-state, that is when [(T)/( t)] = 0, and when heat flow is only in one direction, z, Eqn 3 reduces to the ordinary differential equation:

0 = k  d2 T
d  z2
+ r H
which can be readily integrated with appropriate boundary conditions to yield a geotherm. For example, assuming a constant distribution of heat producing elements with depth, Eqn 6 can be integrated using the boundary conditions appropriate to the continental lithosphere: q = - qs at z = 0 (i.e., the heat flow at the surface of the lithosphere in the direction of increasing depth is the negative of the surface heat flow), and T = Ts at z = 0 (i.e., the temperature at the top of the lithosphere is the surface temperature). The resulting integration yields the expression
Tz = Ts + qs z
- r H z2
2 k
relating the temperature at any depth, z, to the conductivity, surface temperature, Ts, surface heat flow, qs, and heat production, H. The base of the lithosphere is defined by the temperature sensitive rheological properties of peridotite and can therefore be treated as an isotherm, Tl = 1280o C. The heat flow at the surface of the lithosphere, qs, therefore reflects the thickness of the lithosphere as well as the distribution of heat sources within the lithosphere.

5  Natural convection

Natural convection is the motion induced by thermally created density gradients in viscous materials, otherwise known as fluids. In the presence of gravitational field a compositionally homogeneous fluid heated from below will develop an unstable density stratification. Once convective motion is initiated and hotter, less dense material are juxtaposed against cooler more dense material by slight upward displacement (and vice versa in a downwelling zone) the motion will be maintained by buoyancy forces operating on the individual fluid packets (Figure 3.1). Counteracting the tendency to rise will be the viscous drag of the material and the tendency for the rising material to loose heat to its surroundings by thermal diffusion (Figure 3.2 & 3.3). Whether convective motion is initiated depends the material properties (thermal diffusivity, coefficient of thermal expansion, viscosity). The likelihood of convective motion will be enhanced by the fluid having a high coefficient of thermal expansion, subject to a large temperature gradient over a large vertical extent and diminished by a high thermal diffusivity or by high viscosity. Initiation of the convective motion requires a finite amplitude perturbation of the thermal, and hence density, structure in order to localize upwelling and or downwelling.

Figure 1: For the initiation of convection an imbalance of forces is required. A fluid layer subject to an inverted density gradient due, for example, to basal heating and/or surface cooling is capable of undergoing convection. Slight upward displacement of part of the basal layer will produce a horizontal density gradient which will lead to a net upward directed buoyancy force on the displaced fluid parcel, further enhancing the prospect for upward motion and eventually fully developed convective motion. Initiation of the convective motion requires a finite amplitude perturbation of the thermal and hence density structure in order to localize upwelling and or downwelling.

Figure 2: The displacement of a parcel of hot light fluid is resisted by viscous dissipation or drag in the surrounding fluid (strictly the rate at which momentum diffuses into the surrounding medium), a measure of which is given by the kinematic viscosity of the fluid.

The Rayleigh number, RaT, is a measure of the relative magnitude of the destabilizing buoyancy forces (Figure 7.1) and the stabilizing effects of momentum and heat diffusion (Figure 7.2 & 7.3) and is given by:

RaT = a DT L3
k n
Where g is the acceleration due to gravity (m s-2), a is the coefficient of thermal expansion (oK-1), dT is the temperature difference across the fluid layer (oK), L is the vertical length of the system (m), k is the thermal diffusivity (m2 s-1) and n is the kinematic viscosity of the fluid (m2 s-1). The kinematic viscosity is the ratio of the viscosity m to the density r.
n = m
A measure of the relative thickness of the viscous and thermal boundary layers and is given by the Prandtl number, Pr :
Pr = n
Note that both the Rayleigh and Prandtl numbers are dimensionless. When the Prandtl number is large the fluid diffuses momentum relatively faster than it does heat. In such a fluid the slow diffusion of heat from a vertical boundary into the fluid produces a thin layer of hot, low viscosity thermal boundary layer adjacent to the walls. Through the rapid dissipation of momentum, the buoyancy of this layer will drag more fluid into motion with the result that the thickness of the viscous boundary layer is much greater than the thermal boundary layer. In a system in which the temperature gradient is between horizontal surfaces, the critical Rayleigh number for the onset of convection is RaT = 103.

Figure 3: Displaced packets of fluids are no longer in thermal equilibrium with the surrounds and thus gain or lose heat at a rate proportional to the thermal diffusivity of the fluid. For large diffusivities the thermal equilibration takes place quickly thereby annihilating the density contrasts that drive the convective motion.

The geometrical form of convection in 3-dimensions may vary considerably with the physical properties of the fluid. However, in all cases continuity requirements dictate a general self-similar organization. At initiation buoyancy forces induce upward motion of the heated fluid and due to the continuity requirement simultaneously cooling parts of the fluid body must begin to sink. These opposite motions obviously cannot occur at the same place and the convective motion must organize itself into convection cells. Upwelling (and downwelling) zones may either be planar or columnar in form.

Two types of fluid flow in convection can be defined; laminar or turbulent. Conditions for these two modes are described by the Reynolds number, Re (again a dimensionless parameter):

Re = w d
where w is the mean velocity of the fluid (m s-1), d is the width or diameter of the source (m), and n the kinematic viscosity (m2 s-1). If not confined to a pipe then flow becomes turbulent at values of Re > 100. Confined flow requires much higher values (Re = 2000) to become turbulent. The nature of this flow is clearly important during any sort of process of mixing since turbulence leads to much more efficient exchange of heat and material between two fluids. Convective mixing in laminar flow regimes, as appropriate to the mantle, is ultimately a process of stretching and folding. During this mixing process initial compositional heterogeneities must thin with time, eventually to the extent that diffusional exchange is capable of complete homogenisation and equilibration with the host. At temperatures close to those of the liquidi of magmas, convective flow in magma chambers will be vigorous and turbulent and under these circumstances both efficient mixing and wall-rock exchange is to be expected. As a magma body cools towards its solidus, it will change to laminar flow and eventually cease to convect.

Plate tectonics is in part a consequence of thermal convection in the mantle largely driven by radiogenic heat sources as well as residual heat generated due to the release of gravitational energy during core-mantle separation early in earth history. In contrast, convection does not occur in the interior of the geologically inactive planetary bodies like the Moon and, most probably, Mars. The rate and geometry of mantle convection control the rate of heat flow from the earth's interior to its surface, the rate at which crustal extraction takes place and the rate at which material is recycled via subduction. Any fully developed geodynamic hypothesis must therefore address the nature of mantle convective geometry in the present earth. Convection requires viscous flow or deformation. Clearly liquids or gases flow at rates observable on the scale of human patience however the asthenosphere is also capable of viscous flow albeit with very slow strain rates (of the order of 10-13 - 10-15s-1). The very high viscosity of the asthenosphere (1020 - 1022 poise) leads to the possibility of Rayleigh numbers appropriate for convection only when coupled with the very large length scales of appropriate to the mantle thickness (note that the RaT increases as a function of the cube of the length scale, but decreases only to the first power of viscosity). Convection of at least the upper part of the asthenosphere is probable on the grounds of the relative homogeneity of MORB suggesting a quite well mixed source. This conclusion is also suggested by the timescale and efficiency of extraction of the lithophile elements to the crust from the depleted upper mantle.

Systems such as the Earth with substantial density stratification organise themselves into discrete layered systems. These layers may show internal convection with little material, but significant thermal energy, exchange across layer boundaries. In this case, heat transfer between a lower denser system and the base of an overlying layer will be mediated by conduction through a thermal boundary layer. This situation occurs in the Earth at the core-mantle boundary at the so-called D'' (or double prime) layer. Here, the very dense metallic core heats the base of the silicate mantle resulting in periodic, thermally-induced, gravitational instabilities (Rayleigh-Taylor instabilities). When these instabilities achieve critical positive bouyancy they initiate the ascent of a plume of hot mantle material.

File translated from TEX by TTH, version 2.25.
On 7 Oct 2000, 11:26.