5.2 Plate Tectonics

Plate tectonics is the unifying theory of Earth science. The ocean floor is created at mid-ocean ridges, moves laterally as rigid plates, and is recycled at subduction zones. This framework explains earthquakes, volcanism, mountain building, and the evolution of ocean basins through geological time.

Seafloor Spreading

Harry Hess (1962) proposed that new ocean crust forms at mid-ocean ridges and spreads away symmetrically. The Vine–Matthews–Morley hypothesis (1963) provided decisive evidence: magnetic anomaly stripes parallel to ridge axes record reversals of Earth's magnetic field, creating a symmetric "tape-recorder" pattern. The age of the ocean floor increases monotonically away from the ridge, reaching $\sim 200$ Ma in the oldest Pacific crust and $\sim 180$ Ma in the western Atlantic.

Magnetic Stratigraphy

The Geomagnetic Polarity Time Scale (GPTS) is calibrated using radiometric dating of volcanic rocks. Chrons (named intervals of constant polarity) range from $< 0.1$ Myr to $> 10$ Myr. The current normal polarity epoch (Brunhes) began at $0.78$ Ma. The half-spreading rate is computed as:

$$v_{half} = \frac{\Delta x}{\Delta t}$$

where $\Delta x$ is the distance between identified magnetic anomalies and $\Delta t$ is the time interval from the GPTS. The full spreading rate is $v_{full} = 2\,v_{half}$.

Depth–Age Relationship

As oceanic lithosphere cools and thickens while moving away from the ridge, the seafloor deepens following:

$$d(t) = d_0 + 350\sqrt{t}$$

where $d_0 \approx 2500\;\text{m}$ is the ridge crest depth, $t$ is the age in Myr, and $d$ is in metres (Parsons & Sclater, 1977). For ages beyond $\sim 70$ Myr the plate model is preferred:

$$d(t) = d_\infty - (d_\infty - d_0)\,\sum_{n=1}^{\infty} \frac{2}{n\pi}\sin\!\left(\frac{n\pi}{2}\right) e^{-n^2 t/\tau}$$

where $d_\infty \approx 6400\;\text{m}$ is the asymptotic depth and $\tau$ is the thermal time constant of the plate.

Euler Poles and Plate Rotations

On a sphere, the relative motion between two rigid plates is described by a rotation about an Euler pole. Any point on a plate moves along a small circle centred on the Euler pole. The angular velocity vector$\boldsymbol{\omega}$ fully specifies the relative motion:

$$\mathbf{v} = \boldsymbol{\omega} \times \mathbf{r}$$

where $\mathbf{v}$ is the velocity at position $\mathbf{r}$ on Earth's surface. The speed varies as $v = \omega R \sin\theta$, where $\theta$ is the angular distance from the Euler pole and $R = 6371\;\text{km}$ is Earth's radius.

Finite Rotation

The finite rotation of a plate from time $t_1$ to $t_2$ is represented by a rotation matrix $\mathbf{R}$ about the Euler pole:

$$\mathbf{R}(\hat{\mathbf{e}}, \Omega) = \cos\Omega\,\mathbf{I} + (1 - \cos\Omega)\,\hat{\mathbf{e}}\hat{\mathbf{e}}^T + \sin\Omega\,[\hat{\mathbf{e}}]_\times$$

where $\hat{\mathbf{e}}$ is the unit vector along the Euler pole axis, $\Omega$ is the total rotation angle, and $[\hat{\mathbf{e}}]_\times$ is the skew-symmetric cross-product matrix.

Triple Junctions

Three plate boundaries meet at a triple junction. The velocity triangle closure condition requires that the three relative velocity vectors sum to zero:

$$\mathbf{v}_{AB} + \mathbf{v}_{BC} + \mathbf{v}_{CA} = \mathbf{0}$$

This is equivalent to the composition of Euler vectors: $\boldsymbol{\omega}_{AB} + \boldsymbol{\omega}_{BC} + \boldsymbol{\omega}_{CA} = \mathbf{0}$. Ridge–Ridge–Ridge (RRR) junctions are always stable, while other combinations may be stable or unstable depending on geometry (McKenzie & Morgan, 1969).

Hot Spots and the Wilson Cycle

Hot Spots

Mantle plumes rise from the deep mantle (possibly the core–mantle boundary at $\sim 2900\;\text{km}$) and create volcanic chains as the plate moves over the stationary source. The Hawaiian chain extends $\sim 6000\;\text{km}$ from the active volcano (Kilauea) to the Emperor Seamounts, with a $60°$ bend at $\sim 47$ Ma recording a change in Pacific plate motion. The plate velocity over a hotspot is:

$$v_{plate} = \frac{d_{island}}{t_{island}}$$

where $d_{island}$ is the distance between dated volcanic islands.

Wilson Cycle

J. Tuzo Wilson described the cyclic opening and closing of ocean basins over $\sim 200$–$500$ Myr. Stages: (1) Continental rifting (East Africa), (2) Young ocean basin (Red Sea), (3) Mature ocean (Atlantic), (4) Declining ocean (Pacific, with subduction), (5) Terminal ocean (Mediterranean), (6) Continental collision (Himalayas). The cycle is driven by the balance between ridge push and slab pull forces.

Driving Forces

Slab pull ($F_{SP} \sim 3 \times 10^{13}\;\text{N/m}$) is the dominant driving force, arising from the negative buoyancy of cold, dense subducting lithosphere. Ridge push ($F_{RP} \sim 3 \times 10^{12}\;\text{N/m}$) results from the gravitational potential energy of elevated ridges. Basal drag from mantle convection may either drive or resist plate motion depending on the coupling.

The net force balance on a plate is:

$$F_{SP} + F_{RP} + F_{drag} + F_{resist} = 0$$

where $F_{drag}$ is the basal drag from mantle flow and $F_{resist}$ includes fault friction and bending resistance at the trench.

Subduction Zone Processes

Island Arc Volcanism

As the subducting slab descends, it releases water (from hydrated minerals) at $\sim 100$–$150\;\text{km}$ depth, lowering the melting point of the overlying mantle wedge and triggering arc volcanism. Island arcs form when oceanic crust subducts beneath oceanic crust (e.g., Mariana, Tonga). Continental arcs (e.g., Andes, Cascades) form above oceanic–continental subduction. The volcanic front lies $\sim 100$–$200\;\text{km}$ above the Wadati–Benioff seismic zone.

Thermal Parameter

The thermal state of the subducting slab is characterised by the thermal parameter:

$$\Phi = t_{age} \cdot v \cdot \sin\delta$$

where $t_{age}$ is the age of the subducting lithosphere, $v$ is the convergence rate, and $\delta$ is the slab dip angle. Slabs with $\Phi > 5000\;\text{km}$ remain cold enough to transmit deep earthquakes ($> 300\;\text{km}$).

Back-Arc Basins

Extension behind volcanic arcs produces back-arc basins with their own spreading centres (e.g., Lau Basin, Sea of Japan). They form from trench rollback—the hinge of the subducting slab retreats, stretching the overriding plate. The rollback velocity is approximately:

$$v_{rollback} \approx \frac{\Delta\rho \cdot g \cdot h^2}{2\eta}$$

where $\Delta\rho$ is the slab–mantle density contrast, $h$ is the slab length, and $\eta$ is the mantle viscosity.

Seismic Tomography & GPS Geodesy

Seismic tomography images the 3-D velocity structure of the mantle by inverting travel times from thousands of earthquakes. Fast regions (blue in tomographic images) correspond to cold subducted slabs; slow regions (red) indicate hot mantle upwellings. GPS geodesy directly measures plate velocities at $\text{mm/yr}$ precision, confirming the NUVEL-1A and MORVEL plate motion models.

Travel time residual for tomographic inversion:

$$\delta t = \int_{\text{ray}} \frac{\delta v}{v_0^2}\,ds$$

where $\delta v$ is the velocity perturbation and the integral is along the ray path $s$.

The linearised inverse problem is expressed as:

$$\mathbf{d} = \mathbf{G}\,\mathbf{m} + \mathbf{e}$$

where $\mathbf{d}$ is the vector of travel-time residuals, $\mathbf{G}$ is the kernel matrix of ray-path integrals,$\mathbf{m}$ is the model vector of slowness perturbations, and $\mathbf{e}$ represents data errors.

Python: Euler Pole Rotations & Plate Velocity Vectors

Python

!/usr/bin/env python3

script.py88 lines

#!/usr/bin/env python3
"""plate_tectonics.py - Euler poles, plate velocities, ridge geometry."""
import numpy as np
import matplotlib.pyplot as plt

# ── Constants ──
R_earth = 6371.0   # km

# ── 1. Euler Pole Rotation ──
def euler_velocity(lat, lon, pole_lat, pole_lon, omega_deg_myr):
    """
    Compute plate velocity at (lat, lon) given an Euler pole.
    Returns velocity magnitude in cm/yr and azimuth.
    """
    lat_r, lon_r = np.radians(lat), np.radians(lon)
    plat_r, plon_r = np.radians(pole_lat), np.radians(pole_lon)
    omega = np.radians(omega_deg_myr) / 1e6  # rad/yr

# Angular distance from Euler pole (spherical law of cosines)
    cos_theta = (np.sin(plat_r) * np.sin(lat_r) +
                 np.cos(plat_r) * np.cos(lat_r) * np.cos(lon_r - plon_r))
    theta = np.arccos(np.clip(cos_theta, -1, 1))

# Speed = omega * R * sin(theta)
    speed_km_yr = omega * R_earth * np.sin(theta)
    speed_cm_yr = speed_km_yr * 1e5 / 1e3  # km/yr -> cm/yr

# Azimuth (direction of motion)
    sin_az = np.cos(plat_r) * np.sin(lon_r - plon_r) / np.sin(theta)
    cos_az = (np.sin(plat_r) - np.cos(theta) * np.sin(lat_r)) / (
              np.sin(theta) * np.cos(lat_r))
    azimuth = np.degrees(np.arctan2(sin_az, cos_az))
    return speed_cm_yr, azimuth

# Africa-South America Euler pole (NUVEL-1A)
pole_lat, pole_lon, omega = 62.5, -39.4, 0.31   # deg/Myr

# Grid of points
lats = np.arange(-60, 61, 10)
lons = np.arange(-70, 31, 10)
LON, LAT = np.meshgrid(lons, lats)

speeds = np.zeros_like(LAT, dtype=float)
azimuths = np.zeros_like(LAT, dtype=float)
for i in range(LAT.shape[0]):
    for j in range(LAT.shape[1]):
        s, a = euler_velocity(LAT[i,j], LON[i,j], pole_lat, pole_lon, omega)
        speeds[i,j] = s
        azimuths[i,j] = a

# ── 2. Plot plate velocity vectors ──
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

ax = axes[0]
c = ax.contourf(LON, LAT, speeds, levels=15, cmap='hot_r')
plt.colorbar(c, ax=ax, label='Speed (cm/yr)')
U = speeds * np.sin(np.radians(azimuths))
V = speeds * np.cos(np.radians(azimuths))
ax.quiver(LON[::2,::2], LAT[::2,::2], U[::2,::2], V[::2,::2],
          color='blue', scale=30)
ax.plot(pole_lon, pole_lat, 'r*', markersize=15, label='Euler Pole')
ax.set_xlabel('Longitude'); ax.set_ylabel('Latitude')
ax.set_title('Africa-South America Relative Velocity')
ax.legend()

# ── 3. Ridge geometry: transform offsets ──
ax = axes[1]
ridge_x = []
ridge_y = []
for seg in range(6):
    y_start = seg * 200
    offset = 50 * (-1)**seg
    ridge_x.extend([offset, offset])
    ridge_y.extend([y_start, y_start + 180])
    if seg < 5:
        ridge_x.extend([offset, 50 * (-1)**(seg+1)])
        ridge_y.extend([y_start + 180, y_start + 200])

ax.plot(ridge_x, ridge_y, 'r-', lw=3, label='Ridge / Transform')
ax.set_xlabel('Offset (km)')
ax.set_ylabel('Along-axis distance (km)')
ax.set_title('Ridge-Transform Geometry')
ax.legend(); ax.grid(True, alpha=0.3)
ax.set_xlim(-100, 100)

plt.tight_layout()
plt.savefig('plate_tectonics.png', dpi=150)
plt.show()

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Thermal Model of a Subducting Slab

This Fortran program computes the 2-D temperature field within a subducting slab using a finite-difference solution to the steady-state advection–diffusion equation:

$$v_s \frac{\partial T}{\partial s} = \kappa \frac{\partial^2 T}{\partial n^2}$$

where $s$ is the along-slab coordinate, $n$ is the across-slab coordinate,$\kappa = 10^{-6}\;\text{m}^2/\text{s}$ is the thermal diffusivity. The slab enters the mantle at angle $\delta$ with velocity$v_s$ and is heated by the surrounding mantle at $T_m = 1350\;°\text{C}$.

Fortran: Thermal Model of a Subducting Slab

Fortran

============================================================

program.f9082 lines

program subducting_slab_thermal
  ! ============================================================
  ! 2-D thermal model of a subducting slab
  ! Solves steady-state advection-diffusion in slab coordinates
  !   v * dT/ds = kappa * d^2T/dn^2
  ! s = along-slab direction, n = across-slab direction
  ! ============================================================
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: ns = 300, nn = 100
  real(dp) :: kappa, Tm, Ts, vs_cm_yr, vs, slab_thick
  real(dp) :: ds, dn, s_max, n_max
  real(dp) :: T(0:ns, 0:nn)
  real(dp) :: coeff
  integer  :: i, j

! Physical parameters
  kappa      = 1.0d-6            ! thermal diffusivity (m^2/s)
  Tm         = 1350.0_dp         ! mantle temperature (C)
  Ts         = 0.0_dp            ! surface temperature (C)
  vs_cm_yr   = 8.0_dp            ! subduction rate (cm/yr)
  vs         = vs_cm_yr * 1.0d-2 / (365.25d0 * 86400.0_dp)  ! m/s
  slab_thick = 80.0d3            ! slab thickness (m)

s_max = 600.0d3                ! 600 km along-slab distance
  n_max = slab_thick
  ds    = s_max / real(ns, dp)
  dn    = n_max / real(nn, dp)
  coeff = kappa * ds / (vs * dn * dn)

! Initial condition: cold slab with error-function profile
  do j = 0, nn
    do i = 0, ns
      if (i == 0) then
        T(i, j) = Ts + (Tm - Ts) * erf_approx(real(j,dp)*dn / 30.0d3)
      else
        T(i, j) = Tm
      end if
    end do
  end do

! Boundary conditions: slab surfaces in contact with mantle
  do i = 0, ns
    T(i, 0)  = Ts + (Tm - Ts) * min(1.0_dp, real(i,dp)*ds / 200.0d3)
    T(i, nn) = Tm
  end do

! March along slab (downstream integration)
  do i = 1, ns
    do j = 1, nn - 1
      T(i, j) = T(i-1, j) + coeff * (T(i-1,j+1) - 2.0_dp*T(i-1,j) + T(i-1,j-1))
    end do
    T(i, 0)  = Ts + (Tm - Ts) * min(1.0_dp, real(i,dp)*ds / 200.0d3)
    T(i, nn) = Tm
  end do

! Output temperature field
  open(unit=10, file='slab_temperature.dat', status='replace')
  write(10, '(A)') '# s(km)  n(km)  T(C)'
  do i = 0, ns, 5
    do j = 0, nn, 2
      write(10, '(F10.1, 1X, F10.1, 1X, F10.2)') &
           real(i,dp)*ds/1.0d3, real(j,dp)*dn/1.0d3, T(i,j)
    end do
  end do
  close(10)

write(*,'(A)') 'Slab thermal model complete. Output: slab_temperature.dat'
  write(*,'(A,F6.1,A)') 'Centre temp at 300 km depth: ', T(ns/2, nn/2), ' C'

contains

real(dp) function erf_approx(x)
    real(dp), intent(in) :: x
    real(dp) :: t, s
    t = 1.0_dp / (1.0_dp + 0.3275911_dp * abs(x))
    s = t * (0.254829592_dp + t*(-0.284496736_dp + &
        t*(1.421413741_dp + t*(-1.453152027_dp + t*1.061405429_dp))))
    erf_approx = sign(1.0_dp - s * exp(-x*x), x)
  end function erf_approx

end program subducting_slab_thermal

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Plate Boundary Classification

Divergent Boundaries

Plates move apart at mid-ocean ridges. New oceanic lithosphere forms from upwelling mantle material. Spreading rates range from ultra-slow ($< 1\;\text{cm/yr}$, Gakkel Ridge in the Arctic) to super-fast (up to $18\;\text{cm/yr}$, East Pacific Rise near Easter Island). The lithospheric thickness at distance $x$ from the ridge grows as:

$$h(x) \approx 2.32\sqrt{\kappa\,x / v}$$

Continental rifting (East African Rift) represents the early stage of divergence before a new ocean basin forms.

Convergent Boundaries

Three sub-types: (1) Oceanic–oceanic convergence produces island arcs and deep trenches (Mariana, Tonga). (2) Oceanic–continental convergence produces continental volcanic arcs (Andes, Cascades). (3) Continental–continental collision produces orogenic belts (Himalayas, Alps). The subduction angle varies from shallow ($\delta \sim 10°$, Peru) to steep ($\delta \sim 90°$, Mariana), influencing volcanic arc width and the extent of back-arc extension.

Transform Boundaries

Plates slide horizontally past each other. Oceanic transforms connect offset ridge segments and are the most common type. Continental transforms (San Andreas Fault, Alpine Fault) accommodate relative plate motion without creating or destroying lithosphere. The slip rate on a transform fault is given by $v_t = \omega R \sin\theta_t$, where $\theta_t$ is the angular distance from the Euler pole to the fault.

Key Equations Summary

Ridge Push Force

$$F_{RP} = g \rho_m \alpha_v T_m \kappa t \left(1 + \frac{2\rho_m \alpha_v T_m}{\pi(\rho_m - \rho_w)}\right)$$

Slab Pull Force

$$F_{SP} = \Delta\rho \cdot g \cdot h_{slab} \cdot L_{slab}$$

where $\Delta\rho \approx 80\;\text{kg/m}^3$ is the density excess of the cold slab.

Plate Velocity from Euler Pole

$$|\mathbf{v}| = \omega R \sin\theta$$

Maximum velocity at 90° from the pole; zero at the pole itself.

← Seafloor Features Marine Sediments →