Introduction to Geophysics

Geothermics (Heat)
by Jashodhara Chaudhury
IISER Berhampur (Course EES320)
 Worldwide heat flow: total heat loss from the Earth
● The total present-day worldwide rate of heat loss by the Earth is estimated to be (4.2–4.4) × 1013 W.
Table next shows how this heat loss is distributed by area: 71% of this heat loss occurs through the
oceans (which cover 60% of the Earth’s surface). Thus, most of the heat loss results from the creation
and cooling of oceanic lithosphere as it moves away from the mid-ocean ridges. Plate tectonics is a
primary consequence of a cooling Earth. Conversely, it seems clear that the mean rate of plate
generation is determined by some balance between the total rate at which heat is generated within
the Earth and the rate of heat loss at the surface. Some models of the thermal behaviour of the
Earth during the Archaean (before 2500 Ma) suggest that the plates were moving around the surface of
the Earth an order of magnitude faster then than they are today. Other models suggest less marked
differences from the present.
● The heat generated within the Archaean Earth by long-lived radioactive isotopes was probably
three-to-four times greater than that generated now (see Table next to next slide). A large amount of
heat also has been left over from the gravitational energy that was dissipated during accretion of the
Earth and from short-lived but energetic isotopes such as 26Al, which decayed during the first few
million years of the Earth’s history.
● Evidence from Archaean lavas that were derived from the mantle suggests that the Earth has
probably cooled by several hundred degrees since the Archaean as the original inventory of heat has
dissipated. The Earth is gradually cooling, and the plates and rates of plate generation may be
slowing to match. Presumably, after many billion years all plate motion will cease.
Heat loss and heat flow from the Earth
Table: Relative abundances of isotopes and crustal heat generation in the past relative to the present
 Figure: Observed heat flows
for the Atlantic, Indian and
Pacific Oceans. Heat flow
predicted by the plate
models: solid line, GDH1;
and dashed line, PSM. Heat
flow predicted by the
half-space model HS is not
shown – it is almost
coincident with PSM (Table
7.5).

● Measured values of heat flow depend on the age of the underlying

crust, be it oceanic or continental (Figs. above and side).
● Over the oceanic crust the heat flow generally decreases with age: the
highest and very variable measurements occur over the mid-ocean ridges
and young crust, and the lowest values are found over the deep ocean
basins.
● In continental regions the highest heat flows are measured in regions that Figure: Heat flow versus crustal
age for the continents. The
are subject to the most recent tectonic activity, and the lowest values heights of the boxes indicate the
standard deviation about the
occur over the oldest, most stable regions of the Precambrian Shield. mean heat flow, and the widths
indicate the age ranges.
 ● To apply the heat-conduction equation (below left) to the Earth as a whole, we need to use
spherical polar coordinates (r, θ , φ) (below Fig.). If temperature is not a function of θ or φ
but only of radius r, Eq. (7.15) is

● First, let us assume that there is no internal heat generation. The equilibrium temperature
is then the solution to
Figure: In spherical polar coordinates (Fig. A1.4), a point P
is located by specifying r, the radius of the sphere on which
it lies, θ, the colatitude, and φ, the longitude or azimuth,
where r ≥ 0, 0 ≤ φ ≤ 2π, 0 ≤ θ ≤ π . From Fig. A1.4 it can be
● On integrating once, we obtain seen that

● and integrating the second time gives

where c1 and c2 are the constants of integration.

 ● Now impose the boundary conditions for a hollow sphere b < r < a:

(i) zero temperature T = 0 at the surface r = a and

(ii) constant heat flow Q = −k∂ T /∂r = Q b at r = b.

● The temperature in the spherical shell b < r < a is then given by

● An expression such as this could be used to estimate a steady temperature for the lithosphere. However, since the
thickness of the lithosphere is very small compared with the radius of the Earth, (a − b)/a ⪡ 1, this solution is the
same as the solution to the one-dimensional equation (below left) with A = 0. There is no non-zero solution to Eq.
(below middle) for the whole sphere, which has a finite temperature at the origin (r = 0). However, there is a
steady-state solution to Eq. (below right) with constant internal heat generation A within the sphere:

● On integrating twice, the temperature is given by

where c1 and c2 are the constants of integration.

● Let us impose the following two boundary conditions:
(i) T finite at r = 0 and
(ii) T = 0 at r = a.
Then Eq. (side) becomes

and the heat flux is given by

The surface heat flow (at r = a) is therefore equal to Aa/3. If we model the Earth as a solid
sphere with constant thermal properties and uniform heat generation, Eqs. (above to above)
and (above) yield the temperature at the centre of this model Earth, given a value for the
surface heat flow. Assuming values of surface heat flow 80 × 10−3 W m−2 , a = 6370 km
and k = 4 W m−1 ◦ C−1 , we obtain a temperature at the centre of this model solid Earth of
● This temperature is clearly too high for the real Earth because the temperature at the visible
surface of the Sun is only about 5700 ◦ C – the Earth is not a star. The model is unrealistic
since, in fact, heat is not conducted but convected through the mantle, and the
heat-generating elements are concentrated in the upper crust rather than being uniformly
distributed throughout the Earth. These facts mean that the actual temperature at the centre
of the Earth is much lower than this estimate. Convection is important because it allows the
surface heat flow to exploit the entire internal heat of the Earth, instead of just the surface
portions of a conductive Earth.
● Another fact that we have neglected to consider is the decrease of the radioactive heat
generation with time. Equation (slide 6 right) can be solved for an exponential time decay
and a non-uniform distribution of the internal heat generation; the temperature solutions
are rather complicated and still are not applicable to the Earth because heat is convected
through mantle and outer core rather than conducted.
● It is thought that the actual temperature at the centre of the Earth is about 7000 ◦ C, on the
basis of available evidence: thermal and seismic data, laboratory behaviour of solids at high
temperatures and pressures and laboratory melting of iron-rich systems at high pressures.
 Oceanic heat flow- a. Heat flow and the depth of the oceans
● Figure shows the mean heat flow as a function of age for the three major oceans. The average heat flow is higher over young
oceanic crust but exhibits a much greater standard deviation than does that over the older ocean basins. This decrease of heat
flow with increasing age is to be expected if we consider hot volcanic material rising along the axes of the mid-ocean ridges and
plates cooling as they move away from the spreading centres. The very scattered heat-flow values measured over young oceanic
crust are a consequence of the hydrothermal circulation of sea water through the crust. Heat flow is locally high in the vicinity
of hot-water vents and low where cold sea water enters the crust. Water temperatures approaching 400 ◦ C have been measured at
the axes of spreading centres by submersibles, and the presence of hot springs on Iceland (which is located on the Reykjanes Ridge)
and other islands or regions proximal to spreading centres is well known. As will be discussed, the oceanic crust is formed by the
intrusion of basaltic magma from below. Contact with sea water causes rapid cooling of the new crust, and many cracks form in
the lava flows and dykes. Convection of sea water through the cracked crust occurs, and it is probable that this circulation penetrates
through most of the crust, providing an efficient cooling mechanism (unless you drive a Volkswagen, your car’s engine is cooled in
the same manner).

Figure: Observed heat flows

for the Atlantic, Indian and
Pacific Oceans. Heat flow
predicted by the plate
models: solid line, GDH1;
and dashed line, PSM. Heat
flow predicted by the
half-space model HS is not
shown – it is almost
coincident with PSM (Table
7.5).
● As the newly formed plate moves away from the ridge, sedimentation occurs. Deep-sea sediments
have a low permeability (Permeability and porosity are not the same. Sediments have a higher
porosity than that of crustal rocks but lack the connected pore spaces needed for high permeability.)
and, in sufficient thickness, are impermeable to sea water. In well-sedimented and therefore
generally older crust, measurements of conductive heat flow yield reliable estimates of the
actual heat flow. Another important factor affecting cessation of hydrothermal circulation is that,
in older crust, pores and cracks will become plugged with mineral deposits. As a result,
hydrothermal circulation will largely cease. Loss of heat due to hydrothermal circulation is difficult
to measure, so heat-flow estimates for young crust are generally very scattered and also
significantly lower than the theoretical estimates of heat loss (Fig.). That heat-flow measurements
are generally less than the predicted values for oceanic lithosphere younger than 65 ± 10 Ma
indicates that this is the ‘sealing age’.
Figure: Observed heat flows
for the Atlantic, Indian and
Pacific Oceans. Heat flow
predicted by the plate
models: solid line, GDH1;
and dashed line, PSM. Heat
flow predicted by the
half-space model HS is not
shown – it is almost
coincident with PSM (Table
7.5).
 ● As an oceanic plate moves away from the ridge axis and cools, it contracts and thus increases in density. If we
