Listening to Earthquakes: Ray Tracing a Whole Planet

An S wave that reaches the core–mantle boundary is finished. The outer core is liquid, liquid has no shear strength, and a shear wave needs something to shear — so 2,891.5 km down, the phase ends. What that leaves on the far side of a large earthquake is a seismic shadow zone the size of a continent, and it is one of the observations that told us the middle of the planet is molten, decades before anyone could model it.

That shadow is on screen in Tremor whenever you select a large event, and how it gets there is the whole design in miniature. No code draws a shadow. No threshold on angular distance exists anywhere in the renderer. We trace rays through a velocity model, ask each one where it is at the current replay time, and write the answers into a preallocated buffer.

Rays that died at the core write NaN. The line renderer breaks its run at a NaN and starts a new one. The gaps are the shadow zones.

The ak135 travel-time model, in 22 nodes

Everything downstream rests on a table. ak135 is a standard one-dimensional reference model of the Earth: for each depth, a compressional and a shear velocity averaged over the whole planet — approximate enough to describe in a paragraph, good enough that seismologists still locate real earthquakes with it. Our version is 22 rows of [depthKm, vp, vs], interpolated piecewise-linearly in radius rather than depth, within about 0.1 km/s of the published model. The source calls it visual grade and means it: a drawing tool, not an ak135 travel time table you would time a phase pick against.

The trick in the table is how discontinuities are encoded: a depth that appears twice is a discontinuity. The builder walks the node list, splits it into continuous segments wherever a depth repeats, and records the jump. Nothing else in the engine knows the boundaries by name — the Moho, the 410, the 660, the core–mantle boundary and the inner-core boundary are one data structure, and adding another would be adding two rows.

The liquid outer core is not a special case in the engine. It is a zero in a table.

That zero does an enormous amount of work. When an S ray reaches the core–mantle boundary, the refraction step looks up the shear velocity on the far side, finds nothing to refract into, and ends the ray — the tracer never learns it has hit a core. Below the inner-core boundary the shear velocity returns, because the inner core is solid; the tracer would propagate there happily if a shear wave could ever arrive. None does, and an assertion checks that none ever has.

Bending 360 rays through a planet

A spherically symmetric Earth hands you a dimension for free: a ray leaving a source never leaves the plane containing the source and the planet's centre. So the tracer is two-dimensional, marching each ray through that great-circle plane and storing it as a radius, an angular distance and an accumulated time — three floats per stored sample, in a flat Float32Array, one sample kept every six steps.

The integration is classical ray theory: with refractive index n = 1/v, the ray equation d/ds(n·d̂) = ∇n says a ray curves toward higher index — toward slower material. In a radially symmetric medium the gradient points along the radius, so the whole bend collapses into a dozen lines of scalar arithmetic, run once per eight-kilometre step.

if (r > 0.5) {
  const dndr = -velocityGradientAt(r, phase) / (v * v);
  const k = (dndr * DS_KM) / r;
  const n = 1 / v;
  const wx = n * dx + k * x;
  const wz = n * dz + k * z;
  const wl = Math.hypot(wx, wz);
  if (wl > 0) {
    dx = wx / wl;
    dz = wz / wl;
  }
}

Straight stepping is fine until a ray meets a discontinuity, where approximation shows immediately: cross the core–mantle boundary a few kilometres late and the refraction angle is wrong and the ray lands somewhere it has no business being. So crossings are solved exactly. stepToRadius() roots the quadratic that puts the ray precisely on the boundary, the step truncates there, and Snell's law is applied in the local tangential/normal basis: sin θ₂ = sin θ₁ · v₂/v₁. If that exceeds one there is no transmitted ray and the direction mirrors — total internal reflection, which is how the core-reflected phases appear without being asked for. Otherwise the refracted direction is rebuilt and nudged 0.02 km clear, so the next step cannot re-detect the same crossing.

Three hundred and sixty rays fan out from the hypocentre between 0.25° and 179.75° of takeoff angle, source depth clamped to 1–700 km. A full fan costs 10–20 ms per phase — on the main thread, no Web Worker, which is defensible only because it happens once per selection rather than once per frame.

The phases we never wrote

Seismology has a rich phase alphabet — P, S, PKP, PKIKP, PcP, ScS — and the tempting way to render it is to special-case each one: draw this arc for a core phase, hide the direct P past that angle. Tremor has none of that. The tracer knows velocity, gradients, boundaries and Snell's law; every phase in the picture follows from those four things.

The S fan stops at the core because the shear velocity there is zero. The P fan does something more interesting: compressional velocity drops from 13.66 to 8.0 km/s crossing into the outer core, so rays refract sharply downward, direct surface arrivals run out past 104°, and core-transiting energy re-emerges as PKP much further around. Both edges of the P shadow come out of the refraction, not out of an if — which is exactly why they need checking.

So the engine ships with an executable self-check, run in development, in the same spirit as the Orrery's ISS assertion: a physics engine nobody can eyeball needs a test that can.

• Travel time increases monotonically along every ray. (A ray that goes backwards in time means the integrator broke.)
• Direct P at 60° arrives near 600 s, within ±60 s; at 90°, near 780 s within ±80 s.
• Direct S at 60° arrives near 1,090 s, within ±110 s.
• The last direct-P surface arrival falls inside 90°–118° — wide enough to tolerate a visual-grade model, narrow enough to catch a broken refraction.
• PKP arrivals exist beyond 130°.
• No S ray is ever found more than 5 km below the core–mantle boundary.

One isochron per frame, zero allocation

A traced ray is a path through space and time; a wavefront is a slice through it. sampleIsochron() binary-searches each ray's samples for the bracket containing the current replay time and interpolates one position out of it — one radius/distance pair per ray, into a preallocated array. Rays that never reached this time, or that died at the core, write (NaN, NaN).

The fronts are drawn as LineSegments rather than a polyline, and that choice is load-bearing. A polyline would connect the last live sample before a gap to the first one after it, welding a chord straight across the shadow. Segments do not: a NaN ends a run, the next valid pair starts a new one, and the shadow carves itself out of the geometry. The glow ribbon is built from independent quads for the same reason — a shared strip would weld the gap the lines just left.

Nothing here is rebuilt per frame: 720 sample slots, 2,880 line vertices and 1,440 ribbon quads are allocated once, and setDrawRange does the work. Each front is drawn twice, mirrored to ±delta, because a wavefront in a spherically symmetric Earth is a surface of revolution about the source axis and the cutaway exposes two half-planes of it. Fronts hold for 120 s of simulated time and fade as they age — P from near-white through ember, S from pale mint through teal — so a busy replay reads as a sequence rather than a blur.

A cutaway Earth built from geometry, not clipping

To watch a wavefront cross the interior you have to cut the planet open, and the obvious way — clip planes — is the wrong tool for the shells: fragment work on every pixel of every layer, and hollow faces at the end of it. Instead the shells are built incomplete. A SphereGeometry with phiStart = π and phiLength = 1.5π is a 270° shell with a missing quadrant, aimed at the default camera. The hole is in the mesh, so it costs nothing.

One layer cannot use that trick. The coastlines and the graticule have to rotate — that is the point of them — so they are the one thing in the scene that is genuinely GPU-clipped, by two planes with clipIntersection enabled. The coastline data is decoded from the module we packed for the Orrery and re-slerped along great circles at about 2° per subdivision, so long chords never sag below the surface and vanish inside the sphere.

The scene then runs on two reference frames, and this is the piece we are quietly proudest of. The cut frame never moves: wedge, strata, boundary hairlines and polar axis are static, so the world-space clip planes are constants and the cut-face art bakes into one indexed mesh and one draw call. The geographic frame — coastlines, graticule, the whole catalog — is slerped so the selected hypocentre rides the +Y polar axis, the edge shared by both exposed faces. Any plane through the source axis is a valid great-circle section of a spherically symmetric Earth, and here two of them happen to be walls you can see. The cross-section is honest whatever event you pick — and instant, with no animation at all, for visitors who have asked for reduced motion.

Four thousand earthquakes in one mesh

The catalog arrives from two USGS summary feeds — the past day at any magnitude, plus M4.5 and above from the past month — through a proxy that slims each event to a positional tuple and owns its own caching. Route-level revalidation would cache a 503 as happily as data, which on the Orrery once froze an offline state in place for six hours. Here the 90-second TTL is anchored to the older of the two feeds' timestamps, so a failed feed keeps being retried while the healthy one keeps serving. Stale beats blank.

On screen, the whole catalog — up to 4,000 events — is one InstancedMesh of octahedra. Each instance sits at radius (6371 − depth) / 1000 scene units: at its true hypocentral depth rather than pinned to the surface, so a subduction zone reads as a slab leaning into the mantle. Scale is logarithmic in magnitude, 0.012 + 0.02 · 1.9^(mag − 4), clamped at both ends — capped so an M8 stays a bead rather than a moon, floored so the micro-quakes still exist. Colour is age: hot orange through ember across the first 24 hours, then ember to slate across a 30-day tail. The past week's story is visible before you click anything.

Earthquake sonification, from a hand-written miniSEED decoder

Seismograms are public. Anyone can pull raw ground motion from the FDSN dataselect services, and the Global Seismographic Network's IU stations have been recording continuously for decades. So the last thing Tremor does with an event is offer to play it to you — real earthquake sonification, not a synthesizer imitating one.

The browser cannot do this itself, for two independent reasons. The data arrives as Steim-compressed miniSEED, a difference-encoded binary format no browser has heard of, and the service sends no CORS headers, so a page could not fetch the bytes even if it could read them. Both halves live on our server: a proxy fetches, a hand-written decoder unpacks, and a normalized Int16 buffer reaches the client.

Choosing a station is a small ranking problem with a physical constraint inside it. Tremor knows 31 IU broadband stations and wants forty minutes of vertical-component record, but nearest is not best: closer than 2° the instrument clips on a large event, past 90° the core starts shadowing P. So the teleseismic band leads the ranking, distance breaks ties, three candidates are tried in order, and a record younger than 45 minutes is refused outright — archives need time, and the client pre-empts that refusal rather than making you discover it.

The miniSEED decoder is about as low-level as web code gets: a DataView, big-endian reads, and a blockette-1000 header saying which encoding it is looking at. Steim packs samples into 64-byte frames of sixteen 32-bit words, the first of which is a nibble map — sixteen two-bit codes describing how the other fifteen are packed. Steim-1 offers one, two or four differences per word; Steim-2 adds sub-codes fitting five 6-bit, six 5-bit or seven 4-bit differences into one word, plus a 30-bit slot for outliers. Words one and two of frame zero are not differences at all but the integration constants X0 and Xn, which makes the stream self-checking. Both encodings are implemented; where record timing misaligns by more than 1.5 samples we count the gap and concatenate anyway, because for sonification continuity beats exactness.

Which brings us to the part with no code in it. Seismic energy lives around 0.05–10 Hz; hearing starts near 20 Hz. The textbook fix is a pitch shifter, and we wrote exactly none of one, because the Web Audio API lets you declare a buffer's sample rate. Forty samples per second, declared as 40 × 150 = 6,000 Hz, plays 150 times faster and 150 times higher — inside the legal buffer-rate window of 3,000 to 768,000 Hz — landing the seismic band at 7.5–1500 Hz. Forty minutes becomes sixteen seconds. It is not a transformation of the data; it is the data, read quickly.

Dynamic range is the harder problem. A teleseismic P onset can be a hundred times the coda behind it, and played linearly the visitor hears one tick and then silence. So the client compands with an exponent near 0.6 — in the browser, not the server, because the server's job is honest counts and the shaping belongs where the visitor can hear it applied. Before transport the record is demeaned (SNZO, for one, sits about 1,500 counts off zero) and peak-normalized. The sparkline plots peak per bin rather than RMS, for the same reason: the arrival is a spike, and averaging is the one operation that would hide it.

The player itself is deliberately small: one GainNode, 30 ms ramps scheduled one ahead, an audio context created inside the click handler so autoplay policy has nothing to object to, playback that stops when the tab goes away. Mute is not Tremor's own — it reads the sitewide synth mute shared with every game in the Lab. Mute a game, and the planet goes quiet too.

A chip in the corner of the screen says visual-grade travel times, for the same reason the Orrery's says visual-grade propagation. A 1-D model with no station corrections and no three-dimensional structure puts a wavefront where the eye expects it, not where an instrument would record it. Faithful to the eye, not to the meter; never use it for anything but wonder. Saying that plainly costs nothing and buys the rest of the picture its credibility.

The honest gaps: the ray tracing runs on the main thread, and the page has no automated end-to-end coverage yet — only the development-time assertions, which check the physics and not the buttons. Both are on the list.