In this article, I apply basic math to derive some interesting results on the geodynamics of oceanic and continental lithosphere. In section 1, I use the half-space model to derive relations for the geothermal flow and depth of the seafloor. In section 2, I introduce McKenzie’s model of subsidence and include a simple spreadsheet that readers are free to play around with. Finally, section 3 is dedicated to continental lithosphere and the concept of heat flow province, which I illustrate with a solved example.

**1. Heat flow and oceanic depth**

Thermal conditions of an oceanic plate are different from those of a continental plate, because the material cools as it moves away from the ridge and the radiogenic heat in the crust is negligible. The simplest thermal model of oceanic lithosphere is to consider a semi-infinite half-space. If *x *and *z *are the horizontal and vertical coordinates, respectively, and *v *= *dx*/*dt *denotes the plate velocity, assumed to be constant and in practice determined by the age *t *of the crust at a distance *x *from the ridge, we may write

Further, if horizontal heat conduction is negligible compared to heat conduction by the plate’s vertical motion, the Laplacian of temperature can be simplified as

It follows that the Fourier equation simplifies to

where is thermal diffusivity. The problem of finding *T(x,z)* in a moving plate is then equivalent to finding *T(z,t)* in a stationary plate. If one considers that plate cooling is similar to a semi-infinite solid initially at temperature *T _{s}* and with surface maintained at temperature

*T*, the solution to (3) can be shown to be

_{0}where *erf *denotes the Gaussian error function. Equation (4) provides the temperature distribution along the plate with age. One immediate application of this result is the calculation of the geothermal flow *q _{0}*. Rearranging (4) gives

where *u *= *z*/(4*t*)^{1/2}. Differentiating and applying the Leibniz rule, we find the thermal gradient

Noting that geothermal flow equals the product of thermal conductivity and geothermal gradient, it follows that, with *z *= 0,

This result shows that, except for factor , the heat flux into the rock is *k *times the temperature difference (*T _{s}* –

*T*) divided by the thermal diffusion length . Equation (7) indicates that geothermal flow

_{0}*q*decreases with age (or distance from ridge axis), a trend that has been confirmed by experimental data (Figure 1). However, the correlation is somewhat poor in the case of younger material (

_{0}*t*

*40 Myr), for which the observed*

*q*is less than expected. This is due to convective cooling caused by hydrothermal circulation in the crust next to the oceanic ridges. As age increases, the sedimentary cover (less permeable than fractured basalts) prevents the seawater infiltration into the oceanic crust and consequently hinders convective cooling. For

_{0}*t*

*80 Myr, the observed*

*q*reaches an asymptotic value, given by the heat flow from the asthenosphere plus the heat flow produced by the radioactive elements of the oceanic plate.

_{0}*Figure 1. Geothermal heat flow perpendicular to ridge as a function of age. The solid curve is equation (7) and the data points are taken from Stein and Stein (1992). *

Note that an undesirable feature of equation (7) is that it yields an infinite heat flux at *t *= 0. A simple workaround would be to use the total heat *H* instead; *H *is given by integrating *q _{0} *from

*t*= 0 to

*t*,

In the mid-1800s, William Thompson, later Lord Kelvin, used the theory for the conductive cooling of a semi-infinite half-space to estimate the age of the Earth. He hypothesized that the Earth was formed at a uniform temperature *T _{1}* and that its surface was subsequently maintained at a low temperature

*T*. He assumed that a thin near-surface boundary layer developed as the Earth cooled. Since the boundary layer would be thin compared with the radius of the Earth, he reasoned that the one-dimensional model developed in this section could be applied. Using an equation similar to (7), he concluded that the age of the Earth

_{0}*t*was given by