assume the oceanic regions to be compensated, the depth of the oceans should increase with increasing age
(and thus plate density). For any model of the cooling lithosphere, the expected ocean depth can be calculated
simply.
● Figure (a) shows the mean depth of the oceans plotted against age. For ages less than 20 Ma a simple
relation between bathymetric depth d (km) and lithosphere age t (Ma) is observed:

Depth increases linearly with the square root of age.

● For ages greater than 20 Ma this simple relation does not hold; depth increases more slowly with increasing
age and approximates a negative exponential:

Figure: Mean oceanic depth (a) and oceanic heat flow (b) with standard deviations plotted every 10 Ma against age. The data are from the north Pacific and northwest Atlantic. These global depths exclude data
from the hotspot swells. The three model predictions for ocean depth and heat flow are shown as solid and dashed lines. The plate model GDH1 fits both data sets overall better than does either the half-space
model HS or the alternative plate model PSM. Data shown in black were used to determine GDH1. Heat flow data at <50 Ma are shown in grey – these are affected by hydrothermal circulation and were not
used to determine GDH1.(c) Mean oceanic depth plotted against the square root of the age of the √lithosphere ( t). The solid line is the best-fitting half-space model: d = 2.607 + 0.344t1/2 .(After Stein and
Stein (1992) Thermo-mechanical evolution of oceanic lithosphere: implications for the subduction process and deep earthquakes (overview),
 ● Figure (b) shows the measured heat flow plotted against age. A simple relationship, linked to that for ocean
depth, between heat flow Q (10−3 W m−2 ) and lithosphere age t (Ma) is predicted for crust younger than 55
Ma:

Heat flow decreases linearly with the inverse square root of age.
● For ages greater than 55 Ma this simple relation does not hold; heat flow decreases more slowly with
increasing age and follows a negative exponential:

Figure: Mean oceanic depth (a) and oceanic heat flow (b) with standard deviations plotted every 10 Ma against age. The data are from the north Pacific and
northwest Atlantic. These global depths exclude data from the hotspot swells. The three model predictions for ocean depth and heat flow are shown as solid and
dashed lines. The plate model GDH1 fits both data sets overall better than does either the half-space model HS or the alternative plate model PSM. Data shown
in black were used to determine GDH1. Heat flow data at <50 Ma are shown in grey – these are affected by hydrothermal circulation and were not used to
determine GDH1.(c) Mean oceanic depth plotted against the square root of the age of the √lithosphere ( t). The solid line is the best-fitting half-space model: d =
2.607 + 0.344t1/2 .(After Stein and Stein (1992) Thermo-mechanical evolution of oceanic lithosphere: implications for the subduction process and deep
earthquakes (overview),
 Thermal structure of the oceanic lithosphere
● The lithospheric plate is thought to consist of two parts: an upper rigid layer and a lower viscous thermal
boundary layer (Fig.). At about 60 Ma this thermal boundary layer becomes unstable; hence small-scale
convection develops within it, resulting in an increase in heat flow to the base of the rigid layer and a thermal
structure similar to that predicted by the plate model. Very detailed, accurate measurements of heat flow,
bathymetry and the geoid on old oceanic crust and across fracture zones may improve our knowledge of the
thermal structure of the lithosphere.
Figure: A schematic diagram of the oceanic lithosphere,
showing the proposed division of the lithospheric plate.
The base of the mechanical boundary layer is the
isotherm chosen to represent the transition between rigid
and viscous behaviour. The base of the thermal
boundary layer is another isotherm, chosen to represent
correctly the temperature gradient immediately beneath
the base of the rigid plate. In the upper mantle beneath
these boundary layers, the temperature gradient is
approximately adiabatic. At about 60–70 Ma the thermal
boundary layer becomes unstable, and small-scale
convection starts to occur. With a mantle heat flow of
about 38 × 10−3 W m−2 the equilibrium thickness of
the mechanical boundary layer is approximately 90 km.
Continental heat flow- a. The mantle contribution to continental heat flow
● Continental heat flow is harder to understand than oceanic heat flow and harder
to fit into a general theory of thermal evolution of the continents or of the Earth.
Continental heat-flow values are affected by many factors, including erosion,
deposition, glaciation, the length of time since any tectonic events, local
concentrations of heat-generating elements in the crust, the presence or absence
of aquifers and the drilling of the hole in which the measurements were made.
Nevertheless, it is clear that the measured heat-flow values decrease with
increasing age (Fig.). This suggests that, like the oceanic lithosphere, the
continental lithosphere is cooling and slowly thickening with time. The mean
surface heat flow for the continents is ∼65 mW m−2 . The mean surface heat
flow in non-reactivated Archaean cratons is 41 ± 11 mW m−2 , which is Figure: Heat flow versus crustal
significantly lower than the mean value of 55 ± 17 mW m−2 for stable age for the continents. The
Proterozoic crust well away from Archaean craton boundaries. heights of the boxes indicate the
standard deviation about the
● All erosional, depositional, tectonic and magmatic processes occurring in the mean heat flow, and the widths
continental crust affect the measured surface heat-flow values. The particularly indicate the age ranges.
scattered heat-flow values measured at ages less than about 800 Ma are
evidence of strong influence of these transient processes and are therefore very
difficult to interpret in terms of the deeper thermal structure of the
continents.
● In some specific areas known as heat-flow
provinces, there is a linear relationship
between surface heat flow ( Q0) and surface
radioactive heat generation (A0)(Fig.).
Using this relationship, one can make an
approximate estimate of the contribution of
the heat-generating elements in the
continental crust to the surface heat flow. In
these heat-flow provinces, some of which are
listed in Table (next slide), the surface heat
Figure: Measured heat flow Q0 plotted against internal heat
flow Q0 can be expressed in terms of the generation A0 for (a) the eastern-U.S.A. heat-flow province.
measured surface radioactive heat The straight line Q0 = Qr + DA0 that can be fitted to these
generation A0 as measurements has Qr = 33 x 10−3 W m−2 and D = 7.5 km.
(After Roy et al. (1968).) (b) Best-fitting straight lines for
other heat-flow provinces: CA, central Australia; B, Baltic
shield; BR, Basin and Range; C, Atlantic Canada; EW,
where Qr and D are constants for each heat-flow England and Wales; EUS, eastern USA; I, India; S, Superior
Province; SN, Sierra Nevada; and Y, Yilgarn block,
province. Australia.
Some continental heat-flow provinces
● We consider here two extreme models of the distribution of the radioactive heat generation in the crust, both of which
yield a surface heat flow in agreement with this observed linear observation.
1. Heat generation is uniformly concentrated within a slab with thickness D. In this case, using Eq. (steady state
below left), we obtain

Integrating once gives

where c is the constant of integration. At the surface, z = 0, the upward heat flow Q (0) is

Therefore, the constant c is given by

● Thus, in this case, the heat flow Q(D) into the

At depth D, the upward heat flow is base of the uniform slab (and the base of the
crust, since all the heat generation is assumed to
be concentrated in the slab) is the Qr of Eq.
(below).
2. Heat generation is an exponentially decreasing function of depth within a slab of
thickness z*. Equation (below left) then becomes

where
Integrating Eq. (7.79) once gives

where c is the constant of integration. At the surface, z = 0, the heat flow is Q(0)

The constant c is given by

At depth z* (which need not be uniform throughout the heat-flow province), the heat flow is

Thus, by rearranging, we obtain

● Equation (below left) is the same as Eq. (below mid) if we write

as from the basic assumption in slide 19,

● Thus, the linear relation is valid for this model if the heat generation A(z*) at depth z* is constant
throughout the heat-flow province. Unless A(z*)D is small, the observed value of Qr may be very
different from the actual heat flow Q(z*) into the base of the layer of thickness z*.
● However, it can be shown that, for some heat-flow provinces, A(z*)D is small, and thus Qr is a
reasonable estimate of Q(z*). This removes the constraint that A(z*) must be the same throughout
the heat-flow province.
● Additionally, for those provinces in which A(z*)D is small, it can be shown that z* must be
substantially greater than D. Thus, the exponential distribution of heat production satisfies the
observed linear relationship between surface heat flow and heat generation and does so even in
cases of differential erosion.
● In this model, D is a measure of the upward migration of the heat-producing radioactive
isotopes (which can be justified on geochemical grounds), and Qr is approximately the heat flow
into the base of the crust (because z* is probably approximately the thickness of the crust).
● Neither of these models of the distribution of heat
generation within the crust allows for different vertical
distributions among the various radioactive isotopes.
There is some evidence for such variation.
Nevertheless, it is clear that much of the variation in
measured surface heat flow is caused by the
radioactive heat generation in the crust and that the
reduced heat flow Qr is a reasonable estimate of the
heat flow into the base of the crust. Figure shows this
reduced heat flow plotted against age. After about 300
Ma since the last tectonic/thermal event, the reduced
heat flow exhibits no variation and attains a value of
(25 ± 8) × 10−3 W m−2 . This is within experimental
error of the value predicted by the plate model of the
oceanic lithosphere and suggests that there should be Figure: Reduced heat flow Qr vs time since the last tectono-thermal
event for the continental heat flow provinces. The error bars
no significant difference between models of the represent the uncertainties in the data. The solid lines show the
reduced heat flows predicted by the plate model. BR and BR’, Basin
thermal structures of the oceanic and continental and Range; SEA, southeast Appalachians; SN, Sierra Nevada; EUS,
lithospheres. The present-day thermal differences are eastern U.S.A.; SP1 and SP2 , Superior Province; Bz, Brazilian
coastal shield; B, Baltic shield; BM, Bohemian massif; U, the
primarily a consequence of the age disparity between Ukraine; EW, England and Wales; N, Niger; Z and Z’, Zambia; WA,
oceanic and continental lithospheres. western Australia; CA, central Australia; EA, eastern Australia; I1
and I’ 1 Indian shield; and I2 , Archaean Indian shield
Continental heat flow-b. The temperature
structure of the continental lithosphere
● Figure shows two extreme temperature
models of the equilibrium oceanic
lithosphere, O1 and O2 , and two extreme
models of the old stable continental
lithosphere, C1 and C2 . These have been
calculated by using the one-dimensional
heat-conduction equation. The extensive
region of overlap of these four geotherms
indicates that, on the basis of surface
measurements, for depths greater than about
80 km there need be little difference in
equilibrium temperature structure beneath
oceans and continents. All the proposed
oceanic thermal models fall within the shaded
region of overlap.
Figure: (a) Extreme thermal models used to calculate equilibrium geotherms beneath an ocean, O1 and O2 , and beneath an old
stable continent, C1 and C2 . Heat flows Q0 and Qd are in mW m−2 ; heat generation A0 is in W m−3 . (b) Predicted geotherms
for these models. Thin dashed lines, oceanic geotherms; thin solid lines, continental geotherms; heavy solid line,
equilibrium geotherm for the PSM plate model, taking into account the small-scale convection occurring in the thermal
boundary layer (see Fig. slide 14). Grey shading, region of overlap. The heavy dashed line is an error function for the
geotherm of age 70 Ma. The mantle temperature Ta is taken as 1300 ◦ C.
Figure: Thermal models of the lithospheric plates beneath oceans and continents. The
dashed line is the plate thickness predicted by the PSM plate model; k (values of 2.5 and
3.3) is the conductivity in W m−1 ◦ C−1 .