🔧 Fix Natural hazards & risk Data cleaning Natural Hazards & Risk · step 16 Notebook By · Updated

Earthquake counts rise over time in my catalogue

On this page

Problem statement

Plot events per year from a downloaded catalogue and the line climbs. Plot them per decade and the recent decades tower over the old ones. The obvious reading is that earthquakes are becoming more frequent — and it is almost always wrong.

A catalogue records what the network detected. Networks get denser, thresholds fall, magnitude formulas change, and one large earthquake adds thousands of aftershocks to whichever year it happened in. In 55 years of central California, the annual count of M≥2 events goes from 67 in 1973 to 1,944 in 1974, when the regional network came up. Nothing happened to the faults that year.

Quick answer

Count above a magnitude the catalogue is complete for over the whole period, not above the API's minimum:

import io, json, zipfile
import pandas as pd, requests, time

API = "https://earthquake.usgs.gov/fdsnws/event/1/query"
BOX = [-122.5, 35.5, -117.5, 38.5]                   # central California
# each year's answer as ComCat served it on 9 October 2026 (reviews keep revising the catalogue)
SAVED = zipfile.ZipFile(io.BytesIO(requests.get(
    "https://data.spatialworkflow.io/comcat-central-california-1970-2024-m2-2026-10-09.zip", timeout=300).content))

def fetch_year(y, minmag=2.0):
    params = {"format": "geojson", "starttime": f"{y}-01-01",
              "endtime": f"{y}-12-31T23:59:59", "minmagnitude": minmag,
              "minlongitude": BOX[0], "minlatitude": BOX[1],
              "maxlongitude": BOX[2], "maxlatitude": BOX[3],
              "orderby": "time-asc", "limit": 20000}
    # live: r = requests.get(API, params=params, timeout=300); f = r.json()["features"]
    f = json.loads(SAVED.read(f"{y}.json"))["features"]
    if len(f) == 20000:
        raise RuntimeError(f"{y} hit the cap")
    return pd.DataFrame([{"time": pd.to_datetime(x["properties"]["time"], unit="ms", utc=True),
                          "mag": x["properties"]["mag"], "magtype": x["properties"]["magType"],
                          "lat": x["geometry"]["coordinates"][1],
                          "lon": x["geometry"]["coordinates"][0]} for x in f])

eq = pd.concat([fetch_year(y) for y in range(1970, 2025)], ignore_index=True)
eq["year"] = eq.time.dt.year

for m in (2, 3, 4, 5):
    by_dec = eq[eq.mag >= m].groupby(eq.year // 10 * 10).size()
    print(f"M>={m}: " + "  ".join(f"{d}s {n:6,}" for d, n in by_dec.items()))
M>=2: 1970s  9,240  1980s 25,560  1990s 18,255  2000s 15,249  2010s 14,403  2020s  9,358
M>=3: 1970s  2,020  1980s  4,187  1990s  1,934  2000s  1,366  2010s  1,665  2020s  1,045
M>=4: 1970s    173  1980s    470  1990s    158  2000s    105  2010s    148  2020s    108
M>=5: 1970s      7  1980s     52  1990s     14  2000s      9  2010s     11  2020s     10

The 1980s lead every row, and the 2020s (five years so far) trail every row. Neither fact is about seismicity: the 1980s carry two large aftershock sequences, the 1970s are missing events at every magnitude, and the last row is half a decade. As they stand these counts are not comparable at any threshold.

Three causes of rising earthquake counts, with the test and the fix for each. why the count roseWhat is actually happeningThe test, and the fixThe network got denser, so smallevents are detectedRepeat the count above a fixedmagnitude: cut at the worst McA big mainshock filled a decadewith aftershocksDecluster, then recount: the peakmoves or disappearsThe magnitude scale changedbetween erasGroup by magType: a thresholdcatches a different populationIt is a genuine rise (fluidinjection, volcano)It survives all three tests — thenreport it spatialworkflow.io
Three causes, three tests: none of them needs the seismicity to have changed.

Find where the catalogue begins, count above completeness

1. Find where the catalogue starts being a catalogue

Plot annual counts at a low magnitude and look for a step. A step is a network, an operator or a processing change, and everything before it is a different dataset.

2. Estimate the magnitude of completeness, per period

Completeness improves over time. The usable threshold for a trend analysis is the worst Mc across the whole period, not today's.

3. Count only above that threshold

Counts above a time-invariant Mc are comparable between decades. Counts above the API minimum are a record of instrumentation.

4. Decluster before counting

One M7 produces thousands of aftershocks. Whether a decade looks busy usually depends on whether a large mainshock fell in it.

5. Check the magnitude types

A catalogue that was duration magnitude in 1985 and local magnitude in 2015 has a systematic offset between the two eras, of a few tenths, which moves counts around a fixed threshold.

6. Say what period the conclusion covers

"Declustered, central California had between 592 and 757 M≥3 mainshocks per decade from the 1970s to the 2000s, and no more since" is a defensible sentence — it is what Example 2 measures. "Earthquakes are increasing" from raw counts is not.

Annual counts of magnitude 2 and above in central California from 1970 to 1979. M≥2 per year, central California01000200010519707419717419726719731,94419741,91119751,65719761,070197795419781,3841979regional network starts, 197467 events in 1973, 1,944 in 1974 — nothing happened to the faults that year spatialworkflow.io
The regional network comes online: 67 events in 1973, 1,944 in 1974 — a 29-fold jump in one year, with no change in the earth.

The step, declustering and what trend is left

Example 1 — the step, and whether it reaches the large events

print("M>=2 per year, 1970–1979:",
      [int(((eq.year == y) & (eq.mag >= 2)).sum()) for y in range(1970, 1980)])
for m in (2, 4, 5):
    before = eq[(eq.year < 1974) & (eq.mag >= m)]
    after = eq[(eq.year >= 1974) & (eq.mag >= m)]
    print(f"M>={m}: 1970–73 {len(before) / 4:6.1f}/yr    1974–2024 {len(after) / 51:6.1f}/yr")
M>=2 per year, 1970–1979: [105, 74, 74, 67, 1944, 1911, 1657, 1070, 954, 1384]
M>=2: 1970–73   80.0/yr    1974–2024 1798.9/yr
M>=4: 1970–73    5.8/yr    1974–2024   22.3/yr
M>=5: 1970–73    0.5/yr    1974–2024    2.0/yr

The step is not confined to small events: M≥4 is four times more frequent after 1974 and M≥5 four times as well, on small numbers. That rules out detection threshold as the whole story and points at the catalogue's own history — for this region, the events before the regional network are a sparser, differently sourced compilation. The practical consequence is the same either way: start the analysis after the step.

Example 2 — declustering, and what is left of the trend

import numpy as np

def decluster(df):
    """Gardner–Knopoff windows: an event inside another's space-time window is an aftershock."""
    df = df.sort_values("mag", ascending=False).reset_index(drop=True)
    m = df.mag.to_numpy()
    days = df.time.values.astype("datetime64[s]").astype(float) / 86400
    lat, lon = np.radians(df.lat.to_numpy()), np.radians(df.lon.to_numpy())
    radius = 10 ** (0.1238 * m + 0.983)                          # km
    window = np.where(m >= 6.5, 10 ** (0.032 * m + 2.7389), 10 ** (0.5409 * m - 0.547))
    keep = np.ones(len(df), bool)
    for i in range(len(df)):
        if not keep[i]:
            continue
        dist = 6371 * np.hypot(lat - lat[i], (lon - lon[i]) * np.cos(lat[i]))
        dt = days - days[i]
        keep &= ~((dist <= radius[i]) & (dt > 0) & (dt <= window[i]) & (m < m[i]))
    return df[keep]

above = eq[eq.mag >= 3.0]
main = decluster(above)
print(f"M>=3.0: {len(above):,} events, {len(main):,} mainshocks "
      f"({1 - len(main) / len(above):.0%} were aftershocks)")
for name, d in [("all events", above), ("mainshocks", main)]:
    s = d.groupby(d.time.dt.year // 10 * 10).size()
    print(f"  {name:11} " + "  ".join(f"{k}s {v:5,}" for k, v in s.items()))
M>=3.0: 12,217 events, 3,415 mainshocks (72% were aftershocks)
  all events  1970s 2,020  1980s 4,187  1990s 1,934  2000s 1,366  2010s 1,665  2020s 1,045
  mainshocks  1970s   757  1980s   658  1990s   663  2000s   592  2010s   475  2020s   270

The raw series has a 1980s peak of 4,187 against 1,366 in the 2000s — a threefold swing that looks like a trend if only the endpoints are quoted. Declustered, every decade from the 1970s to the 2000s sits between 592 and 757 mainshocks, and the apparent peak is gone: it was the aftershocks of the 1983 Coalinga and 1984 Morgan Hill sequences.

Example 3 — the magnitude types drift too

types = eq.groupby([eq.year // 10 * 10, "magtype"]).size().unstack(fill_value=0)
keep = types.sum().sort_values(ascending=False).head(4).index
print(types[keep].to_string())
print("\nshare of events by duration magnitude (md):")
print("  " + "  ".join(f"{d}s {row.md / row.sum():.0%}" for d, row in types.iterrows()))
magtype     md    ml    mc    mh
year                            
1970      7317  1197     0   648
1980     18312  3156  2716  1300
1990     12257  1362  4257   235
2000     11051  3251   574   123
2010      6518  7260     0    46
2020      3491  5695     0    16

share of events by duration magnitude (md):
  1970s 79%  1980s 72%  1990s 67%  2000s 72%  2010s 45%  2020s 37%

Duration magnitude carried three quarters of the catalogue in the 1980s and about a third by the 2020s, with local magnitude taking over. The two scales differ by a few tenths of a unit in a way that varies with magnitude, so a fixed threshold catches a slightly different population in each era.

Completeness is the whole game

Why the magnitude of completeness is the whole game

Below Mc the catalogue is missing events, and how many it misses depends on the network on that day. Any statistic computed below Mc — a count, a b-value, a rate — measures the network. Above it, the statistic measures the earth. The threshold for comparing periods is the worst Mc in any of them, which usually means throwing away most of the events.

Why aftershocks dominate raw counts

Aftershock sequences follow Omori's law: the rate decays as roughly 1/t, so a large earthquake adds events for years. Because sequence size grows with mainshock magnitude, one M7 can outnumber a decade of background seismicity, and raw decade counts end up reporting where the big ones happened to land.

Why the fix is not a better model

Declustering algorithms disagree with each other, completeness estimators disagree with each other, and both add their own assumptions. The defence is not to pick the best one but to say which was used and show the result with and without it — a conclusion that survives both is worth something, one that only appears in raw counts is not.

Earthquake counts per decade above magnitude 3, raw and after declustering. M≥3 per decade, raw and declustered1970s raw2,020 → 757 mainshocks1980s raw4,187 → 6581990s raw1,934 → 6632000s raw1,366 → 5922010s raw1,665 → 47502100420072% of events were aftershocksThe 1980s peak is Coalinga and Morgan Hill; the mainshock counts barely move spatialworkflow.io
Raw counts swing threefold between decades; the mainshock counts do not.

Induced seismicity, merges and swarms

  • Induced seismicity is a real rise in counts — Oklahoma after 2009 is not an artefact. Test it the same way and it survives.
  • Catalogue merges: global catalogues splice regional ones with different completeness by region and era.
  • Mc varies in space as well as time; a network is denser near cities.
  • Moment magnitude has been retroactively recomputed for many historic events, so old and new catalogue downloads disagree.
  • Foreshocks are removed by most declustering algorithms too, which matters if foreshocks are the subject.
  • Swarms without a clear mainshock are poorly handled by window-based declustering.
  • Partial decades wreck a decade plot: the 2020s row above covers five years, not ten.

FAQ

Are earthquakes becoming more frequent?

Globally, no — the apparent increase in catalogues is detection. Locally, induced seismicity can be a genuine increase.

What is the magnitude of completeness?

The lowest magnitude at which a catalogue records essentially every event. It improves over time and varies by region.

How do I count earthquakes fairly over decades?

Above a threshold the catalogue met throughout, after declustering, and say which declustering you used.

What is declustering?

Removing aftershocks and foreshocks so the remaining events approximate an independent series.

Why do old events have different magnitude types?

Practice changed: duration magnitude was standard for regional networks for decades, local and moment magnitude are now preferred.

Does the API magnitude filter solve this?

No. Filtering at M≥2 returns everything the network detected, which is exactly the quantity that changed.