Fundamental Equations — Derivations

1. Conservation of Momentum (Stokes Flow)

From the Cauchy momentum equation to creeping viscous flow in the mantle

Starting Point: Cauchy Momentum Equation

Newton's second law applied to a continuous fluid element gives the Cauchy momentum equation. For a fluid parcel of density ρ with velocity field v, the material acceleration equals the divergence of the Cauchy stress tensor σ plus body forces:

where D/Dt = ∂/∂t + v·∇ is the material (Lagrangian) derivative.

Step 1: Decompose the Stress Tensor

Split the total Cauchy stress into an isotropic pressure part and a deviatoric (viscous) part:

P is the thermodynamic pressure and τ_ij is the deviatoric stress tensor (trace-free).

Step 2: Newtonian Constitutive Law

For a Newtonian viscous fluid, the deviatoric stress is linearly proportional to the strain rate tensor:

where η is the dynamic viscosity and the strain rate tensor is the symmetric part of the velocity gradient: \(\dot{\varepsilon}_{ij} = \tfrac{1}{2}(\partial_i v_j + \partial_j v_i)\).

Step 3: Substitute into the Cauchy Equation

Inserting the decomposed stress tensor with the Newtonian constitutive law:

This is the full Navier–Stokes equation for an incompressible Newtonian fluid with spatially variable viscosity.

Step 4: The Infinite Prandtl Number Approximation

In the mantle, the Prandtl number Pr = η/(ρκ) ~ 10²³, meaning viscous forces overwhelmingly dominate inertial forces. The Reynolds number Re = ρvL/η ~ 10⁻²⁰. The inertial term ρ Dv/Dt is negligible compared to viscous and pressure terms. Setting the left-hand side to zero:

Step 5: Incompressibility Constraint

The system is supplemented by the continuity equation for an incompressible fluid:

If η is spatially constant, the viscous term simplifies to η∇²v. But the general form retains the divergence of the full viscous stress to handle spatially variable viscosity — critical in the mantle where η varies by orders of magnitude with temperature and pressure.

Physical significance: This is the master equation for geodynamic modeling. Every numerical mantle convection code (CITCOM, ASPECT, STAG3D) solves some form of this equation. The balance is between pressure gradients driving flow, viscous resistance opposing it, and gravity pulling denser material downward.

2. Heat Equation (Thermal Evolution)

From the first law of thermodynamics to the advection-diffusion equation

Starting Point: Energy Conservation

The first law of thermodynamics for a material fluid element states that the rate of change of thermal energy equals heat flux through boundaries plus internal heat generation:

where q is the heat flux vector, H is volumetric internal heating (radioactive decay of U, Th, K), and Φ is viscous dissipation.

Step 1: Fourier's Law of Heat Conduction

The heat flux is proportional to the negative temperature gradient (heat flows from hot to cold):

Substituting into the energy equation:

Step 2: Constant Thermal Conductivity

If k is spatially uniform, \(\nabla \cdot (k\nabla T) = k\nabla^2 T\).

Step 3: Expand the Material Derivative

The material (Lagrangian) derivative splits into a local (Eulerian) rate of change plus advective transport:

Step 4: Neglect Viscous Dissipation

The dissipation number Di = αgd/C_p ~ 0.5 for the mantle, meaning viscous dissipation is not always negligible. However, in the standard form it is dropped, yielding:

Step 5: Thermal Diffusivity Form

Dividing through by ρC_p and introducing the thermal diffusivity κ = k/(ρC_p):

Physical significance: This equation governs the thermal evolution of the Earth. The left-hand side: how temperature changes by local heating plus transport by mantle flow. The right-hand side: conductive smoothing and radioactive heat production. Coupling this to the Stokes equation (through temperature-dependent buoyancy) is what drives mantle convection.

3. Euler's Rotation Theorem (Plate Kinematics)

From rigid-body motion on a sphere to the kinematic description of plate motions

Starting Point: Rigidity of Tectonic Plates

Observationally, plate interiors show negligible internal deformation (< 2% of plate boundary motion). Each plate can therefore be treated as a rigid body constrained to the surface of a sphere of radius R.

Step 1: Euler's Theorem on the Sphere

Theorem (Euler, 1776): Any rigid-body displacement on the surface of a sphere that preserves the origin can be described as a single rotation about an axis passing through the sphere's center.

Proof sketch: Consider two successive infinitesimal displacements of a rigid cap on the sphere. Each defines a rotation. The composition of two rotations about different axes through the origin is itself a rotation about a third axis through the origin (closure property of SO(3), the rotation group). Therefore, any finite displacement reduces to a single rotation.

Step 2: Parameterization

The angular velocity vector ω has direction along the rotation axis (right-hand rule) and magnitude |ω| equal to the angular rate (rad/s, or °/Myr in geophysics). The axis intersects Earth's surface at the Euler pole (λ_E, φ_E).

Step 3: Velocity at a Surface Point

For a point on the plate at position vector r from Earth's center, the instantaneous velocity is the cross product:

The magnitude of the surface velocity is:

where Δ is the angular distance from the Euler pole to the point. Velocity is zero at the Euler pole and maximum 90° away.

Step 4: Relative Motion & Plate Circuits

For plates A and B, the relative angular velocity is simply the difference of their absolute angular velocities:

For a circuit of three plates A, B, C the closure condition holds (McKenzie & Morgan, 1969):

Physical significance: This is the kinematic foundation of plate tectonics. Transform faults trace small circles about the Euler pole; spreading rates vary sinusoidally with angular distance from the pole. The NUVEL-1A and MORVEL global plate motion models parameterize all relative plate motions as Euler vectors.

4. Rayleigh Number (Convective Instability)

From linear stability analysis of a heated fluid layer

Starting Point: Boussinesq Equations

Consider a fluid layer of depth d heated from below with temperature drop ΔT. The quiescent base state has a linear conductive temperature profile\(T_0(z) = T_\text{top} + \Delta T(1 - z/d)\). We linearize the coupled Stokes + heat equations about this base state by introducing perturbations T′, v′, P′.

Step 1: Linearized Equations

Dropping products of perturbation quantities, the linearized momentum equation (vertical component) is:

The linearized energy equation:

Step 2: Non-dimensionalize

Scale lengths by d, time by d²/κ, velocity by κ/d, temperature by ΔT, and pressure by ηκ/d². The single dimensionless group that emerges is:

Step 3: Normal Mode Analysis

Assume perturbations of the form \(v_z', T' \propto f(z)\,e^{i(k_x x + k_y y)}\,e^{\sigma t}\), where σ is the growth rate and k = \(\sqrt{k_x^2 + k_y^2}\) is the horizontal wavenumber. Substituting yields an eigenvalue problem for σ as a function of Ra and k.

Step 4: Critical Rayleigh Number

Marginal stability (σ = 0) defines a neutral curve Ra(k). The minimum of this curve is the critical Rayleigh number:

Step 5: Physical Interpretation

The Rayleigh number is the ratio of the timescale for thermal diffusion across the layer (d²/κ) to the timescale for viscous sinking of a buoyant anomaly (η/(αρgΔTd)):

Physical significance: For Earth's mantle, Ra ~ 10⁶–10⁷, vastly exceeding Ra_cr. The mantle must convect vigorously. Ra controls boundary layer thickness (δ ~ d·Ra^−1/3) and heat flux (Nusselt number Nu ~ Ra^1/3).

5. Half-Space Cooling Model

Solving the 1D heat diffusion equation for oceanic lithosphere cooling

Starting Point: 1D Diffusion Equation

From the heat equation (Eq. 2), eliminate advection (v = 0) and internal heating (H = 0):

Step 1: Initial & Boundary Conditions

Model oceanic lithosphere created at a mid-ocean ridge:

Step 2: Similarity Variable

There is no intrinsic length scale in the problem, so the diffusion equation admits a similarity solution. Define the dimensionless variable:

The quantity \(\sqrt{\kappa t}\) is the thermal diffusion length — the characteristic distance heat penetrates in time t.

Step 3: Transform PDE to ODE

Define the normalized temperature \(\theta(\eta) = (T - T_s)/(T_m - T_s)\). Using the chain rule for the partial derivatives:

Substituting into the diffusion equation and simplifying:

Step 4: Solve the ODE

Let φ = θ′. Then φ′ + 2ηφ = 0, which gives\(\phi = A\,e^{-\eta^2}\). Integrating:

Boundary conditions: θ(0) = 0 gives B = 0; θ(∞) = 1 gives A = 2/√π. Therefore:

Step 5: Final Solution

Step 6: Derived Quantities

Surface heat flux (differentiating and evaluating at z = 0):

Ocean depth from thermal contraction (the classic \(d \propto \sqrt{t}\) subsidence law):

Physical significance: Explains the first-order observation that ocean depth increases as √t (verified to ~80 Ma), predicts surface heat flux decay, and defines lithosphere thickness growth with age.

6. Lithospheric Flexure

From Euler–Bernoulli beam theory to the elastic plate on a fluid foundation

Step 1: Moment–Curvature Relation

For a thin elastic plate bent in one direction (cylindrical bending), the internal bending moment per unit width is related to the curvature:

Step 2: Derive the Flexural Rigidity D

Consider a plate of elastic thickness T_e. Under bending, the neutral surface has zero strain. At distance y from the neutral surface, the strain is \(\varepsilon_{xx} = -y\,d^2w/dx^2\). Using Hooke's law for plane strain: \(\sigma_{xx} = \frac{E}{1-\nu^2}\varepsilon_{xx}\). The bending moment is:

Evaluating the integral:

Step 3: Force Balance on a Plate Element

For an infinitesimal element dx, the vertical force balance and moment balance give:

Step 4: Buoyant Restoring Force

When the plate deflects downward by w, it displaces mantle (ρ_m) and is filled with material of density ρ_fill (water, sediment, or air). The net upward restoring force per unit area is:

Step 5: Combine

Differentiating the moment–curvature relation twice and substituting the equilibrium relations:

Step 6: Flexural Parameter

The homogeneous solution involves exponentially decaying oscillatory functions with characteristic length:

Physical significance: Explains trench–forebulge morphology at subduction zones, foreland basin geometry, moats around volcanic islands (Hawaii), and lithospheric response to glacial loading. The flexural parameter α ~ 100–250 km for oceanic lithosphere.

7. Seismic Wave Equation

From Newton's second law and Hooke's law to P-waves and S-waves

Step 1: Equation of Motion for a Continuum

Newton's second law for a deformable continuum (ignoring body forces for wave propagation):

Step 2: Hooke's Law for an Isotropic Solid

The stress–strain relation for a linear, isotropic elastic solid:

where \(\varepsilon_{ij} = \frac{1}{2}(\partial u_i/\partial x_j + \partial u_j/\partial x_i)\) is the infinitesimal strain tensor, λ is the first Lamé parameter, and μ is the shear modulus.

Step 3: Substitute and Simplify

Computing \(\partial\sigma_{ij}/\partial x_j\) using the constitutive law and simplifying:

Step 4: Apply Vector Identity

Using \(\nabla^2 \mathbf{u} = \nabla(\nabla \cdot \mathbf{u}) - \nabla \times (\nabla \times \mathbf{u})\):

Step 5: Helmholtz Decomposition — P-Waves and S-Waves

Decompose u into irrotational and solenoidal parts via\(\mathbf{u} = \nabla\phi + \nabla \times \boldsymbol{\psi}\)(with \(\nabla \cdot \boldsymbol{\psi} = 0\)). This yields two independent wave equations:

P-waves are compressional (v_P ~ 5–14 km/s), S-waves are shear (v_S ~ 3–7 km/s). The ratio v_P/v_S = \(\sqrt{(\lambda+2\mu)/\mu} \approx \sqrt{3}\) for a Poisson solid (ν = 0.25).

Physical significance: Predicts P-waves and S-waves. The absence of S-waves in the outer core (μ = 0) proved it is liquid. Travel-time tomography from this equation reveals the 3D velocity structure of the mantle.

8. Gutenberg–Richter Frequency–Magnitude Relation

From empirical earthquake statistics to power-law scaling and self-organized criticality

Step 1: Empirical Observation

Gutenberg & Richter (1944, 1954) compiled global earthquake catalogs and plotted the cumulative number of events N with magnitude ≥ M against magnitude on a semi-logarithmic scale. They observed a remarkably linear relationship across many orders of magnitude.

Step 2: Statistical Interpretation

Rewriting in non-logarithmic form:

The probability density function for magnitudes is obtained by differentiating:

Step 3: Connection to Seismic Moment

The moment magnitude scale relates M to seismic moment M₀:

Substituting into the GR relation:

For b = 1, this gives N ∝ M₀^−2/3 — a power-law (Pareto) distribution.

Step 4: Fault Area Scaling

Seismic moment M₀ = μAū, where A is rupture area and ū is mean slip. Empirically, ū ∝ A^1/2 (constant stress drop), so M₀ ∝ A^3/2. The GR law with b = 1 then implies:

This is a fractal size distribution consistent with self-similar fracture of a heterogeneous crust.

Step 5: Physical Basis from Self-Organized Criticality

Models of interacting faults (e.g., Bak–Tang–Wiesenfeld sandpile model, slider-block models) naturally produce power-law frequency–size distributions without parameter tuning. The fault system evolves to a critical state where stress is marginally stable everywhere, and perturbations trigger cascading failures of all sizes — reproducing the GR law.

Physical significance: The a-value measures seismicity rate; the b-value reflects the proportion of small to large events. A b-value of 1 means a tenfold decrease in count per unit increase in magnitude. Foundational for probabilistic seismic hazard analysis (PSHA).

9. Byerlee's Friction Law

From the Mohr–Coulomb criterion to empirical rock friction

Step 1: Mohr–Coulomb Failure Criterion

The general criterion for frictional sliding on a pre-existing surface is:

where τ is the shear stress on the fault plane, σ_n is the normal stress (compressive, positive), c is the cohesion, and μ_f is the coefficient of friction.

Step 2: Amontons' Law

For two surfaces in contact, friction force is proportional to normal force and independent of apparent contact area. At the microscopic level, the real contact area (sum of asperity contacts) is proportional to normal load, explaining the linear τ–σ_n relationship.

Step 3: Byerlee's Compilation (1978)

Byerlee compiled hundreds of laboratory friction experiments on diverse rock types (granite, gabbro, sandstone, limestone) and found the law is remarkably independent of rock type. Two linear regimes emerge:

Step 4: Effective Stress (Terzaghi Principle)

In the presence of pore fluid pressure P_f, the effective normal stress is:

Byerlee's law applies to the effective stress: \(\tau = 0.85(\sigma_n - P_f)\). High pore pressures dramatically weaken faults — critical for explaining the San Andreas paradox.

Step 5: Yield Strength Envelope

Combined with ductile creep laws (Equation 12) at depth, Byerlee's law defines the upper (brittle) part of the lithospheric strength envelope. The brittle–ductile transition occurs where the two curves intersect. This depth depends on temperature gradient, strain rate, and composition.

Physical significance: Sets the strength of the brittle upper crust, constrains maximum differential stress before faulting, controls earthquake stress drops, and defines the upper part of the yield strength envelope.

10. Isostasy (Airy & Pratt Models)

From Archimedes' principle to pressure balance at the compensation depth

Airy Model: Constant Density, Variable Thickness

Step 1: Define Geometry

A mountain of elevation h above sea level, underlain by crust of uniform density ρ_c. Normal crust has thickness T_c. The mountain has a crustal root of extra thickness r below.

Step 2: Pressure Balance at the Compensation Depth

At the base of the deepest root, the pressure in the mountain column equals the reference column:

Step 3: Solve for the Root

Expanding and canceling ρ_cgT_c from both sides:

With ρ_c ~ 2700 kg/m³ and ρ_m ~ 3300 kg/m³, the root is r ≈ 4.5h — about 4.5 times the topographic elevation.

Pratt Model: Constant Thickness, Variable Density

Step 1: Pressure Balance

All columns extend from their surface to a common compensation depth D below sea level. For a column at elevation h with density ρ(h):

Step 2: Solve for Density

Higher elevations correspond to lower densities.

Physical significance: Airy isostasy explains mountain roots (Himalaya: ~70 km crust under 5 km elevation). Pratt isostasy explains mid-ocean ridge topography (hot, less-dense lithosphere stands higher). Together they provide the framework for interpreting free-air and Bouguer gravity anomalies.

11. Omori's Law (Aftershock Decay)

From empirical observation to rate-and-state friction theory

Step 1: Original Omori Law (1894)

Fusakichi Omori observed that aftershock rates following the 1891 Nobi earthquake (Japan) decayed as:

where n(t) is the aftershock rate at time t, K is a productivity constant, and c regularizes the t = 0 singularity.

Step 2: Modified Omori Law (Utsu 1961)

Utsu generalized the exponent:

with p typically in the range 0.7–1.5, mean near 1.0–1.1 globally.

Step 3: Physical Derivation from Rate-and-State Friction (Dieterich, 1994)

Starting from rate-and-state friction theory, the earthquake nucleation rate on a fault population with heterogeneous initial stresses evolves after a sudden stress step Δσ (the mainshock) as:

where r is the background seismicity rate, a is the rate-and-state friction parameter, σ_n is the effective normal stress, and t_a = aσ_n/τ̇ is the aftershock duration.

For times much shorter than t_a and large stress changes (Δσ >> aσ_n):

This recovers the Omori 1/t decay, providing a physical basis for the empirical law.

Parameter Interpretation

Physical significance: The oldest quantitative law in seismology. Governs post-earthquake hazard assessment, informs operational earthquake forecasting (USGS aftershock advisories), and constrains fault-zone rheology.

12. Constitutive Law for Mantle Creep

From Arrhenius rate theory to thermally-activated dislocation and diffusion creep

Step 1: Arrhenius Rate Theory

Any thermally activated process has a rate proportional to the Boltzmann factor. Atoms and dislocations must overcome energy barriers to move; higher temperature provides more thermal energy:

Step 2: Pressure Effect

At high pressure, the activation barrier increases because atoms are more tightly packed. The effective activation energy becomes:

where E* is the activation energy at zero pressure and V* is the activation volume.

Step 3: Stress Dependence

For dislocation creep, dislocations glide and climb through the crystal. Dislocation velocity increases with stress and dislocation density also increases (Taylor hardening), giving a power-law dependence:

For diffusion creep, atoms migrate along grain boundaries or through interiors. The relationship is linear:

Step 4: Grain-Size Dependence

Diffusion creep depends on diffusion path length, which scales with grain size d:

Step 5: Assemble the Complete Law

Step 6: Effective Viscosity

Define effective viscosity as η_eff = σ/(2ε̇):

For diffusion creep (n = 1): viscosity is stress-independent (Newtonian). For dislocation creep (n > 1): viscosity decreases with stress (shear-thinning behavior).

Physical significance: Governs mantle viscosity structure and hence convection rate, postglacial rebound timescales, and plate velocities. The extreme temperature sensitivity (η changes by orders of magnitude for a few hundred K) is why the lithosphere is rigid while the asthenosphere flows.

13. Navier–Stokes with Boussinesq Approximation

From the full compressible equations to the coupled system governing mantle convection

Starting Point: Full Compressible Equations

The complete equations for a compressible, viscous, heat-conducting fluid:

Step 1: Linear Thermal Expansion

For small density variations, the equation of state is:

The key parameter: \(\delta\rho/\rho_0 = \alpha\Delta T \sim 2 \times 10^{-5} \times 2500 \approx 0.05\). Density variations are ~5% — small compared to unity.

Step 2: The Three Boussinesq Statements

Step 3: Subtract the Hydrostatic Reference

The reference state satisfies \(-\nabla P_0 + \rho_0\mathbf{g} = 0\). Defining the dynamic pressure P′ = P − P₀ and subtracting:

Step 4: Simplified Energy Equation

Under the Boussinesq approximation, the pressure-work term αT DP/Dt and viscous dissipation Φ are both O(αΔT) smaller than leading terms and are neglected:

Step 5: Infinite Prandtl Number (Mantle Limit)

For the mantle, Pr = η/(ρ₀κ) ~ 10²³. The inertial term is negligible. For constant viscosity, the viscous stress simplifies to η∇²v:

Step 6: Non-dimensionalization

Scale length by d, time by d²/κ, velocity by κ/d, temperature by ΔT, pressure by η₀κ/d². The non-dimensional momentum equation becomes:

The Rayleigh number (Equation 4) is the single parameter governing the dynamics.

Physical significance: This coupled system is the starting point for virtually all numerical mantle convection simulations. It captures the essential physics: heated fluid becomes buoyant, rises, cools at the surface, becomes dense, and sinks — the convective cycle that drives plate tectonics.