Part II: Geophysics | Chapter 4

Geodynamics

Mantle convection, viscous flow, postglacial rebound, and the geoid

2.1 Mantle Convection and the Rayleigh Number

Earth's mantle convects because internal heating and bottom heating create density inversions that are gravitationally unstable. The key dimensionless number governing the onset and vigor of thermal convection is the Rayleigh number.

Derivation 1: The Rayleigh Number from Scaling Analysis

Consider a fluid layer of depth $d$ heated from below with temperature difference $\Delta T$. The buoyancy force driving convection scales as:

$$F_{\text{buoy}} \sim \rho_0 \alpha g \Delta T$$

where $\alpha$ is thermal expansivity. This buoyancy is opposed by viscous resistance ($\sim \eta v / d^2$) and thermal diffusion ($\sim \kappa / d$). The Rayleigh number is the ratio of buoyancy to these stabilizing effects:

$$Ra = \frac{\rho_0 \alpha g \Delta T d^3}{\eta \kappa}$$

Convection begins when $Ra$ exceeds a critical value $Ra_c$. For a layer heated from below with free-slip boundaries, $Ra_c = 657.5$; for rigid boundaries, $Ra_c = 1707.8$ (Lord Rayleigh, 1916).

For Earth's mantle with $\alpha \approx 2 \times 10^{-5}$ K$^{-1}$,$\Delta T \approx 2500$ K, $d \approx 2900$ km,$\eta \approx 10^{21}$ Pa$\cdot$s, and $\kappa \approx 10^{-6}$ m$^2$/s:

$$Ra_{\text{mantle}} \approx \frac{3300 \times 2 \times 10^{-5} \times 10 \times 2500 \times (2.9 \times 10^6)^3}{10^{21} \times 10^{-6}} \approx 10^7$$

Since $Ra \gg Ra_c$, the mantle convects vigorously. The Nusselt number (ratio of total to conductive heat transport) scales as:

$$Nu \sim Ra^{1/3}$$

giving $Nu \approx 20$ for the mantle, meaning convection transports ~20 times more heat than conduction alone.

Historical Context

Arthur Holmes (1931) first proposed mantle convection as the driving mechanism for continental drift. The mathematical theory was developed by Lord Rayleigh (1916) for the onset problem. Harry Hess (1962) linked convection to seafloor spreading. Modern numerical simulations by Tackley, Zhong, Bercovici, and others now model 3D spherical convection with plates, phase transitions, and compositional heterogeneity.

2.2 Viscous Flow: The Stokes Equations

Derivation 2: The Stokes Equation for Mantle Flow

The mantle deforms as a highly viscous fluid. The Reynolds number is extremely small ($Re \sim 10^{-20}$), so inertial terms are negligible. The governing equations are the Stokes equations (conservation of momentum at zero Reynolds number):

$$-\nabla P + \nabla \cdot [\eta(\nabla \mathbf{v} + \nabla \mathbf{v}^T)] + \rho \mathbf{g} = \mathbf{0}$$

combined with the incompressibility condition:

$$\nabla \cdot \mathbf{v} = 0$$

For constant viscosity, this simplifies to:

$$-\nabla P + \eta \nabla^2 \mathbf{v} + \rho \mathbf{g} = \mathbf{0}$$

Using the Boussinesq approximation, density variations appear only in the buoyancy term: $\rho = \rho_0(1 - \alpha(T - T_0))$. The resulting thermal buoyancy force drives the flow.

The Stream Function Formulation

In 2D, the velocity can be expressed via a stream function $\psi$:$v_x = \partial\psi/\partial z$, $v_z = -\partial\psi/\partial x$. Taking the curl of the Stokes equation eliminates pressure:

$$\eta \nabla^4 \psi = \rho_0 \alpha g \frac{\partial T}{\partial x}$$

This biharmonic equation, coupled with the energy equation, forms the basis of 2D mantle convection models.

2.3 Postglacial Rebound

Derivation 3: Viscous Relaxation Time

When ice sheets melted at the end of the last glaciation (~10,000 years ago), the underlying land surface was depressed. The rate of uplift provides a direct constraint on mantle viscosity. For a surface load of wavelength $\lambda$ on a viscous half-space:

