TREMOR

The Elastic Wave Equation

Rock is an elastic solid. Disturb it — an earthquake rupture, an explosion, a hammer blow — and the disturbance propagates as waves governed by the elastic wave equation. For a homogeneous, isotropic solid with density ρ and Lamé parameters λ and μ, the displacement vector u(x, t) obeys:

ρ ∂²u/∂t² = (λ+2μ) ∇(∇·u) − μ ∇×(∇×u)

where:
  ρ      — density
  λ, μ   — Lamé parameters (μ = shear modulus)
  ∇·u    — divergence of displacement (compression/expansion)
  ∇×u    — curl of displacement (rotation/shear)

In 2D plane strain (displacement in x-y plane only):
  ρ ∂²ux/∂t² = (λ+2μ) ∂²ux/∂x² + μ ∂²ux/∂y² + (λ+μ) ∂²uy/∂x∂y
  ρ ∂²uy/∂t² = (λ+2μ) ∂²uy/∂y² + μ ∂²uy/∂x² + (λ+μ) ∂²ux/∂x∂y

The key insight is the Helmholtz decomposition: any vector field can be split into a curl-free (irrotational) part and a divergence-free (solenoidal) part. For elastic waves, these two parts decouple and propagate independently at different speeds. The irrotational part carries compression — these are P-waves. The solenoidal part carries shear — these are S-waves.

P-Waves and S-Waves

Applying the Helmholtz decomposition u = ∇φ + ∇×ψ (φ a scalar potential, ψ a vector potential) to the elastic wave equation separates it into two independent wave equations:

P-wave (primary, compressional):
  ∂²φ/∂t² = cP² ∇²φ       where cP = √((λ+2μ)/ρ)

S-wave (secondary, shear):
  ∂²ψ/∂t² = cS² ∇²ψ       where cS = √(μ/ρ)

Speed ratio:
  cP/cS = √((λ+2μ)/μ)  >  1  always (since λ ≥ 0)

For a Poisson solid (ν = 1/4, the Earth's crust approximation):
  λ = μ  →  cP/cS = √3 ≈ 1.732

Physical picture:
  P-wave: rock alternately compressed and rarefied along the travel direction
           (like sound in air — but much faster in solid rock)
  S-wave: rock sheared sideways at right angles to travel direction
           (cannot propagate in fluids — μ=0 for liquids)

Typical values in the crust:
  cP ≈ 5–7 km/s
  cS ≈ 3–4 km/s
  ratio ≈ 1.7–1.8 (close to the Poisson solid value √3)

In Tremor, ρ = 1 and the Lamé parameters default to λ = μ = 1, giving c_P = √3 and c_S = 1. You can drag the c_P/c_S ratio slider from 1.2 to 2.8 to see how the gap between the two rings changes — and how the time delay at the seismometer changes with it.

Visualising P and S Separately

The natural observable for P-waves is the divergence of the displacement field — ∇·u = ∂u_x/∂x + ∂u_y/∂y — which measures local expansion and compression. For S-waves, it is the curl — ∂u_y/∂x − ∂u_x/∂y — which measures local rotation. These two fields are computed live from the displacement arrays at every frame:

// Per pixel, from the displacement field ux, uy:

div  = ½(ux[col+1,row] − ux[col−1,row])   // ∂ux/∂x
     + ½(uy[col,row+1] − uy[col,row−1])   // ∂uy/∂y
     → P-wave potential

curl = ½(uy[col+1,row] − uy[col−1,row])   // ∂uy/∂x
     − ½(ux[col,row+1] − ux[col,row−1])   // ∂ux/∂y
     → S-wave potential

Colour mapping:
  div > 0 (compression):  amber-orange   (rock squeezing together)
  div < 0 (rarefaction):  deep blue-teal  (rock pulling apart)
  |curl| > 0 (shear):     teal-green / blue-violet

The coloring makes the physics directly visible: after a click, an amber ring expands rapidly outward (the P-wave compression front), followed a few moments later by a teal ring at a lower speed (the S-wave shear front). Between the two rings is a quiet zone — the elastic medium has been disturbed by the P-wave but the S-wave hasn’t arrived yet.

Locating Earthquakes — The Wadati Method

The seismometer marker (cross symbol on the canvas) records the total displacement at a fixed receiver location and displays it as a scrolling trace below the simulation. After a source fires, the P-wave arrives first, causing a small deflection in the trace. The S-wave arrives later, causing a larger deflection (S-waves typically carry more energy in the horizontal component that surface seismometers measure).

For a source at distance D from the seismometer:
  P arrival time:  tP = D / cP
  S arrival time:  tS = D / cS

  Delay:  Δt = tS − tP = D × (1/cS − 1/cP) = D × (cP − cS) / (cP × cS)

  Rearranging:  D = Δt × cP × cS / (cP − cS)

For a Poisson solid (cP = √3 · cS):
  D = Δt × cP × cS / (cP − cS) = Δt × √3/(√3−1) × cS ≈ Δt × 2.37 × cS

The Wadati diagram (1933): plot tS against tP for multiple stations.
The slope dts/dtP = cP/cS ≈ 1.73; the x-intercept gives the origin time.
With 3+ stations, triangulation gives the epicentre location.

This is how every earthquake is located today — not by GPS or satellite, but by measuring when two kinds of elastic waves arrive at a network of seismometers. The Global Seismographic Network (GSN), IRIS consortium, and regional networks worldwide run on exactly this principle, often locating magnitude-3 quakes to within a few kilometres within minutes of their occurrence.

Source Types: Fault Slip vs Explosion

The character of the waves that emerge depends on the source mechanism. Tremor offers two:

Fault slip (default):
  — Horizontal displacement impulse: ux += Gaussian, uy = 0
  — Drives both div (∂ux/∂x ≠ 0) and curl (−∂ux/∂y ≠ 0)
  — Radiation pattern: P is strongest in ±x, S is strongest in ±y
  — At 45°, both rings visible with similar amplitudes
  — Approximates the double-couple source mechanism of real faults

Explosion:
  — Radially outward displacement: ux += A·(dx/r)·G, uy += A·(dy/r)·G
  — Pure divergence source → primarily P-waves
  — Very weak S-wave (numerical noise only)
  — Approximates underground detonation, volcanic eruption
  — Used by nuclear monitoring networks to detect test blasts
    (no S-wave → "it was an explosion, not a fault rupture")

The distinction between explosion and double-couple sources is how seismologists identify nuclear weapons tests: a fault produces roughly equal P and S radiation; an underground explosion produces almost pure P radiation with a characteristic compressive first motion in all directions. The Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) uses exactly this diagnostic to verify compliance.

