Earthquake Location — Exercises¶

Run in Google Colab

These exercises reinforce the key concepts from the earthquake location lecture:

1. Forward problem: predicting arrival times from a source location
2. Jacobian and linearization: setting up the inverse problem
3. Geiger's method: iterative least-squares location
4. Station geometry: how it controls location uncertainty

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

np.random.seed(42)
plt.rcParams["figure.figsize"] = (10, 6)
plt.rcParams["font.size"] = 12

# Plotting utility — draws an error ellipse from a 2x2 covariance matrix
def plot_error_ellipse(ax, Cm_2d, center, n_std=2, **kwargs):
    eigvals, eigvecs = np.linalg.eigh(Cm_2d)
    angle = np.degrees(np.arctan2(eigvecs[1, 1], eigvecs[0, 1]))
    width  = 2 * n_std * np.sqrt(max(eigvals[1], 0))
    height = 2 * n_std * np.sqrt(max(eigvals[0], 0))
    ell = Ellipse(xy=center, width=width, height=height, angle=angle, **kwargs)
    ax.add_patch(ell)

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Ellipse np.random.seed(42) plt.rcParams["figure.figsize"] = (10, 6) plt.rcParams["font.size"] = 12 # Plotting utility — draws an error ellipse from a 2x2 covariance matrix def plot_error_ellipse(ax, Cm_2d, center, n_std=2, **kwargs): eigvals, eigvecs = np.linalg.eigh(Cm_2d) angle = np.degrees(np.arctan2(eigvecs[1, 1], eigvecs[0, 1])) width = 2 * n_std * np.sqrt(max(eigvals[1], 0)) height = 2 * n_std * np.sqrt(max(eigvals[0], 0)) ell = Ellipse(xy=center, width=width, height=height, angle=angle, **kwargs) ax.add_patch(ell)

---

Exercise 1: The Forward Problem¶

The forward function predicts the arrival time at station $i$ given a source at $(x_0, z_0)$ with origin time $t_0$ in a medium with uniform velocity $v$:

$$\hat{t}^i = \frac{\sqrt{(x^i - x_0)^2 + (z^i - z_0)^2}}{v} + t_0$$

(a) Implement the forward function below.

In [ ]:

def forward(x0, z0, t0, v, sx, sz):
    """Predict arrival time at each station.
    
    Parameters
    ----------
    x0, z0 : float — source position (km)
    t0 : float — origin time (s)
    v : float — velocity (km/s)
    sx, sz : arrays — station positions (km)
    
    Returns
    -------
    t_pred : array — predicted arrival times (s)
    """
    # YOUR CODE HERE
    dist = ???
    t_pred = ???
    return t_pred

def forward(x0, z0, t0, v, sx, sz): """Predict arrival time at each station. Parameters ---------- x0, z0 : float — source position (km) t0 : float — origin time (s) v : float — velocity (km/s) sx, sz : arrays — station positions (km) Returns ------- t_pred : array — predicted arrival times (s) """ # YOUR CODE HERE dist = ??? t_pred = ??? return t_pred

(b) Use your forward function to generate synthetic data. Compute the true arrival times, then add Gaussian noise to simulate picking errors.

In [ ]:

# Physical setup
v_true = 6.0       # km/s, uniform velocity
sigma_true = 0.05  # picking noise standard deviation (s)

# True earthquake location (unknown to us in practice)
x0_true, z0_true, t0_true = 5.0, 8.0, 0.0

# 8 stations at the surface
station_x = np.array([-5, -2, 0, 3, 6, 10, 13, 15], dtype=float)
station_z = np.zeros_like(station_x)
n_stations = len(station_x)

# Compute true arrival times using your forward function
t_true = forward(x0_true, z0_true, t0_true, v_true, station_x, station_z)

# Add random noise to simulate picking errors
t_obs = t_true + np.random.normal(0, sigma_true, n_stations)

# Plot the geometry and travel time curve
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

ax = axes[0]
ax.scatter(station_x, station_z, marker="v", s=150, c="blue", zorder=5, label="Stations")
ax.scatter(x0_true, z0_true, marker="*", s=300, c="red", zorder=5, label="Earthquake")
for i in range(n_stations):
    ax.plot([station_x[i], x0_true], [0, z0_true], "k--", alpha=0.3)
ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)")
ax.set_title("2D Geometry"); ax.legend(loc="upper right")
ax.set_ylim(12, -2); ax.set_aspect("equal")

ax = axes[1]
ax.plot(station_x, t_true, "ro-", label="True arrival times")
ax.errorbar(station_x, t_obs, yerr=sigma_true, fmt="bs", capsize=3, label=f"Observed (sigma={sigma_true}s)")
ax.set_xlabel("Station x (km)"); ax.set_ylabel("Arrival time (s)")
ax.set_title("Travel Time Curve"); ax.legend()

plt.tight_layout(); plt.show()

# Physical setup v_true = 6.0 # km/s, uniform velocity sigma_true = 0.05 # picking noise standard deviation (s) # True earthquake location (unknown to us in practice) x0_true, z0_true, t0_true = 5.0, 8.0, 0.0 # 8 stations at the surface station_x = np.array([-5, -2, 0, 3, 6, 10, 13, 15], dtype=float) station_z = np.zeros_like(station_x) n_stations = len(station_x) # Compute true arrival times using your forward function t_true = forward(x0_true, z0_true, t0_true, v_true, station_x, station_z) # Add random noise to simulate picking errors t_obs = t_true + np.random.normal(0, sigma_true, n_stations) # Plot the geometry and travel time curve fig, axes = plt.subplots(1, 2, figsize=(14, 5)) ax = axes[0] ax.scatter(station_x, station_z, marker="v", s=150, c="blue", zorder=5, label="Stations") ax.scatter(x0_true, z0_true, marker="*", s=300, c="red", zorder=5, label="Earthquake") for i in range(n_stations): ax.plot([station_x[i], x0_true], [0, z0_true], "k--", alpha=0.3) ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)") ax.set_title("2D Geometry"); ax.legend(loc="upper right") ax.set_ylim(12, -2); ax.set_aspect("equal") ax = axes[1] ax.plot(station_x, t_true, "ro-", label="True arrival times") ax.errorbar(station_x, t_obs, yerr=sigma_true, fmt="bs", capsize=3, label=f"Observed (sigma={sigma_true}s)") ax.set_xlabel("Station x (km)"); ax.set_ylabel("Arrival time (s)") ax.set_title("Travel Time Curve"); ax.legend() plt.tight_layout(); plt.show()

Question: The travel time curve has a characteristic "U" shape. Which station has the shortest arrival time, and why? What does the minimum of the curve tell you about the earthquake's horizontal position?

Your answer here:

---

Exercise 2: The Jacobian Matrix¶

To solve the nonlinear location problem with linear algebra, we linearize the forward function around an initial guess $\mathbf{m}_0 = (x_0, z_0, t_0)$:

$$\hat{t}^i(\mathbf{m}_0 + \Delta\mathbf{m}) \approx \hat{t}^i(\mathbf{m}_0) + \frac{\partial \hat{t}^i}{\partial x_0}\Delta x_0 + \frac{\partial \hat{t}^i}{\partial z_0}\Delta z_0 + \frac{\partial \hat{t}^i}{\partial t_0}\Delta t_0$$

The partial derivatives are:

$$\frac{\partial \hat{t}^i}{\partial x_0} = \frac{-(x^i - x_0)}{v \cdot d^i}, \quad \frac{\partial \hat{t}^i}{\partial z_0} = \frac{-(z^i - z_0)}{v \cdot d^i}, \quad \frac{\partial \hat{t}^i}{\partial t_0} = 1$$

where $d^i = \sqrt{(x^i - x_0)^2 + (z^i - z_0)^2}$ is the source-station distance.

(a) Implement the Jacobian matrix $\mathbf{G}$ (size $n_\text{stations} \times 3$).

In [ ]:

def compute_jacobian(x0, z0, v, sx, sz):
    """Jacobian matrix G (n_stations x 3).
    
    Columns: [dt/dx0, dt/dz0, dt/dt0]
    """
    dist = np.sqrt((sx - x0)**2 + (sz - z0)**2)
    dtdx  = ???  # partial derivative with respect to x0
    dtdz  = ???  # partial derivative with respect to z0
    dtdt0 = ???  # partial derivative with respect to t0
    return np.column_stack([dtdx, dtdz, dtdt0])

def compute_jacobian(x0, z0, v, sx, sz): """Jacobian matrix G (n_stations x 3). Columns: [dt/dx0, dt/dz0, dt/dt0] """ dist = np.sqrt((sx - x0)**2 + (sz - z0)**2) dtdx = ??? # partial derivative with respect to x0 dtdz = ??? # partial derivative with respect to z0 dtdt0 = ??? # partial derivative with respect to t0 return np.column_stack([dtdx, dtdz, dtdt0])

(b) Evaluate the Jacobian at an initial guess and inspect its columns. Run the cell below and answer the question.

In [ ]:

# Evaluate at an initial guess (deliberately wrong)
x0_init, z0_init, t0_init = 2.0, 5.0, 0.5
G = compute_jacobian(x0_init, z0_init, v_true, station_x, station_z)

print("Jacobian G (n_stations x 3): columns = [dt/dx0, dt/dz0, dt/dt0]")
print("-" * 65)
for i in range(n_stations):
    print(f"  Station {i} (x={station_x[i]:5.1f}): [{G[i,0]:+.4f}, {G[i,1]:+.4f}, {G[i,2]:+.0f}]")

# Evaluate at an initial guess (deliberately wrong) x0_init, z0_init, t0_init = 2.0, 5.0, 0.5 G = compute_jacobian(x0_init, z0_init, v_true, station_x, station_z) print("Jacobian G (n_stations x 3): columns = [dt/dx0, dt/dz0, dt/dt0]") print("-" * 65) for i in range(n_stations): print(f" Station {i} (x={station_x[i]:5.1f}): [{G[i,0]:+.4f}, {G[i,1]:+.4f}, {G[i,2]:+.0f}]")

Question: Look at column 1 (dt/dx0). Does the sign change across stations? Why does this happen physically?

Your answer here:

Question: Look at column 2 (dt/dz0). Does the sign change across stations? Why or why not?

Your answer here:

Question: Based on these observations, which parameter ($x_0$ or $z_0$) do you expect to be better constrained by the data, and why?

Your answer here:

---

Exercise 3: Geiger's Method¶

Geiger's method iteratively solves the linearized location problem:

1. Compute predicted times at current guess: $\hat{\mathbf{t}} = f(\mathbf{m}_0)$
2. Compute residuals: $\mathbf{r} = \mathbf{t}_\text{obs} - \hat{\mathbf{t}}$
3. Compute Jacobian: $\mathbf{G}$ at $\mathbf{m}_0$
4. Solve for update: $\Delta\mathbf{m} = (\mathbf{G}^T\mathbf{G})^{-1}\mathbf{G}^T\mathbf{r}$
5. Update: $\mathbf{m}_0 \leftarrow \mathbf{m}_0 + \Delta\mathbf{m}$
6. Repeat until $|\Delta\mathbf{m}|$ is small

(a) Complete the Geiger's method function below.

In [ ]:

def geiger_method(x0, z0, t0, v, sx, sz, t_obs, max_iter=20, tol=1e-6):
    """Geiger's method for earthquake location.
    
    Returns: x0, z0, t0, G_final, history
    """
    history = [(x0, z0, t0)]
    for it in range(max_iter):
        # Step 1: predict arrival times at current guess
        t_pred = ???
        
        # Step 2: compute residuals (observed - predicted)
        r = ???
        
        # Step 3: compute Jacobian at current guess
        G = ???
        
        # Step 4: solve least squares for model update
        # Hint: use np.linalg.lstsq(G, r, rcond=None)[0]
        dm = ???
        
        # Step 5: update the guess
        x0 += dm[0]
        z0 += dm[1]
        t0 += dm[2]
        history.append((x0, z0, t0))
        
        rms = np.sqrt(np.mean(r**2))
        print(f"  Iter {it+1}: x={x0:.4f}, z={z0:.4f}, t0={t0:.4f}, RMS={rms:.6f} s")
        
        # Step 6: check convergence
        if np.linalg.norm(dm) < tol:
            print("  Converged!")
            break
    
    return x0, z0, t0, G, np.array(history)

def geiger_method(x0, z0, t0, v, sx, sz, t_obs, max_iter=20, tol=1e-6): """Geiger's method for earthquake location. Returns: x0, z0, t0, G_final, history """ history = [(x0, z0, t0)] for it in range(max_iter): # Step 1: predict arrival times at current guess t_pred = ??? # Step 2: compute residuals (observed - predicted) r = ??? # Step 3: compute Jacobian at current guess G = ??? # Step 4: solve least squares for model update # Hint: use np.linalg.lstsq(G, r, rcond=None)[0] dm = ??? # Step 5: update the guess x0 += dm[0] z0 += dm[1] t0 += dm[2] history.append((x0, z0, t0)) rms = np.sqrt(np.mean(r**2)) print(f" Iter {it+1}: x={x0:.4f}, z={z0:.4f}, t0={t0:.4f}, RMS={rms:.6f} s") # Step 6: check convergence if np.linalg.norm(dm) < tol: print(" Converged!") break return x0, z0, t0, G, np.array(history)

(b) Run Geiger's method starting from the initial guess $(x_0, z_0, t_0) = (2, 5, 0.5)$ and plot the convergence path.

In [ ]:

# Run Geiger's method
x0_est, z0_est, t0_est, G_final, history = geiger_method(
    x0_init, z0_init, t0_init, v_true, station_x, station_z, t_obs)

print(f"\nTrue:      x={x0_true}, z={z0_true}, t0={t0_true}")
print(f"Estimated: x={x0_est:.4f}, z={z0_est:.4f}, t0={t0_est:.4f}")

# Plot convergence path
fig = plt.figure(figsize=(14, 6))
gs = fig.add_gridspec(1, 2, width_ratios=[2, 1])

ax = fig.add_subplot(gs[0])
ax.scatter(station_x, np.zeros(n_stations), marker="v", s=100, c="blue", zorder=5, label="Stations")
ax.scatter(x0_true, z0_true, marker="*", s=300, c="red", zorder=5, label="True")
ax.plot(history[:, 0], history[:, 1], "go-", markersize=8, label="Iterations")
ax.scatter(history[0, 0], history[0, 1], s=150, c="green", edgecolors="k", zorder=6, label="Initial guess")
ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)")
ax.set_title("Geiger's Method: Convergence Path")
ax.legend(fontsize=9); ax.set_ylim(12, -2); ax.set_aspect("equal")

ax2 = fig.add_subplot(gs[1])
rms = [np.sqrt(np.mean((t_obs - forward(h[0], h[1], h[2], v_true, station_x, station_z))**2)) for h in history]
ax2.plot(range(len(rms)), rms, "ko-")
ax2.set_xlabel("Iteration"); ax2.set_ylabel("RMS residual (s)")
ax2.set_title("RMS Convergence")

plt.tight_layout(); plt.show()

# Run Geiger's method x0_est, z0_est, t0_est, G_final, history = geiger_method( x0_init, z0_init, t0_init, v_true, station_x, station_z, t_obs) print(f"\nTrue: x={x0_true}, z={z0_true}, t0={t0_true}") print(f"Estimated: x={x0_est:.4f}, z={z0_est:.4f}, t0={t0_est:.4f}") # Plot convergence path fig = plt.figure(figsize=(14, 6)) gs = fig.add_gridspec(1, 2, width_ratios=[2, 1]) ax = fig.add_subplot(gs[0]) ax.scatter(station_x, np.zeros(n_stations), marker="v", s=100, c="blue", zorder=5, label="Stations") ax.scatter(x0_true, z0_true, marker="*", s=300, c="red", zorder=5, label="True") ax.plot(history[:, 0], history[:, 1], "go-", markersize=8, label="Iterations") ax.scatter(history[0, 0], history[0, 1], s=150, c="green", edgecolors="k", zorder=6, label="Initial guess") ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)") ax.set_title("Geiger's Method: Convergence Path") ax.legend(fontsize=9); ax.set_ylim(12, -2); ax.set_aspect("equal") ax2 = fig.add_subplot(gs[1]) rms = [np.sqrt(np.mean((t_obs - forward(h[0], h[1], h[2], v_true, station_x, station_z))**2)) for h in history] ax2.plot(range(len(rms)), rms, "ko-") ax2.set_xlabel("Iteration"); ax2.set_ylabel("RMS residual (s)") ax2.set_title("RMS Convergence") plt.tight_layout(); plt.show()

Question: How many iterations does Geiger's method need to converge? Does the final RMS residual reach zero? Why or why not?

Your answer here:

---

Exercise 4: Error Ellipse and Station Geometry¶

At the solution, the model covariance matrix quantifies the uncertainty:

$$\mathbf{C}_m = (\mathbf{G}^T \mathbf{G})^{-1} \, \sigma^2$$

The diagonal elements give the variance of each parameter: $\sigma_x^2 = C_{m,11}$, $\sigma_z^2 = C_{m,22}$.

(a) Compute the covariance matrix and the uncertainty in each parameter. Use the estimated sigma from the residuals:

$$\hat{\sigma}^2 = \frac{\sum r_i^2}{n_\text{obs} - n_\text{params}}$$

In [ ]:

# Jacobian at the best-fit location (already computed as G_final)
r_final = t_obs - forward(x0_est, z0_est, t0_est, v_true, station_x, station_z)
n_params = 3  # x, z, t0

# Estimate sigma from the residuals
sigma_est = ???  # sqrt(sum(r^2) / (n_obs - n_params))

# Compute the covariance matrix: C_m = (G^T G)^{-1} * sigma^2
Cm = ???

print(f"Estimated sigma: {sigma_est:.4f} s  (true: {sigma_true} s)")
print(f"sigma_x  = {np.sqrt(Cm[0,0]):.4f} km")
print(f"sigma_z  = {np.sqrt(Cm[1,1]):.4f} km")
print(f"sigma_t0 = {np.sqrt(Cm[2,2]):.4f} s")

# Jacobian at the best-fit location (already computed as G_final) r_final = t_obs - forward(x0_est, z0_est, t0_est, v_true, station_x, station_z) n_params = 3 # x, z, t0 # Estimate sigma from the residuals sigma_est = ??? # sqrt(sum(r^2) / (n_obs - n_params)) # Compute the covariance matrix: C_m = (G^T G)^{-1} * sigma^2 Cm = ??? print(f"Estimated sigma: {sigma_est:.4f} s (true: {sigma_true} s)") print(f"sigma_x = {np.sqrt(Cm[0,0]):.4f} km") print(f"sigma_z = {np.sqrt(Cm[1,1]):.4f} km") print(f"sigma_t0 = {np.sqrt(Cm[2,2]):.4f} s")

(b) Plot the error ellipse at the estimated location. The code below extracts the x-z block of the covariance matrix and uses the provided plot_error_ellipse function.

In [ ]:

# Extract the x-z block of the covariance matrix (ignore t0 for the 2D ellipse)
Cm_xz = Cm[:2, :2].copy()

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

ax = axes[0]
ax.scatter(station_x, np.zeros(n_stations), marker="v", s=100, c="blue", zorder=5, label="Stations")
ax.scatter(x0_true, z0_true, marker="*", s=300, c="red", zorder=5, label="True")
ax.scatter(x0_est, z0_est, marker="x", s=200, c="green", lw=3, zorder=5, label="Estimated")
for ns, alpha in [(1, 0.4), (2, 0.2)]:
    plot_error_ellipse(ax, Cm_xz, (x0_est, z0_est), n_std=ns,
                       fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2, label=f"{ns}σ ellipse")
ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)")
ax.set_title("Error Ellipse — Overview"); ax.legend(fontsize=9); ax.set_ylim(12, -2); ax.set_aspect("equal")

ax = axes[1]
ax.scatter(x0_true, z0_true, marker="*", s=500, c="red", zorder=5, label="True")
ax.scatter(x0_est, z0_est, marker="x", s=300, c="green", lw=3, zorder=5, label="Estimated")
for ns, alpha in [(1, 0.4), (2, 0.2)]:
    plot_error_ellipse(ax, Cm_xz, (x0_est, z0_est), n_std=ns,
                       fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2, label=f"{ns}σ")
margin = max(2 * np.sqrt(Cm_xz[0, 0]), 2 * np.sqrt(Cm_xz[1, 1]), 0.5)
ax.set_xlim(x0_est - margin, x0_est + margin)
ax.set_ylim(z0_est + margin, z0_est - margin)
ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)")
ax.set_title("Error Ellipse — Zoom"); ax.legend(fontsize=9); ax.set_aspect("equal")

plt.tight_layout(); plt.show()

# Extract the x-z block of the covariance matrix (ignore t0 for the 2D ellipse) Cm_xz = Cm[:2, :2].copy() fig, axes = plt.subplots(1, 2, figsize=(14, 6)) ax = axes[0] ax.scatter(station_x, np.zeros(n_stations), marker="v", s=100, c="blue", zorder=5, label="Stations") ax.scatter(x0_true, z0_true, marker="*", s=300, c="red", zorder=5, label="True") ax.scatter(x0_est, z0_est, marker="x", s=200, c="green", lw=3, zorder=5, label="Estimated") for ns, alpha in [(1, 0.4), (2, 0.2)]: plot_error_ellipse(ax, Cm_xz, (x0_est, z0_est), n_std=ns, fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2, label=f"{ns}σ ellipse") ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)") ax.set_title("Error Ellipse — Overview"); ax.legend(fontsize=9); ax.set_ylim(12, -2); ax.set_aspect("equal") ax = axes[1] ax.scatter(x0_true, z0_true, marker="*", s=500, c="red", zorder=5, label="True") ax.scatter(x0_est, z0_est, marker="x", s=300, c="green", lw=3, zorder=5, label="Estimated") for ns, alpha in [(1, 0.4), (2, 0.2)]: plot_error_ellipse(ax, Cm_xz, (x0_est, z0_est), n_std=ns, fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2, label=f"{ns}σ") margin = max(2 * np.sqrt(Cm_xz[0, 0]), 2 * np.sqrt(Cm_xz[1, 1]), 0.5) ax.set_xlim(x0_est - margin, x0_est + margin) ax.set_ylim(z0_est + margin, z0_est - margin) ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)") ax.set_title("Error Ellipse — Zoom"); ax.legend(fontsize=9); ax.set_aspect("equal") plt.tight_layout(); plt.show()

(c) Now explore how station geometry affects the error ellipse. For each configuration below, compute the Jacobian, covariance matrix, and plot the error ellipse. Fill in the ??? to compute the covariance matrix for each configuration.

In [ ]:

configs = [
    ("Good coverage (8 stations)",      np.array([-5.,-2.,0.,3.,6.,10.,13.,15.]), np.zeros(8)),
    ("Azimuthal gap (all one side)",     np.array([10.,12.,14.,16.,18.,20.,22.,24.]), np.zeros(8)),
    ("Few stations (4)",                 np.array([-2.,3.,8.,13.]),                  np.zeros(4)),
    ("With borehole sensor at depth",    np.array([-5.,-2.,0.,3.,6.,10.,13.,5.]),    np.array([0.,0.,0.,0.,0.,0.,0.,10.])),
]

for title, sx, sz in configs:
    # Compute Jacobian at the TRUE location for this station configuration
    G_cfg = compute_jacobian(x0_true, z0_true, v_true, sx, sz)
    
    # Compute the covariance matrix using sigma_true
    # YOUR CODE HERE
    Cm_cfg = ???
    
    Cm_xz_cfg = Cm_cfg[:2, :2]
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax = axes[0]
    ax.scatter(sx, sz, marker="v", s=100, c="blue", zorder=5, label="Stations")
    ax.scatter(x0_true, z0_true, marker="*", s=200, c="red", zorder=5, label="True")
    for ns, alpha in [(1, 0.4), (2, 0.2)]:
        plot_error_ellipse(ax, Cm_xz_cfg, (x0_true, z0_true), n_std=ns,
                           fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2)
    ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)")
    ax.set_title(f"{title} — Overview"); ax.set_ylim(14, -2); ax.set_aspect("equal"); ax.legend(fontsize=9)
    
    ax = axes[1]
    ax.scatter(x0_true, z0_true, marker="*", s=400, c="red", zorder=5, label="True")
    for ns, alpha in [(1, 0.4), (2, 0.2)]:
        plot_error_ellipse(ax, Cm_xz_cfg, (x0_true, z0_true), n_std=ns,
                           fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2, label=f"{ns}σ")
    ax.set_xlim(0, 10); ax.set_ylim(13, 3)
    ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)")
    ax.set_title(f"Zoom: σ_x={np.sqrt(Cm_cfg[0,0]):.3f}, σ_z={np.sqrt(Cm_cfg[1,1]):.3f} km")
    ax.legend(); ax.set_aspect("equal")
    
    plt.tight_layout(); plt.show()

configs = [ ("Good coverage (8 stations)", np.array([-5.,-2.,0.,3.,6.,10.,13.,15.]), np.zeros(8)), ("Azimuthal gap (all one side)", np.array([10.,12.,14.,16.,18.,20.,22.,24.]), np.zeros(8)), ("Few stations (4)", np.array([-2.,3.,8.,13.]), np.zeros(4)), ("With borehole sensor at depth", np.array([-5.,-2.,0.,3.,6.,10.,13.,5.]), np.array([0.,0.,0.,0.,0.,0.,0.,10.])), ] for title, sx, sz in configs: # Compute Jacobian at the TRUE location for this station configuration G_cfg = compute_jacobian(x0_true, z0_true, v_true, sx, sz) # Compute the covariance matrix using sigma_true # YOUR CODE HERE Cm_cfg = ??? Cm_xz_cfg = Cm_cfg[:2, :2] fig, axes = plt.subplots(1, 2, figsize=(14, 5)) ax = axes[0] ax.scatter(sx, sz, marker="v", s=100, c="blue", zorder=5, label="Stations") ax.scatter(x0_true, z0_true, marker="*", s=200, c="red", zorder=5, label="True") for ns, alpha in [(1, 0.4), (2, 0.2)]: plot_error_ellipse(ax, Cm_xz_cfg, (x0_true, z0_true), n_std=ns, fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2) ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)") ax.set_title(f"{title} — Overview"); ax.set_ylim(14, -2); ax.set_aspect("equal"); ax.legend(fontsize=9) ax = axes[1] ax.scatter(x0_true, z0_true, marker="*", s=400, c="red", zorder=5, label="True") for ns, alpha in [(1, 0.4), (2, 0.2)]: plot_error_ellipse(ax, Cm_xz_cfg, (x0_true, z0_true), n_std=ns, fill=True, facecolor="green", alpha=alpha, edgecolor="green", lw=2, label=f"{ns}σ") ax.set_xlim(0, 10); ax.set_ylim(13, 3) ax.set_xlabel("x (km)"); ax.set_ylabel("Depth (km)") ax.set_title(f"Zoom: σ_x={np.sqrt(Cm_cfg[0,0]):.3f}, σ_z={np.sqrt(Cm_cfg[1,1]):.3f} km") ax.legend(); ax.set_aspect("equal") plt.tight_layout(); plt.show()

Question: When all stations are on one side (azimuthal gap), how does the error ellipse change compared to good coverage?

Your answer here:

Question: Why does the borehole sensor dramatically improve the depth uncertainty?

Your answer here:

Question: In real seismology, why is it hard to constrain earthquake depth using only surface stations?

Your answer here: