Earthquake Statistics — Exercises¶

These exercises reinforce the three key laws from the lecture:

1. Gutenberg-Richter: $\log_{10} N = a - bM$
2. Omori: $n(t) = K / (t + c)^p$
3. Poisson recurrence: $P = 1 - e^{-\lambda T}$

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

import numpy as np import matplotlib.pyplot as plt import pandas as pd from scipy.optimize import curve_fit

---

Setup¶

Run the two cells below to load the data. (Same catalogs as the lecture.)

In [ ]:

# Load pre-downloaded catalogs — just run this cell
global_cat = pd.read_csv("https://github.com/AI4EPS/EPS130_Seismology/releases/download/earthquake-statistics-data/global_M4_2015_2025.csv")
global_cat["time"] = pd.to_datetime(global_cat["time"])

ridgecrest = pd.read_csv("https://github.com/AI4EPS/EPS130_Seismology/releases/download/earthquake-statistics-data/ridgecrest_2019.csv")
ridgecrest["time"] = pd.to_datetime(ridgecrest["time"])
ridgecrest = ridgecrest.sort_values("time").reset_index(drop=True)

mainshock = ridgecrest.loc[ridgecrest["mag"].idxmax()]
mainshock_time = mainshock["time"]

print(f"Global: {len(global_cat):,} earthquakes (M4+)")
print(f"Ridgecrest: {len(ridgecrest):,} earthquakes (all magnitudes)")
print(f"Mainshock: M{mainshock['mag']:.1f}")

# Load pre-downloaded catalogs — just run this cell global_cat = pd.read_csv("https://github.com/AI4EPS/EPS130_Seismology/releases/download/earthquake-statistics-data/global_M4_2015_2025.csv") global_cat["time"] = pd.to_datetime(global_cat["time"]) ridgecrest = pd.read_csv("https://github.com/AI4EPS/EPS130_Seismology/releases/download/earthquake-statistics-data/ridgecrest_2019.csv") ridgecrest["time"] = pd.to_datetime(ridgecrest["time"]) ridgecrest = ridgecrest.sort_values("time").reset_index(drop=True) mainshock = ridgecrest.loc[ridgecrest["mag"].idxmax()] mainshock_time = mainshock["time"] print(f"Global: {len(global_cat):,} earthquakes (M4+)") print(f"Ridgecrest: {len(ridgecrest):,} earthquakes (all magnitudes)") print(f"Mainshock: M{mainshock['mag']:.1f}")

---

Exercise 1: The Gutenberg-Richter Law¶

The Gutenberg-Richter law says:

$$\log_{10} N(M) = a - bM \quad \Longleftrightarrow \quad N(M) = 10^{\,a - bM}$$

where $N(M)$ is the number of earthquakes with magnitude $\geq M$.

1a. Explore the data¶

Run the cell below to plot the cumulative G-R distribution and magnitude histogram for the Ridgecrest sequence.

In [ ]:

# Just run this cell
mag_bins = np.arange(-1, 7.5, 0.1)
cum_count = np.array([np.sum(ridgecrest["mag"] >= m) for m in mag_bins])

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

axes[0].semilogy(mag_bins, cum_count, "o-", color="steelblue", markersize=3)
axes[0].set_xlabel("Magnitude (M)")
axes[0].set_ylabel("Cumulative count N(M)")
axes[0].set_title("Cumulative G-R Plot")
axes[0].grid(True, alpha=0.3)

axes[1].hist(ridgecrest["mag"], bins=np.arange(-1.0, 7.5, 0.1),
             color="steelblue", edgecolor="white", linewidth=0.5)
axes[1].set_yscale("log")
axes[1].set_xlabel("Magnitude")
axes[1].set_ylabel("Count")
axes[1].set_title("Magnitude Histogram")

plt.tight_layout()
plt.show()

# Just run this cell mag_bins = np.arange(-1, 7.5, 0.1) cum_count = np.array([np.sum(ridgecrest["mag"] >= m) for m in mag_bins]) fig, axes = plt.subplots(1, 2, figsize=(13, 5)) axes[0].semilogy(mag_bins, cum_count, "o-", color="steelblue", markersize=3) axes[0].set_xlabel("Magnitude (M)") axes[0].set_ylabel("Cumulative count N(M)") axes[0].set_title("Cumulative G-R Plot") axes[0].grid(True, alpha=0.3) axes[1].hist(ridgecrest["mag"], bins=np.arange(-1.0, 7.5, 0.1), color="steelblue", edgecolor="white", linewidth=0.5) axes[1].set_yscale("log") axes[1].set_xlabel("Magnitude") axes[1].set_ylabel("Count") axes[1].set_title("Magnitude Histogram") plt.tight_layout() plt.show()

1b. Estimate $M_c$¶

Look at the histogram above. At what magnitude does the count peak? This is your estimate of $M_c$ (magnitude of completeness) — below this, the catalog is missing events.

In [ ]:

# YOUR CODE HERE: set Mc based on the histogram peak
Mc = ???

# YOUR CODE HERE: set Mc based on the histogram peak Mc = ???

1c. Write the G-R formula¶

Complete the function below. Given $a$, $b$, and $M$, it should return $N(M) = 10^{\,a - bM}$.

In [ ]:

def gutenberg_richter(M, a, b):
    """Predicted number of earthquakes with magnitude >= M."""
    return ???  # Hint: 10 ** (...)

def gutenberg_richter(M, a, b): """Predicted number of earthquakes with magnitude >= M.""" return ??? # Hint: 10 ** (...)

1d. Fit $a$ and $b$¶

The G-R law says $\log_{10} N = a - bM$ — a straight line on a semi-log plot. We fit this line using np.polyfit, which returns [slope, intercept]. Since slope $= -b$, we have $b = -$slope.

In [ ]:

# Fit a line to log10(N) vs M, using only data above Mc
mask = (mag_bins >= Mc) & (cum_count > 0)
slope, intercept = np.polyfit(mag_bins[mask], np.log10(cum_count[mask]), 1)

# YOUR CODE HERE: extract b and a from the fit
# Remember: slope = -b
b = ???
a = ???

print(f"a = {a:.2f}")
print(f"b = {b:.2f}")

# Fit a line to log10(N) vs M, using only data above Mc mask = (mag_bins >= Mc) & (cum_count > 0) slope, intercept = np.polyfit(mag_bins[mask], np.log10(cum_count[mask]), 1) # YOUR CODE HERE: extract b and a from the fit # Remember: slope = -b b = ??? a = ??? print(f"a = {a:.2f}") print(f"b = {b:.2f}")

1e. Predict and verify¶

Use your gutenberg_richter function and the fitted $a$, $b$ to predict how many M3+, M4+, and M5+ earthquakes occurred. Then compare with the actual catalog.

In [ ]:

# YOUR CODE HERE: use gutenberg_richter(M, a, b) to predict, compare with actual
for M in [3, 4, 5]:
    predicted = ???
    actual = np.sum(ridgecrest["mag"] >= M)
    print(f"M{M}+:  predicted = {predicted:.0f},  actual = {actual}")

# YOUR CODE HERE: use gutenberg_richter(M, a, b) to predict, compare with actual for M in [3, 4, 5]: predicted = ??? actual = np.sum(ridgecrest["mag"] >= M) print(f"M{M}+: predicted = {predicted:.0f}, actual = {actual}")

1f. Plot your fit¶

Run the cell below to see your G-R fit overlaid on the data.

In [ ]:

# Just run this cell
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogy(mag_bins, cum_count, "o", color="steelblue", markersize=4, label="Data")
M_plot = np.linspace(Mc, 8, 100)
ax.semilogy(M_plot, gutenberg_richter(M_plot, a, b), "r-", linewidth=2,
            label=f"G-R fit: b = {b:.2f}")
ax.axvline(Mc, color="gray", linestyle="--", alpha=0.7, label=f"$M_c$ = {Mc}")
ax.set_xlabel("Magnitude (M)")
ax.set_ylabel("Cumulative count N(M)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

# Just run this cell fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogy(mag_bins, cum_count, "o", color="steelblue", markersize=4, label="Data") M_plot = np.linspace(Mc, 8, 100) ax.semilogy(M_plot, gutenberg_richter(M_plot, a, b), "r-", linewidth=2, label=f"G-R fit: b = {b:.2f}") ax.axvline(Mc, color="gray", linestyle="--", alpha=0.7, label=f"$M_c$ = {Mc}") ax.set_xlabel("Magnitude (M)") ax.set_ylabel("Cumulative count N(M)") ax.legend() ax.grid(True, alpha=0.3) plt.show()

---

Exercise 2: Omori's Law¶

The modified Omori law describes how the aftershock rate decays with time:

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

2a. Prepare the aftershock rate data¶

Run the cell below to select M2+ aftershocks, bin them in time, and plot the rate.

In [ ]:

# Just run this cell — it selects aftershocks and computes the rate
aftershocks = ridgecrest[(ridgecrest["time"] > mainshock_time) & (ridgecrest["mag"] >= 2.0)].copy()
aftershocks["dt_days"] = (aftershocks["time"] - mainshock_time).dt.total_seconds() / 86400
print(f"M2+ aftershocks: {len(aftershocks):,}")

bin_edges = np.logspace(np.log10(0.01), np.log10(180), 40)
counts_omori, _ = np.histogram(aftershocks["dt_days"], bins=bin_edges)
bin_widths = np.diff(bin_edges)
bin_centers = np.sqrt(bin_edges[:-1] * bin_edges[1:])
rates = counts_omori / bin_widths
mask = counts_omori > 0
bin_centers = bin_centers[mask]
rates = rates[mask]

fig, ax = plt.subplots(figsize=(8, 5))
ax.loglog(bin_centers, rates, "o", color="steelblue", markersize=5)
ax.set_xlabel("Time after mainshock (days)")
ax.set_ylabel("Aftershock rate (events/day)")
ax.grid(True, alpha=0.3, which="both")
plt.show()

# Just run this cell — it selects aftershocks and computes the rate aftershocks = ridgecrest[(ridgecrest["time"] > mainshock_time) & (ridgecrest["mag"] >= 2.0)].copy() aftershocks["dt_days"] = (aftershocks["time"] - mainshock_time).dt.total_seconds() / 86400 print(f"M2+ aftershocks: {len(aftershocks):,}") bin_edges = np.logspace(np.log10(0.01), np.log10(180), 40) counts_omori, _ = np.histogram(aftershocks["dt_days"], bins=bin_edges) bin_widths = np.diff(bin_edges) bin_centers = np.sqrt(bin_edges[:-1] * bin_edges[1:]) rates = counts_omori / bin_widths mask = counts_omori > 0 bin_centers = bin_centers[mask] rates = rates[mask] fig, ax = plt.subplots(figsize=(8, 5)) ax.loglog(bin_centers, rates, "o", color="steelblue", markersize=5) ax.set_xlabel("Time after mainshock (days)") ax.set_ylabel("Aftershock rate (events/day)") ax.grid(True, alpha=0.3, which="both") plt.show()

2b. Write the Omori law¶

Complete the function below: $n(t) = K / (t + c)^p$

In [ ]:

def omori_law(t, K, c, p):
    """Aftershock rate at time t days after mainshock."""
    return ???  # Hint: K / (...)

def omori_law(t, K, c, p): """Aftershock rate at time t days after mainshock.""" return ??? # Hint: K / (...)

2c. Fit and plot¶

Run the cell below to fit your omori_law function to the data.

In [ ]:

# Just run this cell — it fits your omori_law to the rate data
popt, _ = curve_fit(omori_law, bin_centers, rates, p0=[100, 0.1, 1.0], maxfev=10000)
K_fit, c_fit, p_fit = popt
print(f"K = {K_fit:.0f},  c = {c_fit:.3f} days,  p = {p_fit:.2f}")

fig, ax = plt.subplots(figsize=(8, 5))
ax.loglog(bin_centers, rates, "o", color="steelblue", markersize=5, label="Data")
t_plot = np.logspace(-2, np.log10(200), 200)
ax.loglog(t_plot, omori_law(t_plot, *popt), "r-", linewidth=2, label=f"Omori fit: p = {p_fit:.2f}")
ax.set_xlabel("Time after mainshock (days)")
ax.set_ylabel("Aftershock rate (events/day)")
ax.legend()
ax.grid(True, alpha=0.3, which="both")
plt.show()

# Just run this cell — it fits your omori_law to the rate data popt, _ = curve_fit(omori_law, bin_centers, rates, p0=[100, 0.1, 1.0], maxfev=10000) K_fit, c_fit, p_fit = popt print(f"K = {K_fit:.0f}, c = {c_fit:.3f} days, p = {p_fit:.2f}") fig, ax = plt.subplots(figsize=(8, 5)) ax.loglog(bin_centers, rates, "o", color="steelblue", markersize=5, label="Data") t_plot = np.logspace(-2, np.log10(200), 200) ax.loglog(t_plot, omori_law(t_plot, *popt), "r-", linewidth=2, label=f"Omori fit: p = {p_fit:.2f}") ax.set_xlabel("Time after mainshock (days)") ax.set_ylabel("Aftershock rate (events/day)") ax.legend() ax.grid(True, alpha=0.3, which="both") plt.show()

2d. Apply the Omori law¶

Using your fitted parameters, compute the aftershock rate at 1 day, 10 days, and 100 days. What pattern do you see?

In [ ]:

# YOUR CODE HERE: compute rate at t = 1, 10, 100 days
for t in [1, 10, 100]:
    rate = ???
    print(f"Day {t:3d}: {rate:.0f} events/day")

# YOUR CODE HERE: compute rate at t = 1, 10, 100 days for t in [1, 10, 100]: rate = ??? print(f"Day {t:3d}: {rate:.0f} events/day")

---

Exercise 3: Earthquake Hazard¶

If large earthquakes occur randomly with an average rate $\lambda = 1/\tau$ per year, the probability of at least one event in $T$ years is:

$$P = 1 - e^{-\lambda T}$$

The Hayward Fault runs through the UC Berkeley campus. M6.8+ earthquakes occur roughly every $\tau = 150$ years. The last one was in 1868.

3a. Write the Poisson formula¶

Complete the function below.

In [ ]:

def earthquake_probability(tau, T):
    """Probability of at least one event in T years, given recurrence interval tau."""
    lam = 1 / tau
    return ???  # Hint: 1 - np.exp(...)

def earthquake_probability(tau, T): """Probability of at least one event in T years, given recurrence interval tau.""" lam = 1 / tau return ??? # Hint: 1 - np.exp(...)

3b. Apply to the Hayward Fault¶

What is the probability of a M6.8+ earthquake in the next 30 years? The next 100 years?

In [ ]:

tau = 150  # years

# YOUR CODE HERE: compute probability for T = 30, 50, 100, 150 years
for T in [30, 50, 100, 150]:
    P = ???
    print(f"Next {T:3d} years: {P*100:.1f}%")

tau = 150 # years # YOUR CODE HERE: compute probability for T = 30, 50, 100, 150 years for T in [30, 50, 100, 150]: P = ??? print(f"Next {T:3d} years: {P*100:.1f}%")

3c. How long until 50% probability?¶

Solve $1 - e^{-\lambda T} = 0.5$ for $T$. (Work it out on paper first, then compute.)

In [ ]:

# YOUR CODE HERE: solve for T when P = 0.5
T_50 = ???
print(f"50% probability after {T_50:.0f} years")

# YOUR CODE HERE: solve for T when P = 0.5 T_50 = ??? print(f"50% probability after {T_50:.0f} years")

3d. G-R extrapolation¶

The global G-R fit from the lecture gives $a \approx 9.80$ and $b \approx 1.11$ (for a 10-year catalog). Use your gutenberg_richter function to estimate how many M7+ and M8+ earthquakes occur per year globally.

In [ ]:

a_gl = 9.80
b_gl = 1.11
years = 10  # the catalog covers 10 years

# YOUR CODE HERE: compute M5+, M6+, M7+, M8+ per year
for M in [5, 6, 7, 8]:
    N_per_year = ???
    print(f"M{M}+ per year: {N_per_year:.0f}")

a_gl = 9.80 b_gl = 1.11 years = 10 # the catalog covers 10 years # YOUR CODE HERE: compute M5+, M6+, M7+, M8+ per year for M in [5, 6, 7, 8]: N_per_year = ??? print(f"M{M}+ per year: {N_per_year:.0f}")