_{0}where *dT/dz* is the present near-surface thermal gradient. With *dT*/*dz* = 25 K/km, *T _{1}* –

*T*= 2000 K, and = 1 mm

_{0}^{2}/s, the age of the Earth from equation (9) can be calculated to be

or approximately 65 million years. One of the reasons why this approximation misses the real value of 4.54 billion years by such a large margin is that it does not account for radiogenic heat. It was not until radioactivity was discovered at about 1900 that Kelvin’s estimate was seriously questioned.

Since the oceanic lithosphere is approximately in isostatic equilibrium, the average seafloor depth should increase with age as a consequence of the density increase caused by cooling. The principle of isostasy states that there is the same mass per unit area between the surface and the depth of compensation for any vertical column of material. Thus, by denoting with *h _{w}* the depth of the ocean floor,

*h*and the thickness and density of the lithosphere, and the densities of water and astenospheric mantle, we may write

_{l}or

where > because of the thermal contraction of the cooling lithosphere. Densities and temperatures are linearly related through the volume expansion coefficient ,

so that, substituting in (11), we have

Substituting the temperature profile from (7), we write

where we have changed the upper bound of the integral from *z *= *h _{l}* to

*z*= because and

*T*

*T*at the base of the lithosphere. At this point, we introduce the variable

_{s}*u*=

*z/*(4

*t*)

^{1/2}, so that

Here, we appeal to the standard result

so that

This expression shows that, if isostatic conditions prevail, the average depth of the ocean will increase with the square root of age. Figure 2 plots equation (17) for three temperatures *T _{s} *= 500

^{o}C, 1000

^{o}C, and 1500

^{o}C.

*Figure 2. Evolution of oceanic depth given by equation (17) for for three temperatures T _{s}. Values used were = 310^{-5} ^{o}C^{-1}, = 3300 kg/m^{3}, = 1030 kg/m^{3}, T_{0} = 5^{o}C, and = 30 km^{2}/Ma.*

**2. McKenzie’s (1978) subsidence model **

McKenzie (1978) used the observations of the deep structure of rift-type basins together with the exponential nature of their subsidence to suggest a kinematic model for rift-type basin evolution, which predicts subsidence and uplift, crustal thinning, and heat flow for different times since rifting. He showed that, as a consequence of extension, there is an initial subsidence due to crustal thinning followed by a thermal subsidence as the stretched lithosphere cools. Consider some initial thickness of the crust, *t _{c}*, and lithosphere,

*a*, a steady-state geothermal gradient in the lithosphere (which includes the crust) and an isothermal sub-lithospheric mantle. If the temperature at the surface of the lithosphere is 0

^{o}C and the temperature at its base is

*T*, then the

_{1}*average*temperature at the base of the crust is given by

Similarly, the *average *temperature in the sub-crustal mantle is

The pressure at the base of the *unstretched *column is given by

where

and is the coefficient of volume expansion, is the density of the lithosphere at 0^{o}C and is the density of the crust at 0^{o}C. Substituting (21) and (22) into (20) gives for pressure,

The effect of rifting is to thin the crust and lithosphere and to raise the geotherm. We can write the thickness of the thinned crust as *t _{c}*/ and the thickness of the crustal mantle as

*a*/ –

*t*/, where is the

_{c}*stretching factor*. The pressure at the base of the stretched column then becomes

where *S _{i}* is initial subsidence and is the density of water. Combining equations (23) and (24) brings to

This equation gives the initial subsidence for different values of initial thickness, temperature, density of crust and sub-crustal mantle, and stretching factor. This is a lengthy formula, so I’ve set up a simple spreadsheet to automate calculations; download it here. Using the sample values listed in Table 1, the spreadsheet returns densities = 2785 kg/m^{3}, = 3239 kg/m^{3}, and initial subsidence *S _{i}* = 2.279 km. Note that, as the crust is stretched, the densities of the crust and mantle decrease and the lithosphere subsides (

*S*> 0).

_{i}*Table 1. Data for sample subsidence calculations*

The *final subsidence *can be obtained by balancing the pressure at the base of the unstretched column with the final column. The pressure at the base of the unstretched column is given by (23). The pressure at the base of the final column is, in turn,

where is the density of the cooled stretched crust and is the density of the cooled sub-crustal mantle. Equating (23) and (26) and rearranging gives

where is the density of crust before stretching and is the density of the sub-crustal mantle before stretching; these have been defined above. The values of and can be expressed in terms of the coefficient of volume expansion and the new, cooled temperature structure; that is,

Using the same parameters as before, we obtain = 2795 kg/m^{3}, = 3251 kg/m^{3} and *S _{f}* = 4.927 km. Therefore, the densities of the crust and lithosphere increase from initial conditions and the lithosphere subsides.

The thermal subsidence, , is given by the difference between final and initial subsidence:

Substituting the values for *S _{i}* and

*S*above gives = 2.648 km. The initial and thermal subsidence are therefore similar. Watts (2001) notes that we have so far computed the subsidence at two “end member” points: one corresponding to the onset of rifting and the other after a long time has elapsed since rifting. The intermediate values of subsidence can be established by using the Fourier equation, which relates the temperature,

_{f}*T*, at a time since rifting,

*t*, to

*z,*the depth measured downwards, and , the thermal diffusivity of the lithosphere:

Since the temperature distribution determines the density structure, equation (3) can be used to calculate subsidence as a function of time. The elevation *e(t) *above the final depth to which the lithosphere (which includes the crust) sinks is given, for example, by McKenzie (1978), who proposed

where

and

The thermal subsidence, *(t)*, *since *rifting can then be obtained from

where *(0) *is the elevation at *t *= 0. *e(0) *is simply *E _{0}r*; therefore,

A plot of versus (1 – ) will have a slope of *E _{0}r.* This function depends on thermal parameters of the lithosphere as well as , the amount of extension. Since can be estimated from either Airy or flexural backstripping, the slope can be used to estimate

*directly.*This procedure has now been used to quantitatively estimate the amount of extension at a number of rift-type basins in various tectonic settings.

**3. Continental lithosphere and heat flow provinces **

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. This suggests that, like the oceanic lithosphere, continental lithosphere is cooling and slowly thickening with time.

In some specific areas, there is a linear relationship between surface heat flow and surface radioactive heat generation. 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 areas, the surface heat flow *Q _{0}* can be expressed in terms of the measured surface radioactive heat generation

*A*as

_{0}where *Q _{r}* is a constant component of heat flow from the mantle and

*D*represents a depth scale for the vertical distribution of heat-producing elements. A reliable consistency in parameter

*D*at some interval of depths, typically between 10 and 15 km, led workers to define a

*heat flow province*as a region of common tectonothermal history within which consistent values of

*Q*and

_{r}*D*are obtained. Much effort has been expended in identifying and quantifying heat flow provinces around the world. Several regional correlations are given in Figure 3. In each case, the corresponding length scale

*D*is the slope of the best-fit straight line and the reduced heat flow

*Q*is the vertical intercept of the line. These data are summarized in Table 2. Note that in each case a linear correlation appears to fit the data quite well. In all cases the length scale is near 10 km and the reduced heat flow

_{r}*Q*in general appears to be consistent with the range quoted in Pasquale

_{r}*et al.*(2014), who suggest

*Q*values ranging from 18 mW/m

_{r}^{2}for Precambrian areas to 45 mW/m

^{2}for Cenozoic areas.

*Figure 3. Dependence of surface heat flow Q _{0} on radiogenic heat production A_{0} per unit volume in surface rock for selected geological provinces: Sierra Nevada (solid squares and very long dashed line), eastern US (solid circles and intermediate dashed line), Norway and Sweden (open circles and solid line), eastern Canadian shield (open squares and short dashed line). In each case the data are fit with the linear relationship (36).*

*Table 2. Coefficients for linear fit equations shown in Figure 3.*

Although equation (36) provides a convenient, simplified description of the thermal structure of continental lithosphere, several objections have been raised over its validity in specific situations. For one, models involving heat-flow provinces often lump together regions and belts with clearly distinct origins and lithologies. Artemieva (2011) mentions the Paleoproterozoic Trans-Hudson Orogen (the central Canadian shield) as an example. The orogen is formed by markedly different structures: the Flin Flon Belt is made up of arc volcanic and plutonic rocks; the Lynn Lake Belt is formed by island arc volcanic rocks; the Thompson Nickel Belt contains metasedimentary rocks; other domains are made up of Archean gneisses. While a single correlation could be suggested for the entire orogen, heat data gathered for individual regions suggest that at least parameter *D *should be treated as variable from one belt to the next, as illustrated in Figure 4.

