Seismic Tomography — Exercises¶

Complete the following exercises to practice the concepts from the tomography lecture.

Key equations:

Formula	When to use
$\mathbf{m} = (\mathbf{G}^T \mathbf{G})^{-1} \mathbf{G}^T \mathbf{d}$	Least-squares inversion (no regularization)
$\mathbf{m} = (\mathbf{G}^T \mathbf{G} + \alpha^2 \mathbf{L}^T \mathbf{L})^{-1} \mathbf{G}^T \mathbf{d}$	Regularized inversion (smoothing)

In [6]:

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.1.4)
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)

In [7]:

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. Just run it.

In [8]:

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."""
    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."""
    n = nx * nz
    rows = []
    for iz in range(nz):
        for ix in range(nx):
            k = iz * nx + ix
            if ix < nx - 1:
                row = np.zeros(n)
                row[k] = -1
                row[k + 1] = 1
                rows.append(row)
            if iz < nz - 1:
                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."""
    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."""
    n = nx * nz
    rows = []
    for iz in range(nz):
        for ix in range(nx):
            k = iz * nx + ix
            if ix < nx - 1:
                row = np.zeros(n)
                row[k] = -1
                row[k + 1] = 1
                rows.append(row)
            if iz < nz - 1:
                row = np.zeros(n)
                row[k] = -1
                row[k + nx] = 1
                rows.append(row)
    return np.array(rows)

Exercise 1¶

Recap: The forward problem is $\mathbf{d} = \mathbf{G}\mathbf{s}$, where $\mathbf{G}$ is the path-length matrix, $\mathbf{s}$ is slowness, and $\mathbf{d}$ is travel time data. The inverse solution (no noise) is $\mathbf{s} = (\mathbf{G}^T\mathbf{G})^{-1}\mathbf{G}^T\mathbf{d}$.: Forward and Inverse Problem (1D)

We have a 4-layer velocity model with 6 rays. The path-length matrix $\mathbf{G}$ relates the slowness $\mathbf{m}$ to the travel times $\mathbf{d}$:

$$\mathbf{d} = \mathbf{G} \mathbf{m}$$

(a) Compute the synthetic travel times (the forward problem).

(b) Recover the slowness from the travel times (the inverse problem) using:

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

In [9]:

Copied!





# Set up the 1D model (same as lecture)
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)

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

G = build_G(p_values, v_true, dz)
print(f"G matrix shape: {G.shape} ({n_rays} rays x {n_layers} layers)")
# Set up the 1D model (same as lecture)
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)

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

G = build_G(p_values, v_true, dz)
print(f"G matrix shape: {G.shape} ({n_rays} rays x {n_layers} layers)")

G matrix shape: (6, 4) (6 rays x 4 layers)

In [ ]:

Copied!





# (a) Forward problem: compute synthetic travel times
# d = G @ m (matrix-vector multiplication)
d = ???  # d = G @ s_true (matrix-vector multiply)

print("Synthetic travel times (s):")
for i in range(n_rays):
    print(f"  Ray {i+1} (p={p_values[i]:.2f}): t = {d[i]:.3f} s")
# (a) Forward problem: compute synthetic travel times
# d = G @ m (matrix-vector multiplication)
d = ???  # d = G @ s_true (matrix-vector multiply)

print("Synthetic travel times (s):")
for i in range(n_rays):
    print(f"  Ray {i+1} (p={p_values[i]:.2f}): t = {d[i]:.3f} s")

In [ ]:

Copied!





# (b) Inverse problem: recover slowness from travel times
# m = (G^T G)^{-1} G^T d
# Hint: use np.linalg.inv() for matrix inverse, and @ for matrix multiplication
s_recovered = ???  # d = G @ s_true (matrix-vector multiply)

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.6f}")
# (b) Inverse problem: recover slowness from travel times
# m = (G^T G)^{-1} G^T d
# Hint: use np.linalg.inv() for matrix inverse, and @ for matrix multiplication
s_recovered = ???  # d = G @ s_true (matrix-vector multiply)

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.6f}")

Question: Why is the recovery perfect (error ~ 0%) with no noise? Would it still be perfect if we had only 3 rays instead of 6 (fewer rays than layers)?

Your answer here:

Exercise 2¶

Recap: With noisy data, the least-squares solution amplifies noise. Regularization adds a smoothness penalty: $\mathbf{s} = (\mathbf{G}^T\mathbf{G} + \alpha^2\mathbf{L}^T\mathbf{L})^{-1}\mathbf{G}^T\mathbf{d}$. Larger $\alpha$ = smoother model but worse data fit.: Effect of Noise and Smoothing (1D)

Real data always has noise. When we add noise to the travel times, the undamped inversion can produce wild velocity values. Regularization adds a smoothness penalty:

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

(a) Invert the noisy data without smoothing ($\alpha = 0$).

(b) Invert with smoothing for three values of $\alpha$.

In [ ]:

Copied!





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

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

print(f"Noise level: {noise_level} s")
print(f"Smoothing matrix L shape: {L1d.shape}")
print(f"L^T L shape: {LtL_1d.shape}")
# Add noise to the data
np.random.seed(42)
noise_level = 0.05  # seconds
d_noisy = d + noise_level * np.random.randn(n_rays)

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

print(f"Noise level: {noise_level} s")
print(f"Smoothing matrix L shape: {L1d.shape}")
print(f"L^T L shape: {LtL_1d.shape}")

In [ ]:

Copied!





# (a) Invert WITHOUT smoothing (same formula as Exercise 1, but with noisy data)
s_undamped = ???  # d = G @ s_true (matrix-vector multiply)

# (b) Invert WITH smoothing for three alpha values
alphas = [0.05, 0.5, 50.0]
s_smooth = {}
for alpha in alphas:
    # m = (G^T G + alpha^2 L^T L)^{-1} G^T d
    s_smooth[alpha] = ???  # (G^T G + alpha^2 L^T L)^{-1} G^T d

print("Done! Run the next cell to see the results.")
# (a) Invert WITHOUT smoothing (same formula as Exercise 1, but with noisy data)
s_undamped = ???  # d = G @ s_true (matrix-vector multiply)

# (b) Invert WITH smoothing for three alpha values
alphas = [0.05, 0.5, 50.0]
s_smooth = {}
for alpha in alphas:
    # m = (G^T G + alpha^2 L^T L)^{-1} G^T d
    s_smooth[alpha] = ???  # (G^T G + alpha^2 L^T L)^{-1} G^T d

print("Done! Run the next cell to see the results.")

In [ ]:

Copied!





# Plotting (just run this cell)
z_interfaces = np.concatenate([[0], np.cumsum(dz)])
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()
# Plotting (just run this cell)
z_interfaces = np.concatenate([[0], np.cumsum(dz)])
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()

Question: As $\alpha$ increases, the recovered model becomes smoother. But what happens to the data fit (how well $\mathbf{Gm}$ matches $\mathbf{d}$)? Why is there a trade-off between fitting the data and smoothness?

Your answer here:

Exercise 3¶

Recap: In 2D tomography, the same $\mathbf{d} = \mathbf{G}\mathbf{m}$ framework applies but $\mathbf{G}$ is larger and sparser. A checkerboard test reveals which parts of the model are well-resolved by the ray coverage.: 2D Checkerboard Test

In 2D, we divide the model into a grid of blocks. The textbook provides a 20$\times$20 grid with 118 rays. A checkerboard test checks how well the inversion can resolve structure: we create a checkerboard model, generate synthetic data, and see if the inversion can recover the pattern.

(a) Compute synthetic data from the checkerboard model (forward problem).

(b) Invert the data using the regularized formula.

(c) Compare results for three different $\alpha$ values.

In [ ]:

Copied!





# Load textbook G matrix (just run this cell)
nx2, nz2 = 20, 20
n_blocks2 = nx2 * nz2

gmat_raw = np.loadtxt('tomography_data/tomo_gmat.txt')
n_rays2 = int(gmat_raw[:, 0].max())
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

# Build 2D smoothing matrix
L2d = smoothing_matrix_2d(nx2, nz2)
LtL_2d = L2d.T @ L2d

print(f"Grid: {nx2}x{nz2} = {n_blocks2} blocks")
print(f"Rays: {n_rays2}")
print(f"G matrix: {G2.shape}")
print(f"L^T L: {LtL_2d.shape}")
# Load textbook G matrix (just run this cell)
nx2, nz2 = 20, 20
n_blocks2 = nx2 * nz2

gmat_raw = np.loadtxt('tomography_data/tomo_gmat.txt')
n_rays2 = int(gmat_raw[:, 0].max())
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

# Build 2D smoothing matrix
L2d = smoothing_matrix_2d(nx2, nz2)
LtL_2d = L2d.T @ L2d

print(f"Grid: {nx2}x{nz2} = {n_blocks2} blocks")
print(f"Rays: {n_rays2}")
print(f"G matrix: {G2.shape}")
print(f"L^T L: {LtL_2d.shape}")

In [ ]:

Copied!





# Build a checkerboard true model (just run this cell)
m_checker = 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_checker[iz, ix] = 1.0
        else:
            m_checker[iz, ix] = -1.0

plt.figure(figsize=(5, 5))
plt.imshow(m_checker, cmap='RdBu', aspect='equal', extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
plt.colorbar(label='Slowness perturbation')
plt.title('True checkerboard model')
plt.xlabel('Column')
plt.ylabel('Row')
plt.show()
# Build a checkerboard true model (just run this cell)
m_checker = 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_checker[iz, ix] = 1.0
        else:
            m_checker[iz, ix] = -1.0

plt.figure(figsize=(5, 5))
plt.imshow(m_checker, cmap='RdBu', aspect='equal', extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
plt.colorbar(label='Slowness perturbation')
plt.title('True checkerboard model')
plt.xlabel('Column')
plt.ylabel('Row')
plt.show()

In [ ]:

Copied!





# (a) Compute synthetic data: d = G @ m + noise
# Note: m_checker is 2D (20x20), but G expects a 1D vector (400,)
# Use .ravel() to flatten: m_checker.ravel() converts (20,20) → (400,)
noise_level = 5.0
np.random.seed(42)

d_checker = ???  # G2 @ m_checker.ravel() + noise

print(f"Data vector: {len(d_checker)} travel time perturbations")
# (a) Compute synthetic data: d = G @ m + noise
# Note: m_checker is 2D (20x20), but G expects a 1D vector (400,)
# Use .ravel() to flatten: m_checker.ravel() converts (20,20) → (400,)
noise_level = 5.0
np.random.seed(42)

d_checker = ???  # G2 @ m_checker.ravel() + noise

print(f"Data vector: {len(d_checker)} travel time perturbations")

In [ ]:

Copied!





# (b) Invert with a single alpha value
alpha = 1.0

# m = (G^T G + alpha^2 L^T L)^{-1} G^T d
m_recovered = ???  # d = G @ s_true (matrix-vector multiply)

# Plot true vs recovered
fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)

axes[0].imshow(m_checker, cmap='RdBu', aspect='equal',
               extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True checkerboard')

axes[1].imshow(m_recovered.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
               extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[1].set_title(f'Recovered ($\\alpha$={alpha})')

for ax in axes:
    ax.set_xlabel('Column')
    ax.set_ylabel('Row')
plt.suptitle('Checkerboard test', fontsize=13)
plt.show()
# (b) Invert with a single alpha value
alpha = 1.0

# m = (G^T G + alpha^2 L^T L)^{-1} G^T d
m_recovered = ???  # d = G @ s_true (matrix-vector multiply)

# Plot true vs recovered
fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)

axes[0].imshow(m_checker, cmap='RdBu', aspect='equal',
               extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[0].set_title('True checkerboard')

axes[1].imshow(m_recovered.reshape(nz2, nx2), cmap='RdBu', aspect='equal',
               extent=[0, nx2, nz2, 0], vmin=-1, vmax=1)
axes[1].set_title(f'Recovered ($\\alpha$={alpha})')

for ax in axes:
    ax.set_xlabel('Column')
    ax.set_ylabel('Row')
plt.suptitle('Checkerboard test', fontsize=13)
plt.show()

In [ ]:

Copied!





# (c) Compare three different alpha values
alpha_values = [0.005, 5.0, 50.0]

fig, axes = plt.subplots(1, 4, figsize=(18, 4), constrained_layout=True)

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

for ax, alpha in zip(axes[1:], alpha_values):
    # Invert with this alpha
    m_rec = ???  # regularized inversion

    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(f'Effect of smoothing (noise = {noise_level})', fontsize=13)
plt.show()
# (c) Compare three different alpha values
alpha_values = [0.005, 5.0, 50.0]

fig, axes = plt.subplots(1, 4, figsize=(18, 4), constrained_layout=True)

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

for ax, alpha in zip(axes[1:], alpha_values):
    # Invert with this alpha
    m_rec = ???  # regularized inversion

    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(f'Effect of smoothing (noise = {noise_level})', fontsize=13)
plt.show()

Question: Which parts of the 20$\times$20 grid are best resolved (show a clear checkerboard pattern)? Which parts are poorly resolved? Why? (Hint: think about where the rays cross.)

Your answer here:

Exercise 4: Invert Real Data¶

Now we use real travel time observations from the textbook. The same $\mathbf{G}$ matrix applies, but instead of synthetic data, we load measured travel time perturbations $\delta \mathbf{d}$.

(a) Invert the real data with $\alpha = 1.0$.

(b) Try four different $\alpha$ values and compare the results.

In [ ]:

Copied!





# Load real travel time data (just run this cell)
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")
# Load real travel time data (just run this cell)
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")

In [ ]:

Copied!





# (a) Invert the real data with alpha = 1.0
alpha = 1.0

m_real = ???  # same regularized inversion formula

plt.figure(figsize=(6, 5))
plt.imshow(m_real.reshape(nz2, nx2).T, cmap='RdBu', aspect='equal',
           extent=[0, nx2, nz2, 0])
plt.colorbar(label='Slowness perturbation')
plt.title(f'Real data inversion ($\\alpha$ = {alpha})')
plt.xlabel('Column')
plt.ylabel('Row')
plt.show()
# (a) Invert the real data with alpha = 1.0
alpha = 1.0

m_real = ???  # same regularized inversion formula

plt.figure(figsize=(6, 5))
plt.imshow(m_real.reshape(nz2, nx2).T, cmap='RdBu', aspect='equal',
           extent=[0, nx2, nz2, 0])
plt.colorbar(label='Slowness perturbation')
plt.title(f'Real data inversion ($\\alpha$ = {alpha})')
plt.xlabel('Column')
plt.ylabel('Row')
plt.show()

In [ ]:

Copied!





# (b) Compare four different alpha values
alpha_values = [0.1, 1.0, 5.0, 20.0]

fig, axes = plt.subplots(1, 4, figsize=(18, 4), constrained_layout=True)

for ax, alpha in zip(axes, alpha_values):
    m_rec = ???  # regularized inversion

    vmax = np.percentile(np.abs(m_rec), 95)
    im = ax.imshow(m_rec.reshape(nz2, nx2).T, cmap='RdBu', aspect='equal',
                   extent=[0, nx2, nz2, 0], vmin=-vmax, vmax=vmax)
    ax.set_title(f'$\\alpha$ = {alpha}')
    ax.set_xlabel('Column')

axes[0].set_ylabel('Row')
plt.colorbar(im, ax=axes, label='Slowness perturbation', shrink=0.8)
fig.suptitle('Real data inversion: effect of smoothing', fontsize=13)
plt.show()
# (b) Compare four different alpha values
alpha_values = [0.1, 1.0, 5.0, 20.0]

fig, axes = plt.subplots(1, 4, figsize=(18, 4), constrained_layout=True)

for ax, alpha in zip(axes, alpha_values):
    m_rec = ???  # regularized inversion

    vmax = np.percentile(np.abs(m_rec), 95)
    im = ax.imshow(m_rec.reshape(nz2, nx2).T, cmap='RdBu', aspect='equal',
                   extent=[0, nx2, nz2, 0], vmin=-vmax, vmax=vmax)
    ax.set_title(f'$\\alpha$ = {alpha}')
    ax.set_xlabel('Column')

axes[0].set_ylabel('Row')
plt.colorbar(im, ax=axes, label='Slowness perturbation', shrink=0.8)
fig.suptitle('Real data inversion: effect of smoothing', fontsize=13)
plt.show()

Question: Which features appear consistently across different $\alpha$ values? Those are the most robust. Which features change a lot with $\alpha$?

Your answer here:

Exercise 5: Interpreting Global Tomography¶

The cells below plot depth slices from two real tomography models:

S362ANI: global shear-wave velocity ($\delta V_s / V_s$)
MITPS-20: P-wave velocity beneath North America ($\delta V_p / V_p$)

Just run the plotting cells, then answer the interpretation questions below.

Recall: Blue = fast (cold material), Red = slow (hot material).

In [ ]:

Copied!





# Load S362ANI global model (just run this cell)
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)

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):
    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())

print(f"S362ANI: {len(lat)} lat x {len(lon)} lon x {len(depth)} depths")
print(f"Depths (km): {depth}")
# Load S362ANI global model (just run this cell)
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)

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):
    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())

print(f"S362ANI: {len(lat)} lat x {len(lon)} lon x {len(depth)} depths")
print(f"Depths (km): {depth}")

In [ ]:

Copied!





# Plot S362ANI depth slices (just run this cell)
depth_slices = [100, 250, 600, 2800]
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()
# Plot S362ANI depth slices (just run this cell)
depth_slices = [100, 250, 600, 2800]
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()

(a) Look at the S362ANI map at 100 km depth. Where are the fast (blue) anomalies? Where are the slow (red) anomalies? Why are continents generally fast at this depth?

Your answer here:

(b) Look at the S362ANI map at 600 km depth. Can you identify any fast anomalies that might be subducting slabs? Where are they?

Your answer here:

(c) Look at the S362ANI map at 2800 km depth (near the core-mantle boundary). Describe the large-scale pattern. Where are the two large slow regions (LLSVPs)?

Your answer here: