Showcase — the 2023 Türkiye–Syria earthquake sequence (FDSN)¶
On 6 February 2023 an Mw 7.8 earthquake ruptured ~350 km of the East Anatolian Fault in southern Türkiye; nine hours later an Mw 7.5 event broke a second fault ~90 km to the north. Together they killed more than 55,000 people. Seismometer networks publish their event catalogs through the FDSN web-service standard, and earthlens returns them as a GeoDataFrame — enough to reproduce the core statistics seismologists use to characterise a sequence.
What this notebook does¶
- Query the FDSN service for every M ≥ 4 event in the region for Feb–Mar 2023.
- Map the epicentres — they trace the ruptured faults.
- Gutenberg–Richter: the magnitude–frequency distribution and its b-value.
- Omori–Utsu: the aftershock decay rate and its p-value.
- Energy: cumulative seismic moment and event count through time.
- Depth distribution of the sequence.
Live USGS query, no credentials; offline it skips gracefully.
Setup¶
Imports and a bit of noise-suppression. earthlens provides the unified EarthLens entry point; the FDSN backend returns events as a GeoDataFrame, so the rest is plain pandas / numpy / matplotlib.
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from loguru import logger
from earthlens import EarthLens
logger.disable("earthlens")
logger.disable("pyramids")
1 · Query the FDSN service¶
We build the request first — USGS via FDSN, a Feb–Mar 2023 window, a bounding box over the rupture zone, and a magnitude floor of 4.0 (low enough to give enough events for robust statistics). load() returns a GeoDataFrame, one row per event. The whole thing is wrapped in a try so the notebook skips gracefully when run offline.
quakes = None
try:
fdsn = EarthLens(
data_source="fdsn",
variables=["USGS"],
start="2023-02-06",
end="2023-03-31",
aoi=[35.0, 35.0, 39.5, 39.5],
min_magnitude=4.0,
)
quakes = fdsn.load()
quakes["time"] = pd.to_datetime(quakes["time"])
print(f"{len(quakes)} events M>=4.0; largest M{quakes.magnitude.max():.1f}")
except Exception as exc: # noqa: BLE001
print("live FDSN query skipped (offline?):", exc)
2026-06-13 18:20:15 | INFO | pyramids.base.config | Logging is configured.
2026-06-13 18:20:19 | INFO | pyogrio._io | Created 495 records
495 events M>=4.0; largest M7.8
2 · Epicentre map¶
The dots align along the two ruptured fault strands; marker size grows with magnitude and the two M7+ mainshocks are annotated.
if quakes is not None and len(quakes):
q = quakes.sort_values("magnitude")
fig, ax = plt.subplots(figsize=(8, 7))
sc = ax.scatter(
q.longitude,
q.latitude,
c=q.magnitude,
cmap="hot_r",
s=2**q.magnitude,
edgecolor="k",
linewidth=0.3,
alpha=0.8,
)
for _, r in q[q.magnitude >= 7.0].iterrows():
ax.annotate(f"M{r.magnitude:.1f}", (r.longitude, r.latitude), fontweight="bold")
ax.set(
xlabel="longitude",
ylabel="latitude",
title=f"Türkiye–Syria, Feb–Mar 2023 — {len(q)} events M>=4.0",
)
fig.colorbar(sc, ax=ax, label="magnitude")
plt.show()
3 · Gutenberg–Richter relation and b-value¶
Earthquake magnitudes follow log10 N(≥M) = a − bM: each unit of magnitude is ~10× rarer. The b-value (slope) is typically near 1; lower values indicate relatively more large events / higher stress. We estimate it with the maximum-likelihood (Aki) estimator above the completeness magnitude Mc.
if quakes is not None and len(quakes):
m = quakes.magnitude.values
Mc = 4.0 # completeness ~ the query floor
mags = np.sort(m[m >= Mc])
N = np.arange(len(mags), 0, -1) # cumulative count N(>=M)
b = np.log10(np.e) / (mags.mean() - (Mc - 0.05)) # Aki (1965) MLE
a = np.log10(len(mags)) + b * Mc
fig, ax = plt.subplots(figsize=(7, 5))
ax.semilogy(mags, N, "o", ms=4, label="observed")
mm = np.linspace(Mc, mags.max(), 50)
ax.semilogy(mm, 10 ** (a - b * mm), "r-", label=f"GR fit, b={b:.2f}")
ax.set(
xlabel="magnitude M", ylabel="N(≥ M)", title=f"Gutenberg–Richter — b = {b:.2f}"
)
ax.legend()
plt.show()
print(f"b-value = {b:.2f} (≈1 is typical; the sequence is rich in aftershocks)")
b-value = 0.88 (≈1 is typical; the sequence is rich in aftershocks)
4 · Omori–Utsu aftershock decay¶
Aftershock rate falls off with time after the mainshock as n(t) = K/(t+c)^p. Plotting daily rate against days-since-mainshock on log–log axes makes the decay a straight line whose slope is the p-value (typically ~1).
if quakes is not None and len(quakes):
t0 = quakes.loc[quakes.magnitude.idxmax(), "time"]
days = (quakes["time"] - t0).dt.total_seconds() / 86400
days = days[days > 0]
edges = np.arange(0, days.max() + 1, 1)
rate, _ = np.histogram(days, bins=edges)
centre = edges[:-1] + 0.5
ok = rate > 0
p, logK = np.polyfit(np.log10(centre[ok]), np.log10(rate[ok]), 1)
fig, ax = plt.subplots(figsize=(7, 5))
ax.loglog(centre[ok], rate[ok], "o", ms=4, label="events/day")
ax.loglog(centre[ok], 10**logK * centre[ok] ** p, "r-", label=f"Omori p={-p:.2f}")
ax.set(
xlabel="days since mainshock",
ylabel="aftershocks per day",
title="Omori–Utsu aftershock decay",
)
ax.legend()
plt.show()
print(f"Omori p-value = {-p:.2f} (aftershock rate roughly halves as time doubles)")
Omori p-value = 1.00 (aftershock rate roughly halves as time doubles)
5 · Cumulative seismic moment and count¶
Seismic moment M0 = 10^(1.5·Mw + 9.1) N·m measures the energy of each event; it is dominated by the very largest shocks. The step-like cumulative-moment curve shows how the two mainshocks release almost all the energy, while the event count keeps climbing with aftershocks.
if quakes is not None and len(quakes):
q = quakes.sort_values("time")
M0 = 10 ** (1.5 * q.magnitude + 9.1) # scalar moment, N·m
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(q.time, np.cumsum(M0), color="tab:red", label="cumulative moment (N·m)")
ax.set_ylabel("cumulative moment (N·m)", color="tab:red")
ax2 = ax.twinx()
ax2.plot(q.time, np.arange(1, len(q) + 1), color="tab:blue")
ax2.set_ylabel("cumulative event count", color="tab:blue")
ax.set_title("Energy released vs events recorded")
fig.autofmt_xdate()
plt.show()
6 · Depth distribution¶
The sequence is shallow crustal — most events in the top ~20 km, consistent with continental strike-slip faulting.
if quakes is not None and len(quakes):
fig, ax = plt.subplots(figsize=(7, 3.4))
ax.hist(quakes.depth_km.dropna(), bins=30, color="slategray")
ax.set(xlabel="depth (km)", ylabel="events", title="Hypocentre depth distribution")
plt.show()
Recap¶
From a single FDSN query we reproduced the textbook description of an earthquake sequence: the Gutenberg–Richter b-value, the Omori aftershock decay, the cumulative-moment jump at the two mainshocks, and a shallow-crustal depth distribution.
Try it yourself¶
- Lower
min_magnitudefor a denser catalogue (and a more stable b-value). - Change the region/dates to any other sequence, or query a different agency (
variables=['EMSC']). - See the FDSN backend reference.