Earthquake Magnitude — Exercises¶

These exercises use the 2025 Berkeley M4.3 earthquake from the lecture notebook. You'll compute $M_L$ yourself and answer questions about earthquake magnitude and intensity.

In [ ]:

Copied!





!pip install obspy cartopy -q

import numpy as np
import matplotlib.pyplot as plt
import requests
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal.invsim import corn_freq_2_paz
from obspy.geodetics import gps2dist_azimuth
!pip install obspy cartopy -q

import numpy as np
import matplotlib.pyplot as plt
import requests
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal.invsim import corn_freq_2_paz
from obspy.geodetics import gps2dist_azimuth

Setup: Download Data¶

Run the cells below to download the waveform data (same as the lecture).

In [ ]:

Copied!





client = Client("NCEDC")

origin_time = UTCDateTime("2025-09-22T09:56:12.52")
eq_lat, eq_lon, eq_depth = 37.8630, -122.2540, 7.6

station_list = ["MCCM", "FARB", "MHC", "RUSS", "MNRC",
                "HOPS", "CMB", "BUCR", "WINE", "HELL"]

t1 = origin_time - 30
t2 = origin_time + 180

streams = {}
for sta in station_list:
    try:
        st = client.get_waveforms("BK", sta, "*", "HH?", t1, t2,
                                  attach_response=True)
        streams[sta] = st
    except:
        pass

print(f"Downloaded data from {len(streams)} stations")
client = Client("NCEDC")

origin_time = UTCDateTime("2025-09-22T09:56:12.52")
eq_lat, eq_lon, eq_depth = 37.8630, -122.2540, 7.6

station_list = ["MCCM", "FARB", "MHC", "RUSS", "MNRC",
                "HOPS", "CMB", "BUCR", "WINE", "HELL"]

t1 = origin_time - 30
t2 = origin_time + 180

streams = {}
for sta in station_list:
    try:
        st = client.get_waveforms("BK", sta, "*", "HH?", t1, t2,
                                  attach_response=True)
        streams[sta] = st
    except:
        pass

print(f"Downloaded data from {len(streams)} stations")

In [ ]:

Copied!





# Wood-Anderson response (same as lecture)
wa_paz = corn_freq_2_paz(1.0, damp=0.8)
wa_paz["sensitivity"] = 2080
wa_paz["zeros"] = [0j, 0j]
# Wood-Anderson response (same as lecture)
wa_paz = corn_freq_2_paz(1.0, damp=0.8)
wa_paz["sensitivity"] = 2080
wa_paz["zeros"] = [0j, 0j]

Exercise 1: Compute $M_L$ for One Station¶

In the lecture, we computed $M_L$ using code that looped over all stations automatically. Now, do it step by step for a single station (BK.MCCM) to make sure you understand each part of the pipeline.

Recall the pipeline:

$$\text{Raw counts} \xrightarrow{\text{remove response}} \text{Displacement} \xrightarrow{\text{simulate WA}} \text{WA trace} \xrightarrow{\text{half p2p in S window}} A \xrightarrow{\text{formula}} M_L$$

Bakun-Joyner formula: $M_L = \log_{10} A + \log_{10} R + 0.00301\,R + 0.70$, where $A$ is in mm and $R$ is hypocentral distance in km.

1a. Remove instrument response and simulate Wood-Anderson¶

Fill in the ??? to process the waveform.

In [ ]:

Copied!





sta = "MCCM"
st = streams[sta].copy()

# Step 1: Remove instrument response to get displacement
st.detrend("demean")
st.detrend("linear")
st.taper(0.05)
st.remove_response(output=???)  # What output type gives displacement?

# Step 2: Simulate Wood-Anderson
for tr in st:
    tr.stats.response = None
st.simulate(paz_simulate=???)   # What PAZ should we simulate?

# Pick the North component
tr = st.select(component="N")[0]
times = tr.times(reftime=origin_time)

plt.figure(figsize=(10, 3))
plt.plot(times, tr.data * 1e3, "r", lw=0.5)
plt.xlabel("Time after origin (s)")
plt.ylabel("WA amplitude (mm)")
plt.title(f"BK.{sta}.HHN — Wood-Anderson Simulation")
plt.xlim(-10, 150)
plt.tight_layout()
plt.show()
sta = "MCCM"
st = streams[sta].copy()

# Step 1: Remove instrument response to get displacement
st.detrend("demean")
st.detrend("linear")
st.taper(0.05)
st.remove_response(output=???)  # What output type gives displacement?

# Step 2: Simulate Wood-Anderson
for tr in st:
    tr.stats.response = None
st.simulate(paz_simulate=???)   # What PAZ should we simulate?

# Pick the North component
tr = st.select(component="N")[0]
times = tr.times(reftime=origin_time)

plt.figure(figsize=(10, 3))
plt.plot(times, tr.data * 1e3, "r", lw=0.5)
plt.xlabel("Time after origin (s)")
plt.ylabel("WA amplitude (mm)")
plt.title(f"BK.{sta}.HHN — Wood-Anderson Simulation")
plt.xlim(-10, 150)
plt.tight_layout()
plt.show()

1b. Measure the amplitude¶

Compute the epicentral distance, estimate the S-wave arrival, and measure the half peak-to-trough amplitude in the S-wave window.

In [ ]:

Copied!





# Get station location and compute distance
inv = client.get_stations(network="BK", station=sta, level="station")
sta_lat = inv[0][0].latitude
sta_lon = inv[0][0].longitude
dist_m, _, _ = gps2dist_azimuth(eq_lat, eq_lon, sta_lat, sta_lon)
dist_km = dist_m / 1000.0

# Hypocentral distance
R = ???  # How do you compute hypocentral distance from epicentral distance and depth?

# S-wave arrival time (S-wave velocity ≈ 3.5 km/s)
s_arrival = ???  # distance / velocity

# Extract data in the S-wave window
s_mask = (times >= s_arrival - 2) & (times <= s_arrival + 30)
s_data = tr.data[s_mask]

# Half peak-to-trough amplitude in mm
amp_mm = ???  # (max - min) / 2, converted to mm

print(f"Station: BK.{sta}")
print(f"Epicentral distance: {dist_km:.1f} km")
print(f"Hypocentral distance: {R:.1f} km")
print(f"Half peak-to-trough amplitude: {amp_mm:.2f} mm")
# Get station location and compute distance
inv = client.get_stations(network="BK", station=sta, level="station")
sta_lat = inv[0][0].latitude
sta_lon = inv[0][0].longitude
dist_m, _, _ = gps2dist_azimuth(eq_lat, eq_lon, sta_lat, sta_lon)
dist_km = dist_m / 1000.0

# Hypocentral distance
R = ???  # How do you compute hypocentral distance from epicentral distance and depth?

# S-wave arrival time (S-wave velocity ≈ 3.5 km/s)
s_arrival = ???  # distance / velocity

# Extract data in the S-wave window
s_mask = (times >= s_arrival - 2) & (times <= s_arrival + 30)
s_data = tr.data[s_mask]

# Half peak-to-trough amplitude in mm
amp_mm = ???  # (max - min) / 2, converted to mm

print(f"Station: BK.{sta}")
print(f"Epicentral distance: {dist_km:.1f} km")
print(f"Hypocentral distance: {R:.1f} km")
print(f"Half peak-to-trough amplitude: {amp_mm:.2f} mm")

1c. Compute $M_L$¶

Apply the Bakun-Joyner formula.

In [ ]:

Copied!

# Bakun-Joyner (1984): ML = log10(A) + log10(R) + 0.00301*R + 0.70
ml = ???

print(f"Our ML = {ml:.2f}")
print(f"NCEDC catalog ML = 4.55")
# Bakun-Joyner (1984): ML = log10(A) + log10(R) + 0.00301*R + 0.70
ml = ???

print(f"Our ML = {ml:.2f}")
print(f"NCEDC catalog ML = 4.55")

Exercise 2: Why Does $M_L$ Vary Between Stations?¶

In the lecture, we saw that different stations give slightly different $M_L$ values for the same earthquake. Let's quantify this.

In [ ]:

Copied!





# This code computes ML for all stations (from the lecture). Just run it.
station_coords = {}
for s in streams:
    inv = client.get_stations(network="BK", station=s, level="station")
    d, _, _ = gps2dist_azimuth(eq_lat, eq_lon, inv[0][0].latitude, inv[0][0].longitude)
    dk = d / 1000.0
    station_coords[s] = {"dist_km": dk, "hypo_km": np.sqrt(dk**2 + eq_depth**2),
                         "s_arrival": dk / 3.5}

streams_disp = {}
for s, st in streams.items():
    st_d = st.copy()
    st_d.detrend("demean"); st_d.detrend("linear"); st_d.taper(0.05)
    st_d.remove_response(output="DISP")
    streams_disp[s] = st_d

streams_wa = {}
for s, st in streams_disp.items():
    st_w = st.copy()
    for tr in st_w:
        tr.stats.response = None
    st_w.simulate(paz_simulate=wa_paz)
    streams_wa[s] = st_w

ml_results = []
for s, st in streams_wa.items():
    c = station_coords[s]
    for comp in ["N", "E"]:
        trs = st.select(component=comp)
        if not trs: continue
        tr = trs[0]
        t = tr.times(reftime=origin_time)
        mask = (t >= c["s_arrival"] - 2) & (t <= c["s_arrival"] + 30)
        sd = tr.data[mask]
        if len(sd) == 0: continue
        a = (sd.max() - sd.min()) / 2 * 1e3
        m = np.log10(a) + np.log10(c["hypo_km"]) + 0.00301 * c["hypo_km"] + 0.70
        ml_results.append({"station": s, "comp": comp, "dist_km": c["dist_km"], "ml": m})

ml_all = np.array([r["ml"] for r in ml_results])
print(f"ML values across all station-components: median = {np.median(ml_all):.2f}, std = {np.std(ml_all):.2f}")
# This code computes ML for all stations (from the lecture). Just run it.
station_coords = {}
for s in streams:
    inv = client.get_stations(network="BK", station=s, level="station")
    d, _, _ = gps2dist_azimuth(eq_lat, eq_lon, inv[0][0].latitude, inv[0][0].longitude)
    dk = d / 1000.0
    station_coords[s] = {"dist_km": dk, "hypo_km": np.sqrt(dk**2 + eq_depth**2),
                         "s_arrival": dk / 3.5}

streams_disp = {}
for s, st in streams.items():
    st_d = st.copy()
    st_d.detrend("demean"); st_d.detrend("linear"); st_d.taper(0.05)
    st_d.remove_response(output="DISP")
    streams_disp[s] = st_d

streams_wa = {}
for s, st in streams_disp.items():
    st_w = st.copy()
    for tr in st_w:
        tr.stats.response = None
    st_w.simulate(paz_simulate=wa_paz)
    streams_wa[s] = st_w

ml_results = []
for s, st in streams_wa.items():
    c = station_coords[s]
    for comp in ["N", "E"]:
        trs = st.select(component=comp)
        if not trs: continue
        tr = trs[0]
        t = tr.times(reftime=origin_time)
        mask = (t >= c["s_arrival"] - 2) & (t <= c["s_arrival"] + 30)
        sd = tr.data[mask]
        if len(sd) == 0: continue
        a = (sd.max() - sd.min()) / 2 * 1e3
        m = np.log10(a) + np.log10(c["hypo_km"]) + 0.00301 * c["hypo_km"] + 0.70
        ml_results.append({"station": s, "comp": comp, "dist_km": c["dist_km"], "ml": m})

ml_all = np.array([r["ml"] for r in ml_results])
print(f"ML values across all station-components: median = {np.median(ml_all):.2f}, std = {np.std(ml_all):.2f}")

Question 2a: The standard deviation of $M_L$ across stations is about 0.17. Give three physical reasons why different stations give different $M_L$ for the same earthquake.

Your answer here:

Question 2b: NCEDC applies "station corrections" — a fixed offset added to each station's $M_L$ to reduce scatter. Why would one station consistently give higher $M_L$ than another?

Your answer here:

Exercise 3: Magnitude Saturation¶

The NCEDC reported three magnitudes for this earthquake:

$M_w$ = 4.29 (52 stations)
$M_L$ = 4.55 (167 stations)
$M_d$ = 4.63 (191 stations)

Question 3: For very large earthquakes ($M > 7$), $M_L$ and $m_b$ saturate — they stop growing even as the earthquake gets bigger. Explain why, using what you know about the source spectrum and corner frequency.

Your answer here:

Exercise 4: Intensity vs Magnitude¶

Run the cell below to fetch and plot the DYFI intensity data.

In [ ]:

Copied!





import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import BoundaryNorm