*Figure 4. Heat-flow province fits for regions of the central Canadian shield.*

The second main issue in the heat-province linear equation concerns the fact that it is an effectively 1D treatment and that the scatter around data is, in part, caused by 2D and 3D heat flow effects that smooth lateral variations in heat production. The effect of horizontal variations in heat production in the crust on surface heat flow and thus on the *Q*–*A _{0}* relationship was examined by harmonic (Jaupart, 1983) and random (Vasseur and Singh, 1986) lateral distributions of heat sources in the crust with a constant or exponential vertical distribution. Those workers found that, in general, lateral flow smoothes variations in heat production so that only shallow anomalies are well reflected in surface heat flow. In almost all cases, lateral heat flow results in an apparently lower

*D*-value as obtained from the heat flow-heat production plot than the actual thickness of the radioactive layer. The amplitude of this reduction depends on the radial distance over which the (randomly distributed) heat production is correlated. If the thickness of the radioactive

*D*-layer is much smaller than the scale of horizontal fluctuations in radioactivity, vertical heat flow prevails and the effect of lateral heat production heterogeneity on

*Q*and on the linear correlation (36) is insignificant. If, however, the scale of horizontal fluctuations is small compared to

*D*, the contribution of lateral heat flow to surface heat flow variations becomes important (Vasseur and Singh, 1986).

Let’s close this article with a simple applied example.

**Example**

Suppose that, in the eastern United States, a heat flow province can be represented by the usual equation

Data from the US Geological Survey suggest that the local lithosphere can be approximated by a heat flow province equation with constants *Q _{r}* = 34 mW/m

^{2}and

*D*= 7.2 km. Assuming a simple, steady-state, two-layer model of heat production, an average surface temperature of 12

^{o}C, an average thermal conductivity of 2.8 W/mK, and a measured surface heat generation of 3.2 W/m

^{3}, what is the temperature,

*T*, at the base of the heat-generation layer?

The linear relationship tells us that surface heat flow

and the two-layer model tells us that heat flow decreases linearly with depth to *Q _{r}* = 34 mW/m

^{2}at 7.2 km. The thermal gradient (

*T*/

*z*=

*Q*/

*k*) therefore decreases linearly with depth from 57.0/2.8 = 20.36

^{o}C/km at the surface to 34/2.8 = 12.14

^{o}C/km at 7.2 km; fitting these two data points to a linear equation gives

*T*/

*z*= 20.36 – 1.142

*z*. Integrating this relation with respect to

*z*brings to

As a boundary condition, *T(0) *= 12^{o}C, giving

Finally, the temperature profile is given by

and the temperature at *z *= 7.2 km is

**References**

• ARTEMIEVA, I. (2011). *The Lithosphere: An Interdisciplinary Approach*. Cambridge: Cambridge University Press.

• BEARDSMORE, G. and CULL, J. (2001). *Crustal Heat Flow*. Cambridge: Cambridge University Press.

• Jaupart, C. (1983). Horizontal heat transfer due to radioactivity contrasts: causes and consequences of the linear heat flow relation.* Geophys. J. Roy. Astron. Soc.*, 75, 411 – 435.

• McKenzie, D. (1978). Some remarks on the development of sedimentary basins. *Earth Planet. Sci. Lett.*, 40, 25 – 32.

• PASQUALE, V., VERDOYA, M. and CHIOZZI, P. (2014). *Geothermics: Heat Flow in the Lithosphere. *Berlin/Heidelberg: Springer.

• Stein, C. and Stein, S. (1992). A model for the global variation in oceanic depth and heat flow with lithospheric age. *Nature*, 359, 123 – 128.