# Heat Flow

## Contents

# Heat Flow#

## Half-space Cooling Model#

Heat equation is as follow:

where \(T\) is temperature and \(\kappa\) is the thermal diffusivity. The first term represents the lateral diffusion of heat, the second term represents the vertical diffusion of heat, and the third term (on the right side) is the advection of heat by the motion of the plate.

Away from the ridge axis (\(x >> 0\)), one can show that the lateral heat diffusion is much smaller than the vertical heat diffusion. As this is a 2D problem with no heat sources so the heat equation has only diffusive and advective terms:

Dropping this term simplifies the differential equation. The boundary and initial conditions are:

The inifinite half-space has constant thermal diffusivity and an initially constant temperature \(T_{m}\). At times greater than zero, the surface temperature reduced to \(T_{0}\). The temperature will evolve with time. Note for this problem, time also corresponds to the age of the cooling oceanic lithosphere.

The solution of the half-space cooling equation is:

where \(T_{m}\) is the temperature of mantle/asthenosphere, \(T_{0}\) is the surface temperature, \(T_{z}\) is the temperature at depth \(z\), \(t\) is the age of the lithosphere, \(\kappa\) is the thermal diffusivity (~\(10^{-6}\,W\,m^{-1}K^{-1}\))

The definition of the error function \(\text{erf}(x)\):

Error function table

### Practical 3.2#

Calculate the \(60\,\text{Ma}\) geotherm \(T_{z}\) in the oceanic lithosphere for the half-space cooling model. What is the thickness of \(60\,\text{Ma}\) oceanic lithosphere? Use a surface temperature of \(0°C\); an asthenospheric temperature of \(1300°C\); and a temperature of \(1150°C\) for the base of the lithosphere. Assume a thermal conductivity \(\kappa\) of \(10^{-6}\,W\,m^{-1}K^{-1}\)

```
import math
import numpy as np
import matplotlib.pyplot as plt
```

```
# Convert years to seconds
t = 60 * 10**(6) * 365.25 * 24 * 60 * 60
T_m = 1300
T_0 = 0
T_z = 1150
k = 10**(-6)
def erf(T_m, T_0, T_z):
err_f = (T_z - T_0) / (T_m - T_0)
return err_f
# compare err_f with table
# you should obtain 1.115
z = 1.115 * 2 * math.sqrt(k * t)/1000
print("The thickness of the oceanic lithosphere is", round(z,2), "km")
```

```
The thickness of the oceanic lithosphere is 97.04 km
```

# Water Depth and Plate Models#

## Isotasy Model#

The model assumes that the lithosphere is a simply cooled asthenosphere.

The two columns shown are isostatically balanced. They have the same mass per unit area, because the older column contains water to offset the added weight of dense lithospheric rock.

## Half-space Model#

We can now use a heat conduction formulation and isostatic balancing and calculate the masses per unit area in these two columns to derive a formula for the depth of the ocean floor \(w\) as a function of age \(t\) . The solution will be applicable to columns of all ages.

Depth (\(w\)) vs. age (\(t\)) for the ocean floor:

By using the following values,

\({\rho}_{w} = 1000\,kg\,m^{-3}\) (density of water)

\({\rho}_{a} = 3300\,kg\,m^{-3}\) (density of top asthenosphere)

\({\alpha} = 3 \times 10^{-5}\,°C^{-1}\) (coefficient of thermal expansion)

\({\kappa} = 10^{-6}\,m^{2}\,s^{-1}\) (thermal conductivity)

\(T_{a} = 1300\,°C\) (temperature at the base of the plate)

\(w_{0} = 2500\,m\) (depth of ridge)

we can obtain a simplified formula as given by (Parson and Sclater, 1997):

This formula satisties observational constrains quite well for ocean crust younger than \(80\,\text{Ma}\). Note: \(w\) = depth of ocean floor = water depth

The ocean depth with age curve “flattens” for ages older than \(80\)-\(100\) \(\text{Ma}\), i.e. the observed subsidence is smaller than expected from the half-space cooling model. One possibility to explain the observed “flattening” is that by some means there is more heat beneath old lithosphere than predicted by the half-space cooling model.

Possible mechanisms include:

small scale convection,

radioactive heat decay,

phase changes, and

mantle plumes.

```
# Half space model age vs depth plot
x = np.arange(4, 170, 1)
y = (2500+350*np.sqrt(x))/1000
plt.plot(x, y, label='Half-space Model')
plt.yticks(range(2, 9))
plt.xticks(range(0, 200, 50))
plt.gca().invert_yaxis()
plt.title("Half space model age vs depth")
plt.xlabel('Age (Ma)')
plt.ylabel('Depth (km)')
plt.legend()
plt.show()
```

## Stein and Stein Plate Model#

A popular alternative to the half-space cooling model is the so-called plate (PSM) model first proposed by Parsons and Sclater (1977) and refined to GDH1 model by Stein & Stein (1992):

For the plate model beyond a critical age the lithosphere is assumed to be a constant thickness \(L\), with a constant temperature \(T_{a}\) at its base. The proposed mechanism for the required supply of heat is small-scale convection in a lower thermal boundary layer.

The subsidence is obtained by isostatic balancing:

for \(t > 20\,\text{Ma}\) where \(w_{\text{final}}\) is the equilibrium depth (Parsons and Sclater, 1977).

Note: at age \(0\), the depth is \(2.6\,km\).

```
# Half space and plate model age vs depth plot
x = np.arange(4, 170, 1)
y = (2500+350*np.sqrt(x))/1000
x2 = np.arange(20, 170, 1)
y2 = (6400 - 3200 * np.exp(-x2/62.8))/1000
plt.plot(x2, y2, label='Plate Model')
plt.plot(x, y, label='Half-space Model')
plt.yticks(range(2, 9))
plt.xticks(range(0, 200, 50))
plt.gca().invert_yaxis()
plt.title("Half space model age vs depth")
plt.xlabel('Age (Ma)')
plt.ylabel('Depth (km)')
plt.legend()
plt.show()
```

### Practical 3.3#

(a) Using the Stein and Stein half space model, calculate the difference in depth of the seabed at the intersection of a mid-ocean ridge and a transform fault. Assume that the mid-ocean ridge is spreading at \(4\,\text{cm}\,\text{yr}^{-1}\) (full rate) and that the ridge is offset \(300\,km\) by the transform fault.

```
def halfspace_model(age):
depth = (2500+350*np.sqrt(age))/1000
return depth #in km
age_ridge = 0
age_offset = 300 * 1000 * 100 / (4 / 2) / 10**6 #age in Ma
depth_difference = np.abs(halfspace_model(age_ridge) - halfspace_model(age_offset))
print("The difference in depth of the seabed at the intersection of a mid-ocean ridge and a transform fault is",
round(depth_difference, 2), "km." )
```

```
The difference in depth of the seabed at the intersection of a mid-ocean ridge and a transform fault is 1.36 km.
```

(b) Calculate the difference in depth on either side of the same fault at \(1000/1300\,km\) from the ridge axis.

```
def plate_model(age):
depth = (5.65 - 2.47 * np.exp(-age/36))
# depth = (6400 - 3200 * np.exp(-age/62.8))/1000
return depth #in km
age_1000 = 1000 * 1000 * 100 / (4 / 2) / 10**6
age_1300 = 1300 * 1000 * 100 / (4 / 2) / 10**6
depth_difference = np.abs(plate_model(age_1000) - plate_model(age_1300))
print("The difference in depth on either side of the same fault at 1000/1300 km from the ridge axis is",
round(depth_difference, 2), "km." )
```

```
The difference in depth on either side of the same fault at 1000/1300 km from the ridge axis is 0.21 km.
```

## Plate models contraints#

We still can’t discriminate between the different plate models because:

The quality of heat flow data varies substantially

The sediment thickness is often not known well enough to be corrected for

The oceanic crustal thickness is not known well enough in areas where it deviates from normal ocean crust (i.e. oceanic plateaus)

Different ocean basins may have subsided from different ridge elevations at zero age

The thermal effects of hotspots are difficult to constrain, as they varying both in time and space

### Reference#

2022 notes and practical from Lecture 3 of the module ESE 60028 Tectonics of the Ocean