methods

1. Question and scope

After a large earthquake, the rate of aftershocks decays roughly as a power of time and their locations cluster around the rupture. Statistical models of this behaviour — the Omori–Utsu law for the rate, the Gutenberg–Richter law for magnitudes, and the Epidemic-Type Aftershock Sequence (ETAS) model that lets every aftershock trigger its own — are the basis of operational aftershock forecasting. This study asks one question:

Can modern machine learning improve the characterisation of how an aftershock sequence evolves, relative to established statistical baselines?

"Characterisation" is meant literally. QuakeCast forecasts the number and the spatial distribution of aftershocks above a completeness magnitude after a mainshock has already happened, over horizons of one to thirty days. It does not predict the time, place or size of any future mainshock, it is not an operational forecasting system, and nothing in it should be read as one.

Five sub-questions structure the work: (Q1) does ML improve aftershock-rate forecasts, (Q2) does it improve spatial density forecasts, (Q3) how does it compare with the classical laws and models, (Q4) how stable are the results across sequences, and (Q5) does generalisation fail across tectonic regimes?

2. Data

2.1 Catalogs

All events come from authoritative public catalogs through the FDSN event web-service standard; no event, magnitude, location or time is synthesised.

For every event the pipeline keeps the provider event identifier, origin time, latitude, longitude, depth, magnitude, magnitude type, contributing network, review status where the provider reports one, and the catalog of origin. Raw responses are cached compressed and hashed by query so every run is reproducible offline.

Conventions validated in code. ComCat GeoJSON times are epoch milliseconds in UTC and coordinates are ordered [longitude, latitude, depth]; FDSN text times carry no zone designator and are UTC by specification. Validation rejects any event outside the valid latitude and longitude ranges (the signature of a coordinate swap), rejects timezone-naive times, drops events without a magnitude or location, and removes duplicate identifiers. Depths may be negative (above the datum) and are kept as reported. Magnitude types are heterogeneous — SCSN reports mostly coda (mc) and local (ml) magnitudes for small events — and are retained rather than silently homogenised; §7 discusses the consequence.

Why regional catalogs outside the US. An initial survey drew fifteen candidate sequences from ComCat alone. Outside the United States ComCat holds only the largest events: the 2016 Kumamoto sequence has 73 events inside its aftershock zone, the 2016 Norcia sequence 54. Such catalogs cannot support any rate model, and a cross-regime comparison built on them would measure catalog quality rather than tectonics. Italian and New Zealand sequences are therefore taken from INGV and GeoNet. No public FDSN access to the Japanese (JMA), Taiwanese or Chilean national catalogs was used, so those regions are out of scope; the excluded sequences and the reason for each are recorded in sequences.EXCLUDED.

2.2 Sequences

Fourteen mainshocks were resolved against their providers' live services; verify_registry re-checks every identifier, origin time, location and magnitude against the catalog before a run. Magnitudes are each agency's preferred value, so the same earthquake can differ between agencies (Darfield 2010: ComCat Mw 7.0, GeoNet Mw 7.2).

Table 1. The benchmark sequences. Mc from the protocol of §3; events are above Mc, inside the zone and inside the window. Generated from the artifacts by scripts/paper_tables.py.

2.3 Sequence definition

Sequence membership is a modelling choice, so each part of it is explicit and varied in §6:

• Zone. A circle about the mainshock epicentre whose radius is the Gardner & Knopoff (1974) distance window for the mainshock magnitude (77 km for Landers, Mw 7.3). The Wells & Coppersmith (1994) subsurface rupture length is the alternative rule.
• Window. From t0 = 0.5 days to 365 days after the mainshock.
• Truncation rule (pre-registered). The window closes just before the earliest of (a) a later in-zone event with M ≥ M_main − 0.5, which starts a new sequence rather than continuing this one, and (b) the origin of another registered mainshock whose zone intersects this one, so that no two benchmark sequences share target events.
• Trigger history. Events up to 30 days before the mainshock (foreshocks) are kept as ETAS triggers but never as targets. Ridgecrest 2019, whose Mw 6.4 foreshock came 34 hours before the Mw 7.1, is why.
• Local geometry. Spatial models work in kilometres on an azimuthal-equidistant projection centred on the epicentre, which preserves distance from the centre exactly (tested against the haversine distance to machine precision). Kernels in degrees would be silently anisotropic on the ground.

3. Catalog completeness

Aftershock models are only meaningful above the magnitude at which the catalog records essentially every event (Mc). Immediately after a large mainshock, small events are hidden in the coda of larger ones, so Mc is elevated for hours to days — short-term aftershock incompleteness (Kagan 2004; Helmstetter et al. 2006; Page et al. 2016). In the Landers catalog, maximum-curvature Mc falls from 2.7 in the first days to 1.8 after about fifty days, and the b-value estimated at M ≥ 2.8 in the interval 0.2–0.5 days is 0.70 — an incompleteness artefact, against 1.02 at M ≥ 3.2.

Completeness estimators disagree. On synthetic catalogs with a gradual detection ramp (tests/test_completeness.py), maximum curvature (Wiemer & Wyss 2000) and the 90 % goodness-of-fit level underestimate the true threshold, while 95 % goodness-of-fit and b-value stability (Cao & Gao 2002) land near it. A single goodness-of-fit pass can also succeed by chance on a sparse upper tail: the Parkfield catalog, complete near M 1.0–1.2, yields a spurious 2.3 in one time bin.

Mc protocol (version 3, fixed before any model comparison). (1) Split (t0, t_end] into eight log-spaced time bins, merging any bin with fewer than 100 events into the next. (2) In each bin, the per-sample estimate is the goodness-of-fit Mc at R ≥ 95 %, else at R ≥ 90 %, else b-value stability; the bin's Mc is the median of that estimate over 200 bootstrap resamples (Woessner & Wiemer 2005 recommend bootstrapping Mc). (3) The sequence Mc is the largest bin Mc, rounded up to 0.1, so the catalog is complete at that level throughout the window — including its early, worst-recorded part. Mc is a property of the catalog, fixed once per sequence and shared by every model; because estimators still disagree by about 0.5 early in a sequence, raising Mc is a primary robustness axis (§6) rather than something tuned.

b-values use the Aki (1965) maximum-likelihood estimator with the Shi & Bolt (1982) uncertainty.

4. Models

Every model forecasts the number of above-Mc aftershocks in a window (T, T + h] from events with time ≤ T only.

4.1 Classical baselines

• Poisson (mean rate) — a constant rate over (t0, T]. Deliberately naive; reported so that improvement over "no decay at all" can be told apart from improvement over the physical models.

• Persistence — the rate over (T/2, T] carried forward. A stronger naive reference.

• Omori–Utsu — λ(t) = K/(t + c)^p fitted by maximum likelihood on (t0, T], with K profiled out analytically and a multi-start search over (c, p). Recovers known parameters from simulated catalogs to within 1 % in p (tests/test_omori.py).

• ETAS — temporal ETAS (Ogata 1988): λ(t) = μ + Σ_{t_i<t} K e^{α(m_i − Mc)} (t − t_i + c)^{−p}. Maximum likelihood with an analytic gradient (checked against finite differences to 10⁻¹⁰). Forecasts come from Monte Carlo simulation of the branching process (500 catalogs), with magnitudes drawn from the fitted Gutenberg–Richter law truncated at the mainshock magnitude, because future aftershocks trigger further aftershocks and no closed form exists. With the true parameters the forecaster is calibrated on simulated catalogs (probability-integral-transform χ² p = 0.88; 90 % interval coverage 0.90). Parameter recovery on simulated catalogs of about 900 events is within 1–9 % for all five parameters.

  Two ETAS variants are reported. ETAS (primary) fixes μ to the pre-mainshock background rate — the median yearly count above Mc in the zone over the ten years before the mainshock, which uses past data only — and constrains the branching ratio over 365 days to at most 0.99, the condition under which the process is non-explosive. ETAS (free) estimates all five parameters without constraint. Both decisions were made from diagnostics, before any model comparison: early in a sequence the free fit can trade an enormous background rate against a runaway p (at T = 3 days for Landers it returned μ = 68 events/day with p at its upper bound), and a supercritical fit produces forecasts that grow without bound.

4.2 Machine-learned models

All learned models use the same 21 history features computed from events with time ≤ T (ml/features.py, which raises if it is handed a later event): counts and rates in log-spaced trailing windows, an empirical local decay slope, a fit-free 1/t extrapolation, mainshock and largest-aftershock magnitude relative to Mc, time since the largest aftershock, the b-value so far, recent mean magnitude, and spatial spread. Targets are counts; losses are Poisson.

• Poisson GLM — L2-regularised Poisson regression.
• Gradient boosting (GBM) — LightGBM with a Poisson objective.
• GBM-hybrid — the same booster trained with the Omori forecast as an offset, so it learns a multiplicative correction to the physical model rather than the whole rate.
• Poisson MLP — two tanh layers (32 and 16 units), Adam, weight decay, early stopping; the neural representative. Recurrent, convolutional, transformer and neural point-process models carry far more parameters than a pool of a few thousand forecast windows from fourteen sequences can identify, so they were not fitted.

Two training regimes answer different questions:

• Within-sequence. At origin T the model trains only on examples from the same sequence whose target window closed by T — (o′, o′ + h] with o′ + h ≤ T on a fixed grid of 60 log-spaced sub-origins. Early in a sequence there are too few such examples (fewer than 24), and the model abstains.
• Cross-sequence. The model trains on other sequences and forecasts an unseen one. Three training pools: pooled (every other sequence), chronological (only sequences whose mainshock came first — no calendar-future information), and other region (only sequences from a different region group). Pooled + within adds the target's own closed windows (transfer with local adaptation).

5. Evaluation design

Forward-in-time only. Forecasts are issued at 15 log-spaced origins (1–240 days) for horizons of 1, 7 and 30 days, every model refitted at every origin on data up to that origin. No split is random. Leakage guards are enforced in code and tested (tests/test_leakage.py): feature builders reject future events, training windows must close before the origin, and the features at an origin are provably invariant to anything after it.

Scores. Each model is scored against the full predictive distribution it implies. Deterministic-rate models (Poisson, persistence, Omori, the learned models) imply a Poisson count; ETAS forecasts are over-dispersed and are scored with a negative binomial matched to the simulated mean and variance. Scoring ETAS with a Poisson would penalise it for honestly reporting the spread of a cascading process; scoring Omori with a negative binomial it never claimed would flatter it. Reported: the log predictive score (strictly proper; Gneiting & Raftery 2007) and the per-event information gain over Omori; MAE and RMSE on log₁₀(count + 1); 90 % interval coverage; the randomised probability integral transform for calibration (Czado, Gneiting & Held 2009); and the CSEP number test (Schorlemmer et al. 2007; Zechar et al. 2010) per window. Comparisons use only windows where every compared model issued a forecast. Uncertainty on information gain within a sequence comes from a moving-block bootstrap over windows (forecast windows are autocorrelated); across sequences, from a cluster bootstrap over sequences, the unit across which generalisation is asked.

Spatial. Every spatial model outputs a density over the same 2 km grid covering the zone, normalised to one, with an identical 0.1 % uniform floor mixed into every model so one unexpected epicentre cannot score −∞. Models: uniform; Gaussian kernel density with a forward-cross-validated bandwidth; nearest-neighbour adaptive kernel density; an ETAS-style kernel (Omori-time-weighted power-law kernels about past events, temporal parameters from that origin's ETAS fit, kernel shape by forward cross-validation); and a cell-level LightGBM. All tuning uses pseudo- forecasts that closed before the origin. The score is the mean log density at the observed epicentres.

6. Robustness design

Each variant changes one choice: Mc raised by 0.3 or 0.6 above the protocol; t0 of 1 or 2 days; radius scaled by 0.5, 1.5 or 2; the Wells–Coppersmith radius rule; a 180-day window; truncation off; foreshock lookback off. Variants run at lighter settings (eight origins, 200 ETAS simulations) and are compared with the primary design at the same light settings, so settings never confound a design effect.

6a. Changes made after the first Landers run (disclosed)

The design in §§3–6 was fixed before the multi-sequence runs, but the first full Landers run was used as a validation pass, and it changed five things. Each change was justified by a defect found in that run, not by a better score, and is listed here so a reader can discount accordingly.

1. Spatial tuning horizon. Forward cross-validation needs a pseudo- forecast whose window closed before the origin; for 30-day windows none exists until day 30.5, so every early spatial model silently fell back to fixed defaults. Tuning now uses the longest shorter horizon (7, 3 or 1 day) with closed windows.
2. ETAS-kernel shape grid. With γ ≤ 1 the largest mainshock kernel was ≈ 8 km wide, far below the scale of a M7 rupture; γ = 1.5 was added.
3. ETAS optimiser. The subcritical fit moved from SLSQP to a reparametrisation in which the branching ratio is a box constraint (K solved from n = 0.99·sigmoid(z)) with an analytic gradient. On Landers this reached the same log-likelihood at every origin (to 0.01), which established that the degenerate early fits (α and p on their bounds) are genuine maxima of a flat likelihood, not optimiser failures.
4. Dispersion twins. A "+ NB spread" twin was added for every Poisson-predictive model (§5), after Landers showed Omori's Poisson intervals covering 42 % of outcomes at the 90 % level while ETAS won on log score despite worse point accuracy. The twin separates skill in the mean from honesty of the spread.
5. Snapshot maps. A variable-shadowing bug meant no spatial density maps were saved in the batch run (scores were unaffected). The maps were recomputed with the same procedure and are marked as recomputed; a regression test now covers it.