How to Estimate the Gutenberg-Richter b-Value in Python

Problem statement

In almost every region, each step of one magnitude unit up brings roughly ten times fewer earthquakes. The Gutenberg-Richter law writes that as log10 N = a − bM, where N is the number of events of magnitude M or more, and the b-value, the slope, is the single most used number in statistical seismology: it feeds hazard models, recurrence estimates and forecasts. Estimating it well takes two decisions that a straight line through a log plot hides. The first is where the catalogue is complete, because below the magnitude of completeness small events are missing and the slope flattens. The second is how to fit, because a least-squares line on cumulative counts gives the large, rare events far too much weight.

This page downloads five years of the USGS catalogue for the San Jacinto fault zone in southern California, estimates the magnitude of completeness, computes b by maximum likelihood with its uncertainty, checks how b depends on the cut-off, and uses the fit to predict how often larger earthquakes occur.

Quick answer

Find the magnitude of completeness, keep the events at or above it, and use the Aki-Utsu maximum-likelihood estimator with the correction for binned magnitudes:

m = cat.mag.round(1)
counts = m.value_counts().sort_index()
mc = round(counts.idxmax() + 0.2, 1)                        # maximum curvature, plus the usual 0.2
b, n, sigma = b_mle(m, mc)
print(f"maximum curvature at M{counts.idxmax():.1f}; Mc = M{mc:.1f}")
print(f"b = {b:.3f} ± {sigma:.3f} from {n:,} earthquakes of M{mc:.1f} and above")

maximum curvature at M0.6; Mc = M0.8
b = 1.059 ± 0.008 from 18,452 earthquakes of M0.8 and above

The most populated bin is M0.6, so Mc is taken as M0.8. Above it, 18,452 earthquakes give b = 1.059 with a standard error of 0.008: each magnitude unit up brings about 11.5 times fewer events.

Step-by-step solution

1. The frequency-magnitude distribution

Count events per 0.1 magnitude bin and cumulatively. The non-cumulative counts peak where detection starts to fail:

cum = counts[::-1].cumsum()[::-1]
for mag in (0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0):
    print(f"M{mag:.1f}: {counts.get(mag, 0):>5,} in the bin, {cum.loc[mag:].iloc[0]:>6,} at or above")
a = np.log10(n) + b * mc
print(f"a-value for 5 years: {a:.2f}; per year: {a - np.log10(5):.2f}")

M0.5: 3,372 in the bin, 30,156 at or above
M1.0: 2,798 in the bin, 11,444 at or above
M1.5:   630 in the bin,  3,271 at or above
M2.0:   238 in the bin,  1,021 at or above
M2.5:    49 in the bin,    296 at or above
M3.0:    17 in the bin,    105 at or above
M3.5:    11 in the bin,     30 at or above
M4.0:     2 in the bin,      6 at or above
a-value for 5 years: 5.11; per year: 4.41

Each magnitude unit up divides the count by roughly ten: 11,444 events of M1.0 or more, 1,021 of M2.0 or more, 105 of M3.0 or more and 6 of M4.0 or more. The bin counts below M1.0 no longer grow tenfold, which is where detection begins to miss events. The a-value, the logarithm of the expected number of events of M0 or more, is 4.41 per year.

2. How b depends on the cut-off

The estimate is only meaningful at or above completeness. Compute b for a range of cut-offs; below Mc it is too low, above it it should be stable within its uncertainty:

for c in np.arange(0.5, 2.6, 0.3):
    bc, nc, sc = b_mle(m, round(c, 1))
    print(f"cut-off M{c:.1f}: b = {bc:.3f} ± {sc:.3f} ({nc:,} events)")

cut-off M0.5: b = 0.879 ± 0.004 (30,156 events)
cut-off M0.8: b = 1.059 ± 0.008 (18,452 events)
cut-off M1.1: b = 1.034 ± 0.011 (8,646 events)
cut-off M1.4: b = 1.024 ± 0.016 (4,186 events)
cut-off M1.7: b = 1.020 ± 0.022 (2,048 events)
cut-off M2.0: b = 1.045 ± 0.033 (1,021 events)
cut-off M2.3: b = 1.005 ± 0.043 (484 events)
cut-off M2.6: b = 1.035 ± 0.062 (247 events)

Below completeness the estimate is far too low: 0.879 with a cut-off at M0.5. From M0.8 upwards it stays between 1.005 and 1.059, and every value is within about two standard errors of the others, while the error grows from 0.008 to 0.062 as the sample shrinks from 18,452 to 247 events.

3. The uncertainty by bootstrap

The formula's uncertainty assumes the cut-off is right. Resampling the catalogue and re-estimating both Mc and b gives an uncertainty that includes the choice of Mc:

rng = np.random.default_rng(42)
boot = []
for _ in range(500):
    s = m.sample(len(m), replace=True, random_state=rng)
    c = round(s.value_counts().idxmax() + 0.2, 1)
    boot.append((c, b_mle(s, c)[0]))
boot = pd.DataFrame(boot, columns=["mc", "b"])
print(f"bootstrap Mc: {boot.mc.value_counts().sort_index().to_dict()}")
print(f"bootstrap b: mean {boot.b.mean():.3f}, standard deviation {boot.b.std():.3f}, 95% range {boot.b.quantile(0.025):.3f} to {boot.b.quantile(0.975):.3f}")

bootstrap Mc: {0.8: 500}
bootstrap b: mean 1.059, standard deviation 0.008, 95% range 1.045 to 1.075

Every one of the 500 resamples puts Mc at M0.8, so here the choice of cut-off adds nothing and the bootstrap spread, 0.008, equals the formula's. With a smaller catalogue, Mc would move between resamples and the bootstrap range would be wider than the formula.

Code examples

The data used on this page

Every earthquake of magnitude 0 or more from 2020 to 2024 in the USGS ComCat catalogue between 33° and 34°N and 117.5° and 116°W, the San Jacinto fault zone and its surroundings, from the FDSN event service. The service returns at most 20,000 events per request, so the setup asks for one year at a time and checks each year is below the cap.

from pathlib import Path

import numpy as np
import pandas as pd
import requests

Path("data").mkdir(exist_ok=True)
Q = "https://earthquake.usgs.gov/fdsnws/event/1/query"
BOX = dict(minlatitude=33, maxlatitude=34, minlongitude=-117.5, maxlongitude=-116)
def year(y):
    p = Path(f"data/sanjacinto_m0_{y}.csv")
    if not p.exists():
        r = requests.get(Q, params=dict(format="csv", starttime=f"{y}-01-01", endtime=f"{y + 1}-01-01",
                                        minmagnitude=0, orderby="time-asc", **BOX), timeout=600)
        r.raise_for_status()
        p.write_text(r.text)
    d = pd.read_csv(p)
    assert len(d) < 20000, f"{y}: {len(d)} rows, at the cap"
    return d
def b_mle(mags, mc, dm=0.1):
    """Aki-Utsu b-value for magnitudes binned at dm, with the Shi and Bolt (1982) standard error."""
    x = mags[mags >= mc - 1e-9]
    b = np.log10(np.e) / (x.mean() - (mc - dm / 2))
    sigma = 2.3 * b ** 2 * np.sqrt(((x - x.mean()) ** 2).sum() / (len(x) * (len(x) - 1)))
    return b, len(x), sigma
raw = pd.concat([year(y) for y in range(2020, 2025)], ignore_index=True)
cat = raw[raw["type"] == "earthquake"].reset_index(drop=True)
print(f"{len(raw):,} events, {len(cat):,} earthquakes; M{cat.mag.min():.2f} to M{cat.mag.max():.2f}")
print(f"magnitude types: {cat.magType.value_counts().to_dict()}")

38,503 events, 37,620 earthquakes; M0.00 to M4.87
magnitude types: {'ml': 37582, 'mw': 20, 'mh': 16, 'mlr': 2}

The five years hold 38,503 events, of which 37,620 are earthquakes; the other 883 are 878 quarry blasts and 5 chemical explosions, which the setup drops. Almost all, 37,582, have local magnitudes (ml); only the largest few carry a moment magnitude. The largest is M4.87.

Example 1: what the fit predicts for larger earthquakes

Use the a- and b-values to predict the yearly number of earthquakes above a magnitude, and compare with what the five years recorded:

per_year = a - np.log10(5)
for mag in (3.0, 4.0, 5.0):
    pred = 10 ** (per_year - b * mag)
    obs = (m >= mag).sum() / 5
    print(f"M{mag:.0f}+: predicted {pred:.2f} a year (one every {1 / pred:.1f} years), observed {obs:.2f} a year")

M3+: predicted 17.26 a year (one every 0.1 years), observed 21.00 a year
M4+: predicted 1.51 a year (one every 0.7 years), observed 1.20 a year
M5+: predicted 0.13 a year (one every 7.6 years), observed 0.00 a year

The fit predicts 17 earthquakes of M3 or more a year against 21 observed, and 1.5 of M4 or more against 1.2. It predicts one M5 or larger every 7.6 years, and none happened in these five years; five years is too short to test that, which is the point of extrapolating with b.

Explanation

The Gutenberg-Richter relation is empirical, but it holds remarkably widely, and b is close to 1 in most tectonic settings: each magnitude unit brings about ten times fewer events. Departures carry information. Values well above 1 are reported in volcanic and geothermal areas and in swarms; values below 1 have been linked to high stress. Comparing b-values between places or periods is only meaningful when each is estimated at or above its own completeness, with its uncertainty, because a difference of 0.1 is often within the error of a few hundred events.

Maximum likelihood is the right estimator because the magnitudes above Mc follow an exponential distribution, and the Aki (1965) estimator is the maximum-likelihood estimate of its rate; Utsu's correction subtracts half a bin from Mc, because a magnitude reported as 1.5 stands for 1.45 to 1.55. Fitting a straight line by least squares to the cumulative counts is common and wrong: cumulative counts are not independent, and the last few points, single large events, pull the slope. Maximum curvature, the bin with the most events, is a quick estimate of Mc that tends to be too low, which is why 0.2 is commonly added.

Edge cases or notes

• Mixed magnitude types: local magnitudes for small events and moment magnitudes for large ones are not the same scale; a catalogue that switches types at some magnitude can bend the distribution.
• Aftershock sequences: Mc rises for hours to days after a large earthquake, when small events are hidden in the coda of larger ones; a catalogue dominated by one sequence needs a time-varying Mc.
• Binning: the correction assumes magnitudes rounded to 0.1; ComCat reports two decimals for many networks, so round consistently before counting.
• Small samples: with fewer than about 100 events above Mc, b is too uncertain to compare with anything.

Internal links

• Earthquake catalogues explained: magnitudes, completeness and declustering — what a catalogue contains
• Fixing a b-value that is too low — the cut-off below completeness
• The Gutenberg-Richter law explained — what a and b mean

FAQ

How do I calculate the b-value in Python?

Keep the earthquakes at or above the magnitude of completeness and compute b = log10(e) / (mean magnitude − (Mc − Δm/2)), where Δm is the magnitude bin width, usually 0.1.

What is a normal b-value?

Close to 1 in most tectonic regions; values between about 0.8 and 1.2 are common.

How do I find the magnitude of completeness?

The maximum-curvature method takes the magnitude bin with the most events; adding 0.2 corrects its tendency to underestimate. Check that b is stable for cut-offs above it.

Why not fit a line to the log of the cumulative counts?

The cumulative counts are not independent and the few large events dominate the fit; maximum likelihood is unbiased and gives an uncertainty.