Skip to content

Models Classical

Felipe Santibañez-Leal edited this page Jun 17, 2026 · 2 revisions

Models — Classical Theories

The analytical, physics-informed and statistical models that define the field-standard baseline for short-term probabilistic seismic forecasting. Each model is given here with its governing equation(s), its parameters, how it is estimated, its assumptions and failure modes, and a real peer-reviewed reference (with DOI). These are not legacy curiosities: a well-tuned ETAS remains the de-facto operational baseline that any candidate model — including every neural model — must beat in prospective CSEP testing before it can claim forecasting skill.

Honest framing. Every model on this page produces, at best, a conditional rate or probability, never a deterministic prediction. A prediction is a deterministic statement that an event will or will not occur; a forecast gives a probability strictly in $(0,1)$ (Jordan et al., 2011). Where a method is contested, it is flagged explicitly.

Conventions used throughout. Magnitudes are homogenized to moment magnitude $M_w$ where possible; $M_c$ is the magnitude of completeness; $b$ is the Gutenberg–Richter slope; $\beta = b\ln 10$. The history available at time $t$ is written $\mathcal{H}_t = {(t_i, x_i, y_i, m_i): t_i < t}$.


Table of contents

  1. How these models fit together
  2. Gutenberg–Richter — the magnitude–frequency law
  3. Omori–Utsu — aftershock decay
  4. ETAS — Epidemic-Type Aftershock Sequence
  5. Reasenberg–Jones — operational aftershock probabilities
  6. STEP — Short-Term Earthquake Probability
  7. EEPAS — Every Earthquake a Precursor According to Scale
  8. Smoothed seismicity — the spatial background and mandatory null
  9. BPT / renewal — long-term recurrence
  10. Rate-and-state friction + Coulomb stress — the mechanistic layer
  11. Contested / precursory methods — honest assessment
  12. From a conditional intensity to the published probability
  13. References

0. How these models fit together

The product is a conditional estimator: given recent observations $\mathcal{H}_t$, it outputs a bounded probability of one-or-more target events in a region over horizons of ~1 day, ~2 days, ~1 week. The classical models below play three distinct roles — they are the forecasting baselines an enhanced model must beat, the building blocks of the conditional intensity, and the calibration scaffold that keeps the output honest.

Layer Models Forecast role
Magnitude–frequency Gutenberg–Richter, $M_c$, $b$-value Converts a rate of events $\ge M_{\text{ref}}$ into a rate $\ge$ any $M$; the magnitude term of every forecast.
Aftershock decay Omori–Utsu, Reasenberg–Jones The dominant short-horizon (1–7 day) signal; the baseline for "tomorrow's earthquakes."
Self-exciting clustering ETAS (temporal + spatial) State-of-the-art physics-free short-term baseline; the model to beat.
Long-term smoothed rate Helmstetter et al. smoothed seismicity The time-independent spatial background $\mu(x,y)$ that ETAS needs; also the mandatory null.
Medium-term precursory EEPAS Months-to-years scale; outside the 1-week window, useful as a feature/context source.
Operational hybrid STEP Production reference for the "next-day" probabilistic map output shape.
Renewal / recurrence BPT, characteristic EQ Long-term, fault-specific; conditions the background over years.
Physics-based stress Rate-and-state (Dieterich), Coulomb $\Delta\mathrm{CFS}$ Mechanistic priors on where triggering is enhanced or suppressed.
Contested precursors AMR, LURR, RTL, PI Research features at best — never standalone alarms.

A defensible production baseline is ETAS (or Reasenberg–Jones) for the time-dependent term × smoothed-seismicity for the spatial background × Gutenberg–Richter for the magnitude term, all scored against CSEP consistency tests. Everything else is an enhancement or a feature source.

The unifying mathematical object is the conditional intensity function $\lambda(t, x, y \mid \mathcal{H}_t)$ of a marked space–time point process — the instantaneous expected rate of events at time $t$, location $(x,y)$, magnitude $m$, given the history. Every model below is a special case. The model is fit by maximizing the point-process log-likelihood

$$\ln L = \sum_i \ln \lambda(t_i, x_i, y_i \mid \mathcal{H}_{t_i})

  • \int_0^T!!\int_A \lambda(t, x, y \mid \mathcal{H}_t), dx, dy, dt,$$

where the integral ("compensator" / survival) term is what makes the model probabilistic and calibratable rather than a regressor. See Models — Analytical / ML for the point-process framework in full.


1. Gutenberg–Richter — the magnitude–frequency law

The frequency–magnitude distribution of earthquakes is, to first order, a power law. It supplies the magnitude term of every forecast: it converts a forecast rate of events above one reference magnitude into a rate above any target magnitude.

1.1 Governing equation

$$\log_{10} N(\ge M) = a - b,M, \qquad M \ge M_c$$

  • $N(\ge M)$ — number (or rate) of earthquakes with magnitude $\ge M$.
  • $a$ — productivity / total seismicity of the space–time–volume window (depends on catalog length and area).
  • $b$ — the b-value, the slope of the log-frequency vs. magnitude line; globally $b \approx 1.0$.

Equivalently, above completeness $M_c$ the magnitudes follow a (shifted) truncated-exponential density:

$$f(M) = \beta, e^{-\beta (M - M_c)}, \qquad \beta = b \ln 10, \quad M \ge M_c.$$

1.2 b-value estimation — the Aki–Utsu MLE

Aki (1965) showed the maximum-likelihood estimator (for a homogeneous Poisson process with a continuous exponential magnitude distribution) is

$$\hat b = \frac{\log_{10} e}{\bar M - M_c} \approx \frac{0.4343}{\bar M - M_c},$$

where $\bar M$ is the mean magnitude of events with $M \ge M_c$. With magnitudes reported at finite bin resolution $\Delta M$, the Utsu / Tinti–Mulargia binning correction applies:

$$\hat b = \frac{\log_{10} e}{\bar M - \left(M_c - \tfrac{\Delta M}{2}\right)}.$$

This is the Aki–Utsu estimator. A common rule of thumb requires ~50–100 events above $M_c$ for a stable estimate; its standard error is given by Shi & Bolt (1982). The b-value must never be hard-coded to 1 — it is re-estimated on a rolling space–time window and its uncertainty is propagated into the forecast tail.

1.3 Magnitude of completeness $M_c$

$M_c$ is the lowest magnitude at which (almost) all events in a space–time volume are detected. The Aki–Utsu $\hat b$ is strongly biased if $M_c$ is wrong — the single most-cited weakness of the estimator. Common estimators:

  • Maximum-curvature (MAXC): $M_c$ = magnitude bin at the peak of the non-cumulative frequency–magnitude distribution, often with a +0.2–0.3 correction. (The +0.2 correction was calibrated for California and is not established as universal — it must be re-validated per region.)
  • Goodness-of-fit (Wiemer & Wyss, 2000): the smallest $M_c$ for which the GR model fits the observed distribution to a chosen percentage (e.g. 90/95%).
  • b-value stability (Cao & Gao, 2002; Woessner & Wiemer, 2005): $M_c$ where $\hat b$ stabilizes.
  • EMR, MBS, Lilliefors variants as cross-checks.

1.4 Assumptions and failure modes

  • Self-similar (exponential) magnitudes above $M_c$; this breaks down near the maximum magnitude $M_{\max}$ (a taper / corner-magnitude form is then needed: tapered GR / Kagan distribution).
  • Stationary, complete catalog. The b-value can vary with stress, depth and fault maturity (a real signal) or with $M_c$ mis-estimation (an artifact) — the two must be distinguished carefully.
  • Forecasting role: indispensable, the magnitude term of every forecast. A daily pipeline re-estimates $M_c$ and $b$ on a rolling window and carries their uncertainty forward.

References. Gutenberg & Richter (1944), BSSA 34, 185–188; Aki (1965), Bull. Earthq. Res. Inst. 43, 237–239; Tinti & Mulargia (1987), BSSA 77(6), 2125–2134; Wiemer & Wyss (2000), BSSA 90(4), 859–869, doi:10.1785/0119990114; Woessner & Wiemer (2005), BSSA 95(2), 684–698, doi:10.1785/0120040007.


2. Omori–Utsu — aftershock decay

After a mainshock, the aftershock rate decays as a power law in time. This is the single most reliable short-horizon (hours-to-weeks) signal, and the temporal kernel that both ETAS and Reasenberg–Jones embed.

2.1 Modified Omori (Omori–Utsu) law

$$n(t) = \frac{K}{(t + c)^{,p}}$$

  • $n(t)$ — rate of aftershocks per unit time at elapsed time $t$ after the mainshock.
  • $K$ — productivity constant (scales the total number).
  • $c$ — time offset (small, hours-to-day); a contested quantity strongly affected by catalog incompleteness immediately after the mainshock.
  • $p$ — decay exponent, typically $p \in [0.9, 1.5]$, varying sequence to sequence.

The original Omori (1894) law had $p = 1$; Utsu (1957/1961) generalized it to arbitrary $p$. The expected cumulative number between $t_1$ and $t_2$ (for $p \ne 1$) is

$$N(t_1, t_2) = \int_{t_1}^{t_2} \frac{K}{(t+c)^p},dt = \frac{K}{1-p}\Big[(t_2 + c)^{1-p} - (t_1 + c)^{1-p}\Big],$$

and for $p = 1$, $N = K\big[\ln(t_2+c) - \ln(t_1+c)\big]$.

2.2 Estimation

Maximum-likelihood on the point process (Ogata, 1983) jointly fits $(K, c, p)$. Beware early-time incompleteness, which biases $c$ upward and $p$ downward unless the fit start time is delayed past the time of completeness.

2.3 Honest limit — short-term aftershock incompleteness

Immediately after a large mainshock — exactly the highest-stakes, highest-traffic moment for a forecast — the catalog is grossly incomplete: $M_c$ spikes for hours to days while small events are buried in the coda. A naive fit then underestimates productivity precisely when a large aftershock is most likely. The real-time update therefore uses a time-dependent completeness $M_c(t)$ (an incompleteness-aware, detrended likelihood) rather than a flat threshold. This is a decided method, not an open question.

References. Utsu, Ogata & Matsu'ura (1995), J. Phys. Earth 43, 1–33, doi:10.4294/jpe1952.43.1; Ogata (1983), J. Phys. Earth 31, 115–124.


3. ETAS — Epidemic-Type Aftershock Sequence

ETAS (Ogata, 1988; spatio-temporal extension Ogata, 1998) is a self-exciting (Hawkes) point process: every event can trigger its own offspring. It stitches the background rate, the Utsu productivity, the Omori–Utsu time kernel, and a spatial kernel into a single conditional intensity. It is the gold-standard physics-free short-term baseline — the model any enhanced forecaster must beat.

3.1 Space–time conditional intensity

$$\lambda(t, x, y \mid \mathcal{H}_t) = \mu(x,y) + \sum_{i:,t_i < t} k(m_i), g(t - t_i), f!\left(x - x_i,, y - y_i \mid m_i\right)$$

with the canonical Ogata-1998 separable kernels:

Background term — stationary (tectonic) rate, often from smoothed historical seismicity:

$$\mu(x,y) = \mu_0, u(x,y),$$

where $u(x,y)$ is a spatial pdf of background events and $\mu_0$ the total background rate.

Utsu productivity (magnitude):

$$k(m) = K, e^{,\alpha (m - M_0)},$$

with $K$ the aftershock occurrence rate at the reference magnitude $M_0$ and $\alpha$ the productivity parameter (how strongly a larger event seeds more aftershocks).

Omori–Utsu temporal kernel:

$$g(t) = \frac{p-1}{c}\left(1 + \frac{t}{c}\right)^{-p} \quad\Big(\propto (t + c)^{-p}\Big).$$

Spatial kernel (Ogata 1998, inverse-power form):

$$f(x - x_i,, y - y_i \mid m_i) = \frac{q-1}{\pi,\zeta^2} \left(1 + \frac{r^2}{\zeta^2}\right)^{-q}, \qquad r^2 = (x - x_i)^2 + (y - y_i)^2, \qquad \zeta = D, e^{,\gamma (m_i - M_0)},$$

with $D$ the characteristic triggering length at $M_0$, $\gamma$ the magnitude-scaling of the aftershock zone, and $q$ controlling the spatial tail.

The magnitude distribution is independent of history (separability) and follows Gutenberg–Richter: $\log_{10} N(\ge m) = a - b,m$, $f(m) = \beta e^{-\beta(m - M_c)}$.

A common temporal-only ETAS (Ogata 1988) collapses this to $$\lambda(t) = \mu + \sum_{i:,t_i < t} \frac{K, e^{\alpha(m_i - M_0)}}{(t - t_i + c)^{p}}.$$

3.2 Branching ratio and the two stability gates

The branching ratio $n$ is the average number of directly-triggered offspring per event, integrated over time and magnitude:

$$n = \int_{M_0}^{M_{\max}}!!\int_0^{\infty} k(m),g(\tau),\beta e^{-\beta(m - M_0)}, d\tau, dm.$$

Two logically distinct conditions must both hold:

  1. Finite branching — the productivity × magnitude integral converges only if $\alpha < \beta$ (with $\beta = b\ln 10$). If $\alpha \ge \beta$ the largest events dominate and the model is improper. This is a real fitting gotcha.
  2. Subcriticality / stationaritygiven $\alpha < \beta$, the branching ratio must satisfy $n < 1$ (subcritical, the sequence dies out and is integrable). $n = 1$ is critical; $n > 1$ is supercritical (explosive, non-physical for a real catalog) and signals a mis-fit. Any fit with $n \ge 1$ is rejected.

3.3 Estimation

  • MLE of the point-process log-likelihood $\ln L = \sum_i \ln\lambda(t_i) - \int \lambda, dt, dx, dy$ (Ogata, 1988).
  • EM / stochastic declustering (Zhuang, Ogata & Vere-Jones, 2002) assigns each event a probability of being background vs. triggered, recovering $u(x,y)$ and decoupling the background from the clustering.
  • Bayesian / simulation-based inference (e.g. bayesianETAS, INLAbru) for scalable, likelihood-aware fitting with a posterior — not just a point estimate.
  • The fit is performed on the full, un-declustered catalog, because triggering is the predictable signal.

3.4 Assumptions and honest limits

  • Triggering is isotropic and magnitude-independent in shape, whereas real aftershock zones are elongated along the rupture — a known simplification. Anisotropic / finite-fault ETAS variants exist and matter for great subduction earthquakes.
  • The background $\mu$ is stationary in time, which is violated near swarms, slow-slip and induced seismicity.
  • Forecasting role: the floor any candidate must clear. The genuine openings for an enhanced model are exactly ETAS's gaps — fixed functional kernels (cannot learn anisotropy unless hand-coded), single-catalog (hard to fuse geodesy / InSAR / multiple networks / sub-$M_c$ events), and the separability assumption. See Models — Analytical / ML.

References. Ogata (1988), JASA 83(401), 9–27, doi:10.1080/01621459.1988.10478560; Ogata (1998), Ann. Inst. Statist. Math. 50(2), 379–402, doi:10.1023/A:1003403601725; Zhuang, Ogata & Vere-Jones (2002), JASA 97(458), 369–380, doi:10.1198/016214502760046925.


4. Reasenberg–Jones — operational aftershock probabilities

The operational ancestor of USGS aftershock forecasting (now OAF) and the most transparent "tomorrow's earthquakes" model. It is the Gutenberg–Richter magnitude term multiplied by the modified-Omori time decay.

4.1 Rate model

$$\lambda(t, M) = \frac{10^{,a + b,(M_m - M)}}{(t + c)^{,p}}$$

  • $\lambda(t, M)$ — rate of aftershocks $\ge M$ at elapsed time $t$ after a mainshock of magnitude $M_m$.
  • $a$ — regional aftershock productivity.
  • $b$ — Gutenberg–Richter slope.
  • $c, p$ — Omori–Utsu parameters.

4.2 Expected number and probability

The expected number of $\ge M$ events in $[t_1, t_2]$ is

$$N(M; t_1, t_2) = \int_{t_1}^{t_2} \lambda(t, M), dt = 10^{,a + b(M_m - M)} \int_{t_1}^{t_2} (t+c)^{-p}, dt,$$

and, assuming a non-homogeneous Poisson process, the probability of one or more $\ge M$ aftershocks is

$$P(M; t_1, t_2) = 1 - e^{-N(M; t_1, t_2)}.$$

4.3 Three parameter regimes

  1. Generic — fixed $(a, b, c, p)$ from past sequences in a similar tectonic setting. (The classic Reasenberg–Jones 1989 California generic values do not transfer to Chilean subduction — do not reuse them.)
  2. Sequence-specific — refit to the ongoing sequence.
  3. Bayesian — a generic prior updated by the sequence (the modern OAF default; Page et al., 2016 extended the scheme to a full global, tectonic-regime-aware form).

4.4 Role and limits

  • No spatial term — it gives "how many", not "where". STEP and ETAS supply the spatial component.
  • An excellent, transparent fallback and sanity-check alongside ETAS. The USGS OAF system runs both Reasenberg–Jones and ETAS.

References. Reasenberg & Jones (1989), Science 243(4895), 1173–1176, doi:10.1126/science.243.4895.1173; Reasenberg & Jones (1994), Science 265, 1251–1252, doi:10.1126/science.265.5176.1251; Page et al. (2016), BSSA 106(5), 2290–2301, doi:10.1785/0120160073.


5. STEP — Short-Term Earthquake Probability

STEP is the production reference for the product's output shape: a near-real-time system producing gridded short-interval probabilistic maps — i.e. exactly the "one-inference-per-short- interval" probabilistic regional map this product emits daily.

5.1 What it is

A near-real-time California system (USGS + SCEC + ETH) producing hourly maps of the probability of strong shaking (Modified Mercalli $\ge$ VI) on a ~10-km grid.

5.2 Structure

STEP combines a time-independent background rate with a clustering (aftershock) component built on the Reasenberg–Jones rate model. It blends generic, sequence-specific, and spatially-varying Reasenberg–Jones parameter sets — picking the most informative available for each grid cell — then maps the resulting short-term rate to a ground-motion probability via a GMPE. The conditional rate at a cell is the background rate plus the summed Reasenberg–Jones contribution of recent events.

5.3 Role and limits

The canonical production template for "tomorrow's earthquakes" probabilistic maps, exactly the daily gridded output shape of this product. Limits: California-tuned generic parameters; it is aftershock-dominated by design (weak for mainshock/foreshock anticipation).

Reference. Gerstenberger, Wiemer, Jones & Reasenberg (2005), Nature 435, 328–331, doi:10.1038/nature03622.


6. EEPAS — Every Earthquake a Precursor According to Scale

A space–time–magnitude point-process model in which every earthquake is treated as a long-term precursor, according to scale, of larger events to follow. It is built on the empirical precursory scale increase ($\Psi$) phenomenon. It is a medium-term (months-to-years) model — outside the 1-week primary window — used here as a feature/context source, not a short-term core model.

6.1 Predictive scaling relations

For a precursor of magnitude $M_p$, the expected mainshock magnitude $M_m$, precursor time $T_P$, and precursory area $A$ scale as

$$M_m = a_M + b_M, M_p \quad (b_M \text{ fixed to } 1)$$ $$\log_{10} T_P = a_T + b_T, M_p$$ $$\log_{10} A = a_A + b_A, M_p$$

with associated spreads $\sigma_M, \sigma_T, \sigma_A$.

6.2 Rate density

The forecast rate density is a background term plus a weighted sum of contributions from prior earthquakes, each a product of three densities — normal in magnitude, lognormal in elapsed time, bivariate-normal in location:

$$\lambda(t, m, x, y) = \mu, \lambda_0(m, x, y)

  • \sum_{i:,t_i < t} w_i; f(m \mid M_i); g(t - t_i \mid M_i); h(x, y \mid x_i, y_i, M_i),$$
  • $f$ — normal density, mean $a_M + b_M M_i$, sd $\sigma_M$.
  • $g$ — lognormal density in elapsed time, median $10^{,a_T + b_T M_i}$, log-sd $\sigma_T$.
  • $h$ — bivariate-normal in location, with variance scaling as $A \propto 10^{,a_A + b_A M_i}$.
  • $w_i$ — weighting that down-weights events likely to be aftershocks; $\mu$ mixes EEPAS with a baseline (PPE) model.

Honest caveat. The published EEPAS density constants contain known typos across papers, acknowledged by the authors. If EEPAS is used, pin the constants to a reference implementation (pyCSEP / floatCSEP) rather than transcribing them; the structure above is the load-bearing part. EEPAS has genuine, CSEP-tested forecasting skill at medium term in New Zealand, Japan, California and Italy.

References. Rhoades & Evison (2004), Pure Appl. Geophys. 161, 47–72, doi:10.1007/s00024-003-2434-9; Rhoades (2007), SRL 78, 110–115.


7. Smoothed seismicity — the spatial background and mandatory null

A stationary, time-independent estimate of where earthquakes occur, obtained by smoothing a declustered catalog with an adaptive kernel. Large earthquakes nucleate where small ones cluster, so the smoothed density of past small events predicts the spatial distribution of future large ones. It plays a dual role: it is the spatial background $\mu(x,y)$ that ETAS needs, and it is the stationary Poisson reference — the mandatory null any time-dependent model must beat.

7.1 Adaptive power-law kernel

$$\mu(x,y) = \sum_{i} K_{d_i}!\big(\lVert (x,y) - (x_i,y_i)\rVert\big), \qquad K_{d}(r) = \frac{C(d)}{\big(r^2 + d^2\big)^{s}},$$

with normalization $C(d)$ chosen so each kernel integrates to 1 over the plane, and an adaptive bandwidth $d_i$ set to the horizontal distance from event $i$ to its $n$-th nearest neighbor ($n \approx 6$, optimized retrospectively). The expected number per cell is $\mu(\text{cell})$ × (time scaling) × (fraction $\ge M$ from GR). The catalog is declustered first, so the map reflects the long-term average rather than transient bursts.

Implementation note (do not hard-code the exponent). The kernel exponent $s$ and normalization $C(d)$ vary across the Helmstetter–Kagan–Jackson family of papers — forms with $s = 1$ ($\propto 1/(r^2 + d^2)$) and $s = 3/2$ ($\propto 1/(r^2 + d^2)^{3/2}$) both appear, depending on the specific kernel and normalization. Pin $s$, $C(d)$, and the neighbor count to a specific reference implementation (Helmstetter–Kagan–Jackson 2007 / Werner et al. 2011 or the pyCSEP/floatCSEP code); do not treat $s = 3/2$ as a verified universal constant.

7.2 Role and limits

The best-performing time-independent spatial forecast in the RELM/CSEP California experiments, and the natural $\mu(x,y)$ background for ETAS. Werner et al. (2011) extended it to time-dependent adaptive smoothing. Its limit: it assumes the future resembles the smoothed past locations, so it misses unruptured faults with no recent microseismicity.

References. Helmstetter, Kagan & Jackson (2007), SRL 78(1), 78–86, doi:10.1785/gssrl.78.1.78; Werner, Helmstetter, Jackson & Kagan (2011), BSSA 101, 1630–1648, doi:10.1785/0120090340.


8. BPT / renewal — long-term recurrence

Where paleoseismic or historical data genuinely constrain a fault's mean recurrence interval, a Brownian Passage Time (BPT) renewal model conditions the long-term background over years to decades. Tectonic loading plus Brownian noise drives a state variable to a failure threshold, so inter-event times follow the BPT (inverse-Gaussian / Wald) distribution.

8.1 Recurrence-time density

$$f(t; \mu, \alpha) = \sqrt{\frac{\mu}{2\pi,\alpha^2, t^3}}; \exp!\left(-\frac{(t-\mu)^2}{2,\mu,\alpha^2, t}\right), \qquad t > 0,$$

  • $\mu$ — mean recurrence interval.
  • $\alpha$aperiodicity (coefficient of variation of recurrence times); $\alpha \to 0$ gives quasi-periodic recurrence, $\alpha \approx 1$ gives near-Poisson behaviour.

The key property is that the hazard rate rises from ~0 right after an event, peaks, then plateaus — unlike Poisson's constant hazard and unlike the lognormal's eventually-decreasing hazard. The conditional probability of an event in $[T_e, T_e + \Delta T]$ given $T_e$ years elapsed since the last is

$$P = \frac{\int_{T_e}^{T_e + \Delta T} f(t),dt}{\int_{T_e}^{\infty} f(t),dt}.$$

8.2 Characteristic earthquakes and BPT vs. Poisson

The characteristic-earthquake hypothesis (Schwartz & Coppersmith, 1984) holds that a given fault tends to rupture in similar-size "characteristic" events, producing a magnitude distribution that departs from GR at large $M$ (a bump). It is contested against pure GR.

Poisson (memoryless, constant hazard) is the null. BPT is preferred where paleoseismic / historical recurrence data justify time-dependence (used in UCERF and in the Italian and Japanese national hazard models).

Honest caveat. With only a few observed cycles, $\alpha$ is poorly constrained and the gain over a plain Poisson background is often marginal. Do not claim renewal skill where the data do not support it.

References. Matthews, Ellsworth & Reasenberg (2002), BSSA 92, 2233–2250, doi:10.1785/0120010267; Schwartz & Coppersmith (1984), JGR 89, 5681–5698; Field et al. (2015), UCERF3-TD, BSSA 105.


9. Rate-and-state friction + Coulomb stress — the mechanistic layer

The physics-based layer. Static stress transfer from a fault slip changes the Coulomb failure stress on neighbouring faults; rate-and-state friction then predicts how the seismicity rate responds. These supply mechanistic spatial priors (Coulomb lobes promote or suppress triggering) — optional, feature-flagged covariates, never a standalone forecaster.

9.1 Coulomb failure stress change (King, Stein & Lin, 1994)

$$\Delta\mathrm{CFS} = \Delta\tau - \mu',\Delta\sigma_n$$

(also written $\Delta\mathrm{CFF}$):

  • $\Delta\tau$ — shear-stress change on the target fault plane, positive in the slip direction.
  • $\Delta\sigma_n$ — normal-stress change, positive in compression.
  • $\mu'$ — apparent (effective) friction coefficient absorbing pore-pressure (Skempton) effects, typically $\mu' \approx 0.4$.

$\Delta\mathrm{CFS} &gt; 0$ promotes failure on the target fault; $\Delta\mathrm{CFS} &lt; 0$ inhibits it. The pattern correlates with aftershock locations (the Coulomb "lobes"). It is independent of the regional stress and depends only on source geometry, target geometry, slip sense and $\mu'$.

9.2 Rate-and-state seismicity (Dieterich, 1994)

Dieterich derived how a population of rate-and-state faults responds to a stress change, giving a seismicity rate $R$ relative to background $r$. The state evolves as

$$d\gamma = \frac{1}{A\sigma_n}\big(dt - \gamma, d\tau\big), \qquad R = \frac{r}{\gamma,\dot\tau_r},$$

  • $\gamma$ — state variable (inverse of an instantaneous rate).
  • $A$ — rate-and-state constitutive parameter (lab values $\approx 0.005$–$0.02$).
  • $\sigma_n$ — effective normal stress; $\dot\tau_r$ — reference (background) stressing rate.
  • $A\sigma_n$ — the characteristic stress scale of the response.

For a stress step $\Delta\tau$ (e.g. a Coulomb step), the seismicity-rate response has the closed form

$$R(t) = \frac{r}{\left[\exp!\left(-\dfrac{\Delta\tau}{A\sigma_n}\right) - 1\right] \exp!\left(-\dfrac{t}{t_a}\right) + 1}, \qquad t_a = \frac{A\sigma_n}{\dot\tau_r},$$

where $t_a$ is the aftershock-duration timescale. To first order, the seismicity rate responds exponentially to a Coulomb stress step,

$$\frac{R}{r} = \exp!\left(\frac{\Delta\mathrm{CFS}}{A\sigma}\right).$$

This produces an Omori-like $1/t$ decay of triggered seismicity from first principles — the mechanistic bridge between Coulomb stress transfer and the empirical Omori law. The $\Delta\mathrm{CFS} \rightarrow$ rate-and-state chain is the basis of physics-based CSEP forecasts.

References. King, Stein & Lin (1994), BSSA 84, 935–953; Stein (1999), Nature 402, 605–609, doi:10.1038/45144; Dieterich (1994), JGR 99(B2), 2601–2618, doi:10.1029/93JB02581; Heimisson & Segall (2018), JGR 123, doi:10.1029/2018JB015656.


10. Contested / precursory methods — honest assessment

These are listed because they tempt feature engineering. None is a credible standalone alarm. At most they are weak features whose contribution must be demonstrated in prospective CSEP testing before touching a public number.

Method Core idea Honest verdict
Accelerating Moment Release (AMR) Cumulative Benioff strain accelerates as a power-law toward a critical point before a large event. Largely discredited as a predictor. Hardebeck, Felzer & Michael (2008) showed apparent AMR arises from selection bias (window, area and magnitude range optimized per mainshock) plus ordinary clustering; statistically insignificant under unbiased tests.
LURR (Load–Unload Response Ratio) System response to loading vs. unloading diverges before instability. Plausible lab analogy; field claims of intermediate-term skill exist but are not validated in transparent prospective CSEP tests. Shares the AMR selection-bias risk. Contested.
RTL (Region–Time–Length) Weighted seismicity anomaly (quiescence/activation) around a point. Retrospective successes published; no robust prospective validation. Contested.
Pattern Informatics (PI) Maps regions of anomalous change in seismicity rate as a proxy for stress change. The most testable of this group; entered RELM/CSEP. Some retrospective skill; prospective performance modest and debated. Usable as a feature, not an alarm.

Bottom line. Anything in this table sits behind a feature flag and must beat ETAS + smoothed-seismicity in prospective testing before it touches a public number. The central cautionary lesson — from Hardebeck et al. (2008) — is that optimizing parameters per target manufactures false precursors. (The deep-learning analogue of this trap is told in Models — Analytical / ML.)

References. Hardebeck, Felzer & Michael (2008), JGR 113, B08310, doi:10.1029/2007JB005410; Rundle, Tiampo et al. (2002), PNAS 99(suppl. 1), 2514–2521.


11. From a conditional intensity to the published probability

The single number rendered to a user combines the conditional intensity (from ETAS / R-J / STEP), the Gutenberg–Richter magnitude tail, and the chosen horizon. The expected number of events above a target magnitude $M^*$ in region $A$ over the horizon $[0, T]$ is

$$N_{\ge M^_} = \int_0^T!!\int_A \lambda(t, x, y \mid \mathcal{H}_t),\Phi(M^_), dx, dy, dt, \qquad \Phi(M^_) = 10^{-b(M^_ - M_c)},$$

and, treating events as a non-homogeneous Poisson process over the window, the published probability is

$$P(\ge 1 \text{ event} \ge M^_) = 1 - e^{-N_{\ge M^_}}.$$

The public formula $P = 1 - e^{-N}$ never changes. Only the quality of $\lambda$ improves as the model improves. Quiet days correctly read near climatology, and the number is always shown next to its long-term baseline, so a user reads "X% vs. Y% baseline," never an unanchored figure.

A single threshold vs. the full distribution. A single fixed $M^*$ is brittle (a Chilean $M \ge 6$ is routine; a UK $M \ge 6$ is unheard-of) and discards Gutenberg–Richter information the model already computes. The product instead forecasts the full conditional magnitude distribution (via the GR / $b$-value term) and derives the public scalar as an exceedance probability at one or more region-appropriate thresholds. $M_{\max}$ is specified per region — it bounds the exceedance integral and sets the tail probability of the rare, high-impact events. See Models — Employed for the full target definition and how the bounds are computed, and Evaluation for how every model on this page is scored.


References

  1. Gutenberg, B. & Richter, C.F. (1944). Frequency of earthquakes in California. BSSA 34, 185–188.
  2. Aki, K. (1965). Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits. Bull. Earthq. Res. Inst. 43, 237–239.
  3. Tinti, S. & Mulargia, F. (1987). Confidence intervals of b-values for grouped magnitudes. BSSA 77(6), 2125–2134.
  4. Wiemer, S. & Wyss, M. (2000). Minimum magnitude of completeness in earthquake catalogs. BSSA 90(4), 859–869. doi:10.1785/0119990114
  5. Woessner, J. & Wiemer, S. (2005). Assessing the quality of earthquake catalogues: estimating the magnitude of completeness and its uncertainty. BSSA 95(2), 684–698. doi:10.1785/0120040007
  6. Utsu, T., Ogata, Y. & Matsu'ura, R.S. (1995). The centenary of the Omori formula for a decay law of aftershock activity. J. Phys. Earth 43, 1–33. doi:10.4294/jpe1952.43.1
  7. Ogata, Y. (1983). Estimation of the parameters in the modified Omori formula for aftershock frequency by the maximum likelihood procedure. J. Phys. Earth 31, 115–124.
  8. Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes. JASA 83(401), 9–27. doi:10.1080/01621459.1988.10478560
  9. Ogata, Y. (1998). Space–time point-process models for earthquake occurrences. Ann. Inst. Statist. Math. 50(2), 379–402. doi:10.1023/A:1003403601725
  10. Zhuang, J., Ogata, Y. & Vere-Jones, D. (2002). Stochastic declustering of space–time earthquake occurrences. JASA 97(458), 369–380. doi:10.1198/016214502760046925
  11. Reasenberg, P.A. & Jones, L.M. (1989). Earthquake hazard after a mainshock in California. Science 243(4895), 1173–1176. doi:10.1126/science.243.4895.1173
  12. Reasenberg, P.A. & Jones, L.M. (1994). Earthquake aftershocks: update. Science 265, 1251–1252. doi:10.1126/science.265.5176.1251
  13. Page, M.T. et al. (2016). Three ingredients for improved global aftershock forecasts. BSSA 106(5), 2290–2301. doi:10.1785/0120160073
  14. Gerstenberger, M.C., Wiemer, S., Jones, L.M. & Reasenberg, P.A. (2005). Real-time forecasts of tomorrow's earthquakes in California. Nature 435, 328–331. doi:10.1038/nature03622
  15. Rhoades, D.A. & Evison, F.F. (2004). Long-range earthquake forecasting with every earthquake a precursor according to scale. Pure Appl. Geophys. 161, 47–72. doi:10.1007/s00024-003-2434-9
  16. Rhoades, D.A. (2007). Application of the EEPAS model to forecasting earthquakes of moderate magnitude in Southern California. SRL 78, 110–115.
  17. Helmstetter, A., Kagan, Y.Y. & Jackson, D.D. (2007). High-resolution time-independent grid-based forecast for M ≥ 5 earthquakes in California. SRL 78(1), 78–86. doi:10.1785/gssrl.78.1.78
  18. Werner, M.J., Helmstetter, A., Jackson, D.D. & Kagan, Y.Y. (2011). High-resolution long-term and short-term earthquake forecasts for California. BSSA 101, 1630–1648. doi:10.1785/0120090340
  19. Matthews, M.V., Ellsworth, W.L. & Reasenberg, P.A. (2002). A Brownian model for recurrent earthquakes. BSSA 92, 2233–2250. doi:10.1785/0120010267
  20. Schwartz, D.P. & Coppersmith, K.J. (1984). Fault behavior and characteristic earthquakes. JGR 89, 5681–5698.
  21. Field, E.H. et al. (2015). UCERF3-TD: Time-dependent earthquake rupture forecast. BSSA 105.
  22. King, G.C.P., Stein, R.S. & Lin, J. (1994). Static stress changes and the triggering of earthquakes. BSSA 84, 935–953.
  23. Stein, R.S. (1999). The role of stress transfer in earthquake occurrence. Nature 402, 605–609. doi:10.1038/45144
  24. Dieterich, J. (1994). A constitutive law for rate of earthquake production and its application to earthquake clustering. JGR 99(B2), 2601–2618. doi:10.1029/93JB02581
  25. Heimisson, E.R. & Segall, P. (2018). Constitutive law for earthquake production based on rate-and-state friction: Dieterich 1994 revisited. JGR 123. doi:10.1029/2018JB015656
  26. Hardebeck, J.L., Felzer, K.R. & Michael, A.J. (2008). Improved tests reveal that the accelerating moment release hypothesis is statistically insignificant. JGR 113, B08310. doi:10.1029/2007JB005410
  27. Rundle, J.B., Tiampo, K.F. et al. (2002). Self-organization, criticality and seismic hazard. PNAS 99(suppl. 1), 2514–2521.
  28. Jordan, T.H. et al. (2011). Operational Earthquake Forecasting: state of knowledge and guidelines for utilization (ICEF Report). Annals of Geophysics 54(4), 315–391. doi:10.4401/ag-5350

See also: Models — Analytical / ML · Models — Employed · Evaluation · Methodology. All equations are transcribed from the primary/authoritative sources cited above.

Clone this wiki locally