# Fetch DYFI data
detail = requests.get("https://earthquake.usgs.gov/earthquakes/feed/v1.0/detail/nc75240492.geojson").json()
dyfi_products = detail["properties"]["products"].get("dyfi", [])
dyfi_url = dyfi_products[0]["contents"]["dyfi_zip.geojson"]["url"] if dyfi_products else None

lats, lons, intensities, nresps = [], [], [], []
if dyfi_url:
    for f in requests.get(dyfi_url).json()["features"]:
        if f["geometry"]["type"] == "Point":
            lons.append(f["geometry"]["coordinates"][0])
            lats.append(f["geometry"]["coordinates"][1])
            intensities.append(f["properties"].get("cdi", 0))
            nresps.append(f["properties"].get("nresp", 0))
lats, lons, intensities = np.array(lats), np.array(lons), np.array(intensities)
nresps = np.array(nresps)

# Intensity vs distance
distances_km = np.array([gps2dist_azimuth(eq_lat, eq_lon, la, lo)[0] / 1000
                         for la, lo in zip(lats, lons)])

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(distances_km, intensities, alpha=0.4, s=20, c="steelblue", edgecolors="none")

dist_bins = np.logspace(np.log10(3), np.log10(800), 15)
for i in range(len(dist_bins) - 1):
    mask = (distances_km >= dist_bins[i]) & (distances_km < dist_bins[i+1])
    if mask.sum() > 3:
        ax.plot(np.sqrt(dist_bins[i] * dist_bins[i+1]), np.median(intensities[mask]),
                "ro", markersize=7, zorder=4)

ax.set_xscale("log")
ax.set_xlabel("Distance from epicenter (km)")
ax.set_ylabel("DYFI Intensity (MMI)")
ax.set_title("Intensity vs Distance")
ax.grid(alpha=0.3, which="both")
ax.set_ylim(0.5, 6)
ax.set_xlim(2, 1000)
plt.tight_layout()
plt.show()

print(f"{len(lats)} zip codes, {nresps.sum():.0f} total responses")
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import BoundaryNorm

# Fetch DYFI data
detail = requests.get("https://earthquake.usgs.gov/earthquakes/feed/v1.0/detail/nc75240492.geojson").json()
dyfi_products = detail["properties"]["products"].get("dyfi", [])
dyfi_url = dyfi_products[0]["contents"]["dyfi_zip.geojson"]["url"] if dyfi_products else None

lats, lons, intensities, nresps = [], [], [], []
if dyfi_url:
    for f in requests.get(dyfi_url).json()["features"]:
        if f["geometry"]["type"] == "Point":
            lons.append(f["geometry"]["coordinates"][0])
            lats.append(f["geometry"]["coordinates"][1])
            intensities.append(f["properties"].get("cdi", 0))
            nresps.append(f["properties"].get("nresp", 0))
lats, lons, intensities = np.array(lats), np.array(lons), np.array(intensities)
nresps = np.array(nresps)

# Intensity vs distance
distances_km = np.array([gps2dist_azimuth(eq_lat, eq_lon, la, lo)[0] / 1000
                         for la, lo in zip(lats, lons)])

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(distances_km, intensities, alpha=0.4, s=20, c="steelblue", edgecolors="none")

dist_bins = np.logspace(np.log10(3), np.log10(800), 15)
for i in range(len(dist_bins) - 1):
    mask = (distances_km >= dist_bins[i]) & (distances_km < dist_bins[i+1])
    if mask.sum() > 3:
        ax.plot(np.sqrt(dist_bins[i] * dist_bins[i+1]), np.median(intensities[mask]),
                "ro", markersize=7, zorder=4)

ax.set_xscale("log")
ax.set_xlabel("Distance from epicenter (km)")
ax.set_ylabel("DYFI Intensity (MMI)")
ax.set_title("Intensity vs Distance")
ax.grid(alpha=0.3, which="both")
ax.set_ylim(0.5, 6)
ax.set_xlim(2, 1000)
plt.tight_layout()
plt.show()

print(f"{len(lats)} zip codes, {nresps.sum():.0f} total responses")

Question 4a: At a distance of ~10 km from the epicenter, the DYFI intensity reports range from about MMI 2 to MMI 5. Why is there so much scatter at the same distance?

Your answer here:

Question 4b: What is the key difference between magnitude and intensity? Could you determine the magnitude of an earthquake from intensity data alone?

Your answer here: