Seismic Tomography¶

In the travel times lecture, we solved the forward problem: given a velocity model, compute travel times. Now we flip the question — given observed travel times, can we recover the velocity structure? This is the inverse problem, and it's the foundation of seismic tomography.

In [1]:

Copied!





# Install dependencies
!pip install netCDF4 cartopy

# Download data files
import os, tarfile, urllib.request, json
data_url = "https://github.com/AI4EPS/EPS130_Seismology/releases/download/tomography-data/tomography_data.tar.gz"
if not os.path.exists("tomography_data"):
    print("Downloading tomography data...")
    urllib.request.urlretrieve(data_url, "tomography_data.tar.gz")
    with tarfile.open("tomography_data.tar.gz") as tf:
        tf.extractall()
    os.remove("tomography_data.tar.gz")
    print("Done.")

# Download plate boundary data (Bird, 2003)
pb_file = "tomography_data/PB2002_boundaries.json"
if not os.path.exists(pb_file):
    print("Downloading plate boundaries...")
    pb_url = "https://raw.githubusercontent.com/fraxen/tectonicplates/master/GeoJSON/PB2002_boundaries.json"
    urllib.request.urlretrieve(pb_url, pb_file)
    print("Done.")
# Install dependencies
!pip install netCDF4 cartopy

# Download data files
import os, tarfile, urllib.request, json
data_url = "https://github.com/AI4EPS/EPS130_Seismology/releases/download/tomography-data/tomography_data.tar.gz"
if not os.path.exists("tomography_data"):
    print("Downloading tomography data...")
    urllib.request.urlretrieve(data_url, "tomography_data.tar.gz")
    with tarfile.open("tomography_data.tar.gz") as tf:
        tf.extractall()
    os.remove("tomography_data.tar.gz")
    print("Done.")

# Download plate boundary data (Bird, 2003)
pb_file = "tomography_data/PB2002_boundaries.json"
if not os.path.exists(pb_file):
    print("Downloading plate boundaries...")
    pb_url = "https://raw.githubusercontent.com/fraxen/tectonicplates/master/GeoJSON/PB2002_boundaries.json"
    urllib.request.urlretrieve(pb_url, pb_file)
    print("Done.")

Requirement already satisfied: netCDF4 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (1.7.4)
Requirement already satisfied: cartopy in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (0.25.0)
Requirement already satisfied: cftime in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from netCDF4) (1.6.5)
Requirement already satisfied: certifi in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from netCDF4) (2026.2.25)
Requirement already satisfied: numpy>=1.21.2 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from netCDF4) (2.3.5)
Requirement already satisfied: matplotlib>=3.6 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from cartopy) (3.10.8)
Requirement already satisfied: shapely>=2.0 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from cartopy) (2.1.2)
Requirement already satisfied: packaging>=21 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from cartopy) (25.0)
Requirement already satisfied: pyshp>=2.3 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from cartopy) (3.0.3)
Requirement already satisfied: pyproj>=3.3.1 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from cartopy) (3.7.2)
Requirement already satisfied: contourpy>=1.0.1 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (1.3.3)
Requirement already satisfied: cycler>=0.10 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (4.61.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (1.4.9)
Requirement already satisfied: pillow>=8 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (12.1.0)
Requirement already satisfied: pyparsing>=3 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (3.2.5)
Requirement already satisfied: python-dateutil>=2.7 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from matplotlib>=3.6->cartopy) (2.9.0.post0)
Requirement already satisfied: six>=1.5 in /Users/weiqiang/.local/miniconda3/lib/python3.12/site-packages (from python-dateutil>=2.7->matplotlib>=3.6->cartopy) (1.17.0)

Downloading tomography data...

Done.
Downloading plate boundaries...
Done.

/var/folders/gc/lpnp82h92tv35c_7v97w97cm0000gn/T/ipykernel_58523/2020279363.py:11: DeprecationWarning: Python 3.14 will, by default, filter extracted tar archives and reject files or modify their metadata. Use the filter argument to control this behavior.
  tf.extractall()

In [2]:

Copied!

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt

The cell below defines helper functions used throughout this notebook. You don't need to read this code — just run it. Here's what each function does:

build_G(...) — Computes the path-length matrix $\mathbf{G}$ for 1D layered models (used in Section 1)
smoothing_matrix_1d(n) — Builds the smoothing matrix $\mathbf{L}$ for 1D (differences between adjacent layers)
smoothing_matrix_2d(nx, nz) — Builds the smoothing matrix $\mathbf{L}$ for 2D grids (differences between adjacent blocks)

In [3]:

Copied!





def build_G(p_values, v_true, dz):
    """Build the path-length matrix G for 1D layers.
    G[i,j] = path length of ray i in layer j."""
    n_rays, n_layers = len(p_values), len(v_true)
    G = np.zeros((n_rays, n_layers))
    for i, p in enumerate(p_values):
        for j in range(n_layers):
            if p >= 1.0 / v_true[j]:
                break
            cos_theta = np.sqrt(1 - (p * v_true[j])**2)
            G[i, j] = 2 * dz[j] / cos_theta
    return G

def smoothing_matrix_1d(n):
    """Build 1D smoothing matrix L: differences between adjacent layers.
    For 4 layers: L is 3x4, each row computes (layer[j+1] - layer[j])."""
    L = np.zeros((n - 1, n))
    for i in range(n - 1):
        L[i, i] = -1
        L[i, i + 1] = 1
    return L

def smoothing_matrix_2d(nx, nz):
    """Build 2D smoothing matrix L: differences between adjacent blocks.
    Penalizes differences between horizontal and vertical neighbors."""
    n = nx * nz
    rows = []
    for iz in range(nz):
        for ix in range(nx):
            k = iz * nx + ix
            if ix < nx - 1:  # horizontal neighbor
                row = np.zeros(n)
                row[k] = -1
                row[k + 1] = 1
                rows.append(row)
            if iz < nz - 1:  # vertical neighbor
                row = np.zeros(n)
                row[k] = -1
                row[k + nx] = 1
                rows.append(row)
    return np.array(rows)
def build_G(p_values, v_true, dz):
    """Build the path-length matrix G for 1D layers.
    G[i,j] = path length of ray i in layer j."""
    n_rays, n_layers = len(p_values), len(v_true)
    G = np.zeros((n_rays, n_layers))
    for i, p in enumerate(p_values):
        for j in range(n_layers):
            if p >= 1.0 / v_true[j]:
                break
            cos_theta = np.sqrt(1 - (p * v_true[j])**2)
            G[i, j] = 2 * dz[j] / cos_theta
    return G

def smoothing_matrix_1d(n):
    """Build 1D smoothing matrix L: differences between adjacent layers.
    For 4 layers: L is 3x4, each row computes (layer[j+1] - layer[j])."""
    L = np.zeros((n - 1, n))
    for i in range(n - 1):
        L[i, i] = -1
        L[i, i + 1] = 1
    return L

def smoothing_matrix_2d(nx, nz):
    """Build 2D smoothing matrix L: differences between adjacent blocks.
    Penalizes differences between horizontal and vertical neighbors."""
    n = nx * nz
    rows = []
    for iz in range(nz):
        for ix in range(nx):
            k = iz * nx + ix
            if ix < nx - 1:  # horizontal neighbor
                row = np.zeros(n)
                row[k] = -1
                row[k + 1] = 1
                rows.append(row)
            if iz < nz - 1:  # vertical neighbor
                row = np.zeros(n)
                row[k] = -1
                row[k + nx] = 1
                rows.append(row)
    return np.array(rows)

1. The Inverse Problem: 1D Inversion¶

Consider rays traveling through horizontal layers. The travel time of ray $i$ is:

$$t_i = \sum_{j=1}^{M} G_{ij} \, s_j$$

where $s_j = 1/v_j$ is the slowness (what we want to find) and $G_{ij}$ is the path length of ray $i$ in layer $j$. In matrix form: $\mathbf{Gm} = \mathbf{d}$.

Symbol	Meaning	Size
$\mathbf{G}$	Path-length matrix	$N_{\text{rays}} \times M_{\text{layers}}$
$\mathbf{m}$	Slowness in each layer	$M_{\text{layers}} \times 1$
$\mathbf{d}$	Observed travel times	$N_{\text{rays}} \times 1$

We set up a 4-layer model with 6 rays at different ray parameters $p$. The path length through each layer is:

$$\ell_{ij} = \frac{2 \, \Delta z_j}{\sqrt{1 - p_i^2 v_j^2}}$$

In [4]:

Copied!





# True model (unknown to the "observer")
v_true = np.array([4.0, 5.5, 7.0, 8.0])  # km/s
s_true = 1.0 / v_true                      # slowness (s/km)
dz = np.array([3.0, 4.0, 5.0, 6.0])       # layer thicknesses (km)
n_layers = len(v_true)

# Ray parameters for 6 rays (s/km)
p_values = np.array([0.02, 0.05, 0.08, 0.10, 0.12, 0.14])
n_rays = len(p_values)

G = build_G(p_values, v_true, dz)

print("G matrix (path lengths in km):")
print(f"{'p (s/km)':<12}", ''.join(f'Layer {j+1:>5}' for j in range(n_layers)))
for i in range(n_rays):
    print(f"{p_values[i]:<12.2f}", ''.join(f'{G[i,j]:>9.2f}' for j in range(n_layers)))
# True model (unknown to the "observer")
v_true = np.array([4.0, 5.5, 7.0, 8.0])  # km/s
s_true = 1.0 / v_true                      # slowness (s/km)
dz = np.array([3.0, 4.0, 5.0, 6.0])       # layer thicknesses (km)
n_layers = len(v_true)

# Ray parameters for 6 rays (s/km)
p_values = np.array([0.02, 0.05, 0.08, 0.10, 0.12, 0.14])
n_rays = len(p_values)

G = build_G(p_values, v_true, dz)

print("G matrix (path lengths in km):")
print(f"{'p (s/km)':<12}", ''.join(f'Layer {j+1:>5}' for j in range(n_layers)))
for i in range(n_rays):
    print(f"{p_values[i]:<12.2f}", ''.join(f'{G[i,j]:>9.2f}' for j in range(n_layers)))

G matrix (path lengths in km):
p (s/km)     Layer     1Layer     2Layer     3Layer     4
0.02              6.02     8.05    10.10    12.16
0.05              6.12     8.32    10.68    13.09
0.08              6.33     8.91    12.07    15.62
0.10              6.55     9.58    14.00    20.00
0.12              6.84    10.65    18.43    42.86
0.14              7.24    12.54    50.25     0.00

We set up a 4-layer velocity model and 6 rays with different ray parameters. The build_G function computes how far each ray travels in each layer — this fills the $\mathbf{G}$ matrix.

In [5]:

Copied!





# Visualize the G matrix and ray paths side by side
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Left: G matrix as image
im = axes[0].imshow(G, aspect='auto', cmap='YlOrRd')
axes[0].set_xlabel('Layer index')
axes[0].set_ylabel('Ray index')
axes[0].set_xticks(range(n_layers))
axes[0].set_xticklabels([f'Layer {j+1}' for j in range(n_layers)])
axes[0].set_yticks(range(n_rays))
axes[0].set_yticklabels([f'p={p:.2f}' for p in p_values])
axes[0].set_title('G matrix (path length, km)')
for i in range(n_rays):
    for j in range(n_layers):
        if G[i, j] > 0:
            axes[0].text(j, i, f'{G[i,j]:.1f}', ha='center', va='center', fontsize=8)
plt.colorbar(im, ax=axes[0], shrink=0.8)

# Right: ray paths
ax = axes[1]
colors = plt.cm.tab10(np.linspace(0, 0.6, n_rays))
z_interfaces = np.concatenate([[0], np.cumsum(dz)])
for z_int in z_interfaces:
    ax.axhline(z_int, color='k', lw=0.5, ls='--', alpha=0.5)
for j in range(n_layers):
    z_mid = z_interfaces[j] + dz[j] / 2
    ax.text(18, z_mid, f'v={v_true[j]} km/s', va='center', fontsize=9,
            bbox=dict(boxstyle='round', fc='wheat', alpha=0.8))

for i, p in enumerate(p_values):
    x, z = 0.0, 0.0
    x_pts, z_pts = [0.0], [0.0]
    for j in range(n_layers):
        if p >= 1.0 / v_true[j]:
            break
        cos_theta = np.sqrt(1 - (p * v_true[j])**2)
        dx = dz[j] * p * v_true[j] / cos_theta
        x += dx
        z += dz[j]
        x_pts.append(x)
        z_pts.append(z)
    x_up = [2 * x_pts[-1] - xi for xi in reversed(x_pts)]
    z_up = list(reversed(z_pts))
    ax.plot(x_pts + x_up[1:], z_pts + z_up[1:], '-', color=colors[i], lw=2, label=f'p={p:.2f}')

ax.set_ylim(sum(dz) + 1, -1)
ax.set_xlim(-0.5, 22)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Depth (km)')
ax.set_title('Ray paths through the model')
ax.legend(fontsize=8, loc='lower right')
plt.tight_layout()
plt.show()
# Visualize the G matrix and ray paths side by side
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Left: G matrix as image
im = axes[0].imshow(G, aspect='auto', cmap='YlOrRd')
axes[0].set_xlabel('Layer index')
axes[0].set_ylabel('Ray index')
axes[0].set_xticks(range(n_layers))
axes[0].set_xticklabels([f'Layer {j+1}' for j in range(n_layers)])
axes[0].set_yticks(range(n_rays))
axes[0].set_yticklabels([f'p={p:.2f}' for p in p_values])
axes[0].set_title('G matrix (path length, km)')
for i in range(n_rays):
    for j in range(n_layers):
        if G[i, j] > 0:
            axes[0].text(j, i, f'{G[i,j]:.1f}', ha='center', va='center', fontsize=8)
plt.colorbar(im, ax=axes[0], shrink=0.8)

# Right: ray paths
ax = axes[1]
colors = plt.cm.tab10(np.linspace(0, 0.6, n_rays))
z_interfaces = np.concatenate([[0], np.cumsum(dz)])
for z_int in z_interfaces:
    ax.axhline(z_int, color='k', lw=0.5, ls='--', alpha=0.5)
for j in range(n_layers):
    z_mid = z_interfaces[j] + dz[j] / 2
    ax.text(18, z_mid, f'v={v_true[j]} km/s', va='center', fontsize=9,
            bbox=dict(boxstyle='round', fc='wheat', alpha=0.8))

for i, p in enumerate(p_values):
    x, z = 0.0, 0.0
    x_pts, z_pts = [0.0], [0.0]
    for j in range(n_layers):
        if p >= 1.0 / v_true[j]:
            break
        cos_theta = np.sqrt(1 - (p * v_true[j])**2)
        dx = dz[j] * p * v_true[j] / cos_theta
        x += dx
        z += dz[j]
        x_pts.append(x)
        z_pts.append(z)
    x_up = [2 * x_pts[-1] - xi for xi in reversed(x_pts)]
    z_up = list(reversed(z_pts))
    ax.plot(x_pts + x_up[1:], z_pts + z_up[1:], '-', color=colors[i], lw=2, label=f'p={p:.2f}')

ax.set_ylim(sum(dz) + 1, -1)
ax.set_xlim(-0.5, 22)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Depth (km)')
ax.set_title('Ray paths through the model')
ax.legend(fontsize=8, loc='lower right')
plt.tight_layout()
plt.show()

The left panel shows the $\mathbf{G}$ matrix as an image — each entry is a path length (km). The right panel shows the ray paths through the layers. Notice that steeper rays (small $p$) penetrate deeper, while shallow rays (large $p$) turn before reaching the bottom layer.

Solving the inverse problem¶

With more rays than layers ($N > M$), the system is overdetermined. The least-squares solution minimizes $||\mathbf{Gm} - \mathbf{d}||^2$:

$$\mathbf{m} = (\mathbf{G}^T \mathbf{G})^{-1} \mathbf{G}^T \mathbf{d}$$

In [6]:

Copied!





# Step 1: Forward problem — compute synthetic data (no noise)
d = G @ s_true

# ===== KEY STEP: Inversion =====
# m = (G^T G)^{-1} G^T d
s_recovered = np.linalg.inv(G.T @ G) @ G.T @ d
# ================================

v_recovered = 1.0 / s_recovered

print(f"{'Layer':<8} {'v_true (km/s)':<16} {'v_recovered (km/s)':<20} {'Error (%)':<10}")
for j in range(n_layers):
    err = 100 * abs(v_recovered[j] - v_true[j]) / v_true[j]
    print(f"{j+1:<8} {v_true[j]:<16.2f} {v_recovered[j]:<20.4f} {err:<10.4f}")
# Step 1: Forward problem — compute synthetic data (no noise)
d = G @ s_true

# ===== KEY STEP: Inversion =====
# m = (G^T G)^{-1} G^T d
s_recovered = np.linalg.inv(G.T @ G) @ G.T @ d
# ================================

v_recovered = 1.0 / s_recovered

print(f"{'Layer':<8} {'v_true (km/s)':<16} {'v_recovered (km/s)':<20} {'Error (%)':<10}")
for j in range(n_layers):
    err = 100 * abs(v_recovered[j] - v_true[j]) / v_true[j]
    print(f"{j+1:<8} {v_true[j]:<16.2f} {v_recovered[j]:<20.4f} {err:<10.4f}")

Layer    v_true (km/s)    v_recovered (km/s)   Error (%) 
1        4.00             4.0000               0.0000    
2        5.50             5.5000               0.0000    
3        7.00             7.0000               0.0000    
4        8.00             8.0000               0.0000

With perfect data, we recover the model exactly. But real data always has noise.

Effect of noise and smoothing¶

With noisy data, the least-squares solution can become unstable. Regularized least squares adds a penalty that prefers smooth models — neighboring layers should have similar velocities unless the data strongly require a jump:

$$\mathbf{m} = (\mathbf{G}^T \mathbf{G} + \alpha^2 \mathbf{L}^T \mathbf{L})^{-1} \mathbf{G}^T \mathbf{d}$$

where $\mathbf{L}$ is a smoothing matrix (finite-difference operator between adjacent layers) and $\alpha$ controls the strength:

$\alpha \to 0$: pure least-squares (fits data closely, but noisy model)
$\alpha$ large: heavily smoothed (smooth model, but poor data fit)

In [7]:

Copied!





# Add noise to the data
np.random.seed(3)
noise_level = 0.05  # seconds
d_noisy = d + noise_level * np.random.randn(n_rays)

# Build 1D smoothing matrix
L1d = smoothing_matrix_1d(n_layers)
LtL_1d = L1d.T @ L1d

# ===== KEY STEP: Inversion with different smoothing strengths =====
# Undamped: m = (G^T G)^{-1} G^T d
s_undamped = np.linalg.inv(G.T @ G) @ G.T @ d_noisy

# Smoothed: m = (G^T G + alpha^2 L^T L)^{-1} G^T d
alphas = [0.05, 0.5, 50.0]
s_smooth = {}
for alpha in alphas:
    s_smooth[alpha] = np.linalg.inv(G.T @ G + alpha**2 * LtL_1d) @ G.T @ d_noisy
# ==================================================================

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 4, figsize=(16, 5), sharey=True)

cases = [('Undamped', s_undamped)] + [(f'$\\alpha$ = {a}', s_smooth[a]) for a in alphas]
for ax, (label, s_inv) in zip(axes, cases):
    v_inv = 1.0 / np.clip(s_inv, 1e-6, None)
    ax.step(np.concatenate([[v_true[0]], v_true]), z_interfaces, 'k-', lw=2, label='True', where='pre')
    ax.step(np.concatenate([[v_inv[0]], v_inv]), z_interfaces, 'r--', lw=2, label='Recovered', where='pre')
    ax.set_xlabel('Velocity (km/s)')
    ax.set_title(label)
    ax.legend(fontsize=8)
    ax.set_xlim(0, 12)

axes[0].set_ylabel('Depth (km)')
axes[0].invert_yaxis()
fig.suptitle(f'Effect of smoothing (noise = {noise_level} s)', fontsize=13)
plt.tight_layout()
plt.show()
# Add noise to the data
np.random.seed(3)
noise_level = 0.05  # seconds
d_noisy = d + noise_level * np.random.randn(n_rays)

# Build 1D smoothing matrix
L1d = smoothing_matrix_1d(n_layers)
LtL_1d = L1d.T @ L1d

# ===== KEY STEP: Inversion with different smoothing strengths =====
# Undamped: m = (G^T G)^{-1} G^T d
s_undamped = np.linalg.inv(G.T @ G) @ G.T @ d_noisy

# Smoothed: m = (G^T G + alpha^2 L^T L)^{-1} G^T d
alphas = [0.05, 0.5, 50.0]
s_smooth = {}
for alpha in alphas:
    s_smooth[alpha] = np.linalg.inv(G.T @ G + alpha**2 * LtL_1d) @ G.T @ d_noisy
# ==================================================================

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 4, figsize=(16, 5), sharey=True)

cases = [('Undamped', s_undamped)] + [(f'$\\alpha$ = {a}', s_smooth[a]) for a in alphas]
for ax, (label, s_inv) in zip(axes, cases):
    v_inv = 1.0 / np.clip(s_inv, 1e-6, None)
    ax.step(np.concatenate([[v_true[0]], v_true]), z_interfaces, 'k-', lw=2, label='True', where='pre')
    ax.step(np.concatenate([[v_inv[0]], v_inv]), z_interfaces, 'r--', lw=2, label='Recovered', where='pre')
    ax.set_xlabel('Velocity (km/s)')
    ax.set_title(label)
    ax.legend(fontsize=8)
    ax.set_xlim(0, 12)

axes[0].set_ylabel('Depth (km)')
axes[0].invert_yaxis()
fig.suptitle(f'Effect of smoothing (noise = {noise_level} s)', fontsize=13)
plt.tight_layout()
plt.show()

Try it yourself¶

Change noise_level and alphas in the cell above. What happens when noise is very large? What smoothing strength gives the best recovery?

2. 2D Block Tomography¶

The 1D inversion assumed the Earth varies only with depth. In reality, velocity varies laterally too. We divide the Earth into a 2D grid of $N_x \times N_z$ blocks, each with unknown slowness $s_k$:

$$t_i = \sum_{k=1}^{M} G_{ik} \, s_k$$

where $G_{ik}$ is the length of ray $i$ inside block $k$. For simplicity, we trace straight rays between sources and receivers.

Textbook example: 20x20 grid (Shearer Exercise 2)¶

The textbook provides a pre-computed $\mathbf{G}$ matrix and travel time data for a 20x20 block model (400 unknowns) with 118 rays:

20 horizontal rays (one per row)
20 vertical rays (one per column)
39 diagonal rays at +45 degrees
39 diagonal rays at -45 degrees

Having four ray directions gives good angular coverage, which is critical for resolving 2D structure. This is the same setup you will use in the homework.

We invert for slowness perturbations using smoothing regularization:

$$\delta \mathbf{s} = (\mathbf{G}^T \mathbf{G} + \alpha^2 \mathbf{L}^T \mathbf{L})^{-1} \mathbf{G}^T \, \delta \mathbf{d}$$

where $\mathbf{L}$ penalizes differences between adjacent blocks (prefers smooth models).

The textbook provides a pre-computed $\mathbf{G}$ matrix stored in a sparse format (only non-zero entries). The cell below loads it into a full matrix. Just run it — the file I/O details aren't important.

In [8]:

Copied!





# Load textbook G matrix (sparse format: ray_index, block_index, path_length)
nx2, nz2 = 20, 20
n_blocks2 = nx2 * nz2  # 400 unknowns

# Each row of the file is: [ray_number, block_number, path_length_in_that_block]
gmat_raw = np.loadtxt('tomography_data/tomo_gmat.txt')
n_rays2 = int(gmat_raw[:, 0].max())  # 118 rays total

# Build full G matrix from sparse entries
G2 = np.zeros((n_rays2, n_blocks2))
for ray_i, mod_j, path_len in gmat_raw:
    G2[int(ray_i) - 1, int(mod_j) - 1] = path_len  # -1 for 0-based indexing

print(f"Grid: {nx2}x{nz2} = {n_blocks2} blocks")
print(f"Rays: {n_rays2} (20 horizontal + 20 vertical + 39+39 diagonal)")
print(f"G matrix: {G2.shape}, non-zeros: {np.count_nonzero(G2)}")
# Load textbook G matrix (sparse format: ray_index, block_index, path_length)
nx2, nz2 = 20, 20
n_blocks2 = nx2 * nz2  # 400 unknowns

# Each row of the file is: [ray_number, block_number, path_length_in_that_block]
gmat_raw = np.loadtxt('tomography_data/tomo_gmat.txt')
n_rays2 = int(gmat_raw[:, 0].max())  # 118 rays total

# Build full G matrix from sparse entries
G2 = np.zeros((n_rays2, n_blocks2))
for ray_i, mod_j, path_len in gmat_raw:
    G2[int(ray_i) - 1, int(mod_j) - 1] = path_len  # -1 for 0-based indexing

print(f"Grid: {nx2}x{nz2} = {n_blocks2} blocks")
print(f"Rays: {n_rays2} (20 horizontal + 20 vertical + 39+39 diagonal)")
print(f"G matrix: {G2.shape}, non-zeros: {np.count_nonzero(G2)}")

Grid: 20x20 = 400 blocks
Rays: 118 (20 horizontal + 20 vertical + 39+39 diagonal)
G matrix: (118, 400), non-zeros: 1600

In [9]:

Copied!





# Visualize the four ray types
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
ray_groups = [
    (range(0, 20), 'Horizontal (20)'),
    (range(20, 40), 'Vertical (20)'),
    (range(40, 79), 'Diagonal +45\u00b0 (39)'),
    (range(79, 118), 'Diagonal \u221245\u00b0 (39)'),
]

for ax, (rays, title) in zip(axes, ray_groups):
    for i in range(nx2 + 1):
        ax.axvline(i, color='gray', lw=0.3)
    for j in range(nz2 + 1):
        ax.axhline(j, color='gray', lw=0.3)

    for r in rays:
        blocks = np.where(G2[r] > 0)[0]
        if len(blocks) == 0:
            continue
        rows, cols = blocks // nx2, blocks % nx2
        if title.startswith('Horizontal'):
            order = np.argsort(cols)
            ax.plot([cols[order[0]] + 0.5, cols[order[-1]] + 0.5],
                    [rows[0] + 0.5, rows[0] + 0.5], 'b-', alpha=0.4, lw=0.8)
        elif title.startswith('Vertical'):
            order = np.argsort(rows)
            ax.plot([cols[0] + 0.5, cols[0] + 0.5],
                    [rows[order[0]] + 0.5, rows[order[-1]] + 0.5], 'r-', alpha=0.4, lw=0.8)
        else:
            order = np.argsort(cols)
            ax.plot([cols[order[0]] + 0.5, cols[order[-1]] + 0.5],
                    [rows[order[0]] + 0.5, rows[order[-1]] + 0.5], 'g-', alpha=0.4, lw=0.8)

    ax.set_xlim(0, nx2)
    ax.set_ylim(nz2, 0)
    ax.set_aspect('equal')
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('Column')

axes[0].set_ylabel('Row')
plt.suptitle('Textbook ray geometry (Shearer Fig. 5.8)', fontsize=13)
plt.tight_layout()
plt.show()
# Visualize the four ray types
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
ray_groups = [
    (range(0, 20), 'Horizontal (20)'),
    (range(20, 40), 'Vertical (20)'),
    (range(40, 79), 'Diagonal +45\u00b0 (39)'),
    (range(79, 118), 'Diagonal \u221245\u00b0 (39)'),
]

for ax, (rays, title) in zip(axes, ray_groups):
    for i in range(nx2 + 1):
        ax.axvline(i, color='gray', lw=0.3)
    for j in range(nz2 + 1):
        ax.axhline(j, color='gray', lw=0.3)

    for r in rays:
        blocks = np.where(G2[r] > 0)[0]
        if len(blocks) == 0:
            continue
        rows, cols = blocks // nx2, blocks % nx2
        if title.startswith('Horizontal'):
            order = np.argsort(cols)
            ax.plot([cols[order[0]] + 0.5, cols[order[-1]] + 0.5],
                    [rows[0] + 0.5, rows[0] + 0.5], 'b-', alpha=0.4, lw=0.8)
        elif title.startswith('Vertical'):
            order = np.argsort(rows)
            ax.plot([cols[0] + 0.5, cols[0] + 0.5],
                    [rows[order[0]] + 0.5, rows[order[-1]] + 0.5], 'r-', alpha=0.4, lw=0.8)
        else:
            order = np.argsort(cols)
            ax.plot([cols[order[0]] + 0.5, cols[order[-1]] + 0.5],
                    [rows[order[0]] + 0.5, rows[order[-1]] + 0.5], 'g-', alpha=0.4, lw=0.8)

    ax.set_xlim(0, nx2)
    ax.set_ylim(nz2, 0)
    ax.set_aspect('equal')
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('Column')

axes[0].set_ylabel('Row')
plt.suptitle('Textbook ray geometry (Shearer Fig. 5.8)', fontsize=13)
plt.tight_layout()
plt.show()

The plot below shows the four types of rays. Notice that horizontal and vertical rays each cover one row/column of blocks, while diagonal rays cut across at 45°. Having rays at multiple angles is key to resolving 2D structure — a single direction can only see averages along that direction.

Spike test (point-spread function)¶

Place a single anomalous block in the model and see how the inversion smears it. This reveals the point-spread function — if we can't recover a single point perfectly, we can't resolve fine details either.

In [10]:

Copied!





# Build 2D smoothing matrix for the 20x20 grid
L2d_20 = smoothing_matrix_2d(nx2, nz2)
LtL_20 = L2d_20.T @ L2d_20

# Spike test: single anomalous block at center
alpha = 0.0
noise_level = 0.0

m_spike = np.zeros(n_blocks2)
m_spike[10 * nx2 + 10] = 1.0  # single anomaly at row 10, col 10

np.random.seed(42)
d_spike = G2 @ m_spike + noise_level * np.random.randn(n_rays2)

# ===== KEY STEP: Inversion =====
# m = (G^T G + alpha^2 L^T L)^{-1} G^T d
m_spike_rec = np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d_spike
# ================================

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)

im0 = axes[0].imshow(m_spike.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True spike model')
plt.colorbar(im0, ax=axes[0], shrink=0.8, label='Slowness perturbation')

im1 = axes[1].imshow(m_spike_rec.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-0.3, vmax=0.3)
axes[1].set_title(f'Recovered ($\\alpha$={alpha}, noise={noise_level})')
plt.colorbar(im1, ax=axes[1], shrink=0.8, label='Slowness perturbation')

for ax in axes:
    ax.set_xlabel('Column')
    ax.set_ylabel('Row')
plt.suptitle('Spike test: how a point anomaly gets smeared', fontsize=13)
plt.show()
# Build 2D smoothing matrix for the 20x20 grid
L2d_20 = smoothing_matrix_2d(nx2, nz2)
LtL_20 = L2d_20.T @ L2d_20

# Spike test: single anomalous block at center
alpha = 0.0
noise_level = 0.0

m_spike = np.zeros(n_blocks2)
m_spike[10 * nx2 + 10] = 1.0  # single anomaly at row 10, col 10

np.random.seed(42)
d_spike = G2 @ m_spike + noise_level * np.random.randn(n_rays2)

# ===== KEY STEP: Inversion =====
# m = (G^T G + alpha^2 L^T L)^{-1} G^T d
m_spike_rec = np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d_spike
# ================================

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)

im0 = axes[0].imshow(m_spike.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True spike model')
plt.colorbar(im0, ax=axes[0], shrink=0.8, label='Slowness perturbation')

im1 = axes[1].imshow(m_spike_rec.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-0.3, vmax=0.3)
axes[1].set_title(f'Recovered ($\\alpha$={alpha}, noise={noise_level})')
plt.colorbar(im1, ax=axes[1], shrink=0.8, label='Slowness perturbation')

for ax in axes:
    ax.set_xlabel('Column')
    ax.set_ylabel('Row')
plt.suptitle('Spike test: how a point anomaly gets smeared', fontsize=13)
plt.show()

Try it yourself¶

Change alpha and noise_level in the cell above. What happens when alpha = 0? Why does it fail? Try moving the spike to a corner block (e.g., m_spike[0] = 1.0) — how does the recovery change for blocks at the edge vs. the center?

Checkerboard test¶

A more realistic resolution test: alternate fast and slow blocks in a checkerboard pattern, generate synthetic data (with noise), and see how well the inversion recovers the pattern. Regions that are well-resolved will show a clear checkerboard; poorly-resolved regions will be smeared.

In [11]:

Copied!





# Checkerboard test
alpha = 0.001
noise_level = 5.0

# Build checkerboard: alternating +1/-1 in 4x4 blocks
m_checker2 = np.zeros((nz2, nx2))
block_size = 4
for iz in range(nz2):
    for ix in range(nx2):
        if ((iz // block_size) + (ix // block_size)) % 2 == 0:
            m_checker2[iz, ix] = 1.0
        else:
            m_checker2[iz, ix] = -1.0

np.random.seed(42)
d_checker2 = G2 @ m_checker2.ravel() + noise_level * np.random.randn(n_rays2)

# ===== KEY STEP: Inversion =====
# m = (G^T G + alpha^2 L^T L)^{-1} G^T d
m_checker_rec = np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d_checker2
# ================================

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)

im0 = axes[0].imshow(m_checker2, cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True checkerboard')
plt.colorbar(im0, ax=axes[0], shrink=0.8, label='Slowness perturbation')

im1 = axes[1].imshow(m_checker_rec.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[1].set_title(f'Recovered ($\\alpha$={alpha}, noise={noise_level})')
plt.colorbar(im1, ax=axes[1], shrink=0.8, label='Slowness perturbation')

for ax in axes:
    ax.set_xlabel('Column')
    ax.set_ylabel('Row')
plt.suptitle('Checkerboard test: can the inversion resolve the pattern?', fontsize=13)
plt.show()
# Checkerboard test
alpha = 0.001
noise_level = 5.0

# Build checkerboard: alternating +1/-1 in 4x4 blocks
m_checker2 = np.zeros((nz2, nx2))
block_size = 4
for iz in range(nz2):
    for ix in range(nx2):
        if ((iz // block_size) + (ix // block_size)) % 2 == 0:
            m_checker2[iz, ix] = 1.0
        else:
            m_checker2[iz, ix] = -1.0

np.random.seed(42)
d_checker2 = G2 @ m_checker2.ravel() + noise_level * np.random.randn(n_rays2)

# ===== KEY STEP: Inversion =====
# m = (G^T G + alpha^2 L^T L)^{-1} G^T d
m_checker_rec = np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d_checker2
# ================================

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)

im0 = axes[0].imshow(m_checker2, cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True checkerboard')
plt.colorbar(im0, ax=axes[0], shrink=0.8, label='Slowness perturbation')

im1 = axes[1].imshow(m_checker_rec.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[1].set_title(f'Recovered ($\\alpha$={alpha}, noise={noise_level})')
plt.colorbar(im1, ax=axes[1], shrink=0.8, label='Slowness perturbation')

for ax in axes:
    ax.set_xlabel('Column')
    ax.set_ylabel('Row')
plt.suptitle('Checkerboard test: can the inversion resolve the pattern?', fontsize=13)
plt.show()

Effect of smoothing strength¶

The smoothing parameter $\alpha$ controls the balance between fitting the data and keeping the model smooth. Too little smoothing amplifies noise; too much smoothing smears out the structure.

In [12]:

Copied!





# Effect of smoothing on checkerboard recovery
alpha_values = [0.005, 5.0, 50.0]

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 4, figsize=(18, 4), constrained_layout=True)

# True model
im0 = axes[0].imshow(m_checker2, cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True checkerboard')
axes[0].set_ylabel('Row')

for ax, alpha in zip(axes[1:], alpha_values):
    # ===== KEY STEP: Inversion with different alpha =====
    m_rec = np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d_checker2
    # ====================================================
    im = ax.imshow(m_rec.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                   extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
    ax.set_title(f'$\\alpha$ = {alpha}')

for ax in axes:
    ax.set_xlabel('Column')
plt.colorbar(im, ax=axes, label='Slowness perturbation', shrink=0.8)
fig.suptitle('Effect of smoothing on checkerboard recovery (with noise)', fontsize=13)
plt.show()
# Effect of smoothing on checkerboard recovery
alpha_values = [0.005, 5.0, 50.0]

# ----- Plotting (you can ignore the details below) -----
fig, axes = plt.subplots(1, 4, figsize=(18, 4), constrained_layout=True)

# True model
im0 = axes[0].imshow(m_checker2, cmap='RdBu', aspect='equal',
                      extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True checkerboard')
axes[0].set_ylabel('Row')

for ax, alpha in zip(axes[1:], alpha_values):
    # ===== KEY STEP: Inversion with different alpha =====
    m_rec = np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d_checker2
    # ====================================================
    im = ax.imshow(m_rec.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
                   extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
    ax.set_title(f'$\\alpha$ = {alpha}')

for ax in axes:
    ax.set_xlabel('Column')
plt.colorbar(im, ax=axes, label='Slowness perturbation', shrink=0.8)
fig.suptitle('Effect of smoothing on checkerboard recovery (with noise)', fontsize=13)
plt.show()

Try it yourself¶

Change noise_level, alpha, and block_size in the checkerboard cell above, and alpha_values in the smoothing comparison cell. What happens with smaller blocks (e.g., block_size=2)? Can the inversion still resolve the pattern?

Inverting real data¶

In the synthetic tests above, we computed the data: $\mathbf{d} = \mathbf{G} \mathbf{m}_{\text{true}}$ (+ noise). Now we load real travel time observations — the same $\mathbf{G}$ matrix, but with measured $\delta \mathbf{d}$ from the textbook. Unlike synthetic tests, we don't know the true model.

Invert the data using the same formula and try different $\alpha$ values:

$$\delta \mathbf{s} = (\mathbf{G}^T \mathbf{G} + \alpha^2 \mathbf{L}^T \mathbf{L})^{-1} \mathbf{G}^T \, \delta \mathbf{d}$$

Think about:

Which features appear consistently across different $\alpha$ values? Those are likely real.
Which features change dramatically? Those may be artifacts.
What happens when $\alpha = 0$? Why?

In [13]:

Copied!





# Load real travel time data
data_raw = np.loadtxt('tomography_data/tomo_data.txt')
d2 = data_raw[:, 1]
print(f"Data: {len(d2)} travel time perturbations, {np.sum(d2 != 0)} non-zero")

# YOUR CODE HERE: invert d2 using np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d2
# Try several alpha values and plot the results with imshow
alpha = 1.0
# Load real travel time data
data_raw = np.loadtxt('tomography_data/tomo_data.txt')
d2 = data_raw[:, 1]
print(f"Data: {len(d2)} travel time perturbations, {np.sum(d2 != 0)} non-zero")

# YOUR CODE HERE: invert d2 using np.linalg.inv(G2.T @ G2 + alpha**2 * LtL_20) @ G2.T @ d2
# Try several alpha values and plot the results with imshow
alpha = 1.0

Data: 118 travel time perturbations, 58 non-zero

3. Global Tomography Models¶

The toy inversions used synthetic data on small grids. Real seismic tomography applies the same principles to millions of travel times. The result: 3D images of Earth's interior.

S362ANI: Global shear-wave velocity¶

S362ANI (Kustowski et al., 2008) is a global model of shear-wave velocity perturbations from 25 km to the core-mantle boundary (2890 km), available from the IRIS Earth Model Collaboration.

In [14]:

Copied!





import cartopy.crs as ccrs
import cartopy.feature as cfeature
import netCDF4 as nc

ds = nc.Dataset('tomography_data/S362ANI_percent.nc')
lat = ds.variables['latitude'][:]
lon = ds.variables['longitude'][:]
depth = ds.variables['depth'][:]
dvs = np.ma.filled(ds.variables['dvs'][:], fill_value=np.nan)

print(f"Model grid: {len(lat)} lat x {len(lon)} lon x {len(depth)} depths")
print(f"Depth levels (km): {depth}")
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import netCDF4 as nc

ds = nc.Dataset('tomography_data/S362ANI_percent.nc')
lat = ds.variables['latitude'][:]
lon = ds.variables['longitude'][:]
depth = ds.variables['depth'][:]
dvs = np.ma.filled(ds.variables['dvs'][:], fill_value=np.nan)

print(f"Model grid: {len(lat)} lat x {len(lon)} lon x {len(depth)} depths")
print(f"Depth levels (km): {depth}")

Model grid: 91 lat x 181 lon x 25 depths
Depth levels (km): [  25.   50.   75.  100.  125.  150.  200.  250.  300.  350.  400.  500.
  600.  700.  800. 1000. 1250. 1500. 1750. 2000. 2250. 2500. 2750. 2800.
 2890.]

Each map shows $\delta V_s / V_s$ (%) at a fixed depth. Blue = fast (cold, dense material), Red = slow (hot, buoyant material).

The plotting code below uses cartopy for map projections. Just run it and focus on interpreting the maps — the plotting details aren't important.

In [15]:

Copied!





# Load plate boundaries (Bird, 2003)
with open('tomography_data/PB2002_boundaries.json') as f:
    pb_data = json.load(f)

def plot_plate_boundaries(ax, color='k', lw=0.8, alpha=0.6):
    """Overlay PB2002 plate boundaries on a cartopy axis."""
    for feature in pb_data['features']:
        coords = feature['geometry']['coordinates']
        if feature['geometry']['type'] == 'LineString':
            coords = [coords]
        for segment in coords:
            lons, lats = zip(*segment)
            ax.plot(lons, lats, '-', color=color, lw=lw, alpha=alpha,
                    transform=ccrs.PlateCarree())

depth_slices = [100, 250, 600, 2800]  # km
LON, LAT = np.meshgrid(lon, lat)

fig, axes = plt.subplots(2, 2, figsize=(14, 10),
                         subplot_kw={'projection': ccrs.Robinson()})

for ax, target_depth in zip(axes.ravel(), depth_slices):
    iz = np.argmin(np.abs(depth - target_depth))
    data = dvs[iz]
    vmax = np.nanpercentile(np.abs(data), 98)

    im = ax.pcolormesh(LON, LAT, data, cmap='RdBu', vmin=-vmax, vmax=vmax,
                       transform=ccrs.PlateCarree(), shading='auto')
    ax.coastlines(lw=0.5)
    plot_plate_boundaries(ax)
    ax.set_title(f'Depth = {depth[iz]:.0f} km', fontsize=12)
    plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05,
                 label='$\\delta V_s / V_s$ (%)', shrink=0.7)

fig.suptitle('S362ANI: Global shear-wave velocity perturbations', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
# Load plate boundaries (Bird, 2003)
with open('tomography_data/PB2002_boundaries.json') as f:
    pb_data = json.load(f)

def plot_plate_boundaries(ax, color='k', lw=0.8, alpha=0.6):
    """Overlay PB2002 plate boundaries on a cartopy axis."""
    for feature in pb_data['features']:
        coords = feature['geometry']['coordinates']
        if feature['geometry']['type'] == 'LineString':
            coords = [coords]
        for segment in coords:
            lons, lats = zip(*segment)
            ax.plot(lons, lats, '-', color=color, lw=lw, alpha=alpha,
                    transform=ccrs.PlateCarree())

depth_slices = [100, 250, 600, 2800]  # km
LON, LAT = np.meshgrid(lon, lat)

fig, axes = plt.subplots(2, 2, figsize=(14, 10),
                         subplot_kw={'projection': ccrs.Robinson()})

for ax, target_depth in zip(axes.ravel(), depth_slices):
    iz = np.argmin(np.abs(depth - target_depth))
    data = dvs[iz]
    vmax = np.nanpercentile(np.abs(data), 98)

    im = ax.pcolormesh(LON, LAT, data, cmap='RdBu', vmin=-vmax, vmax=vmax,
                       transform=ccrs.PlateCarree(), shading='auto')
    ax.coastlines(lw=0.5)
    plot_plate_boundaries(ax)
    ax.set_title(f'Depth = {depth[iz]:.0f} km', fontsize=12)
    plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05,
                 label='$\\delta V_s / V_s$ (%)', shrink=0.7)

fig.suptitle('S362ANI: Global shear-wave velocity perturbations', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

What to notice:

100 km: Continents are fast (thick, cold lithosphere); ocean ridges are slow (hot upwelling)
250 km: Slow anomalies under East Africa (hotspot), fast under cratonic continents
600 km: Subducting slabs appear as fast anomalies (e.g., western Pacific)
2800 km: Near the core-mantle boundary — large low-velocity provinces under Africa and the Pacific ("LLSVPs")

Cross-sections through the Japan subduction zone¶

A vertical slice reveals how a subducting slab descends into the mantle.

The code below plots vertical cross-sections through the Japan subduction zone. Just run it — the cartopy/matplotlib details are complex but just produce the map and slices.

In [16]:

Copied!





# Cross-sections through the Japan subduction zone at three latitudes
lat_slices = [44, 36, 28]  # Hokkaido, central Honshu, Ryukyu (north to south)
lon_range = (125, 155)
lon_mask = (lon >= lon_range[0]) & (lon <= lon_range[1])
lon_sub = lon[lon_mask]
depth_mask = depth <= 700
depth_sub = depth[depth_mask]

colors = ['C0', 'C1', 'C2']
fig = plt.figure(figsize=(14, 10))

# Map panel (left)
ax_map = fig.add_subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax_map.set_extent([122, 158, 22, 50], crs=ccrs.PlateCarree())
ax_map.coastlines(lw=0.8)
ax_map.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5)
ax_map.add_feature(cfeature.OCEAN, facecolor='lightskyblue', alpha=0.2)
for lat_s, c in zip(lat_slices, colors):
    ax_map.plot([lon_range[0], lon_range[1]], [lat_s, lat_s],
                '-', color=c, lw=2.5, transform=ccrs.PlateCarree(), zorder=5)
    ax_map.text(lon_range[1]+0.5, lat_s, f'{lat_s}\u00b0N', fontsize=10, va='center',
                transform=ccrs.PlateCarree(), color=c, fontweight='bold')
gl = ax_map.gridlines(draw_labels=True, lw=0.5, alpha=0.5)
gl.top_labels = False
gl.right_labels = False
ax_map.set_title('Cross-section locations')

# Cross-section panels (right, stacked north to south)
LON_xs, DEPTH_xs = np.meshgrid(lon_sub, depth_sub)
gs = fig.add_gridspec(len(lat_slices), 2, width_ratios=[1, 1.5], wspace=0.35, hspace=0.35)

for i, (lat_s, c) in enumerate(zip(lat_slices, colors)):
    ilat = np.argmin(np.abs(lat - lat_s))
    dvs_slice = dvs[:, ilat, lon_mask][depth_mask, :]
    vmax = np.nanpercentile(np.abs(dvs_slice), 95)

    ax = fig.add_subplot(gs[i, 1])
    im = ax.pcolormesh(LON_xs, DEPTH_xs, dvs_slice, cmap='RdBu',
                       vmin=-vmax, vmax=vmax, shading='auto')
    ax.invert_yaxis()
    ax.set_ylabel('Depth (km)')
    ax.set_title(f'{lat_s}\u00b0N', color=c, fontweight='bold', fontsize=11)
    plt.colorbar(im, ax=ax, label='$\\delta V_s / V_s$ (%)', shrink=0.9)
    if i == len(lat_slices) - 1:
        ax.set_xlabel('Longitude (\u00b0E)')
    else:
        ax.set_xticklabels([])

plt.suptitle('S362ANI: Cross-sections through the Japan subduction zone', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
# Cross-sections through the Japan subduction zone at three latitudes
lat_slices = [44, 36, 28]  # Hokkaido, central Honshu, Ryukyu (north to south)
lon_range = (125, 155)
lon_mask = (lon >= lon_range[0]) & (lon <= lon_range[1])
lon_sub = lon[lon_mask]
depth_mask = depth <= 700
depth_sub = depth[depth_mask]

colors = ['C0', 'C1', 'C2']
fig = plt.figure(figsize=(14, 10))

# Map panel (left)
ax_map = fig.add_subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax_map.set_extent([122, 158, 22, 50], crs=ccrs.PlateCarree())
ax_map.coastlines(lw=0.8)
ax_map.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5)
ax_map.add_feature(cfeature.OCEAN, facecolor='lightskyblue', alpha=0.2)
for lat_s, c in zip(lat_slices, colors):
    ax_map.plot([lon_range[0], lon_range[1]], [lat_s, lat_s],
                '-', color=c, lw=2.5, transform=ccrs.PlateCarree(), zorder=5)
    ax_map.text(lon_range[1]+0.5, lat_s, f'{lat_s}\u00b0N', fontsize=10, va='center',
                transform=ccrs.PlateCarree(), color=c, fontweight='bold')
gl = ax_map.gridlines(draw_labels=True, lw=0.5, alpha=0.5)
gl.top_labels = False
gl.right_labels = False
ax_map.set_title('Cross-section locations')

# Cross-section panels (right, stacked north to south)
LON_xs, DEPTH_xs = np.meshgrid(lon_sub, depth_sub)
gs = fig.add_gridspec(len(lat_slices), 2, width_ratios=[1, 1.5], wspace=0.35, hspace=0.35)

for i, (lat_s, c) in enumerate(zip(lat_slices, colors)):
    ilat = np.argmin(np.abs(lat - lat_s))
    dvs_slice = dvs[:, ilat, lon_mask][depth_mask, :]
    vmax = np.nanpercentile(np.abs(dvs_slice), 95)

    ax = fig.add_subplot(gs[i, 1])
    im = ax.pcolormesh(LON_xs, DEPTH_xs, dvs_slice, cmap='RdBu',
                       vmin=-vmax, vmax=vmax, shading='auto')
    ax.invert_yaxis()
    ax.set_ylabel('Depth (km)')
    ax.set_title(f'{lat_s}\u00b0N', color=c, fontweight='bold', fontsize=11)
    plt.colorbar(im, ax=ax, label='$\\delta V_s / V_s$ (%)', shrink=0.9)
    if i == len(lat_slices) - 1:
        ax.set_xlabel('Longitude (\u00b0E)')
    else:
        ax.set_xticklabels([])

plt.suptitle('S362ANI: Cross-sections through the Japan subduction zone', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

/var/folders/gc/lpnp82h92tv35c_7v97w97cm0000gn/T/ipykernel_58523/3229784002.py:50: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

MITPS-20: Regional P-wave velocity beneath North America¶

MITPS-20 (Golos et al., 2020) is a higher-resolution model of P-wave velocity at depths from 20 to 300 km.

The code below plots depth slices of P-wave velocity beneath North America. Just run it and focus on interpreting the velocity patterns.

In [17]:

Copied!





ds2 = nc.Dataset('tomography_data/MITPS-20.nc')
lat2 = ds2.variables['latitude'][:]
lon2 = ds2.variables['longitude'][:]
depth2 = ds2.variables['depth'][:]
dvp = np.ma.filled(ds2.variables['dVp'][:], fill_value=np.nan)

depth_slices = [60, 100, 160, 200]
LON2, LAT2 = np.meshgrid(lon2, lat2)

fig, axes = plt.subplots(2, 2, figsize=(14, 10),
                         subplot_kw={'projection': ccrs.AlbersEqualArea(
                             central_longitude=-95, central_latitude=40)})

for ax, target_depth in zip(axes.ravel(), depth_slices):
    iz = np.argmin(np.abs(depth2 - target_depth))
    data = dvp[iz]
    vmax = np.nanpercentile(np.abs(data), 95)

    im = ax.pcolormesh(LON2, LAT2, data, cmap='RdBu', vmin=-vmax, vmax=vmax,
                       transform=ccrs.PlateCarree(), shading='auto')
    ax.coastlines(lw=0.5)
    ax.add_feature(cfeature.BORDERS, lw=0.3)
    ax.add_feature(cfeature.STATES, lw=0.2)
    ax.plot(-122.27, 37.87, 'k^', ms=8, transform=ccrs.PlateCarree())
    ax.set_title(f'Depth = {depth2[iz]:.0f} km', fontsize=12)
    plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05,
                 label='$\\delta V_p / V_p$ (%)', shrink=0.7)

fig.suptitle('MITPS-20: P-wave velocity perturbations beneath North America\n(\u25b2 = Berkeley)',
             fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
ds2 = nc.Dataset('tomography_data/MITPS-20.nc')
lat2 = ds2.variables['latitude'][:]
lon2 = ds2.variables['longitude'][:]
depth2 = ds2.variables['depth'][:]
dvp = np.ma.filled(ds2.variables['dVp'][:], fill_value=np.nan)

depth_slices = [60, 100, 160, 200]
LON2, LAT2 = np.meshgrid(lon2, lat2)

fig, axes = plt.subplots(2, 2, figsize=(14, 10),
                         subplot_kw={'projection': ccrs.AlbersEqualArea(
                             central_longitude=-95, central_latitude=40)})

for ax, target_depth in zip(axes.ravel(), depth_slices):
    iz = np.argmin(np.abs(depth2 - target_depth))
    data = dvp[iz]
    vmax = np.nanpercentile(np.abs(data), 95)

    im = ax.pcolormesh(LON2, LAT2, data, cmap='RdBu', vmin=-vmax, vmax=vmax,
                       transform=ccrs.PlateCarree(), shading='auto')
    ax.coastlines(lw=0.5)
    ax.add_feature(cfeature.BORDERS, lw=0.3)
    ax.add_feature(cfeature.STATES, lw=0.2)
    ax.plot(-122.27, 37.87, 'k^', ms=8, transform=ccrs.PlateCarree())
    ax.set_title(f'Depth = {depth2[iz]:.0f} km', fontsize=12)
    plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05,
                 label='$\\delta V_p / V_p$ (%)', shrink=0.7)

fig.suptitle('MITPS-20: P-wave velocity perturbations beneath North America\n(\u25b2 = Berkeley)',
             fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

What to notice:

The western US is broadly slow (hot, actively deforming: Basin & Range, Yellowstone)
The eastern US / Canadian Shield is fast (cold, thick, ancient cratonic lithosphere)
The Juan de Fuca slab appears as a fast anomaly under the Pacific Northwest (Cascadia subduction)
The Yellowstone hotspot shows as a prominent slow anomaly

Physical interpretation¶

Seismic velocity perturbations are primarily controlled by temperature:

Anomaly	Velocity	Temperature	Physical interpretation
Blue (fast)	$\delta V > 0$	Cold	Subducting slabs, cratonic lithosphere
Red (slow)	$\delta V < 0$	Hot	Mantle plumes, mid-ocean ridges, active tectonics

Typical velocity perturbations are 1-5%, corresponding to temperature variations of ~100-500 K.

Limitations of tomographic models:

Uneven ray coverage: oceans have fewer stations, so lower resolution in oceanic mantle
Regularization: damping/smoothing means we can never recover the true Earth perfectly
Trade-offs: velocity vs. anisotropy, temperature vs. composition

Summary¶

The big picture¶

Seismic tomography turns travel time data into images of Earth's interior. The same linear algebra — $\mathbf{Gm} = \mathbf{d}$ — works from a 4-layer toy model all the way up to 3D global imaging with millions of rays.

What tomography has revealed:

Subducting slabs: cold oceanic plates sinking into the mantle (fast anomalies dipping from trenches)
Mantle plumes / hotspots: hot upwellings from the deep mantle (slow anomalies beneath Yellowstone, Hawaii, East Africa)
LLSVPs: two continent-sized slow structures at the base of the mantle beneath Africa and the Pacific
Cratonic roots: thick, cold, ancient lithosphere beneath stable continents (Canadian Shield, Siberia)

Key concepts¶

Concept	Key idea
Forward vs. inverse	Forward: given model, compute data. Inverse: given data, recover model
$\mathbf{Gm} = \mathbf{d}$	Path-length matrix $\times$ slowness = travel times
Smoothing regularization	$\alpha^2 \mathbf{L}^T \mathbf{L}$ penalizes rough models — prefers smooth solutions
Checkerboard test	Synthetic test showing which parts of the model can be resolved
Spike test	Point-spread function — how a point anomaly gets smeared
Ray coverage	Blocks sampled by many rays at different angles are best resolved

Always check resolution before interpreting!¶

A tomographic image is only as good as the data behind it. Before interpreting any feature:

Check ray coverage — are there enough rays crossing that region at different angles?
Look at checkerboard/spike tests — can the inversion resolve features of this size at this location?
Consider smoothing — strong smoothing smears structure; weak smoothing creates artifacts

If a feature doesn't pass these checks, it might be an artifact of poor coverage or regularization, not real Earth structure.