$$\tau = \frac{4\pi\eta}{\rho g \lambda}$$

where $\tau$ is the exponential relaxation time. The surface displacement decays as:

$$w(t) = w_0 \exp(-t/\tau)$$

For Fennoscandia (ice sheet diameter ~1000 km), the observed relaxation time of ~4500 years implies a mantle viscosity of:

$$\eta = \frac{\tau \rho g \lambda}{4\pi} = \frac{4500 \times 3.156 \times 10^7 \times 3300 \times 10 \times 10^6}{4\pi} \approx 10^{21} \text{ Pa}\cdot\text{s}$$

This was one of the first quantitative estimates of mantle viscosity, by Niskanen (1939) and later refined by Haskell (1935). Modern GPS measurements in Scandinavia show uplift rates of ~10 mm/yr, with ~100 m of remaining uplift. Hudson Bay is still depressed by ~150 m and rising at ~12 mm/yr.

2.4 The Geoid and Dynamic Topography

Derivation 4: Geoid Anomaly from Internal Density Perturbation

The geoid is the equipotential surface of Earth's gravity field that coincides with mean sea level. Geoid anomalies arise from internal density variations. The gravitational potential anomaly at the surface from a density perturbation $\delta\rho$ at depth $z$ is:

$$\delta U = -\frac{G}{r} \int_V \delta\rho(\mathbf{r'}) dV'$$

The geoid height anomaly is related to the potential anomaly by Bruns' formula:

$$\delta N = \frac{\delta U}{g}$$

For a thin layer of density anomaly $\delta\rho$ with thickness $\delta h$ at depth $z$, the geoid anomaly of spherical harmonic degree $l$ is:

$$\delta N_l = \frac{4\pi G R}{g(2l+1)} \left(\frac{r}{R}\right)^{l+2} \delta\rho \cdot \delta h$$

A key observation is that subducting slabs (positive density anomaly) produce negativegeoid anomalies at long wavelengths. This paradox is explained by dynamic topography: the dense slab pulls down the surface, creating a depression that more than compensates the gravitational attraction of the slab itself.

Derivation 5: Convective Velocity Scaling

The characteristic convective velocity in the mantle can be estimated from scaling analysis. Balancing buoyancy with viscous drag in the Stokes equation:

$$\rho_0 \alpha g \Delta T \sim \frac{\eta v}{d^2}$$

Solving for velocity:

$$v \sim \frac{\rho_0 \alpha g \Delta T d^2}{\eta} = \frac{\kappa}{d} Ra$$

For $Ra \sim 10^7$ and $\kappa/d \sim 3 \times 10^{-13}$ m/s, this gives$v \sim 3$ cm/yr, consistent with observed plate velocities. More refined scaling gives $v \sim (\kappa/d) Ra^{2/3}$ for internally heated convection.

2.5 Mantle Rheology

Mantle viscosity is not constant but depends strongly on temperature, pressure, grain size, and water content. The dominant creep mechanism in the mantle is diffusion creep at low stress and dislocation creep at high stress. The effective viscosity follows an Arrhenius law:

$$\eta = A^{-1} d^m \sigma^{1-n} \exp\left(\frac{E^* + PV^*}{RT}\right)$$

where $E^*$ is activation energy, $V^*$ is activation volume, $d$ is grain size, $n$ is the stress exponent (1 for diffusion, ~3.5 for dislocation creep), and $m$ is the grain size exponent. The transition from diffusion to dislocation creep depends on stress and grain size, with typical mantle grain sizes of 1-10 mm.

The viscosity profile from postglacial rebound and geoid modeling shows: upper mantle $\eta \approx 10^{20}$-$10^{21}$ Pa$\cdot$s, a viscosity increase of 10-100x across the 660 km discontinuity, and lower mantle$\eta \approx 10^{22}$-$10^{23}$ Pa$\cdot$s. The low-viscosity asthenosphere ($\eta \approx 10^{18}$-$10^{19}$ Pa$\cdot$s) beneath the lithosphere is likely due to small amounts of partial melt and/or water.

Thermal Boundary Layers

Convective heat transport creates thin thermal boundary layers at the top (lithosphere) and bottom (D$''$ layer) of the mantle. The boundary layer thickness scales as:

$$\delta \sim d \cdot Ra^{-1/3}$$

For $Ra \sim 10^7$ and $d = 2900$ km, this gives $\delta \sim 70$ km, consistent with the ~100 km thickness of the thermal lithosphere. The bottom boundary layer at the CMB is ~200 km thick (the D$''$ layer), where temperature increases by ~1500 K across a thin zone, feeding thermal instabilities that rise as mantle plumes.

The relationship between surface heat flow and boundary layer thickness is:

$$q = \frac{k\Delta T}{\delta} = \frac{k\Delta T}{d} \cdot Ra^{1/3} = \frac{k\Delta T}{d} \cdot Nu$$

For typical mantle values, this predicts a surface heat flow of ~70 mW/m$^2$, close to the observed global average of ~87 mW/m$^2$ (the difference comes from radioactive heat production in the crust). The cooling rate of the Earth can be estimated: with total heat loss ~44 TW and internal radiogenic heating ~20 TW, the mantle is cooling at approximately 50-100 K per billion years.

2.6 Plate Driving Forces and Mantle Plumes

The forces driving plate tectonics can be estimated from the torque balance on each plate. Forsyth and Uyeda (1975) showed that slab pull is the dominant driving force, correlating with the observation that plates attached to subducting slabs (Pacific, Nazca, Indian) move fastest.

Ridge Push

Ridge push arises from the gravitational potential energy difference between the elevated ridge and the adjacent seafloor. The force per unit ridge length is:

$$F_{\text{RP}} = g \rho_m \alpha \Delta T \kappa t \approx 2-3 \times 10^{12} \text{ N/m}$$

This is 5-10 times smaller than slab pull, confirming that slab pull is the primary driver.

Mantle Plumes

W. Jason Morgan (1971) proposed that mantle plumes are narrow upwellings from the deep mantle, producing hotspot volcanism (e.g., Hawaii, Iceland, Yellowstone). The buoyancy flux of a plume is:

$$B = \rho_m \alpha \Delta T_p Q_p$$

where $\Delta T_p$ is the plume temperature excess (~200-300 K) and $Q_p$ is the volume flux. Hawaii has a buoyancy flux of ~8.7 Mg/s, the largest of any plume. The age progression along the Hawaiian-Emperor chain (0-80 Ma) records the motion of the Pacific plate over the relatively fixed plume at ~7-10 cm/yr, with the famous bend at 47 Ma recording a change in Pacific plate motion direction.

Topographic swells surrounding hotspots rise ~1-2 km above the surrounding seafloor over a width of ~1000 km. The swell height constrains the plume buoyancy flux, and the excess temperature can be estimated from the depth of the basalt solidus. Seismic tomography has imaged low-velocity conduits beneath some hotspots extending to the core-mantle boundary, supporting the deep plume hypothesis.

Large Low Shear Velocity Provinces

Two continent-sized structures at the base of the mantle, known as Large Low Shear Velocity Provinces (LLSVPs), are among the most enigmatic features of Earth's interior. Located beneath Africa and the Pacific, they are characterized by $V_S$ reductions of 2-3% and extend 500-1000 km above the CMB. Their composition is debated: they may be thermochemical piles of primordial mantle material, accumulations of subducted oceanic crust, or purely thermal anomalies. Most hotspots originate from the margins of LLSVPs, suggesting these structures organize mantle flow patterns and control the locations of volcanic activity on Earth's surface over hundreds of millions of years.

The long-term stability of LLSVPs (~300+ Myr based on paleogeographic reconstructions) implies they have a compositional component that prevents entrainment into the general circulation. The buoyancy ratio $B = \Delta\rho_c / (\rho_0 \alpha \Delta T)$, where $\Delta\rho_c$ is the intrinsic density excess, must be close to unity: large enough to prevent complete entrainment but small enough to allow partial interaction with surrounding flow. This delicate balance makes LLSVPs sensitive recorders of mantle evolution since early Earth history.

Computational Lab: Geodynamics

Mantle Convection, Postglacial Rebound, and Geoid Modeling

Python

script.py290 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# GEODYNAMICS SIMULATION
# ============================================================
print("GEODYNAMICS SIMULATION")
print("=" * 70)

# ============================================================
# 1. Rayleigh Number and Convection Parameters
# ============================================================
print()
print("1. Rayleigh Number for Different Planetary Bodies")
print("-" * 70)

bodies = [
    ("Earth mantle", 3300, 2e-5, 10.0, 2500, 2.9e6, 1e21, 1e-6),
    ("Earth upper", 3400, 3e-5, 10.0, 1000, 6.7e5, 1e20, 1e-6),
    ("Moon mantle", 3400, 2e-5, 1.6, 500, 1.0e6, 1e21, 1e-6),
    ("Mars mantle", 3500, 2e-5, 3.7, 1500, 1.8e6, 1e21, 1e-6),
    ("Venus mantle", 3300, 2e-5, 8.9, 2500, 2.9e6, 1e21, 1e-6),
]

print(f"  {'Body':18} {'Ra':14} {'Ra/Ra_c':12} {'Nu':10} {'v (cm/yr)'}")
print("  " + "-" * 58)

Ra_c = 1000.0  # critical Rayleigh number
yr_to_s = 3.156e7

for name, rho, alpha, g, dT, d, eta, kappa in bodies:
    Ra = rho * alpha * g * dT * d**3 / (eta * kappa)
    Nu = max(1.0, 0.28 * Ra**(1.0/3.0)) if Ra > Ra_c else 1.0
    v = kappa / d * Ra**(2.0/3.0) * 100 * yr_to_s if Ra > Ra_c else 0.0  # cm/yr
    print(f"  {name:16} {Ra:12.2e} {Ra/Ra_c:10.1f} {Nu:8.1f} {v:8.2f}")

# ============================================================
# 2. Postglacial Rebound: Relaxation Times
# ============================================================
print()
print()
print("2. Postglacial Rebound: Relaxation Times")
print("-" * 65)

rho_m = 3300.0
g = 10.0
eta_values = [1e19, 1e20, 1e21, 1e22]

print(f"  {'Wavelength (km)':18}", end="")
for eta in eta_values:
    print(f"{'eta=' + f'{eta:.0e}':16}", end="")
print()
print("  " + "-" * 70)

for lam_km in [100, 200, 500, 1000, 2000, 3000, 5000]:
    lam_m = lam_km * 1000
    print(f"  {lam_km:16}", end="")
    for eta in eta_values:
        tau = 4 * np.pi * eta / (rho_m * g * lam_m)
        tau_kyr = tau / (1000 * yr_to_s)
        print(f"  {tau_kyr:14.1f}", end="")
    print(" kyr")

# ============================================================
# 3. Viscous Relaxation Simulation
# ============================================================
print()
print()
print("3. Fennoscandian Rebound Simulation")
print("-" * 55)

w0 = 300.0  # initial depression (m)
eta_best = 1e21  # Pa.s
lam_fenn = 1000e3  # m
tau_fenn = 4 * np.pi * eta_best / (rho_m * g * lam_fenn)
tau_kyr_fenn = tau_fenn / (1000 * yr_to_s)

print(f"  Initial depression: {w0} m")
print(f"  Viscosity: {eta_best:.0e} Pa.s")
print(f"  Relaxation time: {tau_kyr_fenn:.1f} kyr")
print()
print(f"  {'Time (kyr ago)':18} {'Remaining (m)':16} {'Uplift rate (mm/yr)'}")
print("  " + "-" * 50)

for t_kyr in [20, 15, 12, 10, 8, 6, 4, 2, 1, 0]:
    t_s = t_kyr * 1000 * yr_to_s
    # Time since deglaciation started ~18 kyr ago, rapid from ~12 kyr
    t_since = (12 - t_kyr) * 1000 * yr_to_s
    if t_since > 0:
        w = w0 * np.exp(-t_since / tau_fenn)
        dw_dt = w0 / tau_fenn * np.exp(-t_since / tau_fenn)
        rate_mm = dw_dt * yr_to_s * 1000
    else:
        w = w0
        rate_mm = 0.0
    print(f"  {t_kyr:16} {w:14.1f} {rate_mm:10.2f}")

# ============================================================
# 4. Convective Heat Transport
# ============================================================
print()
print()
print("4. Convective Heat Transport: Nusselt-Rayleigh Scaling")
print("-" * 60)

print(f"  {'Ra':14} {'Nu':10} {'q_conv (mW/m2)':18} {'BL thick (km)':16} {'v (cm/yr)'}")
print("  " + "-" * 65)

d_mantle = 2900e3  # m
kappa_m = 1e-6
k_cond = 4.0  # W/(m.K)
dT = 2500.0
alpha_m = 2e-5

for Ra in [1e4, 1e5, 1e6, 1e7, 1e8, 1e9]:
    Nu = 0.28 * Ra**(1.0/3.0)
    q_cond = k_cond * dT / d_mantle * 1000  # mW/m^2
    q_conv = Nu * q_cond
    delta_BL = d_mantle / (2 * Nu) / 1000  # km
    eta_for_Ra = rho_m * alpha_m * g * dT * d_mantle**3 / (Ra * kappa_m)
    v = kappa_m / d_mantle * Ra**(2.0/3.0) * 100 * yr_to_s
    print(f"  {Ra:12.0e} {Nu:8.1f} {q_conv:16.1f} {delta_BL:14.1f} {v:8.2f}")

# ============================================================
# 5. Geoid Anomaly from Density Perturbation
# ============================================================
print()
print()
print("5. Geoid Anomaly from Subducted Slab")
print("-" * 55)

G = 6.674e-11
R_earth = 6371e3
g_surf = 9.81

# Slab parameters
drho_slab = 80.0   # kg/m^3 density excess
h_slab = 100e3     # thickness
L_slab = 500e3     # length

print(f"  Slab: drho = {drho_slab} kg/m^3, thickness = {h_slab/1e3:.0f} km")
print()
print(f"  {'Depth (km)':14} {'dN (m)':12} {'Notes'}")
print("  " + "-" * 45)

for z_km in [100, 200, 300, 400, 500, 660, 1000, 1500, 2000]:
    z = z_km * 1000
    r = R_earth - z
    # Simplified geoid from slab (degree 2 approximation)
    # dN ~ 2*pi*G*drho*h_slab*r / (g*R)
    dN = 2 * np.pi * G * drho_slab * h_slab / g_surf * (r / R_earth)
    note = ""
    if z_km == 660:
        note = "<- 660 km discontinuity"
    print(f"  {z_km:12} {dN:10.2f} {note}")

# ============================================================
# PLOTS
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Plot 1: 2D Convection Cell (Analytical Benard)
ax = axes[0, 0]
nx, nz = 80, 40
x = np.linspace(0, 2, nx)
z = np.linspace(0, 1, nz)
X, Z = np.meshgrid(x, z)

# Stream function for single convection cell
psi = np.sin(np.pi * X) * np.sin(np.pi * Z)
# Temperature perturbation
T_pert = -np.cos(np.pi * X) * np.sin(np.pi * Z)
T_total = 1 - Z + 0.3 * T_pert

ax.contourf(X, Z, T_total, levels=20, cmap='hot')
ax.contour(X, Z, psi, levels=10, colors='white', linewidths=0.8)
# Velocity arrows
skip = 4
vx = np.cos(np.pi * X) * np.sin(np.pi * Z)
vz = -np.sin(np.pi * X) * np.cos(np.pi * Z)
ax.quiver(X[::skip, ::skip], Z[::skip, ::skip],
          vx[::skip, ::skip], vz[::skip, ::skip],
          color='cyan', alpha=0.7, scale=15)
ax.set_xlabel('Horizontal distance (d)')
ax.set_ylabel('Depth (d)')
ax.set_title('Benard Convection Cell (Temperature + Flow)')
ax.set_aspect('equal')
ax.invert_yaxis()

# Plot 2: Postglacial Rebound
ax = axes[0, 1]
t_arr = np.linspace(0, 20, 200)  # kyr since deglaciation

for eta_val, label in [(5e20, '5e20'), (1e21, '1e21'), (3e21, '3e21'), (1e22, '1e22')]:
    tau_val = 4 * np.pi * eta_val / (rho_m * g * 1000e3) / (1000 * yr_to_s)
    w_arr = w0 * np.exp(-t_arr / tau_val)
    ax.plot(t_arr, w_arr, linewidth=2, label=f'eta = {label} Pa.s')

# Observed data points (approximate)
obs_t = [0, 2, 4, 6, 8, 10, 12]
obs_w = [300, 250, 180, 120, 70, 35, 15]
ax.plot(obs_t, obs_w, 'ko', markersize=8, label='Approx observations')

ax.set_xlabel('Time since deglaciation (kyr)')
ax.set_ylabel('Remaining depression (m)')
ax.set_title('Fennoscandian Postglacial Rebound')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# Plot 3: Nusselt-Rayleigh scaling
ax = axes[1, 0]
Ra_range = np.logspace(3, 10, 100)
Nu_range = np.where(Ra_range > Ra_c, 0.28 * Ra_range**(1.0/3.0), 1.0)

ax.loglog(Ra_range, Nu_range, 'b-', linewidth=2, label='Nu ~ 0.28 Ra^(1/3)')
ax.axvline(x=Ra_c, color='red', linestyle='--', label=f'Ra_c = {Ra_c:.0f}')
ax.axvline(x=1e7, color='green', linestyle=':', label='Earth mantle')

ax.fill_between([1e6, 1e8], [1, 1], [1e4, 1e4], alpha=0.1, color='green')
ax.set_xlabel('Rayleigh Number')
ax.set_ylabel('Nusselt Number')
ax.set_title('Heat Transport Scaling')
ax.legend()
ax.grid(True, alpha=0.3, which='both')

# Plot 4: Mantle viscosity profile
ax = axes[1, 1]
z_visc = [0, 80, 80, 220, 220, 410, 410, 660, 660, 1000, 1500, 2000, 2500, 2891]
eta_visc = [1e25, 1e25, 1e19, 1e19, 1e20, 1e20, 1e21, 1e21, 3e21, 5e21, 1e22, 2e22, 5e22, 5e22]

ax.semilogx(eta_visc, z_visc, 'b-', linewidth=2)
ax.axhline(y=80, color='red', linestyle=':', alpha=0.5, label='LAB')
ax.axhline(y=410, color='orange', linestyle=':', alpha=0.5, label='410 km')
ax.axhline(y=660, color='green', linestyle=':', alpha=0.5, label='660 km')
ax.fill_betweenx([80, 220], 1e18, 1e24, alpha=0.1, color='yellow', label='Asthenosphere')

ax.set_xlabel('Viscosity (Pa.s)')
ax.set_ylabel('Depth (km)')
ax.set_title('Mantle Viscosity Profile')
ax.invert_yaxis()
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim(1e18, 1e24)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Mantle Plume Parameters")
print("-" * 60)

plumes = [
    ("Hawaii", 8.7, 300, 80, 7.0),
    ("Iceland", 3.0, 200, 20, 2.0),
    ("Reunion", 2.0, 250, 65, 4.0),
    ("Yellowstone", 1.5, 150, 16, 2.5),
    ("Kerguelen", 0.5, 200, 120, 1.0),
    ("Galapagos", 1.0, 150, 15, 5.0),
    ("Canary", 0.3, 100, 70, 1.0),
    ("Tristan", 1.7, 200, 130, 2.0),
]

print(f"  {'Plume':16} {'B (Mg/s)':12} {'dT (K)':10} {'Age (Ma)':12} {'v_plate (cm/yr)'}")
print("  " + "-" * 58)

for name, B, dT, age, v in plumes:
    print(f"  {name:14} {B:10.1f} {dT:8} {age:10} {v:10.1f}")

# Swell height calculation
print()
print("  Predicted swell heights (Crough, 1983 scaling):")
print("  " + "-" * 45)
for name, B, dT, age, v in plumes[:5]:
    # h_swell ~ B / (rho * alpha * g * v * w)
    rho_m = 3300
    alpha_m = 2e-5
    g_val = 10
    v_ms = v / 100 / 3.156e7  # to m/s
    w_plume = 200e3  # plume width
    h_swell = B * 1e3 / (rho_m * v_ms * w_plume) / 1000  # rough scaling, km
    h_swell = min(h_swell, 2.5)
    print(f"    {name}: swell height ~ {h_swell:.1f} km")

print()
print()
print("[Chart saved: convection cell, rebound, Nu-Ra scaling, viscosity profile]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Lithosphere Strength Envelopes and Plate Force Balance

Python

script.py91 lines

import numpy as np

# ============================================================
# LITHOSPHERE STRENGTH AND PLATE FORCES
# ============================================================
print("LITHOSPHERE STRENGTH AND PLATE FORCE ANALYSIS")
print("=" * 70)

g = 9.81
rho_c = 2700  # crust density
rho_m = 3300  # mantle density

# ============================================================
# 1. Brittle-Ductile Transition Depth
# ============================================================
print()
print("1. Continental Lithosphere Strength Envelope")
print("-" * 65)

# Byerlee's law for extension: delta_sigma = 0.75 * rho*g*z
# Power-law creep: delta_sigma = (edot/A)^(1/n) * exp(E/(nRT))

edot = 1e-15  # typical geological strain rate (s^-1)

# Quartz (wet): A=1e-3, n=4, E=135 kJ/mol
# Olivine (dry): A=7e-4, n=3.5, E=510 kJ/mol
R = 8.314

print(f"  Strain rate = {edot:.0e} s^-1")
print()
print(f"  {'Depth (km)':14} {'T (C)':10} {'Brittle (MPa)':16} {'Ductile (MPa)':16} {'Strength (MPa)'}")
print("  " + "-" * 65)

geotherm = 20.0  # C/km in upper lithosphere
T_surface = 15.0

for z_km in range(0, 121, 10):
    if z_km == 0:
        z_km = 1
    z = z_km * 1000
    T_C = T_surface + geotherm * z_km
    T_K = T_C + 273.15

# Brittle (extension regime)
    if z_km <= 35:
        sigma_brittle = 0.75 * rho_c * g * z / 1e6
    else:
        sigma_brittle = 0.75 * rho_m * g * z / 1e6

# Ductile (quartz for z<35km, olivine for z>35km)
    if z_km <= 35:
        A_q = 1e-3
        n_q = 4.0
        E_q = 135e3
        sigma_ductile = (edot / A_q)**(1/n_q) * np.exp(E_q / (n_q * R * T_K)) / 1e6
    else:
        A_o = 7e-4
        n_o = 3.5
        E_o = 510e3
        sigma_ductile = (edot / A_o)**(1/n_o) * np.exp(E_o / (n_o * R * T_K)) / 1e6

sigma_ductile = min(sigma_ductile, 5000)
    strength = min(sigma_brittle, sigma_ductile)

print(f"  {z_km:12} {T_C:8.0f} {sigma_brittle:14.0f} {sigma_ductile:14.0f} {strength:10.0f}")

# ============================================================
# 2. Integrated Lithosphere Strength
# ============================================================
print()
print()
print("2. Integrated Lithosphere Strength by Tectonic Setting")
print("-" * 55)

settings = [
    ("Old oceanic (100 Ma)", 3.0, 50, 1e13),
    ("Young oceanic (10 Ma)", 3.0, 16, 2e12),
    ("Continental shield", 20.0, 200, 5e13),
    ("Continental rift", 35.0, 60, 5e11),
    ("Orogenic belt", 25.0, 100, 2e12),
]

print(f"  {'Setting':28} {'Geotherm (C/km)':18} {'Mech. thick (km)':18} {'F_int (N/m)'}")
print("  " + "-" * 65)

for name, geo, thick, F_int in settings:
    print(f"  {name:26} {geo:16.1f} {thick:16} {F_int:12.1e}")

print()
print("Complete: strength envelope analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Earth's Interior Next: Climate Science →