Tectonic Simulation Programs

Interactive Python simulations running in your browser via Pyodide, plus Fortran numerical models compiled and executed on the server. Edit the code, tweak parameters, and explore the physics of plate tectonics computationally.

Oceanic lithosphere thermal evolution using the half-space cooling and plate models. Computes geotherms at different plate ages and predicts ocean depth vs age (Parsons & Sclater, 1977).

Lithosphere Cooling

Oceanic lithosphere thermal evolution using the half-space cooling and plate models. Computes geotherms at different plate ages and predicts ocean depth vs age (Parsons & Sclater, 1977).

import numpy as np
import matplotlib
matplotlib.use('AGG')
import matplotlib.pyplot as plt
from math import erf, sqrt, pi

# ==============================================================
# Oceanic Lithosphere Cooling: Half-Space vs Plate Model
# ==============================================================
# Compares two fundamental models for the thermal structure
# of oceanic lithosphere as it ages away from a mid-ocean ridge.
#
# Half-space model: T(z,t) = Ts + (Tm-Ts) * erf(z / (2*sqrt(kappa*t)))
# Plate model:      T(z,t) = Ts + (Tm-Ts) * [z/a + (2/pi) * sum(...)]
#
# The half-space model works well for young lithosphere (<80 Ma)
# but over-predicts subsidence for old lithosphere. The plate model
# (with a finite plate thickness a) flattens at large ages.

# Physical parameters
Ts = 0.0          # Surface temperature (deg C)
Tm = 1350.0       # Mantle temperature (deg C)
kappa = 1.0e-6    # Thermal diffusivity (m^2/s)
a = 125e3         # Plate thickness for plate model (m)
rho_m = 3300.0    # Mantle density (kg/m^3)
rho_w = 1030.0    # Water density (kg/m^3)
alpha = 3.1e-5    # Thermal expansivity (1/K)
d_ridge = 2500.0  # Ridge depth (m)

# Age range
ages_Ma = np.array([1, 5, 10, 20, 40, 80, 120, 180])
ages_s = ages_Ma * 1e6 * 365.25 * 24 * 3600  # Convert to seconds

# Depth array for geotherms
z = np.linspace(0, 150e3, 500)  # 0 to 150 km

print("Oceanic Lithosphere Cooling Models")
print("=" * 55)
print(f"Surface temperature:   Ts = {Ts:.0f} deg C")
print(f"Mantle temperature:    Tm = {Tm:.0f} deg C")
print(f"Thermal diffusivity:   kappa = {kappa:.1e} m^2/s")
print(f"Plate thickness (PM):  a = {a/1e3:.0f} km")
print()

# --- Figure 1: Temperature profiles ---
fig, axes = plt.subplots(1, 3, figsize=(16, 7))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.tick_params(colors='#cbd5e1')

colors = plt.cm.hot(np.linspace(0.2, 0.9, len(ages_Ma)))

# Panel 1: Half-space geotherms
ax = axes[0]
for i, (age_Ma, age_s) in enumerate(zip(ages_Ma, ages_s)):
    T_hs = np.array([Ts + (Tm - Ts) * erf(zi / (2.0 * sqrt(kappa * age_s))) for zi in z])
    ax.plot(T_hs, z / 1e3, color=colors[i], linewidth=1.5, label=f'{age_Ma} Ma')

ax.set_xlabel('Temperature (deg C)', color='white', fontsize=11)
ax.set_ylabel('Depth (km)', color='white', fontsize=11)
ax.set_title('Half-Space Model Geotherms', color='#fbbf24', fontsize=13)
ax.invert_yaxis()
ax.set_xlim(0, 1400)
ax.set_ylim(150, 0)
ax.legend(fontsize=7, facecolor='#1e293b', edgecolor='gray', labelcolor='white', ncol=2)
ax.grid(True, alpha=0.15, color='#334155')

# Panel 2: Plate model geotherms
ax = axes[1]
for i, (age_Ma, age_s) in enumerate(zip(ages_Ma, ages_s)):
    T_pm = np.zeros_like(z)
    for j, zi in enumerate(z):
        if zi <= a:
            series_sum = 0.0
            for n in range(1, 30):
                series_sum += (1.0/n) * np.sin(n * pi * zi / a) * np.exp(-n**2 * pi**2 * kappa * age_s / a**2)
            T_pm[j] = Ts + (Tm - Ts) * (zi / a + (2.0 / pi) * series_sum)
        else:
            T_pm[j] = Tm
    ax.plot(T_pm, z / 1e3, color=colors[i], linewidth=1.5, label=f'{age_Ma} Ma')

ax.set_xlabel('Temperature (deg C)', color='white', fontsize=11)
ax.set_title('Plate Model Geotherms', color='#fbbf24', fontsize=13)
ax.invert_yaxis()
ax.set_xlim(0, 1400)
ax.set_ylim(150, 0)
ax.legend(fontsize=7, facecolor='#1e293b', edgecolor='gray', labelcolor='white', ncol=2)
ax.grid(True, alpha=0.15, color='#334155')

# Panel 3: Bathymetry prediction
ax = axes[2]
ages_cont = np.linspace(0.5, 200, 500)
ages_cont_s = ages_cont * 1e6 * 365.25 * 24 * 3600

# Half-space prediction
depth_hs = d_ridge + 2 * rho_m * alpha * (Tm - Ts) * np.sqrt(kappa * ages_cont_s / pi) / (rho_m - rho_w)

# Plate model prediction (integrated thermal contraction)
depth_pm = np.zeros_like(ages_cont)
for i, t_s in enumerate(ages_cont_s):
    integral = 0.0
    dz = a / 200
    for j in range(200):
        zi = (j + 0.5) * dz
        series_sum = 0.0
        for n in range(1, 20):
            series_sum += (1.0/n) * np.sin(n * pi * zi / a) * np.exp(-n**2 * pi**2 * kappa * t_s / a**2)
        T = Ts + (Tm - Ts) * (zi / a + (2.0 / pi) * series_sum)
        integral += alpha * (Tm - T) * dz
    depth_pm[i] = d_ridge + rho_m * integral / (rho_m - rho_w)

ax.plot(ages_cont, depth_hs / 1e3, color='#f97316', linewidth=2, label='Half-space')
ax.plot(ages_cont, depth_pm / 1e3, color='#22d3ee', linewidth=2, label='Plate model')
ax.set_xlabel('Age (Ma)', color='white', fontsize=11)
ax.set_ylabel('Ocean depth (km)', color='white', fontsize=11)
ax.set_title('Predicted Bathymetry vs Age', color='#fbbf24', fontsize=13)
ax.invert_yaxis()
ax.legend(fontsize=9, facecolor='#1e293b', edgecolor='gray', labelcolor='white')
ax.grid(True, alpha=0.15, color='#334155')

# Lithosphere thickness
z_L_hs = 2.32 * np.sqrt(kappa * ages_cont_s) / 1e3
print("Lithospheric thickness (half-space, T=1100C isotherm):")
for age_check in [10, 50, 100, 180]:
    t_s = age_check * 1e6 * 365.25 * 24 * 3600
    zL = 2.32 * sqrt(kappa * t_s) / 1e3
    print(f"  {age_check:>3d} Ma: {zL:.1f} km")

plt.tight_layout()
plt.savefig('lithosphere_cooling.png', dpi=150, bbox_inches='tight',
            facecolor='#0f172a', edgecolor='none')
plt.show()
print("\nPlot generated: geotherms and bathymetry predictions.")

Click Run to execute the Python code

First run will download Python environment (~15MB)

Theoretical Background

Lithospheric Thermal Evolution

The thermal structure of oceanic lithosphere is governed by the 1D diffusion equation (neglecting advection for a static plate):

\[\frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial z^2}\]

With initial condition T(z,0) = T_m (hot mantle) and boundary conditions T(0,t) = T_s (cold surface), T(∞,t) = T_m, the solution is the half-space cooling model:

\[T(z,t) = T_s + (T_m - T_s)\,\text{erf}\!\left(\frac{z}{2\sqrt{\kappa t}}\right)\]

The predicted depth of the ocean floor increases as √t, matching observations for ages < 80 Ma. The plate model (McKenzie, 1967; Parsons & Sclater, 1977) adds a finite plate thickness that causes the bathymetry to flatten at ~6.4 km for old lithosphere. The surface heat flow decays as q(t) = k(T_m − T_s)/√(πκt), from ~200 mW/m² at ridges to ~50 mW/m² for 100 Ma lithosphere.

Plate Kinematics & Euler Poles

On a spherical Earth, the relative motion between two rigid plates is described by a rotation about an axis through Earth's center (Euler's theorem). The velocity at any point is:

\[\mathbf{v} = \boldsymbol{\omega} \times \mathbf{r}, \qquad |\mathbf{v}| = |\boldsymbol{\omega}| R \sin\Delta\]

where Δ is the angular distance from the Euler pole. Transform faults trace small circles about the Euler pole, and spreading rates vary sinusoidally with distance from it. The MORVEL model (DeMets et al., 2010) provides angular velocity vectors for 25 plates, derived from GPS, transform fault azimuths, and spreading rates at mid-ocean ridges. Plate circuits relate all plates to a common reference frame.

Mantle Convection & the Rayleigh Number

Thermal convection occurs when the buoyancy driving force exceeds viscous resistance. The Rayleigh number is the ratio of these:

\[\text{Ra} = \frac{\alpha \rho g \Delta T d^3}{\kappa \eta}, \qquad \text{Ra}_{\text{cr}} \approx 10^3 \text{ (free-slip)}\]

For the Earth's mantle, Ra ~ 10⁷, indicating vigorous convection. The non-dimensional equations in the Boussinesq approximation with infinite Prandtl number (Pr = η/ρκ → ∞) are:

\[\nabla^2 \psi = -\omega, \qquad \nabla^2 \omega = -\text{Ra}\frac{\partial T}{\partial x}, \qquad \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T = \nabla^2 T\]

The stream function ψ relates to velocity as v_x = ∂ψ/∂z, v_z = −∂ψ/∂x. The Nusselt number Nu = q_total/q_conductive measures convective efficiency; scaling laws predict Nu ~ Ra^(1/3) for vigorous convection (boundary layer theory).

Seismic Wave Propagation & Ray Theory

The elastic wave equation for an isotropic medium separates into P-waves (compressional) and S-waves (shear):

\[v_P = \sqrt{\frac{\lambda + 2\mu}{\rho}}, \qquad v_S = \sqrt{\frac{\mu}{\rho}}\]

In the high-frequency limit, seismic energy propagates along rays governed by Snell's law:

\[p = \frac{\sin\theta_i}{v_i} = \frac{\sin\theta_j}{v_j} = \text{const (ray parameter)}\]

The ray parameter p is conserved along a ray path. Total internal reflection occurs when p·v ≥ 1 in a layer, generating head waves and triplication in travel time curves. The IASP91 reference model defines the 1D velocity structure of the Earth from seismic travel time observations, revealing the Moho (35 km), transition zone (410–660 km), and core-mantle boundary (2891 km).

Lithospheric Flexure Theory

The lithosphere responds to loads as a thin elastic plate floating on a viscous fluid. For a 2D plate under a line load V₀, the deflection follows the flexure equation:

\[D\frac{d^4 w}{dx^4} + \Delta\rho\,g\,w = q(x), \qquad D = \frac{Eh^3}{12(1-\nu^2)}\]

The solution for a concentrated line load at x = 0 shows exponentially decaying oscillations:

\[w(x) = \frac{V_0 \alpha^3}{8D} e^{-|x|/\alpha}\left(\cos\frac{|x|}{\alpha} + \sin\frac{|x|}{\alpha}\right), \quad \alpha = \left(\frac{4D}{\Delta\rho\,g}\right)^{1/4}\]

The flexural parameter α determines the wavelength of bending: ~70 km for 10 Ma oceanic lithosphere (Te ~ 10 km) to ~250 km for old cratonic lithosphere (Te ~ 80 km). The forebulge— a topographic rise at ~πα from the load — is a key prediction verified by observations of moats around volcanic islands and foreland basin geometry.

Earthquake Statistics & Scaling Laws

Earthquake populations obey remarkably simple statistical laws spanning many orders of magnitude:

\[\text{Gutenberg-Richter: } \log_{10} N = a - bM, \quad b \approx 1.0\]

\[\text{Omori (aftershocks): } n(t) = \frac{K}{(t+c)^p}, \quad p \approx 1.0\text{--}1.2\]

\[\text{Moment-magnitude: } M_w = \frac{2}{3}\log_{10}(M_0) - 6.07, \quad M_0 = \mu A \bar{D}\]

The b-value near 1 implies scale-invariant (self-similar) faulting: for each magnitude unit increase, events become 10× rarer, but seismic moment (energy) increases ~31.6× (= 10^1.5). The maximum likelihood estimator (Aki, 1965) gives b = log₁₀(e)/(M̄ − M_min). Spatial b-value mapping reveals stress heterogeneity: b < 1 in locked fault patches (stress concentrations), b > 1 in volcanic regions.