Stick-Slip Earthquakes: Gutenberg-Richter

drag Y for loading rate · click to nudge stress

A 1D Burridge-Knopoff chain of $N = 30$ blocks rests on a frictional surface. Each block is coupled to its two neighbors by springs of stiffness $k_{c}$ , and pulled forward by a slowly-creeping tectonic plate through a loading spring of stiffness $k_{ℓ}$ . The instantaneous force on block $i$ is

F_{i} = k_{ℓ} (X_{plate} - x_{i}) + k_{c} (x_{i - 1} - 2 x_{i} + x_{i + 1}) .

Each block carries an independent static-friction threshold

F_{s}^{(i)}

drawn from a narrow distribution — the in-built **heterogeneity** of a real fault zone. Blocks **stick** as long as

∣ F_{i} ∣ \leq F_{s}^{(i)}

and **slip** when the loading plate has wound the springs tight enough to exceed it. A slip releases stored elastic energy by jumping the block until its force drops to

- sign (F) η F_{s}^{(i)}

(dynamic-friction undershoot,

η < 1

); that sudden displacement perturbs the neighbors via the coupling springs and can push *them* over threshold too. The cascade of slips that propagates from one nucleation site **is** the earthquake; its size

s

— number of blocks that slipped before the chain quiets — is a proxy for seismic moment. Bin events by magnitude

M = lo g_{10} (s)

and accumulate a histogram of counts

N (M)

. The hallmark observation of seismology, the **Gutenberg-Richter law**, is that

lo g_{10} N = a - b M,

with

b \approx 1

across most tectonic regions. The same line emerges here from nothing but threshold heterogeneity + nearest-neighbor coupling + slow uniform loading; the running least-squares fit is drawn live in yellow over the histogram bars as the catalog grows. The block strip up top is colored by current

∣ F_{i} ∣/ F_{s}^{(i)}

— cool teal is well below threshold, hot red is right at the brink — and lights up white at the moment of slip so you can watch ruptures propagate along the chain. **Interaction:** dragging the mouse from top to bottom of the canvas scrubs the tectonic loading rate

v_{ℓ}

from glacially slow (rare, system-spanning quakes) to fast (frequent small ones, smaller

b

). Clicking a block injects a localized stress perturbation — manual nucleation of a quake at the cursor.

idle

331 lines · vanilla

view source

// Burridge-Knopoff stick-slip earthquake model.
// A 1D chain of N blocks sits on a frictional surface, connected to neighbors
// by springs of stiffness kC. Each block is also tied to a "loading plate"
// that drifts forward at velocity vLoad via a loading spring of stiffness kL.
// Force on block i: F_i = kL*(plate - x_i) + kC*(x_{i-1} - 2 x_i + x_{i+1}).
// Each block has a per-block static-friction threshold Fs_i (lognormal-ish).
// When |F_i| > Fs_i the block slips: we release the stress instantaneously
// by setting x_i so the net force drops to a fraction of Fs_i (overshoot),
// and that displacement perturbs the neighbors' forces — potentially
// triggering them too. The cascade of slips in one timestep is one quake;
// its size = number of blocks that slipped (a proxy for seismic moment).
//
// Magnitude M = log10(size). We bin events and accumulate a histogram;
// plotted log10(N) vs M it forms the Gutenberg-Richter line  log N = a - b M
// with b ~ 1 over geological time. We also fit a / b by least-squares to
// the populated bins and draw the line live.

const N = 30;              // number of blocks
const KC = 1.0;            // coupling spring stiffness (block <-> block)
const KL = 0.05;           // loading spring stiffness (plate -> block)
const FS_MEAN = 1.0;       // mean static friction threshold
const FS_SPREAD = 0.18;    // gaussian spread (block-to-block heterogeneity)
const OVERSHOOT = 0.25;    // fraction of Fs released past zero on slip (dynamic friction undershoot)
const MAX_CASCADE = 4000;  // safety cap on cascade iterations per frame

// Histogram: log10(size) bins. size in [1, N], so M in [0, ~1.5].
const N_BINS = 14;
const M_MIN = 0.0;
const M_MAX = 1.55;

let W, H;
let x;            // Float32Array(N): block displacements
let fs;           // Float32Array(N): static friction thresholds
let stress;       // Float32Array(N): last computed |F_i| for color
let plate;        // scalar: plate displacement
let vLoad;        // current loading rate (set from mouseY)
let timeAcc;      // total sim time for stats
let hist;         // Int32Array(N_BINS): event count per magnitude bin
let totalEvents;
let totalSlips;
let lastQuakeSize; // last cascade size (for HUD)
let lastQuakeTime; // when (sim time) it happened, for the flash
let recentQuakes;  // ring buffer of [time, size] for the timeline
const RECENT_LEN = 60;
let recentHead;
let pendingNudge; // {idx, amount} or null — set by clicks

