Reasenberg Jones Model

Reasenberg–Jones Model — Operational Aftershock Probabilities

The transparent, closed-form workhorse of operational aftershock forecasting. Given a felt mainshock, the Reasenberg–Jones (R–J) model answers the single most-asked question in the days after a damaging earthquake: what is the probability that a still-larger event follows in the next hours, days, or week? It is the analytical core of the USGS Operational Aftershock Forecasting (OAF) service and a building block of STEP and the OEF-Italy ensemble. This page is a self-contained, deep treatment of that one model: intuition and history, the governing equations with derivation, parameter estimation, assumptions, strengths and limitations, its operational role, a worked example, and references.

Honest framing. R–J produces a conditional probability, never a deterministic prediction. A prediction is a yes/no statement that an event will occur; a forecast gives a probability strictly in $(0,1)$ (Jordan et al., 2011). R–J answers "how many, and how likely is one or more" — it does not say where on a fault, exactly when, or how large. A single outcome neither confirms nor refutes a probabilistic forecast. See Honest-Limits.

---

Table of contents

1. Intuition and history
2. Where it sits among the models
3. The governing rate model
4. From rate to expected number to probability
5. The three parameter regimes
6. Parameter estimation
7. Assumptions and failure modes
8. Strengths and limitations
9. Role in operational earthquake forecasting
10. Worked illustration
11. How R–J is used in this product
12. References

---

1. Intuition and history

Two robust empirical laws of seismology, discovered independently in the 19th and 20th centuries, together pin down almost everything we can say about an aftershock sequence:

• Omori–Utsu decay — the rate of aftershocks falls off roughly as a power of elapsed time. A sequence is most dangerous in the first hours and decays predictably thereafter.
• Gutenberg–Richter (GR) magnitude–frequency — for every magnitude-6 aftershock there are about ten magnitude-5s and a hundred magnitude-4s, with a log-linear slope $b \approx 1$.

Reasenberg and Jones (1989) made the decisive observation that these two laws can simply be multiplied: a separable rate that is exponential in magnitude (GR) and power-law in time (Omori–Utsu). Crucially, they showed that the productivity of an aftershock sequence scales with the mainshock magnitude in a way that is roughly constant across California sequences — so a single set of "generic" parameters gives a useful forecast the instant a mainshock is located, before any of its own aftershocks have even been recorded. This is what made operational, real-time aftershock advisories possible.

The model was first deployed by the USGS and the California Office of Emergency Services to issue public aftershock advisories after the 1989 Loma Prieta earthquake, and it remained the backbone of US aftershock statements for two decades. The 1994 update (Reasenberg & Jones, 1994) refined the California generic parameters. Page et al. (2016) generalized the scheme to a global, tectonic-regime-aware Bayesian form and folded it into the modern USGS OAF service, which now runs R–J alongside a full ETAS model.

Why this matters operationally. The single most reliably forecastable thing in all of seismology is the decay of an aftershock sequence. R–J is the smallest, most transparent model that captures it and turns it into a number a civil-protection officer can act on.

---

2. Where it sits among the models

R–J is a temporal–magnitude model: it forecasts how many aftershocks of a given size will follow, but has no spatial term — it does not say where within the aftershock zone. It is the direct ancestor of two richer models:

flowchart TD
    GR["Gutenberg–Richter<br/>magnitude term<br/>10^{b(M_m − M)}"] --> RJ
    OU["Omori–Utsu<br/>temporal decay<br/>(t + c)^{−p}"] --> RJ
    RJ["Reasenberg–Jones<br/>λ(t, M) = magnitude × time<br/>(temporal–magnitude, no space)"]
    RJ -->|"add spatial grid +<br/>map to shaking via GMPE"| STEP["STEP<br/>gridded shaking probability"]
    RJ -->|"let every event trigger<br/>its own offspring (self-exciting)"| ETAS["ETAS<br/>full space–time clustering"]
    RJ --> OAF["USGS OAF<br/>operational service<br/>(R–J + ETAS)"]

Loading

• Add a spatial grid and a ground-motion conversion and you get STEP.
• Let every event (not just the mainshock) trigger its own aftershocks, and you get the self-exciting ETAS process. R–J is the special case of ETAS in which only the mainshock is treated as a trigger and the rate is integrated over space.

In this product R–J is the transparent fallback and sanity-check that runs alongside ETAS — see Models-Employed.

---

3. The governing rate model

The Reasenberg–Jones rate of aftershocks of magnitude $\ge M$ at elapsed time $t$ after a mainshock of magnitude $M_m$ is

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

The numerator is the Gutenberg–Richter magnitude term and the denominator is the modified-Omori temporal decay. The parameters are:

Symbol	Name	Meaning	Typical range
$a$	Productivity	Aftershock abundance of the sequence (more negative ⇒ fewer aftershocks).	$-2$ to $-1$ (regime-dependent)
$b$	GR slope	Magnitude–frequency slope; globally $\approx 1$.	$0.8$–$1.2$
$c$	Omori offset	Small time offset (days) regularizing the singularity at $t = 0$; sensitive to early incompleteness.	$0.01$–$0.3$ d
$p$	Omori exponent	Decay rate of the sequence.	$0.9$–$1.5$

Reading the structure. Hold $t$ fixed and the rate is log-linear in $M$ with slope $-b$ — pure Gutenberg–Richter. Hold $M$ fixed and the rate decays as $(t+c)^{-p}$ — pure Omori–Utsu. The separability (the magnitude and time factors multiply, never interact) is the model's defining simplification and the source of both its transparency and its limits.

The productivity term is often written so that $a$ absorbs the mainshock magnitude scaling: the combination $10^{a + b M_m}$ sets the overall sequence size, and a larger mainshock produces proportionally more aftershocks at every magnitude. This mainshock-magnitude scaling of productivity is the empirical regularity that lets generic parameters work.

---

4. From rate to expected number to probability

A rate is not yet a probability. Two steps convert $\lambda(t, M)$ into the number the public sees.

4.1 Expected number (integrate the rate over the horizon)

The expected number of aftershocks $\ge M$ in the time window $[t_1, t_2]$ is the integral of the rate. Because the magnitude factor is constant in time, it pulls out of the integral:

$$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.$$

The remaining time integral is elementary. For $p \neq 1$,

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

and for the special case $p = 1$,

$$\int_{t_1}^{t_2} (t+c)^{-1}, dt = \ln(t_2 + c) - \ln(t_1 + c).$$

Note that for $p > 1$ the integral converges as $t_2 \to \infty$, so a sequence has a finite expected total number of aftershocks — a physically necessary property.

4.2 Probability of one or more (Poisson assumption)

R–J treats aftershock occurrence as a non-homogeneous Poisson process with the time-varying rate $\lambda(t, M)$. For a Poisson process, the probability of observing zero events in a window with expected count $N$ is $e^{-N}$, so the probability of one or more $\ge M$ aftershocks is

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

This is the canonical "probability of a $\ge M$ aftershock in the next $N$ days" that operational advisories report. Two limiting checks: when $N \ll 1$, $P \approx N$ (rare-event regime, the probability is just the expected count); when $N \gg 1$, $P \to 1$ (a large aftershock is essentially certain).

Derivation note — why Poisson. Conditioning on the deterministic rate function $\lambda(t)$, aftershock times are an inhomogeneous Poisson process, for which the count in any window is Poisson with mean $N = \int \lambda,dt$. Treating the parameters as fixed (rather than random) is the approximation; the Bayesian regime (§5) relaxes it by integrating over parameter uncertainty, which fattens the tail of $P$ and is the honest choice for a public number.

---

5. The three parameter regimes

R–J is run in three regimes that trade immediacy against specificity:

flowchart LR
    A["Mainshock located<br/>(t = 0)"] --> G["Generic<br/>regional (a,b,c,p)<br/>available instantly"]
    G -->|"aftershocks accumulate"| B["Bayesian<br/>generic prior updated<br/>by the sequence"]
    B -->|"enough data"| S["Sequence-specific<br/>refit to this sequence"]
    G -.->|"if data sparse,<br/>stay generic"| G

Loading

1. Generic. Fixed $(a, b, c, p)$ taken from a catalog of past sequences in the same tectonic setting. Available the instant the mainshock is located — this is what makes a forecast possible before any aftershocks occur. The classic 1989/1994 California generic values exist, but they do not transfer to other tectonic regimes (subduction megathrust productivity and decay differ markedly from California strike-slip — do not reuse California numbers elsewhere).

2. Sequence-specific. Once enough aftershocks have accumulated, refit $(a, c, p)$ (and possibly $b$) to the ongoing sequence by maximum likelihood. More accurate, but only after a delay and only if the sequence is productive enough to constrain the fit.

3. Bayesian. A generic prior is updated by the observed aftershocks via Bayes' rule. This is the modern operational default: it behaves like the generic model when data are sparse and converges to the sequence-specific fit as data accumulate, with calibrated uncertainty throughout. Page et al. (2016) built the global, tectonic-regime-aware version that the USGS OAF uses. Integrating over the posterior — rather than plugging in point estimates — is what keeps the published probability honest.

---

6. Parameter estimation

6.1 Likelihood

Given aftershock times ${t_i}$ above a completeness magnitude $M_c$, with magnitudes contributing through the GR term, the point-process log-likelihood for $(a, c, p)$ is

$$\ln L = \sum_i \ln \lambda(t_i) ;-; \int_{t_{\text{start}}}^{t_{\text{end}}} \lambda(t), dt,$$

the sum over observed events minus the compensator (the integrated rate). The compensator term is exactly the expected-number integral of §4.1; it is what makes the fit probabilistic and calibratable rather than a naive curve-fit to the binned decay. Maximizing $\ln L$ gives the sequence-specific MLE (the Omori–Utsu estimation of Ogata, 1983).

6.2 Practical estimation rules

• $b$-value: estimate with the Aki–Utsu MLE on events above $M_c$ (binning-corrected); never hard-code $b = 1$. See Models-Classical §1 and Glossary.
• Early incompleteness: immediately after a large mainshock, small aftershocks are missed (the seismograms are saturated and overlapping), which biases $c$ upward and $p$ downward if ignored. The standard remedy is to start the fit after a delay $t_{\text{start}}$, past the time-of-completeness, or to model the time-varying $M_c(t)$ explicitly.
• Generic-prior fallback: with few aftershocks, the sequence-specific MLE is unstable — fall back to the Bayesian regime so the regional prior dominates.
• Tooling: the USGS OAF (aftershock/AftershockStatistics) reference codes, and the Omori–Utsu MLE in R ETAS / Python pyetas.

---

7. Assumptions and failure modes

Assumption	What breaks it
Separable magnitude × time (GR slope constant over the sequence).	Real $b$ can drift through a sequence; foreshock/swarm episodes violate separability.
Single trigger — only the mainshock seeds aftershocks.	Large aftershocks have their own aftershocks (secondary triggering); ETAS captures this, R–J does not.
Non-homogeneous Poisson counts.	Spatial-temporal clustering makes counts over-dispersed relative to Poisson; the Bayesian regime partly compensates.
Stationary parameters over the forecast window.	A large secondary event resets the clock; the model must be re-run.
Catalog complete above $M_c$ for the whole window.	Early incompleteness (the most common operational failure) biases $c$, $p$.
Generic parameters transfer within a tectonic regime.	Cross-regime transfer (e.g. California → subduction) is invalid.

The dominant practical failure is under-forecasting secondary sequences: because only the mainshock is a trigger, R–J systematically misses the burst of new aftershocks that follows a large aftershock. This is precisely the gap ETAS fills, and why operational systems run both.

---

8. Strengths and limitations

Strengths

• Transparent and closed-form — every number traces to four interpretable parameters; auditable for public communication.
• Immediate — the generic regime produces a forecast the instant the mainshock is located.
• Cheap and robust — no spatial grid, trivial to compute, hard to break; an ideal sanity-check baseline and fallback.
• Battle-tested — three decades of operational use; the reference any short-term model is compared against.

Limitations

• No spatial resolution — answers "how many", not "where". STEP and ETAS add space.
• No secondary triggering — systematically under-forecasts the aftershocks of large aftershocks.
• Aftershock-only — it is a post-mainshock model; it says nothing about background or mainshock/foreshock anticipation.
• Parameter portability — generic values are regime-specific; misuse across regimes is a real pitfall.

---

9. Role in operational earthquake forecasting

R–J is not a research curiosity — it is live operational infrastructure:

• USGS Operational Aftershock Forecasting (OAF). The OAF service issues public aftershock forecasts after significant earthquakes worldwide. Its analytical core is the Reasenberg–Jones model in the global Bayesian form of Page et al. (2016), run alongside ETAS. Forecasts are published as probabilities of one-or-more $\ge M$ aftershocks over 1-day, 1-week, 1-month and 1-year windows, with uncertainty bounds.
• STEP (Gerstenberger et al., 2005) wraps the R–J clustering term, plus a background, into gridded hourly shaking-probability maps — see STEP-Model.
• OEF-Italy runs an ensemble in which a STEP-type component (built on R–J clustering) is one of three models blended with ETAS and ETES.

The reason R–J persists despite ETAS being more complete is communication and robustness: a four-parameter, closed-form model whose every term has a physical meaning is far easier to explain, audit, and defend to the public and to civil-protection authorities than a black-box. In CSEP-style retrospective tests across a full decade of California next-day forecasting, no single model dominates — aftershock-sequence-tuned models excel during sequences, ETAS is the consistent generalist — which is exactly why operational systems run several models, R–J among them, side by side. See Evaluation-and-Tests.

---

10. Worked illustration

Illustrative numbers only — for intuition, not a forecast.

Suppose an $M_m = 6.0$ mainshock occurs and we adopt generic parameters $a = -1.7$, $b = 0.9$, $c = 0.05$ d, $p = 1.1$. What is the probability of one or more $\ge 5.0$ aftershocks in the first 7 days ($t_1 = 0$, $t_2 = 7$ d)?

Step 1 — magnitude term. $10^{a + b(M_m - M)} = 10^{-1.7 + 0.9(6.0 - 5.0)} = 10^{-0.8} \approx 0.158$.

Step 2 — time integral ($p = 1.1 \neq 1$): $$\int_0^7 (t + 0.05)^{-1.1}, dt = \frac{(7.05)^{-0.1} - (0.05)^{-0.1}}{-0.1} = \frac{0.819 - 1.349}{-0.1} \approx 5.30.$$

Step 3 — expected number. $N = 0.158 \times 5.30 \approx 0.84$.

Step 4 — probability. $P = 1 - e^{-0.84} \approx 0.57$ — about a 57 % chance of at least one $\ge 5.0$ aftershock in the first week.

The same machinery, with $M = 6.0$ (a $\ge$-mainshock event), gives $10^{-1.7} \approx 0.020$, $N \approx 0.020 \times 5.30 \approx 0.106$, and $P = 1 - e^{-0.106} \approx 0.10$ — roughly a 10 % chance the sequence produces an event as large as, or larger than, the mainshock in the first week. That order of magnitude is consistent with the well-known empirical rule that a mainshock is followed by an equal-or-larger event a small but non-negligible fraction of the time.

Read it honestly. A 57 % forecast that does not produce a $\ge 5$ aftershock is not wrong, and a 10 % forecast that does is not a failure — these are probabilities, scored over many sequences (see Evaluation-and-Tests), not promises about one sequence.

---

11. How R–J is used in this product

In CAOS_SEISMIC, Reasenberg–Jones is the transparent fallback and sanity-check that runs beside the ETAS conditional intensity (see Models-Employed and Methodology-History):

• Its closed-form probability is computed every cycle and cross-checked against the ETAS forecast; large divergences flag a data or fitting problem.
• It supplies an interpretable explanation of the short-horizon aftershock signal that accompanies the published number.
• Generic parameters are taken from, and re-estimated for, the operating tectonic regime — the California 1989/1994 values are never reused outside California.
• The R–J probability feeds the same exceedance-probability machinery, Gutenberg–Richter magnitude tail, and CSEP scoring as every other model in the stack (Pipeline, Evaluation-and-Tests).

It is never an alarm. Its output is a bounded, calibrated probability conditioned on the present sequence — see Honest-Limits.

---

References

1. 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
2. Reasenberg, P. A. & Jones, L. M. (1994). Earthquake aftershocks: update. Science 265(5176), 1251–1252. doi:10.1126/science.265.5176.1251
3. Page, M. T., van der Elst, N., Hardebeck, J., Felzer, K. & Michael, A. J. (2016). Three ingredients for improved global aftershock forecasts: tectonic region, time-dependent catalog incompleteness, and intersequence variability. Bulletin of the Seismological Society of America 106(5), 2290–2301. doi:10.1785/0120160073
4. Utsu, T., Ogata, Y. & Matsu'ura, R. S. (1995). The centenary of the Omori formula for a decay law of aftershock activity. Journal of Physics of the Earth 43(1), 1–33. doi:10.4294/jpe1952.43.1
5. Ogata, Y. (1983). Estimation of the parameters in the modified Omori formula for aftershock frequencies by the maximum likelihood procedure. Journal of Physics of the Earth 31(2), 115–124. doi:10.4294/jpe1952.31.115
6. 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
7. Aki, K. (1965). Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits. Bulletin of the Earthquake Research Institute 43, 237–239.
8. Jordan, T. H., Chen, Y.-T., Gasparini, P., Madariaga, R., Main, I., Marzocchi, W., Papadopoulos, G., Sobolev, G., Yamaoka, K. & Zschau, J. (2011). Operational earthquake forecasting: state of knowledge and guidelines for utilization. Annals of Geophysics 54(4), 315–391. doi:10.4401/ag-5350

---

Related pages: Models-Classical · STEP-Model · EEPAS-Model · Models-Employed · Methodology-History · Evaluation-and-Tests · Honest-Limits · Glossary.