Real-World Applications

The physics in Tremor underlies a remarkable breadth of technology:

Earthquake early warning systems — P-waves arrive before S-waves and cause less shaking. Networks that detect P-waves can send automated alerts in the seconds before the more destructive S-waves and surface waves reach a city. Japan’s J-ALERT system, Mexico’s SASMEX, and ShakeAlert in the western United States all exploit this lag.

Seismic exploration — the oil and gas industry uses controlled explosions or vibrator trucks on the surface and records reflected P-waves to image subsurface rock layers. The reflectance at each interface depends on the acoustic impedance ρ·c_P. Multi-billion-dollar reservoirs are found from wave travel time differences measured in milliseconds.

Medical ultrasound — ultrasound in soft tissue is a P-wave (shear waves don’t propagate in fluids). Frequencies of 2–18 MHz instead of 1–100 Hz, but the same physics: pulse-echo ranging from ∂x = ½·c·Δt, where Δt is the round-trip travel time to a reflecting interface. Doppler ultrasound measures the frequency shift of waves reflected from moving blood cells.

Nuclear treaty monitoring — S-waves cannot propagate in the Earth’s outer core (liquid iron, μ=0). This “S-wave shadow zone” between 103° and 143° from an earthquake was how seismologists first proved the outer core is liquid, in 1936 (Inge Lehmann). The absence of S-waves in certain directions is a diagnostic tool at planetary scale.

Structural health monitoring — guided elastic waves (Lamb waves) travel through metal plates and can reveal cracks, delaminations, and voids by changes in wave speed and reflection patterns. Installed on aircraft structures and bridges, they form a “nervous system” that continuously monitors structural integrity.

Implementation

Tremor solves the 2D plane-strain elastic wave equation on a 200×200 grid using a second-order leapfrog finite differencescheme with proper double-buffering. The update is:

// Double-buffer leapfrog (correctness requires reading only old values):
//
// Step 1: compute ux_new and uy_new from ux, uy (current) and uxP, uyP (previous)
//         at all interior nodes simultaneously — no in-place update, no asymmetry.
//
// Step 2: swap with sponge:
//   uxP[n] = ux[n] × damp[n]
//   ux[n]  = ux_new[n] × damp[n]
//
// Stability (CFL condition):
//   dt < dx / (cP × √2)  →  dt = 0.88 / (cP × √2) ≈ 0.36 for default cP=√3
//
// Absorbing sponge: 26-cell boundary zone
//   damp[n] = 1 − 0.085 × r²  where r = 1−(d/SPONGE), d = distance to border
//   Per-step factor: ~0.915 at outer edge → total absorption over ~30 steps ≈ 97%
//
// Divergence and curl computed at render time from ux, uy by centred differences.
// Canvas: 200×200 offscreen scaled to 580×580 main canvas via drawImage().
// Seismometer: records √(ux²+uy²) at a fixed cell; scrolling Float32Array trace.

The six Float32Array fields (ux, uxP, uxN, uy, uyP, uyN) total 6 × 200 × 200 × 4 bytes = ~960 KB. The inner loop runs at 200 × 200 × 3 substeps × ~25 ops/cell ≈ 30M float ops per frame — well within browser performance at 60fps. No WebGL, no libraries, no backend. The colour LUT maps (div, curl) pairs directly to (R, G, B) with gains tuned so the wave rings are bright against the dark background without saturation.