Ray Tracing in Cartesian Coordinates#
Colab note: This notebook is designed to run on Google Colab. The first code cell installs dependencies.
Learning Objectives:
Implement Snell’s Law in 1D layered media
Compute ray parameters and takeoff angles
Calculate travel times for direct and refracted waves
Visualize ray paths in flat-layered models
Prerequisites: Basic Python, Trigonometry
Reference: Shearer, Chapter 4 (Ray Theory)
Notebook Outline:
This notebook implements the II. practicum: computing travel-time fields and ray paths in 2-D cross sections using the fast marching method (FMM) solver in pykonal, then exploring how layering, curvature, and low-velocity zones affect:
# Test pykonal installation
try:
import h5py
print(f"✓ h5py version: {h5py.__version__}")
except ImportError as e:
print(f"✗ h5py not found: {e}")
print("Run: pixi install")
try:
import pykonal
print(f"✓ pykonal version: {pykonal.__version__}")
print("✓ pykonal successfully imported!")
except ImportError as e:
print(f"✗ pykonal not found: {e}")
print("Run: pixi install")
np.infty = np.inf
0. Install dependencies (Colab)#
We install pykonal plus the usual scientific stack.
If the PyPI wheel install fails on Colab (rare, but can happen if compilation is needed), the cell includes a fallback to install from GitHub.
# # Colab install
# !pip -q install --upgrade pip
# !pip -q install numpy scipy matplotlib
# # Try PyPI first (fastest)
# try:
# import pykonal # noqa: F401
# print("pykonal already available")
# except Exception:
# !pip -q install pykonal
# #
# # Fallback (if needed): install a specific tag from GitHub
# try:
# import pykonal
# print("pykonal version:", getattr(pykonal, "__version__", "unknown"))
# except Exception as e:
# print("PyPI install failed, trying GitHub install. Error:", e)
# !pip -q install cython
# !pip -q install "git+https://github.com/malcolmw/pykonal@0.4.1"
# import pykonal
# print("pykonal version:", getattr(pykonal, "__version__", "unknown"))
# !pip install pykonal
1. Imports and plotting helpers#
import numpy as np
import matplotlib.pyplot as plt
import pykonal
from dataclasses import dataclass
from typing import Tuple, Dict, Optional
plt.rcParams.update({
"figure.dpi": 140,
"savefig.dpi": 200,
"axes.spines.top": False,
"axes.spines.right": False,
})
2. Model setup#
We work on a 2-D grid in ((x,z)) with depth (z) positive downward.
PyKonal is formally 3-D, so we represent 2-D by setting ny=1 and embedding ((x,z)) into ((x,y,z)) with (y=0).
Coordinate conventions#
We store velocity on a regular grid with spacing (\Delta x, \Delta z) (km).
The source is a point at ((x_s, z_s)).
Receivers lie on the surface (z=0) for travel-time curve experiments.
What the solver returns#
solver.traveltime.values: a 3-D array of (T(x,y,z)) on the gridsolver.trace_ray(end): ray path as an array of points ((x,y,z)), found by steepest descent on (T)
@dataclass
class Grid2D:
x_min: float
z_min: float
dx: float
dz: float
nx: int
nz: int
@property
def x(self) -> np.ndarray:
return self.x_min + self.dx * np.arange(self.nx)
@property
def z(self) -> np.ndarray:
return self.z_min + self.dz * np.arange(self.nz)
def mesh(self) -> Tuple[np.ndarray, np.ndarray]:
X, Z = np.meshgrid(self.x, self.z, indexing="ij") # (nx,nz)
return X, Z
def make_solver_from_velocity(grid: Grid2D, v_xz: np.ndarray) -> pykonal.EikonalSolver:
"""
Create a PyKonal EikonalSolver for a 2-D (x,z) model embedded in 3-D (x,y,z) with ny=1.
v_xz must be shape (nx, nz) in km/s.
"""
assert v_xz.shape == (grid.nx, grid.nz)
solver = pykonal.EikonalSolver(coord_sys="cartesian")
# PyKonal uses (x, y, z) order for cartesian.
solver.velocity.min_coords = (grid.x_min, 0.0, grid.z_min)
solver.velocity.node_intervals = (grid.dx, 1.0, grid.dz)
solver.velocity.npts = (grid.nx, 1, grid.nz)
# Embed v(x,z) into v(x,y,z) with ny=1
vv = np.zeros((grid.nx, 1, grid.nz), dtype=float)
vv[:, 0, :] = v_xz
solver.velocity.values = vv
return solver
def initialize_point_source(solver: pykonal.EikonalSolver, grid: Grid2D, xs: float, zs: float) -> Tuple[int, int, int]:
"""
Place a point source at the nearest grid node to (xs,zs).
Returns source index (ix, iy, iz).
"""
ix = int(np.round((xs - grid.x_min) / grid.dx))
iz = int(np.round((zs - grid.z_min) / grid.dz))
ix = np.clip(ix, 0, grid.nx - 1)
iz = np.clip(iz, 0, grid.nz - 1)
src_idx = (ix, 0, iz)
solver.traveltime.values[src_idx] = 0.0
solver.unknown[src_idx] = False
solver.trial.push(*src_idx)
return src_idx
def solve_traveltime(grid: Grid2D, v_xz: np.ndarray, xs: float, zs: float) -> Tuple[pykonal.EikonalSolver, np.ndarray, Tuple[int,int,int]]:
"""
Solve T(x,z) for a given velocity model and point source.
Returns solver, T_xz (nx,nz), and src_idx.
"""
solver = make_solver_from_velocity(grid, v_xz)
src_idx = initialize_point_source(solver, grid, xs, zs)
ok = solver.solve()
if not ok:
raise RuntimeError("PyKonal solver did not return success.")
T = solver.traveltime.values[:, 0, :] # (nx,nz)
return solver, T, src_idx
def trace_rays_to_surface_receivers(
solver: pykonal.EikonalSolver,
grid: Grid2D,
receiver_x: np.ndarray,
receiver_z: float = 0.0,
y_fixed: float = 0.0,
) -> Dict[float, np.ndarray]:
"""
Trace rays from each receiver location (x, z=receiver_z) back to the source using solver.trace_ray(end).
Returns dict mapping receiver x -> ray path array (N,3) with columns (x,y,z).
"""
rays = {}
for xr in receiver_x:
end = np.array([float(xr), float(y_fixed), float(receiver_z)], dtype=float)
ray = solver.trace_ray(end) # coords, not indices
rays[float(xr)] = ray
return rays
def plot_velocity_and_rays(grid: Grid2D, v_xz: np.ndarray, rays: Optional[Dict[float, np.ndarray]] = None,
xs: Optional[float]=None, zs: Optional[float]=None, title: str="") -> None:
X, Z = grid.mesh()
fig, ax = plt.subplots(figsize=(10, 4.3))
pm = ax.pcolormesh(X, Z, v_xz, shading="auto")
ax.invert_yaxis()
cbar = fig.colorbar(pm, ax=ax)
cbar.set_label("Velocity (km/s)")
if rays:
for _, ray in rays.items():
ax.plot(ray[:, 0], ray[:, 2], lw=1.8) # x vs z
if xs is not None and zs is not None:
ax.scatter([xs], [zs], marker="*", s=120, edgecolor="k", zorder=10, label="Source")
ax.legend(loc="upper right")
ax.set_xlabel("x (km)")
ax.set_ylabel("z (km)")
ax.set_title(title)
plt.show()
def plot_traveltime_contours(grid: Grid2D, T_xz: np.ndarray, title: str="") -> None:
X, Z = grid.mesh()
fig, ax = plt.subplots(figsize=(10, 4.3))
cs = ax.contour(X, Z, T_xz, levels=25)
ax.clabel(cs, inline=True, fontsize=8, fmt="%.1f")
ax.invert_yaxis()
ax.set_xlabel("x (km)")
ax.set_ylabel("z (km)")
ax.set_title(title)
plt.show()
def travel_time_curve(T_xz: np.ndarray, grid: Grid2D, receiver_x: np.ndarray, receiver_z: float = 0.0) -> np.ndarray:
"""
Extract T(x, z=receiver_z) at receiver_x (nearest node).
"""
iz = int(np.round((receiver_z - grid.z_min)/grid.dz))
iz = np.clip(iz, 0, grid.nz-1)
ix = np.clip(np.round((receiver_x - grid.x_min)/grid.dx).astype(int), 0, grid.nx-1)
return T_xz[ix, iz]
2.1 Velocity model builders#
We will use simple, controllable models so students can form hypotheses:
Depth-dependent gradient: (v(z) = v_0 + g z)
Layered: piecewise constant or piecewise linear
Curved layers: warp the interface depth as a function of (x)
Low-velocity zone (LVZ): Gaussian slow anomaly centered at ((x_c, z_c))
In all experiments, keep the grid spacing coarse enough to run fast in Colab, but fine enough to see ray curvature.
def v_gradient(grid: Grid2D, v0: float=5.5, grad: float=0.02) -> np.ndarray:
"""
Simple depth gradient: v(x,z) = v0 + grad*z
"""
X, Z = grid.mesh()
return v0 + grad * Z
def v_layered(grid: Grid2D, layers: Tuple[Tuple[float, float], ...]) -> np.ndarray:
"""
Piecewise constant layers.
layers: tuple of (z_top, velocity) sorted by increasing z_top
Example: ((0, 5.5), (15, 6.3), (35, 7.0))
"""
_, Z = grid.mesh()
v = np.zeros((grid.nx, grid.nz), dtype=float)
ztops = [lz[0] for lz in layers]
vels = [lz[1] for lz in layers]
# For each depth node, choose last layer whose z_top <= z
for iz, z in enumerate(grid.z):
idx = max([i for i, zt in enumerate(ztops) if z >= zt])
v[:, iz] = vels[idx]
return v
def warp_layers(v_xz: np.ndarray, grid: Grid2D, amp_km: float=5.0, wavelength_km: float=120.0) -> np.ndarray:
"""
Warp a layered model by shifting depth index as a sinusoidal function of x.
This creates curved/undulating interfaces while preserving layer velocities.
"""
x = grid.x
shift = amp_km * np.sin(2*np.pi*x / wavelength_km) # km
shift_idx = np.round(shift / grid.dz).astype(int)
v_new = np.zeros_like(v_xz)
for ix in range(grid.nx):
v_new[ix, :] = np.roll(v_xz[ix, :], shift_idx[ix])
return v_new
def add_lvz(v_xz: np.ndarray, grid: Grid2D, xc: float, zc: float, dv: float=-0.8, sigx: float=25.0, sigz: float=8.0) -> np.ndarray:
"""
Add a Gaussian low-velocity zone: v -> v + dv * exp(...)
dv should be negative for a LVZ.
"""
X, Z = grid.mesh()
anomaly = np.exp(-0.5*((X-xc)/sigx)**2 - 0.5*((Z-zc)/sigz)**2)
v2 = v_xz + dv*anomaly
# prevent non-physical values
v2 = np.maximum(v2, 0.5)
return v2
3. Experiments#
We will compute and compare travel-time fields, ray paths, and surface travel-time curves.
Experiment A — Flat depth-dependent velocity (baseline)#
Model: smooth increase with depth (v(z) = v_0 + g z)
What to look for
Travel-time contours should be roughly elliptical near the source and become distorted with depth gradient.
Rays should curve and turn more strongly for larger offsets.
Prediction prompt
If you increase the gradient (g), do rays turn shallower or deeper for the same offset?
vA = v_gradient(grid, v0=5.8, grad=0.03)
solverA, TA, _ = solve_traveltime(grid, vA, xs, zs)
raysA = trace_rays_to_surface_receivers(solverA, grid, rx, receiver_z=0.0)
plot_velocity_and_rays(grid, vA, raysA, xs, zs, title="Experiment A: v(z) gradient + traced rays")
plot_traveltime_contours(grid, TA, title="Experiment A: travel-time contours T(x,z)")
# Travel-time curve on surface
tA = travel_time_curve(TA, grid, rx, receiver_z=0.0)
plt.figure(figsize=(7.5,3.2))
plt.plot(rx, tA, marker="o")
plt.xlabel("x (km)")
plt.ylabel("T (s) at z=0")
plt.title("Experiment A: surface travel-time curve T(x, z=0)")
plt.show()
Experiment B — From flat to curved layers#
We now introduce layered structure and then warp it to make curved interfaces.
Model B1: flat layers
Model B2: same layers, but warped sinusoidally in (x)
What to look for
In B1, rays should bend only due to velocity contrasts and gradient-like behavior.
In B2, you should see lateral deflection and mild focusing/defocusing tied to interface curvature.
Prediction prompt
Where do you expect shorter travel times: beneath a syncline-like depression of a fast layer, or beneath an anticline-like uplift?
# Flat layered baseline
layers = (
(0.0, 5.8),
(12.0, 6.4),
(28.0, 7.2),
(45.0, 7.8),
)
vB1 = v_layered(grid, layers)
# Warped layers
vB2 = warp_layers(vB1, grid, amp_km=10.0, wavelength_km=40.0)
# Solve + trace for B1
solverB1, TB1, _ = solve_traveltime(grid, vB1, xs, zs)
raysB1 = trace_rays_to_surface_receivers(solverB1, grid, rx, receiver_z=0.0)
plot_velocity_and_rays(grid, vB1, raysB1, xs, zs, title="Experiment B1: flat layers + rays")
plot_traveltime_contours(grid, TB1, title="Experiment B1: travel-time contours")
tB1 = travel_time_curve(TB1, grid, rx, receiver_z=0.0)
# Solve + trace for B2
solverB2, TB2, _ = solve_traveltime(grid, vB2, xs, zs)
raysB2 = trace_rays_to_surface_receivers(solverB2, grid, rx, receiver_z=0.0)
plot_velocity_and_rays(grid, vB2, raysB2, xs, zs, title="Experiment B2: curved layers + rays")
plot_traveltime_contours(grid, TB2, title="Experiment B2: travel-time contours")
tB2 = travel_time_curve(TB2, grid, rx, receiver_z=0.0)
# Compare travel-time curves
plt.figure(figsize=(7.5,3.2))
plt.plot(rx, tB1, marker="o", label="Flat layers")
plt.plot(rx, tB2, marker="s", label="Curved layers")
plt.plot(rx, 10*(tB2-tB1), marker="+", label="Difference")
plt.xlabel("x (km)")
plt.ylabel("T (s) at z=0")
plt.title("Experiment B: surface travel-time curve comparison")
plt.legend()
plt.show()
Experiment C — Add a Low-Velocity Zone (LVZ)#
We add a mid-crustal LVZ (Gaussian slow anomaly) to the layered model.
What to look for
Rays can bend toward slower material.
Travel times can develop subtle changes in curvature; for strong LVZs you may see travel-time curve kinks or gaps in first arrivals (geometry-driven).
Prediction prompt
If the LVZ is strong enough, do you expect a shadow zone on the surface? Where?
# Start from the layered model and add LVZ
vC = add_lvz(vB1, grid, xc=110.0, zc=22.0, dv=-1, sigx=22.0, sigz=7.0)
solverC, TC, _ = solve_traveltime(grid, vC, xs, zs)
raysC = trace_rays_to_surface_receivers(solverC, grid, rx, receiver_z=0.0)
plot_velocity_and_rays(grid, vC, raysC, xs, zs, title="Experiment C: layered model + LVZ + rays")
plot_traveltime_contours(grid, TC, title="Experiment C: travel-time contours")
tC = travel_time_curve(TC, grid, rx, receiver_z=0.0)
plt.figure(figsize=(7.5,3.2))
plt.plot(rx, tB1, marker="o", label="No LVZ (B1)")
plt.plot(rx, tC, marker="s", label="With LVZ (C)")
plt.xlabel("x (km)")
plt.ylabel("T (s) at z=0")
plt.title("Experiment C: LVZ effect on surface travel-time curve")
plt.legend()
plt.show()
3.5 Residuals: source location vs. structure#
This notebook also illustrates two core research approaches in seismology:
Image the subsurface from travel-time residuals.
Locate sources (earthquakes) by minimizing residuals between predicted and observed arrivals.
What are travel-time residuals?#
We define residuals as:
where \(\mathbf{x}_i\) are receiver locations and \(\mathbf{m}\) denotes the model (velocity structure + source location).
Two causes of residual patterns#
1) Unknown velocity structure (model error)
If the Earth is slower than the model along a path, predicted times are too early, so \(r_i > 0\).
If the Earth is faster than the model, predicted times are too late, so \(r_i < 0\).
Spatially coherent residuals (e.g., a broad positive bulge in \(r(x)\)) often indicate missing structure.
2) Source location error
Moving the source changes travel times systematically with distance and azimuth.
A source that is too shallow or too close to the array often yields a monotonic trend in \(r(x)\).
Location errors tend to produce simple, smooth patterns tied to geometry.
Building intuition from \(r(x)\)#
Broad, localized residual highs/lows → likely structure (unknown velocity anomalies).
Smooth, near-linear residual trend with offset → likely source location (mislocated event).
Mixed patterns suggest both effects; this is why imaging and relocation are coupled problems.
In practice, we iteratively update both velocity and source parameters to reduce residuals across the whole network.
4. Toy problem (experiential learning)#
You introduced a mid-crustal LVZ and noticed the surface travel-time curve changed.
Task#
Quantify the LVZ signature:
Compute (\Delta T(x) = T_{\mathrm{LVZ}}(x) - T_{\mathrm{base}}(x)).
Diagnose whether the change is local or broad:
Where is (\Delta T(x)) maximal?
Hypothesis test (write before running):
If you move the LVZ deeper (increase (z_c)), does (\Delta T) get smaller or larger at the surface? Why?
Stretch goal#
Try two LVZ strengths (dv=-0.5 and dv=-1.5) and see whether the ray geometry changes (e.g., paths bend toward the LVZ).
# 1) Quantify ΔT(x)
dT = tC - tB1
plt.figure(figsize=(7.5,3.2))
plt.plot(rx, dT, marker="o")
plt.axhline(0, lw=1)
plt.xlabel("x (km)")
plt.ylabel("ΔT (s) = T_LVZ - T_base")
plt.title("Toy problem: LVZ travel-time perturbation on the surface")
plt.show()
# 2) Where is ΔT maximal?
imax = np.argmax(dT)
print(f"Max ΔT ≈ {dT[imax]:.3f} s at x ≈ {rx[imax]:.1f} km")
# 3) Hypothesis test: move LVZ deeper and recompute quickly
vC_deeper = add_lvz(vB1, grid, xc=110.0, zc=30.0, dv=-1.0, sigx=22.0, sigz=7.0)
solverCd, TCd, _ = solve_traveltime(grid, vC_deeper, xs, zs)
tCd = travel_time_curve(TCd, grid, rx, receiver_z=0.0)
dT_deeper = tCd - tB1
plt.figure(figsize=(7.5,3.2))
plt.plot(rx, dT, marker="o", label="LVZ at zc=22 km")
plt.plot(rx, dT_deeper, marker="s", label="LVZ at zc=30 km")
plt.axhline(0, lw=1)
plt.xlabel("x (km)")
plt.ylabel("ΔT (s)")
plt.title("Toy problem: effect of LVZ depth on surface ΔT(x)")
plt.legend()
plt.show()