function gauss() {
  // Box-Muller, clamped.
  let u = 0;
  while (u === 0) u = Math.random();
  const v = Math.random();
  const z = Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
  return Math.max(-2.5, Math.min(2.5, z));
}

function init({ width, height, ctx }) {
  W = width; H = height;
  x = new Float32Array(N);
  fs = new Float32Array(N);
  stress = new Float32Array(N);
  for (let i = 0; i < N; i++) {
    x[i] = 0;
    fs[i] = Math.max(0.35, FS_MEAN + FS_SPREAD * gauss());
  }
  plate = 0;
  vLoad = 0.02;
  timeAcc = 0;
  hist = new Int32Array(N_BINS);
  totalEvents = 0;
  totalSlips = 0;
  lastQuakeSize = 0;
  lastQuakeTime = -10;
  recentQuakes = new Float32Array(RECENT_LEN * 2); // (time, size) pairs
  recentHead = 0;
  pendingNudge = null;

  ctx.fillStyle = '#06080d';
  ctx.fillRect(0, 0, W, H);
}

// Force on block i given current state.
function forceOn(i) {
  let f = KL * (plate - x[i]);
  if (i > 0) f += KC * (x[i - 1] - x[i]);
  if (i < N - 1) f += KC * (x[i + 1] - x[i]);
  return f;
}

// Run one global cascade of slips until quiescent. Returns slip count.
// On slip we release displacement so the block's net force drops to
// -sign(F)*OVERSHOOT*Fs (a small dynamic undershoot). This re-distributes
// force into the neighbors and can chain.
function relax() {
  let count = 0;
  for (let iter = 0; iter < MAX_CASCADE; iter++) {
    // Find the most overstressed block (largest excess over Fs).
    let worst = -1;
    let worstExcess = 0;
    for (let i = 0; i < N; i++) {
      const f = forceOn(i);
      const excess = Math.abs(f) - fs[i];
      if (excess > worstExcess) {
        worstExcess = excess;
        worst = i;
      }
    }
    if (worst < 0) break; // quiescent

    // Slip block `worst`. Solve for new x[worst] s.t. force on it
    // becomes -sign(F)*OVERSHOOT*Fs[worst].
    // F(x) = KL*(plate - x) + KC*(neighbors' x - 2x) but with x_neighbors
    //        fixed, F is linear in x[worst] with slope -(KL + KC*deg).
    const deg = (worst > 0 ? 1 : 0) + (worst < N - 1 ? 1 : 0);
    const slope = -(KL + KC * deg);
    const fNow = forceOn(worst);
    const fTarget = -Math.sign(fNow) * OVERSHOOT * fs[worst];
    // fTarget = fNow + slope * dx  =>  dx = (fTarget - fNow) / slope
    const dx = (fTarget - fNow) / slope;
    x[worst] += dx;
    count++;
  }
  return count;
}

function magnitudeOf(size) {
  return Math.log10(size);
}

function binOf(mag) {
  const t = (mag - M_MIN) / (M_MAX - M_MIN);
  if (t < 0) return 0;
  if (t >= 1) return N_BINS - 1;
  return Math.floor(t * N_BINS);
}

// Least-squares fit of  log10(N) = a - b * M  to populated bins only.
// Returns {a, b} or null.
function fitGR() {
  // Single-pass least-squares; no per-frame arrays allocated.
  let n = 0, sx = 0, sy = 0, sxx = 0, sxy = 0;
  for (let k = 0; k < N_BINS; k++) {
    if (hist[k] > 0) {
      const M = M_MIN + (k + 0.5) * (M_MAX - M_MIN) / N_BINS;
      const y = Math.log10(hist[k]);
      n++; sx += M; sy += y; sxx += M * M; sxy += M * y;
    }
  }
  if (n < 3) return null;
  const denom = n * sxx - sx * sx;
  if (Math.abs(denom) < 1e-9) return null;
  const slope = (n * sxy - sx * sy) / denom; // -b
  const intercept = (sy - slope * sx) / n;   // a
  return { a: intercept, b: -slope };
}

function tick({ ctx, dt, width, height, input }) {
  if (width !== W || height !== H) { W = width; H = height; }

  // --- input: mouseY -> loading rate (slow at top, fast at bottom) ---
  if (input) {
    // Default to mid range if mouse hasn't moved.
    const ym = (typeof input.mouseY === 'number' && input.mouseY >= 0) ? input.mouseY : H * 0.5;
    const t = Math.max(0, Math.min(1, ym / H));
    // Exponential map: 0.003 .. 0.16 (sim-units per sim-second).
    vLoad = 0.003 * Math.pow(0.16 / 0.003, t);

    if (input.consumeClicks) {
      const clicks = input.consumeClicks();
      if (clicks.length) {
        // Map click x to block index, queue a stress nudge.
        const c = clicks[clicks.length - 1];
        const bx = Math.max(0, Math.min(N - 1, Math.floor((c.x / W) * N)));
        pendingNudge = { idx: bx, amount: 1.2 };
      }
    }
  }

  // --- advance loading plate ---
  const dtClamped = Math.min(0.05, dt);
  timeAcc += dtClamped;
  plate += vLoad * dtClamped;

  // --- inject manual perturbation if pending ---
  if (pendingNudge) {
    const i = pendingNudge.idx;
    // Push the block backward (against the plate) to increase tension on its
    // loading + coupling springs — effectively raise the stress at that site.
    x[i] -= pendingNudge.amount;
    pendingNudge = null;
  }

  // --- relax: cascade of slips until quiescent ---
  const slips = relax();
  if (slips > 0) {
    totalSlips += slips;
    lastQuakeSize = slips;
    lastQuakeTime = timeAcc;
    totalEvents++;
    const m = magnitudeOf(slips);
    hist[binOf(m)]++;
    recentQuakes[recentHead * 2 + 0] = timeAcc;
    recentQuakes[recentHead * 2 + 1] = slips;
    recentHead = (recentHead + 1) % RECENT_LEN;
  }

  // Update displayed stress per block.
  for (let i = 0; i < N; i++) {
    stress[i] = Math.abs(forceOn(i));
  }

  // --- render ---
  ctx.fillStyle = '#06080d';
  ctx.fillRect(0, 0, W, H);

  // Layout: top ~40% chain, bottom ~60% histogram with HUD strip in between.
  const chainTop = 30;
  const chainBot = Math.max(110, H * 0.42);
  const chainH = chainBot - chainTop;
  const stripTop = chainBot + 6;
  const stripH = 16;
  const histTop = stripTop + stripH + 24;
  const histBot = H - 26;

  // --- title strip ---
  ctx.fillStyle = 'rgba(220, 225, 240, 0.92)';
  ctx.font = 'bold 13px sans-serif';
  ctx.textAlign = 'left';
  ctx.fillText('Burridge-Knopoff stick-slip chain', 12, 18);
  ctx.font = '11px sans-serif';
  ctx.fillStyle = 'rgba(180, 190, 215, 0.78)';
  ctx.textAlign = 'right';
  const vTxt = vLoad < 0.01 ? 'slow tectonic' : (vLoad < 0.04 ? 'moderate' : (vLoad < 0.09 ? 'fast' : 'very fast'));
  ctx.fillText(`loading: ${vTxt}  (v = ${vLoad.toFixed(3)})`, W - 12, 18);

  // --- the chain of blocks ---
  // Each block drawn as a square. Color from stress (cool blue -> hot red).
  // x offset shows displacement (visually).
  const slotW = (W - 24) / N;
  const blockSize = Math.min(slotW * 0.78, chainH * 0.55, 26);
  const baselineY = chainTop + chainH * 0.55;

  // Driving plate above blocks.
  ctx.fillStyle = 'rgba(70, 80, 110, 0.5)';
  ctx.fillRect(12, chainTop, W - 24, 6);
  ctx.fillStyle = 'rgba(120, 140, 180, 0.85)';
  ctx.font = '10px sans-serif';
  ctx.textAlign = 'left';
  ctx.fillText('driving plate  →', 16, chainTop - 4);

  // Ground line beneath blocks.
  ctx.strokeStyle = 'rgba(110, 90, 70, 0.7)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(12, baselineY + blockSize * 0.5 + 4);
  ctx.lineTo(W - 12, baselineY + blockSize * 0.5 + 4);
  ctx.stroke();
  // Hatch the ground.
  ctx.strokeStyle = 'rgba(140, 110, 80, 0.45)';
  ctx.lineWidth = 1;
  for (let xH = 14; xH < W - 12; xH += 8) {
    ctx.beginPath();
    ctx.moveTo(xH, baselineY + blockSize * 0.5 + 4);
    ctx.lineTo(xH - 5, baselineY + blockSize * 0.5 + 11);
    ctx.stroke();
  }

  // Springs from plate to each block: zig-zag vertical.
  for (let i = 0; i < N; i++) {
    const cx = 12 + (i + 0.5) * slotW;
    const px = cx + Math.max(-blockSize * 0.6, Math.min(blockSize * 0.6, x[i] * 18));
    ctx.strokeStyle = 'rgba(120, 130, 160, 0.55)';
    ctx.lineWidth = 1;
    ctx.beginPath();
    ctx.moveTo(cx, chainTop + 6);
    const segs = 6;
    for (let s = 1; s < segs; s++) {
      const t = s / segs;
      const sy = chainTop + 6 + t * (baselineY - chainTop - 6 - blockSize * 0.5);
      const sx = (s % 2 === 0 ? cx - 2 : cx + 2);
      // Bend the spring toward the block's actual position as it gets close.
      const w = t * t;
      const fx = sx * (1 - w) + px * w;
      ctx.lineTo(fx, sy);
    }
    ctx.lineTo(px, baselineY - blockSize * 0.5);
    ctx.stroke();
  }

  // Coupling springs between blocks (horizontal).
  for (let i = 0; i < N - 1; i++) {
    const cx1 = 12 + (i + 0.5) * slotW + Math.max(-blockSize * 0.6, Math.min(blockSize * 0.6, x[i] * 18));
    const cx2 = 12 + (i + 1.5) * slotW + Math.max(-blockSize * 0.6, Math.min(blockSize * 0.6, x[i + 1] * 18));
    const y0 = baselineY;
    ctx.strokeStyle = 'rgba(100, 110, 140, 0.55)';
    ctx.lineWidth = 1;
    ctx.beginPath();
    ctx.moveTo(cx1 + blockSize * 0.5, y0);
    const segs = 5;
    for (let s = 1; s < segs; s++) {
      const t = s / segs;
      const sx = cx1 + blockSize * 0.5 + t * (cx2 - cx1 - blockSize);
      const sy = y0 + (s % 2 === 0 ? -3 : 3);
      ctx.lineTo(sx, sy);
    }
    ctx.lineTo(cx2 - blockSize * 0.5, y0);
    ctx.stroke();
  }

  // Blocks themselves.
  for (let i = 0; i < N; i++) {
    const cx = 12 + (i + 0.5) * slotW + Math.max(-blockSize * 0.6, Math.min(blockSize * 0.6, x[i] * 18));
    const sFrac = Math.max(0, Math.min(1, stress[i] / fs[i]));
    // Hue from cool teal (170) at low stress to hot red (5) at threshold.
    const hue = 170 - 165 * sFrac;
    const lum = 38 + 26 * sFrac;
    // Recently slipped: white flash overlay.
    const sinceQuake = timeAcc - lastQuakeTime;
    ctx.fillStyle = `hsl(${hue}, 78%, ${lum}%)`;
    ctx.fillRect(cx - blockSize * 0.5, baselineY - blockSize * 0.5, blockSize, blockSize);
    if (sinceQuake < 0.4 && lastQuakeSize > 0) {
      ctx.fillStyle = `rgba(255, 250, 230, ${(0.4 - sinceQuake) * 1.2})`;
      ctx.fillRect(cx - blockSize * 0.5, baselineY - blockSize * 0.5, blockSize, blockSize);
    }
    // Outline.
    ctx.strokeStyle = 'rgba(10, 12, 18, 0.85)';
    ctx.lineWidth = 1;
    ctx.strokeRect(cx - blockSize * 0.5 + 0.5, baselineY - blockSize * 0.5 + 0.5, blockSize - 1, blockSize - 1);
  }

  // --- HUD strip: stats line ---
  ctx.fillStyle = 'rgba(180, 190, 215, 0.85)';
  ctx.font = '11px sans-serif';
  ctx.textAlign = 'left';
  const lastMag = lastQuakeSize > 0 ? magnitudeOf(lastQuakeSize).toFixed(2) : '—';
  ctx.fillText(`events: ${totalEvents}    last: size ${lastQuakeSize || '—'} (M ${lastMag})`, 12, stripTop + 12);
  ctx.textAlign = 'right';
  ctx.fillText(`sim time: ${timeAcc.toFixed(1)}s`, W - 12, stripTop + 12);

  // --- Gutenberg-Richter histogram ---
  ctx.fillStyle = 'rgba(140, 150, 180, 0.95)';
  ctx.font = 'bold 12px sans-serif';
  ctx.textAlign = 'left';
  ctx.fillText('Gutenberg-Richter:  log₁₀ N  vs  magnitude M', 12, histTop - 6);

  // Plot area frame.
  const padL = 44, padR = 16;
  const px0 = padL, px1 = W - padR;
  const py0 = histTop, py1 = histBot;
  ctx.strokeStyle = 'rgba(110, 120, 150, 0.6)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(px0, py0);
  ctx.lineTo(px0, py1);
  ctx.lineTo(px1, py1);
  ctx.stroke();

  // Determine plot Y range from log10(max count).
  let maxCount = 1;
  for (let k = 0; k < N_BINS; k++) if (hist[k] > maxCount) maxCount = hist[k];
  const yMaxLog = Math.max(1, Math.ceil(Math.log10(maxCount + 1)));

  // Y gridlines at integer log10 levels.
  for (let g = 0; g <= yMaxLog; g++) {
    const y = py1 - (g / yMaxLog) * (py1 - py0);
    ctx.strokeStyle = 'rgba(80, 90, 120, 0.35)';
    ctx.beginPath(); ctx.moveTo(px0, y); ctx.lineTo(px1, y); ctx.stroke();
    ctx.fillStyle = 'rgba(160, 170, 195, 0.8)';
    ctx.font = '10px sans-serif';
    ctx.textAlign = 'right';
    ctx.fillText(`10${g === 0 ? '⁰' : g === 1 ? '¹' : g === 2 ? '²' : g === 3 ? '³' : g === 4 ? '⁴' : '?'}`, px0 - 4, y + 3);
  }

  // X-axis labels.
  ctx.fillStyle = 'rgba(160, 170, 195, 0.85)';
  ctx.font = '10px sans-serif';
  ctx.textAlign = 'center';
  for (let m = 0; m <= 3; m++) {
    const M = m * 0.5;
    if (M < M_MIN || M > M_MAX) continue;
    const xpx = px0 + ((M - M_MIN) / (M_MAX - M_MIN)) * (px1 - px0);
    ctx.fillText(`M ${M.toFixed(1)}`, xpx, py1 + 14);
    ctx.strokeStyle = 'rgba(80, 90, 120, 0.25)';
    ctx.beginPath(); ctx.moveTo(xpx, py0); ctx.lineTo(xpx, py1); ctx.stroke();
  }

  // Bars.
  const barW = (px1 - px0) / N_BINS;
  for (let k = 0; k < N_BINS; k++) {
    const c = hist[k];
    if (c === 0) continue;
    const y = py1 - (Math.log10(c + 1) / yMaxLog) * (py1 - py0);
    const xL = px0 + k * barW + 1;
    const M = M_MIN + (k + 0.5) * (M_MAX - M_MIN) / N_BINS;
    // Color the bar by magnitude (small = teal, big = red).
    const t = Math.max(0, Math.min(1, M / M_MAX));
    const hue = 170 - 165 * t;
    ctx.fillStyle = `hsla(${hue}, 80%, 55%, 0.85)`;
    ctx.fillRect(xL, y, barW - 2, py1 - y);
  }

  // Fitted GR line.
  const fit = fitGR();
  if (fit) {
    ctx.strokeStyle = 'rgba(255, 240, 160, 0.9)';
    ctx.lineWidth = 2;
    ctx.beginPath();
    let first = true;
    for (let s = 0; s < 64; s++) {
      const M = M_MIN + (s / 63) * (M_MAX - M_MIN);
      const yLog = fit.a - fit.b * M;
      if (yLog < 0) continue;
      const xpx = px0 + ((M - M_MIN) / (M_MAX - M_MIN)) * (px1 - px0);
      const ypx = py1 - (yLog / yMaxLog) * (py1 - py0);
      if (ypx < py0 - 4 || ypx > py1 + 4) { first = true; continue; }
      if (first) { ctx.moveTo(xpx, ypx); first = false; }
      else ctx.lineTo(xpx, ypx);
    }
    ctx.stroke();
    // b-value label.
    ctx.fillStyle = 'rgba(255, 240, 160, 0.95)';
    ctx.font = 'bold 11px sans-serif';
    ctx.textAlign = 'right';
    ctx.fillText(`fit:  log₁₀ N = ${fit.a.toFixed(2)} − ${fit.b.toFixed(2)} M`, px1 - 4, py0 + 12);
  } else {
    ctx.fillStyle = 'rgba(200, 210, 230, 0.7)';
    ctx.font = 'italic 11px sans-serif';
    ctx.textAlign = 'right';
    ctx.fillText('accumulating events…', px1 - 4, py0 + 12);
  }

  // Axis title bottom.
  ctx.fillStyle = 'rgba(160, 170, 195, 0.85)';
  ctx.font = '10px sans-serif';
  ctx.textAlign = 'center';
  ctx.fillText('magnitude  M = log₁₀(slip count)', (px0 + px1) / 2, py1 + 24);
}

Stick-Slip Earthquakes: Gutenberg-Richter

Comments (0)