Travel Times Through the Earth¶
Run in Google Colab
Seismic waves traveling through the Earth encounter layers with different velocities. At each interface the ray bends — governed by Snell's law and the conserved ray parameter. In this lecture, we derive the travel time formulas, build travel time curves, explore how velocity structure creates distinctive signatures, and connect everything to the real Earth.
In [1]:
Copied!
!pip install obspy -q
!pip install obspy -q
In [2]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import sympy as sp
from obspy.taup import TauPyModel
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from obspy.geodetics import gps2dist_azimuth
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import sympy as sp
from obspy.taup import TauPyModel
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from obspy.geodetics import gps2dist_azimuth
In [3]:
Copied!
def trace_ray(velocities, thicknesses, p):
"""Trace a ray through horizontal layers. Returns (x, z) of the downgoing path."""
slownesses = 1.0 / np.array(velocities)
x, z = [0.0], [0.0]
for u, dz in zip(slownesses, thicknesses):
if p >= u:
break
dx = p * dz / np.sqrt(u**2 - p**2)
x.append(x[-1] + dx)
z.append(z[-1] + dz)
return np.array(x), np.array(z)
def draw_layers(ax, velocities, thicknesses, xlim=12):
"""Draw colored horizontal layers with velocity labels."""
colors = plt.cm.YlOrBr(np.linspace(0.15, 0.55, len(velocities)))
z_top = 0
for i, (v, dz, c) in enumerate(zip(velocities, thicknesses, colors)):
ax.add_patch(patches.Rectangle((0, z_top), xlim, dz, fc=c, ec='k', lw=1))
ax.text(xlim - 0.3, z_top + dz / 2, f'$v_{i+1}$ = {v} km/s',
va='center', ha='right', fontsize=10,
bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.8))
z_top += dz
ax.set_xlim(0, xlim)
ax.set_ylim(z_top + 0.3, -0.5)
ax.set_xlabel('Horizontal distance (km)')
ax.set_ylabel('Depth (km)')
def compute_XT(velocities, thicknesses, p):
"""Compute total distance X and travel time T for ray parameter p (two-way path)."""
slownesses = 1.0 / np.array(velocities)
X, T = 0.0, 0.0
for u, dz in zip(slownesses, thicknesses):
if p >= u:
break
eta = np.sqrt(u**2 - p**2)
X += 2 * p * dz / eta
T += 2 * u**2 * dz / eta
return X, T
def trace_ray(velocities, thicknesses, p):
"""Trace a ray through horizontal layers. Returns (x, z) of the downgoing path."""
slownesses = 1.0 / np.array(velocities)
x, z = [0.0], [0.0]
for u, dz in zip(slownesses, thicknesses):
if p >= u:
break
dx = p * dz / np.sqrt(u**2 - p**2)
x.append(x[-1] + dx)
z.append(z[-1] + dz)
return np.array(x), np.array(z)
def draw_layers(ax, velocities, thicknesses, xlim=12):
"""Draw colored horizontal layers with velocity labels."""
colors = plt.cm.YlOrBr(np.linspace(0.15, 0.55, len(velocities)))
z_top = 0
for i, (v, dz, c) in enumerate(zip(velocities, thicknesses, colors)):
ax.add_patch(patches.Rectangle((0, z_top), xlim, dz, fc=c, ec='k', lw=1))
ax.text(xlim - 0.3, z_top + dz / 2, f'$v_{i+1}$ = {v} km/s',
va='center', ha='right', fontsize=10,
bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.8))
z_top += dz
ax.set_xlim(0, xlim)
ax.set_ylim(z_top + 0.3, -0.5)
ax.set_xlabel('Horizontal distance (km)')
ax.set_ylabel('Depth (km)')
def compute_XT(velocities, thicknesses, p):
"""Compute total distance X and travel time T for ray parameter p (two-way path)."""
slownesses = 1.0 / np.array(velocities)
X, T = 0.0, 0.0
for u, dz in zip(slownesses, thicknesses):
if p >= u:
break
eta = np.sqrt(u**2 - p**2)
X += 2 * p * dz / eta
T += 2 * u**2 * dz / eta
return X, T
1. Snell's Law and the Ray Parameter¶
At each interface between layers:
$$\frac{\sin\theta_1}{v_1} = \frac{\sin\theta_2}{v_2} = p$$
where $p$ is the ray parameter (units: s/km), constant along the entire ray path.
Defining the slowness $u = 1/v$, we can write $p = u \sin\theta$. The slowness $u$ is the maximum possible ray parameter for a given layer (when $\theta = 90°$).
| $p$ value | Ray behavior |
|---|---|
| Small $p$ | Steep ray (nearly vertical) |
| Large $p$ | Shallow ray (nearly horizontal) |
| $p = u_i$ | Ray goes horizontal in layer $i$ (critical angle) |
| $p > u_i$ | Ray cannot enter layer $i$ (total reflection) |
Critical angle and total reflection¶
When velocity increases across an interface, the ray bends away from the vertical. At the critical angle $\theta_c = \arcsin(v_1/v_2)$, the transmitted ray goes horizontal. Beyond this, the ray is totally reflected.
In [4]:
Copied!
v1, v2 = 4, 6 # km/s
theta_crit = np.degrees(np.arcsin(v1 / v2))
p_crit = 1 / v2
print(f"v₁ = {v1} km/s, v₂ = {v2} km/s")
print(f"Critical angle: θ_c = arcsin(v₁/v₂) = {theta_crit:.1f}°")
print(f"Critical ray parameter: p_c = 1/v₂ = {p_crit:.4f} s/km")
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
cases = [
(0.12, f'Below critical (θ < {theta_crit:.0f}°)\nRay passes through'),
(p_crit, f'At critical angle (θ = {theta_crit:.0f}°)\nRay goes horizontal'),
(0.20, f'Above critical (θ > {theta_crit:.0f}°)\nTotal reflection'),
]
for ax, (p, title) in zip(axes, cases):
for j, (dz, c) in enumerate(zip([3, 3], ['#f0d9b5', '#c4a35a'])):
ax.add_patch(patches.Rectangle((0, j*3), 10, dz, fc=c, ec='k', lw=1))
ax.text(9, 1.5, f'v₁ = {v1}', va='center', ha='right', fontsize=10,
bbox=dict(boxstyle='round', fc='white', alpha=0.8))
ax.text(9, 4.5, f'v₂ = {v2}', va='center', ha='right', fontsize=10,
bbox=dict(boxstyle='round', fc='white', alpha=0.8))
u1 = 1/v1
theta1 = np.arcsin(min(p * v1, 1.0))
dx1 = 3 * np.tan(theta1)
ax.plot([1, 1 + dx1], [0, 3], 'r-', lw=2.5)
ax.plot(1, 0, 'r*', ms=12)
u2 = 1/v2
if p < u2: # passes through
theta2 = np.arcsin(p * v2)
dx2 = 3 * np.tan(theta2)
ax.plot([1+dx1, 1+dx1+dx2], [3, 6], 'r-', lw=2.5)
elif abs(p - u2) < 1e-10: # critical — head wave
ax.plot([1+dx1, 9], [3, 3], 'r-', lw=2.5)
ax.text(5, 2.5, 'head wave', fontsize=10, color='red', ha='center')
else: # total reflection
ax.plot([1+dx1, 1+2*dx1], [3, 0], 'r--', lw=2.5, alpha=0.7)
ax.set_xlim(0, 10)
ax.set_ylim(6.5, -0.5)
ax.set_title(f'{title}\np = {p:.4f} s/km', fontsize=10)
ax.set_xlabel('Distance (km)')
axes[0].set_ylabel('Depth (km)')
plt.tight_layout()
plt.show()
v1, v2 = 4, 6 # km/s
theta_crit = np.degrees(np.arcsin(v1 / v2))
p_crit = 1 / v2
print(f"v₁ = {v1} km/s, v₂ = {v2} km/s")
print(f"Critical angle: θ_c = arcsin(v₁/v₂) = {theta_crit:.1f}°")
print(f"Critical ray parameter: p_c = 1/v₂ = {p_crit:.4f} s/km")
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
cases = [
(0.12, f'Below critical (θ < {theta_crit:.0f}°)\nRay passes through'),
(p_crit, f'At critical angle (θ = {theta_crit:.0f}°)\nRay goes horizontal'),
(0.20, f'Above critical (θ > {theta_crit:.0f}°)\nTotal reflection'),
]
for ax, (p, title) in zip(axes, cases):
for j, (dz, c) in enumerate(zip([3, 3], ['#f0d9b5', '#c4a35a'])):
ax.add_patch(patches.Rectangle((0, j*3), 10, dz, fc=c, ec='k', lw=1))
ax.text(9, 1.5, f'v₁ = {v1}', va='center', ha='right', fontsize=10,
bbox=dict(boxstyle='round', fc='white', alpha=0.8))
ax.text(9, 4.5, f'v₂ = {v2}', va='center', ha='right', fontsize=10,
bbox=dict(boxstyle='round', fc='white', alpha=0.8))
u1 = 1/v1
theta1 = np.arcsin(min(p * v1, 1.0))
dx1 = 3 * np.tan(theta1)
ax.plot([1, 1 + dx1], [0, 3], 'r-', lw=2.5)
ax.plot(1, 0, 'r*', ms=12)
u2 = 1/v2
if p < u2: # passes through
theta2 = np.arcsin(p * v2)
dx2 = 3 * np.tan(theta2)
ax.plot([1+dx1, 1+dx1+dx2], [3, 6], 'r-', lw=2.5)
elif abs(p - u2) < 1e-10: # critical — head wave
ax.plot([1+dx1, 9], [3, 3], 'r-', lw=2.5)
ax.text(5, 2.5, 'head wave', fontsize=10, color='red', ha='center')
else: # total reflection
ax.plot([1+dx1, 1+2*dx1], [3, 0], 'r--', lw=2.5, alpha=0.7)
ax.set_xlim(0, 10)
ax.set_ylim(6.5, -0.5)
ax.set_title(f'{title}\np = {p:.4f} s/km', fontsize=10)
ax.set_xlabel('Distance (km)')
axes[0].set_ylabel('Depth (km)')
plt.tight_layout()
plt.show()
v₁ = 4 km/s, v₂ = 6 km/s Critical angle: θ_c = arcsin(v₁/v₂) = 41.8° Critical ray parameter: p_c = 1/v₂ = 0.1667 s/km
2. Travel Time and Distance Formulas¶
For a ray with parameter $p$ crossing a single layer of thickness $\Delta z$ and slowness $u = 1/v$:
$$\Delta x = \frac{p \, \Delta z}{\sqrt{u^2 - p^2}}, \qquad \Delta t = \frac{u^2 \, \Delta z}{\sqrt{u^2 - p^2}}$$
For a ray going down and back up through $N$ layers, the total distance and travel time are:
$$X(p) = 2 \sum_{i=1}^{N} \frac{p \, \Delta z_i}{\sqrt{u_i^2 - p^2}}, \qquad T(p) = 2 \sum_{i=1}^{N} \frac{u_i^2 \, \Delta z_i}{\sqrt{u_i^2 - p^2}}$$
The factor of 2 accounts for the symmetric downgoing + upgoing path.
Worked example: 3-layer model¶
Model: $v_1 = 4$, $v_2 = 6$, $v_3 = 8$ km/s, all layers 3 km thick. Ray parameter: $p = 0.15$ s/km.
In [5]:
Copied!
velocities = np.array([4.0, 6.0, 8.0]) # km/s
thicknesses = np.array([3.0, 3.0, 3.0]) # km
slownesses = 1.0 / velocities
p = 0.15 # s/km
print("Step 1: Convert to slownesses")
for i, (v, u) in enumerate(zip(velocities, slownesses)):
print(f" Layer {i+1}: v = {v} km/s → u = {u:.4f} s/km")
print(f"\nStep 2: Check penetration (p = {p} s/km)")
X_total, T_total = 0, 0
for i, (u, dz) in enumerate(zip(slownesses, thicknesses)):
if p >= u:
print(f" Layer {i+1}: p = {p} ≥ u = {u:.4f} → ray REFLECTS, does not enter")
else:
eta = np.sqrt(u**2 - p**2)
dx = p * dz / eta
dt = u**2 * dz / eta
X_total += 2 * dx
T_total += 2 * dt
print(f" Layer {i+1}: p = {p} < u = {u:.4f} → passes through, Δx = {dx:.3f} km, Δt = {dt:.3f} s")
print(f"\nStep 3: Total (two-way path)")
print(f" X(p) = {X_total:.2f} km")
print(f" T(p) = {T_total:.2f} s")
velocities = np.array([4.0, 6.0, 8.0]) # km/s
thicknesses = np.array([3.0, 3.0, 3.0]) # km
slownesses = 1.0 / velocities
p = 0.15 # s/km
print("Step 1: Convert to slownesses")
for i, (v, u) in enumerate(zip(velocities, slownesses)):
print(f" Layer {i+1}: v = {v} km/s → u = {u:.4f} s/km")
print(f"\nStep 2: Check penetration (p = {p} s/km)")
X_total, T_total = 0, 0
for i, (u, dz) in enumerate(zip(slownesses, thicknesses)):
if p >= u:
print(f" Layer {i+1}: p = {p} ≥ u = {u:.4f} → ray REFLECTS, does not enter")
else:
eta = np.sqrt(u**2 - p**2)
dx = p * dz / eta
dt = u**2 * dz / eta
X_total += 2 * dx
T_total += 2 * dt
print(f" Layer {i+1}: p = {p} < u = {u:.4f} → passes through, Δx = {dx:.3f} km, Δt = {dt:.3f} s")
print(f"\nStep 3: Total (two-way path)")
print(f" X(p) = {X_total:.2f} km")
print(f" T(p) = {T_total:.2f} s")
Step 1: Convert to slownesses Layer 1: v = 4.0 km/s → u = 0.2500 s/km Layer 2: v = 6.0 km/s → u = 0.1667 s/km Layer 3: v = 8.0 km/s → u = 0.1250 s/km Step 2: Check penetration (p = 0.15 s/km) Layer 1: p = 0.15 < u = 0.2500 → passes through, Δx = 2.250 km, Δt = 0.938 s Layer 2: p = 0.15 < u = 0.1667 → passes through, Δx = 6.194 km, Δt = 1.147 s Layer 3: p = 0.15 ≥ u = 0.1250 → ray REFLECTS, does not enter Step 3: Total (two-way path) X(p) = 16.89 km T(p) = 4.17 s
Visualizing the ray path¶
In [6]:
Copied!
p = 0.15
x_down, z_down = trace_ray(velocities, thicknesses, p)
X_total = 2 * x_down[-1]
x_full = np.concatenate([x_down, X_total - x_down[::-1]])
z_full = np.concatenate([z_down, z_down[::-1]])
fig, ax = plt.subplots(figsize=(10, 6))
draw_layers(ax, velocities, thicknesses, xlim=20)
ax.plot(x_full, z_full, 'r-', lw=2.5)
ax.plot(x_full[0], 0, 'r*', ms=14, label='Source')
ax.plot(x_full[-1], 0, 'bv', ms=10, label=f'Receiver (X = {X_total:.1f} km)')
X, T = compute_XT(velocities, thicknesses, p)
ax.set_title(f'Two-way ray path: p = {p} s/km → X = {X:.1f} km, T = {T:.2f} s', fontsize=12)
ax.legend(loc='lower right')
plt.tight_layout()
plt.show()
p = 0.15
x_down, z_down = trace_ray(velocities, thicknesses, p)
X_total = 2 * x_down[-1]
x_full = np.concatenate([x_down, X_total - x_down[::-1]])
z_full = np.concatenate([z_down, z_down[::-1]])
fig, ax = plt.subplots(figsize=(10, 6))
draw_layers(ax, velocities, thicknesses, xlim=20)
ax.plot(x_full, z_full, 'r-', lw=2.5)
ax.plot(x_full[0], 0, 'r*', ms=14, label='Source')
ax.plot(x_full[-1], 0, 'bv', ms=10, label=f'Receiver (X = {X_total:.1f} km)')
X, T = compute_XT(velocities, thicknesses, p)
ax.set_title(f'Two-way ray path: p = {p} s/km → X = {X:.1f} km, T = {T:.2f} s', fontsize=12)
ax.legend(loc='lower right')
plt.tight_layout()
plt.show()
Try it yourself¶
Change the ray parameter p in the cell below and re-run. How do the distance and travel time change?
In [7]:
Copied!
p_values = [0.05, 0.10, 0.15, 0.20]
fig, axes = plt.subplots(1, len(p_values), figsize=(16, 6), sharey=True)
for ax, p in zip(axes, p_values):
draw_layers(ax, velocities, thicknesses, xlim=20)
x_down, z_down = trace_ray(velocities, thicknesses, p)
X_total = 2 * x_down[-1]
x_full = np.concatenate([x_down, X_total - x_down[::-1]])
z_full = np.concatenate([z_down, z_down[::-1]])
ax.plot(x_full, z_full, 'r-', lw=2.5)
ax.plot(0, 0, 'r*', ms=12)
ax.plot(X_total, 0, 'bv', ms=8)
X, T = compute_XT(velocities, thicknesses, p)
ax.set_title(f'p = {p} s/km\nX = {X:.1f} km, T = {T:.2f} s', fontsize=10)
if ax != axes[0]:
ax.set_ylabel('')
fig.suptitle('How the ray parameter p controls ray geometry', fontsize=13, y=1.02)
plt.tight_layout()
plt.show()
p_values = [0.05, 0.10, 0.15, 0.20]
fig, axes = plt.subplots(1, len(p_values), figsize=(16, 6), sharey=True)
for ax, p in zip(axes, p_values):
draw_layers(ax, velocities, thicknesses, xlim=20)
x_down, z_down = trace_ray(velocities, thicknesses, p)
X_total = 2 * x_down[-1]
x_full = np.concatenate([x_down, X_total - x_down[::-1]])
z_full = np.concatenate([z_down, z_down[::-1]])
ax.plot(x_full, z_full, 'r-', lw=2.5)
ax.plot(0, 0, 'r*', ms=12)
ax.plot(X_total, 0, 'bv', ms=8)
X, T = compute_XT(velocities, thicknesses, p)
ax.set_title(f'p = {p} s/km\nX = {X:.1f} km, T = {T:.2f} s', fontsize=10)
if ax != axes[0]:
ax.set_ylabel('')
fig.suptitle('How the ray parameter p controls ray geometry', fontsize=13, y=1.02)
plt.tight_layout()
plt.show()
3. Building Travel Time Curves¶
Each ray parameter $p$ gives one point $(X, T)$ on the travel time curve. By sweeping through many values of $p$, we build the complete $T(X)$ curve — the fundamental observable in seismology.
In [8]:
Copied!
p_list = [0.05, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
draw_layers(ax1, velocities, thicknesses, xlim=25)
cmap = plt.cm.viridis(np.linspace(0, 0.9, len(p_list)))
for p, color in zip(p_list, cmap):
X, T = compute_XT(velocities, thicknesses, p)
x_d, z_d = trace_ray(velocities, thicknesses, p)
x_full = np.concatenate([x_d, X - x_d[::-1]])
z_full = np.concatenate([z_d, z_d[::-1]])
ax1.plot(x_full, z_full, '-', color=color, lw=1.5)
ax2.plot(X, T, 'o', color=color, ms=8)
ax1.set_title('Ray paths')
ax2.set_xlabel('Distance X (km)')
ax2.set_ylabel('Travel time T (s)')
ax2.set_title('Travel time curve (one point per ray parameter)')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
p_list = [0.05, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
draw_layers(ax1, velocities, thicknesses, xlim=25)
cmap = plt.cm.viridis(np.linspace(0, 0.9, len(p_list)))
for p, color in zip(p_list, cmap):
X, T = compute_XT(velocities, thicknesses, p)
x_d, z_d = trace_ray(velocities, thicknesses, p)
x_full = np.concatenate([x_d, X - x_d[::-1]])
z_full = np.concatenate([z_d, z_d[::-1]])
ax1.plot(x_full, z_full, '-', color=color, lw=1.5)
ax2.plot(X, T, 'o', color=color, ms=8)
ax1.set_title('Ray paths')
ax2.set_xlabel('Distance X (km)')
ax2.set_ylabel('Travel time T (s)')
ax2.set_title('Travel time curve (one point per ray parameter)')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
4. Continuous Velocity Models¶
Real Earth velocity varies continuously with depth, not in discrete jumps. As layer thicknesses become infinitely thin, the summation formulas become integrals:
$$X(p) = 2 \int_0^{z_p} \frac{p}{\sqrt{u(z)^2 - p^2}} \, dz, \qquad T(p) = 2 \int_0^{z_p} \frac{u(z)^2}{\sqrt{u(z)^2 - p^2}} \, dz$$
where $z_p$ is the turning depth — the depth where $u(z_p) = p$ and the ray goes horizontal.
Linear velocity gradient: $v(z) = v_0 + kz$¶
This is the simplest continuous model: velocity increases linearly with depth. A remarkable property is that ray paths are circular arcs in a linear gradient.
In [9]:
Copied!
# Linear gradient parameters
v0 = 4.0 # surface velocity (km/s)
k = 0.1 # gradient (1/s)
z_max = 50 # maximum depth (km)
# Sweep ray parameters (needed for evaluating the analytical expressions)
p_min = 1.0 / (v0 + k * z_max) # slowness at bottom
p_max = 1.0 / v0 - 1e-6 # just below surface slowness
p_arr = np.linspace(p_min + 1e-4, p_max, 300)
# Linear gradient parameters
v0 = 4.0 # surface velocity (km/s)
k = 0.1 # gradient (1/s)
z_max = 50 # maximum depth (km)
# Sweep ray parameters (needed for evaluating the analytical expressions)
p_min = 1.0 / (v0 + k * z_max) # slowness at bottom
p_max = 1.0 / v0 - 1e-6 # just below surface slowness
p_arr = np.linspace(p_min + 1e-4, p_max, 300)
Analytical solution using SymPy¶
For $v(z) = v_0 + kz$, the integrals have closed-form solutions. Let's derive them symbolically.
In [10]:
Copied!
# Symbolic derivation
z, p_sym, v0_sym, k_sym = sp.symbols('z p v_0 k', positive=True)
v_z = v0_sym + k_sym * z
u_z = 1 / v_z
# Turning depth: u(z_p) = p → z_p = (1/p - v0) / k
z_turn = sp.solve(u_z - p_sym, z)[0]
print("Turning depth: z_p =", z_turn)
# X(p) integrand
integrand_X = p_sym / sp.sqrt(u_z**2 - p_sym**2)
X_half = sp.integrate(integrand_X, (z, 0, z_turn))
X_analytic = sp.simplify(2 * X_half)
print("\nX(p) =", X_analytic)
# T(p) integrand
integrand_T = u_z**2 / sp.sqrt(u_z**2 - p_sym**2)
T_half = sp.integrate(integrand_T, (z, 0, z_turn))
T_analytic = sp.simplify(2 * T_half)
print("\nT(p) =", T_analytic)
# Symbolic derivation
z, p_sym, v0_sym, k_sym = sp.symbols('z p v_0 k', positive=True)
v_z = v0_sym + k_sym * z
u_z = 1 / v_z
# Turning depth: u(z_p) = p → z_p = (1/p - v0) / k
z_turn = sp.solve(u_z - p_sym, z)[0]
print("Turning depth: z_p =", z_turn)
# X(p) integrand
integrand_X = p_sym / sp.sqrt(u_z**2 - p_sym**2)
X_half = sp.integrate(integrand_X, (z, 0, z_turn))
X_analytic = sp.simplify(2 * X_half)
print("\nX(p) =", X_analytic)
# T(p) integrand
integrand_T = u_z**2 / sp.sqrt(u_z**2 - p_sym**2)
T_half = sp.integrate(integrand_T, (z, 0, z_turn))
T_analytic = sp.simplify(2 * T_half)
print("\nT(p) =", T_analytic)
Turning depth: z_p = (-p*v_0 + 1)/(k*p)
X(p) = 2*sqrt(-p**2*v_0**2 + 1)/(k*p)
T(p) = 2*Integral(1/((k*z + v_0)*sqrt(-k*p*z - p*v_0 + 1)*sqrt(k*p*z + p*v_0 + 1)), (z, 0, (-p*v_0 + 1)/(k*p)))
In [11]:
Copied!
# Evaluate the analytical expressions numerically
X_func = sp.lambdify((p_sym, v0_sym, k_sym), X_analytic)
T_func = sp.lambdify((p_sym, v0_sym, k_sym), T_analytic)
X_ana = np.array([X_func(p, v0, k) for p in p_arr])
T_ana = np.array([T_func(p, v0, k) for p in p_arr])
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Velocity model
z_plot = np.linspace(0, z_max, 200)
axes[0].plot(v0 + k * z_plot, z_plot, 'b-', lw=2)
axes[0].set_xlabel('Velocity (km/s)')
axes[0].set_ylabel('Depth (km)')
axes[0].set_title(f'v(z) = {v0} + {k}z km/s')
axes[0].invert_yaxis()
axes[0].grid(True, alpha=0.3)
# T(X) curve
axes[1].plot(X_ana, T_ana, 'b-', lw=2)
axes[1].set_xlabel('Distance X (km)')
axes[1].set_ylabel('Travel time T (s)')
axes[1].set_title('T(X) — analytical solution')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Evaluate the analytical expressions numerically
X_func = sp.lambdify((p_sym, v0_sym, k_sym), X_analytic)
T_func = sp.lambdify((p_sym, v0_sym, k_sym), T_analytic)
X_ana = np.array([X_func(p, v0, k) for p in p_arr])
T_ana = np.array([T_func(p, v0, k) for p in p_arr])
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Velocity model
z_plot = np.linspace(0, z_max, 200)
axes[0].plot(v0 + k * z_plot, z_plot, 'b-', lw=2)
axes[0].set_xlabel('Velocity (km/s)')
axes[0].set_ylabel('Depth (km)')
axes[0].set_title(f'v(z) = {v0} + {k}z km/s')
axes[0].invert_yaxis()
axes[0].grid(True, alpha=0.3)
# T(X) curve
axes[1].plot(X_ana, T_ana, 'b-', lw=2)
axes[1].set_xlabel('Distance X (km)')
axes[1].set_ylabel('Travel time T (s)')
axes[1].set_title('T(X) — analytical solution')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Circular ray paths¶
For a linear velocity gradient, ray paths are circular arcs. The center of curvature is at depth $z_c = -v_0/k$ (above the surface) and the radius is $R = 1/(pk)$.
In [12]:
Copied!
# Discretize the linear gradient into thin layers for ray tracing
n_layers = 500
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u_layers = 1.0 / (v0 + k * z_mid)
fig, ax = plt.subplots(figsize=(12, 6))
# Background: shade velocity
z_bg = np.linspace(0, z_max, 100)
for i in range(len(z_bg) - 1):
ax.axhspan(z_bg[i], z_bg[i+1], color=plt.cm.YlOrBr(0.15 + 0.4 * z_bg[i] / z_max), alpha=0.5)
# Trace rays through thin layers
p_rays = np.linspace(p_min + 0.005, p_max - 0.005, 8)
cmap = plt.cm.viridis(np.linspace(0, 0.9, len(p_rays)))
for p, color in zip(p_rays, cmap):
x, z = [0.0], [0.0]
for u in u_layers:
if p >= u:
break
dx = p * dz / np.sqrt(u**2 - p**2)
x.append(x[-1] + dx)
z.append(z[-1] + dz)
x, z = np.array(x), np.array(z)
X_total = 2 * x[-1]
x_full = np.concatenate([x, X_total - x[::-1]])
z_full = np.concatenate([z, z[::-1]])
ax.plot(x_full, z_full, '-', color=color, lw=2)
ax.set_xlim(0, 250)
ax.set_ylim(z_max, -1)
ax.set_xlabel('Distance (km)', fontsize=12)
ax.set_ylabel('Depth (km)', fontsize=12)
ax.set_title(f'Ray paths in a linear gradient: v(z) = {v0} + {k}z km/s\n(circular arcs)', fontsize=13)
plt.tight_layout()
plt.show()
# Discretize the linear gradient into thin layers for ray tracing
n_layers = 500
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u_layers = 1.0 / (v0 + k * z_mid)
fig, ax = plt.subplots(figsize=(12, 6))
# Background: shade velocity
z_bg = np.linspace(0, z_max, 100)
for i in range(len(z_bg) - 1):
ax.axhspan(z_bg[i], z_bg[i+1], color=plt.cm.YlOrBr(0.15 + 0.4 * z_bg[i] / z_max), alpha=0.5)
# Trace rays through thin layers
p_rays = np.linspace(p_min + 0.005, p_max - 0.005, 8)
cmap = plt.cm.viridis(np.linspace(0, 0.9, len(p_rays)))
for p, color in zip(p_rays, cmap):
x, z = [0.0], [0.0]
for u in u_layers:
if p >= u:
break
dx = p * dz / np.sqrt(u**2 - p**2)
x.append(x[-1] + dx)
z.append(z[-1] + dz)
x, z = np.array(x), np.array(z)
X_total = 2 * x[-1]
x_full = np.concatenate([x, X_total - x[::-1]])
z_full = np.concatenate([z, z[::-1]])
ax.plot(x_full, z_full, '-', color=color, lw=2)
ax.set_xlim(0, 250)
ax.set_ylim(z_max, -1)
ax.set_xlabel('Distance (km)', fontsize=12)
ax.set_ylabel('Depth (km)', fontsize=12)
ax.set_title(f'Ray paths in a linear gradient: v(z) = {v0} + {k}z km/s\n(circular arcs)', fontsize=13)
plt.tight_layout()
plt.show()
5. Interpreting Velocity Structure from Travel Time Curves¶
The shape of the $T(X)$ curve encodes information about the velocity structure. Three types of velocity profiles produce distinctive signatures — and seismologists used these signatures to discover major Earth structures.
In [13]:
Copied!
def compute_XT_continuous(v_func, z_max, p, n_layers=1000):
"""Compute X and T for a continuous velocity model v(z) by discretization."""
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u = 1.0 / v_func(z_mid)
X, T = 0.0, 0.0
for ui, _ in zip(u, z_mid):
if p >= ui:
break
eta = np.sqrt(ui**2 - p**2)
X += 2 * p * dz / eta
T += 2 * ui**2 * dz / eta
return X, T
def trace_ray_continuous(v_func, z_max, p, n_layers=1000):
"""Trace a ray through a continuous model. Returns full (down+up) path."""
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u = 1.0 / v_func(z_mid)
x, z = [0.0], [0.0]
for ui, ddz in zip(u, np.full(n_layers, dz)):
if p >= ui:
break
dx = p * ddz / np.sqrt(ui**2 - p**2)
x.append(x[-1] + dx)
z.append(z[-1] + ddz)
x, z = np.array(x), np.array(z)
X_total = 2 * x[-1]
return np.concatenate([x, X_total - x[::-1]]), np.concatenate([z, z[::-1]])
def sweep_TX(v_func, z_max, n_p=400, n_layers=1000):
"""Sweep ray parameters and return X, T arrays."""
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u = 1.0 / v_func(z_mid)
p_min, p_max = u.min() + 1e-6, u.max() - 1e-6
p_arr = np.linspace(p_min, p_max, n_p)
results = np.array([compute_XT_continuous(v_func, z_max, p, n_layers) for p in p_arr])
return results[:, 0], results[:, 1], p_arr
def plot_model_rays_TX(v_func, z_max, title, n_rays=8):
"""Three-panel plot: velocity model, ray paths, T(X) curve."""
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Velocity model
z_plot = np.linspace(0, z_max, 500)
axes[0].plot(v_func(z_plot), z_plot, 'b-', lw=2)
axes[0].set_xlabel('Velocity (km/s)')
axes[0].set_ylabel('Depth (km)')
axes[0].set_title('v(z)')
axes[0].invert_yaxis()
axes[0].grid(True, alpha=0.3)
# Ray paths
dz = z_max / 1000
z_mid = np.arange(1000) * dz + dz / 2
u = 1.0 / v_func(z_mid)
p_lo, p_hi = u.min() + 1e-4, u.max() - 1e-4
p_rays = np.linspace(p_lo, p_hi, n_rays)
cmap = plt.cm.viridis(np.linspace(0, 0.9, n_rays))
for p, color in zip(p_rays, cmap):
rx, rz = trace_ray_continuous(v_func, z_max, p)
axes[1].plot(rx, rz, '-', color=color, lw=1.5)
axes[1].set_xlabel('Distance (km)')
axes[1].set_ylabel('Depth (km)')
axes[1].set_title('Ray paths')
axes[1].invert_yaxis()
# T(X) curve — use scatter so gaps (shadow zones) are visible
X, T, _ = sweep_TX(v_func, z_max)
axes[2].scatter(X, T, c='b', s=4)
axes[2].set_xlabel('Distance X (km)')
axes[2].set_ylabel('Travel time T (s)')
axes[2].set_title('T(X)')
axes[2].grid(True, alpha=0.3)
fig.suptitle(title, fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
return X, T
def compute_XT_continuous(v_func, z_max, p, n_layers=1000):
"""Compute X and T for a continuous velocity model v(z) by discretization."""
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u = 1.0 / v_func(z_mid)
X, T = 0.0, 0.0
for ui, _ in zip(u, z_mid):
if p >= ui:
break
eta = np.sqrt(ui**2 - p**2)
X += 2 * p * dz / eta
T += 2 * ui**2 * dz / eta
return X, T
def trace_ray_continuous(v_func, z_max, p, n_layers=1000):
"""Trace a ray through a continuous model. Returns full (down+up) path."""
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u = 1.0 / v_func(z_mid)
x, z = [0.0], [0.0]
for ui, ddz in zip(u, np.full(n_layers, dz)):
if p >= ui:
break
dx = p * ddz / np.sqrt(ui**2 - p**2)
x.append(x[-1] + dx)
z.append(z[-1] + ddz)
x, z = np.array(x), np.array(z)
X_total = 2 * x[-1]
return np.concatenate([x, X_total - x[::-1]]), np.concatenate([z, z[::-1]])
def sweep_TX(v_func, z_max, n_p=400, n_layers=1000):
"""Sweep ray parameters and return X, T arrays."""
dz = z_max / n_layers
z_mid = np.arange(n_layers) * dz + dz / 2
u = 1.0 / v_func(z_mid)
p_min, p_max = u.min() + 1e-6, u.max() - 1e-6
p_arr = np.linspace(p_min, p_max, n_p)
results = np.array([compute_XT_continuous(v_func, z_max, p, n_layers) for p in p_arr])
return results[:, 0], results[:, 1], p_arr
def plot_model_rays_TX(v_func, z_max, title, n_rays=8):
"""Three-panel plot: velocity model, ray paths, T(X) curve."""
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# Velocity model
z_plot = np.linspace(0, z_max, 500)
axes[0].plot(v_func(z_plot), z_plot, 'b-', lw=2)
axes[0].set_xlabel('Velocity (km/s)')
axes[0].set_ylabel('Depth (km)')
axes[0].set_title('v(z)')
axes[0].invert_yaxis()
axes[0].grid(True, alpha=0.3)
# Ray paths
dz = z_max / 1000
z_mid = np.arange(1000) * dz + dz / 2
u = 1.0 / v_func(z_mid)
p_lo, p_hi = u.min() + 1e-4, u.max() - 1e-4
p_rays = np.linspace(p_lo, p_hi, n_rays)
cmap = plt.cm.viridis(np.linspace(0, 0.9, n_rays))
for p, color in zip(p_rays, cmap):
rx, rz = trace_ray_continuous(v_func, z_max, p)
axes[1].plot(rx, rz, '-', color=color, lw=1.5)
axes[1].set_xlabel('Distance (km)')
axes[1].set_ylabel('Depth (km)')
axes[1].set_title('Ray paths')
axes[1].invert_yaxis()
# T(X) curve — use scatter so gaps (shadow zones) are visible
X, T, _ = sweep_TX(v_func, z_max)
axes[2].scatter(X, T, c='b', s=4)
axes[2].set_xlabel('Distance X (km)')
axes[2].set_ylabel('Travel time T (s)')
axes[2].set_title('T(X)')
axes[2].grid(True, alpha=0.3)
fig.suptitle(title, fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
return X, T
Smoothly increasing velocity (baseline)¶
A smooth velocity increase with depth is the "normal" case for most of the Earth's mantle. It produces a smooth, concave-down T(X) curve — no surprises.
In [14]:
Copied!
v_smooth = lambda z: 4.0 + 0.1 * z # linear gradient
X1, T1 = plot_model_rays_TX(v_smooth, 50, 'Baseline: Smooth velocity increase')
v_smooth = lambda z: 4.0 + 0.1 * z # linear gradient
X1, T1 = plot_model_rays_TX(v_smooth, 50, 'Baseline: Smooth velocity increase')
Sharp velocity jump → Triplication¶
When there is a large, abrupt velocity increase at a boundary, three different ray paths can arrive at the same distance:
- A direct ray that turns above the boundary
- A refracted ray that dives below the boundary
- A head wave traveling along the boundary
This produces a triplication — a fold in the T(X) curve where there are multiple arrivals.
Real Earth example: The Moho (crust-mantle boundary) has a velocity jump from ~6.5 to ~8 km/s and produces a classic triplication at regional distances.
In [15]:
Copied!
def v_jump(z):
"""Velocity model with a sharp jump at z=20 km (like the Moho)."""
z = np.atleast_1d(z)
v = np.where(z < 20, 4.0 + 0.05 * z, 7.5 + 0.02 * (z - 20))
return v
X2, T2 = plot_model_rays_TX(v_jump, 60, 'Sharp velocity jump → Triplication', n_rays=12)
def v_jump(z):
"""Velocity model with a sharp jump at z=20 km (like the Moho)."""
z = np.atleast_1d(z)
v = np.where(z < 20, 4.0 + 0.05 * z, 7.5 + 0.02 * (z - 20))
return v
X2, T2 = plot_model_rays_TX(v_jump, 60, 'Sharp velocity jump → Triplication', n_rays=12)
Low velocity layer → Shadow zone¶
When a low velocity zone (LVZ) is sandwiched between faster layers, rays bend toward the vertical when entering it. Rays that would normally sample a range of distances instead skip over the LVZ, creating a shadow zone — a range of distances where no direct rays arrive.
Real Earth examples:
- The asthenosphere (~100-200 km depth): a slight velocity decrease causing a minor shadow zone
- The outer core (liquid): S-waves cannot enter, creating the S-wave shadow zone
In [16]:
Copied!
def v_lvz(z):
"""Velocity model with a low velocity zone at 15-25 km depth."""
z = np.atleast_1d(z)
v = np.piecewise(z,
[z < 15, (z >= 15) & (z < 25), z >= 25],
[lambda z: 4.0 + 0.1 * z, # above LVZ: increasing
lambda z: 5.5 - 0.06 * (z - 15), # LVZ: decreasing
lambda z: 5.0 + 0.08 * (z - 25)] # below LVZ: increasing again
)
return v
X3, T3 = plot_model_rays_TX(v_lvz, 60, 'Low velocity zone → Shadow zone', n_rays=15)
def v_lvz(z):
"""Velocity model with a low velocity zone at 15-25 km depth."""
z = np.atleast_1d(z)
v = np.piecewise(z,
[z < 15, (z >= 15) & (z < 25), z >= 25],
[lambda z: 4.0 + 0.1 * z, # above LVZ: increasing
lambda z: 5.5 - 0.06 * (z - 15), # LVZ: decreasing
lambda z: 5.0 + 0.08 * (z - 25)] # below LVZ: increasing again
)
return v
X3, T3 = plot_model_rays_TX(v_lvz, 60, 'Low velocity zone → Shadow zone', n_rays=15)
Side-by-side comparison¶
Comparing all three cases together makes the connection between velocity structure and T(X) shape clear.
In [17]:
Copied!
models = [
(v_smooth, 50, 'Smooth increase'),
(v_jump, 60, 'Sharp velocity jump'),
(v_lvz, 60, 'Low velocity zone'),
]
fig, axes = plt.subplots(2, 3, figsize=(16, 9))
for j, (v_func, z_max, label) in enumerate(models):
# Top row: velocity models
z_plot = np.linspace(0, z_max, 500)
axes[0, j].plot(v_func(z_plot), z_plot, 'b-', lw=2)
axes[0, j].set_xlabel('Velocity (km/s)')
axes[0, j].set_title(label, fontsize=12)
axes[0, j].invert_yaxis()
axes[0, j].grid(True, alpha=0.3)
if j == 0:
axes[0, j].set_ylabel('Depth (km)')
# Bottom row: T(X) curves — scatter to show gaps
X, T, _ = sweep_TX(v_func, z_max)
axes[1, j].scatter(X, T, c='b', s=4)
axes[1, j].set_xlabel('Distance X (km)')
axes[1, j].grid(True, alpha=0.3)
if j == 0:
axes[1, j].set_ylabel('Travel time T (s)')
# Annotate
axes[1, 0].set_title('Smooth T(X)', fontsize=11)
axes[1, 1].set_title('Triplication', fontsize=11)
axes[1, 2].set_title('Shadow zone', fontsize=11)
plt.tight_layout()
plt.show()
models = [
(v_smooth, 50, 'Smooth increase'),
(v_jump, 60, 'Sharp velocity jump'),
(v_lvz, 60, 'Low velocity zone'),
]
fig, axes = plt.subplots(2, 3, figsize=(16, 9))
for j, (v_func, z_max, label) in enumerate(models):
# Top row: velocity models
z_plot = np.linspace(0, z_max, 500)
axes[0, j].plot(v_func(z_plot), z_plot, 'b-', lw=2)
axes[0, j].set_xlabel('Velocity (km/s)')
axes[0, j].set_title(label, fontsize=12)
axes[0, j].invert_yaxis()
axes[0, j].grid(True, alpha=0.3)
if j == 0:
axes[0, j].set_ylabel('Depth (km)')
# Bottom row: T(X) curves — scatter to show gaps
X, T, _ = sweep_TX(v_func, z_max)
axes[1, j].scatter(X, T, c='b', s=4)
axes[1, j].set_xlabel('Distance X (km)')
axes[1, j].grid(True, alpha=0.3)
if j == 0:
axes[1, j].set_ylabel('Travel time T (s)')
# Annotate
axes[1, 0].set_title('Smooth T(X)', fontsize=11)
axes[1, 1].set_title('Triplication', fontsize=11)
axes[1, 2].set_title('Shadow zone', fontsize=11)
plt.tight_layout()
plt.show()
6. Travel Times Through the Real Earth¶
We've built all the tools for understanding travel times using simple toy models. Now we connect to the real Earth — using the PREM velocity model and ObsPy's TauP to compute travel times for actual seismic phases, and compare with real earthquake data.
The PREM velocity model¶
The Preliminary Reference Earth Model (PREM, Dziewonski & Anderson, 1981) is the standard 1D model of the Earth's interior. Let's load it and see the velocity structure we've been learning about.
In [18]:
Copied!
# Load PREM through TauP
model = TauPyModel(model="prem")
# Extract velocity model layers
tau_model = model.model
s_mod = tau_model.s_mod
v_mod = s_mod.v_mod
# Get layer properties
depths = []
vp_vals = []
vs_vals = []
for layer in v_mod.layers:
depths.extend([layer['top_depth'], layer['bot_depth']])
vp_vals.extend([layer['top_p_velocity'], layer['bot_p_velocity']])
vs_vals.extend([layer['top_s_velocity'], layer['bot_s_velocity']])
depths = np.array(depths)
vp_vals = np.array(vp_vals)
vs_vals = np.array(vs_vals)
# Load PREM through TauP
model = TauPyModel(model="prem")
# Extract velocity model layers
tau_model = model.model
s_mod = tau_model.s_mod
v_mod = s_mod.v_mod
# Get layer properties
depths = []
vp_vals = []
vs_vals = []
for layer in v_mod.layers:
depths.extend([layer['top_depth'], layer['bot_depth']])
vp_vals.extend([layer['top_p_velocity'], layer['bot_p_velocity']])
vs_vals.extend([layer['top_s_velocity'], layer['bot_s_velocity']])
depths = np.array(depths)
vp_vals = np.array(vp_vals)
vs_vals = np.array(vs_vals)
In [19]:
Copied!
fig, ax = plt.subplots(figsize=(8, 10))
ax.plot(vp_vals, depths, 'b-', lw=2, label='$V_P$')
ax.plot(vs_vals, depths, 'r-', lw=2, label='$V_S$')
# Mark major discontinuities
boundaries = [
(35, 'Moho'),
(2891, 'Core-Mantle\nBoundary'),
(5150, 'Inner Core\nBoundary'),
]
for depth, name in boundaries:
ax.axhline(depth, color='gray', ls='--', alpha=0.5)
ax.text(0.3, depth, f' {name} ({depth} km)', fontsize=10, va='bottom')
# Highlight the low velocity zone
ax.axhspan(80, 220, alpha=0.1, color='orange')
ax.text(13.5, 150, 'LVZ', fontsize=10, ha='right', color='orange')
ax.set_xlabel('Velocity (km/s)', fontsize=12)
ax.set_ylabel('Depth (km)', fontsize=12)
ax.set_title('PREM Velocity Model', fontsize=14)
ax.invert_yaxis()
ax.legend(fontsize=12, loc='lower left')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 14.5)
plt.tight_layout()
plt.show()
print("Key features to recognize from our toy models:")
print(" • Moho (~35 km): Sharp velocity jump → triplication")
print(" • Low velocity zone (~80-220 km): Slight decrease → minor shadow zone")
print(" • Core-mantle boundary (~2891 km): Huge Vp drop, Vs → 0 → P-wave shadow zone")
print(" • Inner core boundary (~5150 km): Vp increase, Vs reappears")
fig, ax = plt.subplots(figsize=(8, 10))
ax.plot(vp_vals, depths, 'b-', lw=2, label='$V_P$')
ax.plot(vs_vals, depths, 'r-', lw=2, label='$V_S$')
# Mark major discontinuities
boundaries = [
(35, 'Moho'),
(2891, 'Core-Mantle\nBoundary'),
(5150, 'Inner Core\nBoundary'),
]
for depth, name in boundaries:
ax.axhline(depth, color='gray', ls='--', alpha=0.5)
ax.text(0.3, depth, f' {name} ({depth} km)', fontsize=10, va='bottom')
# Highlight the low velocity zone
ax.axhspan(80, 220, alpha=0.1, color='orange')
ax.text(13.5, 150, 'LVZ', fontsize=10, ha='right', color='orange')
ax.set_xlabel('Velocity (km/s)', fontsize=12)
ax.set_ylabel('Depth (km)', fontsize=12)
ax.set_title('PREM Velocity Model', fontsize=14)
ax.invert_yaxis()
ax.legend(fontsize=12, loc='lower left')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 14.5)
plt.tight_layout()
plt.show()
print("Key features to recognize from our toy models:")
print(" • Moho (~35 km): Sharp velocity jump → triplication")
print(" • Low velocity zone (~80-220 km): Slight decrease → minor shadow zone")
print(" • Core-mantle boundary (~2891 km): Huge Vp drop, Vs → 0 → P-wave shadow zone")
print(" • Inner core boundary (~5150 km): Vp increase, Vs reappears")
Key features to recognize from our toy models: • Moho (~35 km): Sharp velocity jump → triplication • Low velocity zone (~80-220 km): Slight decrease → minor shadow zone • Core-mantle boundary (~2891 km): Huge Vp drop, Vs → 0 → P-wave shadow zone • Inner core boundary (~5150 km): Vp increase, Vs reappears
Computing travel times with TauP¶
Our flat-layer code doesn't account for Earth's spherical geometry. For global-scale ray tracing, we use ObsPy's TauP — it solves the same physics (Snell's law, ray parameter) but in spherical coordinates.
Major seismic phases¶
| Phase | Path description |
|---|---|
| P, S | Direct waves through the mantle |
| PcP, ScS | Reflected off the core-mantle boundary |
| PKP | P-wave through the outer core |
| PKIKP | P-wave through the inner core |
| PP, SS | Waves that bounce once off the surface |
In [20]:
Copied!
# Compute travel times for all major phases at many distances
source_depth = 100 # km (moderately deep earthquake)
distances = np.arange(0, 181, 1) # epicentral distance in degrees
phases_to_plot = [
('P', 'blue', '-', 2.0),
('S', 'red', '-', 2.0),
('PcP', 'cyan', '--', 1.5),
('ScS', 'orange', '--', 1.5),
('PP', 'blue', ':', 1.5),
('SS', 'red', ':', 1.5),
('PKP', 'green', '-', 2.0),
('PKIKP', 'purple', '-', 2.0),
]
fig, ax = plt.subplots(figsize=(14, 10))
for phase_name, color, ls, lw in phases_to_plot:
phase_dist = []
phase_time = []
for dist in distances:
arrivals = model.get_travel_times(source_depth_in_km=source_depth,
distance_in_degree=dist,
phase_list=[phase_name])
for arr in arrivals:
phase_dist.append(dist)
phase_time.append(arr.time / 60) # convert to minutes
if phase_dist:
ax.scatter(phase_dist, phase_time, c=color, s=3, zorder=3)
# Add label at the last point
ax.text(phase_dist[-1] + 1, phase_time[-1], phase_name,
color=color, fontsize=10, fontweight='bold', va='center')
# Highlight P-wave shadow zone
ax.axvspan(100, 140, alpha=0.1, color='blue')
ax.text(120, 2, 'P-wave\nshadow zone', ha='center', fontsize=11, color='blue', alpha=0.7)
ax.set_xlabel('Epicentral distance (degrees)', fontsize=13)
ax.set_ylabel('Travel time (minutes)', fontsize=13)
ax.set_title(f'Travel time curves for major seismic phases (source depth = {source_depth} km)', fontsize=14)
ax.set_xlim(0, 185)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Compute travel times for all major phases at many distances
source_depth = 100 # km (moderately deep earthquake)
distances = np.arange(0, 181, 1) # epicentral distance in degrees
phases_to_plot = [
('P', 'blue', '-', 2.0),
('S', 'red', '-', 2.0),
('PcP', 'cyan', '--', 1.5),
('ScS', 'orange', '--', 1.5),
('PP', 'blue', ':', 1.5),
('SS', 'red', ':', 1.5),
('PKP', 'green', '-', 2.0),
('PKIKP', 'purple', '-', 2.0),
]
fig, ax = plt.subplots(figsize=(14, 10))
for phase_name, color, ls, lw in phases_to_plot:
phase_dist = []
phase_time = []
for dist in distances:
arrivals = model.get_travel_times(source_depth_in_km=source_depth,
distance_in_degree=dist,
phase_list=[phase_name])
for arr in arrivals:
phase_dist.append(dist)
phase_time.append(arr.time / 60) # convert to minutes
if phase_dist:
ax.scatter(phase_dist, phase_time, c=color, s=3, zorder=3)
# Add label at the last point
ax.text(phase_dist[-1] + 1, phase_time[-1], phase_name,
color=color, fontsize=10, fontweight='bold', va='center')
# Highlight P-wave shadow zone
ax.axvspan(100, 140, alpha=0.1, color='blue')
ax.text(120, 2, 'P-wave\nshadow zone', ha='center', fontsize=11, color='blue', alpha=0.7)
ax.set_xlabel('Epicentral distance (degrees)', fontsize=13)
ax.set_ylabel('Travel time (minutes)', fontsize=13)
ax.set_title(f'Travel time curves for major seismic phases (source depth = {source_depth} km)', fontsize=14)
ax.set_xlim(0, 185)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
The P-wave shadow zone¶
Notice the gap in P-wave arrivals between ~100° and ~140° — this is the P-wave shadow zone caused by the liquid outer core. The core refracts P-waves so strongly that no direct P arrivals reach this distance range. The only arrivals in this zone come from waves that travel through the core (PKP, PKIKP).
This is exactly the low-velocity layer effect we saw in our toy model — but at a planetary scale!
Visualizing ray paths through the Earth¶
TauP can also show us the ray paths — how the waves travel through the Earth's interior.
In [21]:
Copied!
examples = [
('P', 30, 'P at 30°'),
('P', 90, 'P at 90°'),
('PcP', 40, 'PcP at 40°'),
('PKP', 150, 'PKP at 150°'),
('PKIKP',150, 'PKIKP at 150°'),
('PP', 90, 'PP at 90°'),
]
fig, axes = plt.subplots(2, 3, figsize=(16, 10), subplot_kw={'polar': True})
for ax, (phase, dist, title) in zip(axes.flatten(), examples):
arrivals = model.get_ray_paths(source_depth_in_km=source_depth,
distance_in_degree=dist,
phase_list=[phase])
if len(arrivals) > 0:
arrivals.plot_rays(plot_type='spherical', ax=ax, show=False, legend=False)
ax.set_title(title, fontsize=12, pad=15)
plt.tight_layout()
plt.show()
examples = [
('P', 30, 'P at 30°'),
('P', 90, 'P at 90°'),
('PcP', 40, 'PcP at 40°'),
('PKP', 150, 'PKP at 150°'),
('PKIKP',150, 'PKIKP at 150°'),
('PP', 90, 'PP at 90°'),
]
fig, axes = plt.subplots(2, 3, figsize=(16, 10), subplot_kw={'polar': True})
for ax, (phase, dist, title) in zip(axes.flatten(), examples):
arrivals = model.get_ray_paths(source_depth_in_km=source_depth,
distance_in_degree=dist,
phase_list=[phase])
if len(arrivals) > 0:
arrivals.plot_rays(plot_type='spherical', ax=ax, show=False, legend=False)
ax.set_title(title, fontsize=12, pad=15)
plt.tight_layout()
plt.show()
All ray paths together¶
Plotting many ray paths for each phase on a single Earth cross-section reveals the overall geometry — where each phase samples, the shadow zones, and how the core bends and reflects waves.
In [22]:
Copied!
R_earth = 6371 # km
def plot_earth_layers(ax):
"""Draw Earth's major boundaries on a polar axis."""
for r in [R_earth, R_earth - 35, R_earth - 2891, R_earth - 5150]:
theta = np.linspace(0, 2 * np.pi, 300)
ax.plot(theta, np.full_like(theta, r), 'k-', lw=0.8, alpha=0.4)
# Shade outer core and inner core
ax.fill_between(np.linspace(0, 2*np.pi, 300),
0, R_earth - 5150, alpha=0.15, color='gold', label='Inner core')
ax.fill_between(np.linspace(0, 2*np.pi, 300),
R_earth - 5150, R_earth - 2891, alpha=0.10, color='orange', label='Outer core')
ax.set_ylim(0, R_earth)
ax.set_yticks([])
ax.set_xticks([])
def plot_rays_on_ax(ax, model, phase_list, distances, source_depth, colors):
"""Plot ray paths for multiple phases on a single polar axis."""
for phase_name, color in zip(phase_list, colors):
for dist in distances:
arrivals = model.get_ray_paths(source_depth_in_km=source_depth,
distance_in_degree=dist,
phase_list=[phase_name])
for arr in arrivals:
path = arr.path
theta = np.array([p[2] for p in path]) # distance in radians
r = R_earth - np.array([p[3] for p in path]) # radius
ax.plot(theta, r, color=color, lw=0.8, alpha=0.6)
# Dummy line for legend
ax.plot([], [], color=color, lw=2, label=phase_name)
# Same color palette for both panels (separate plots, no ambiguity)
phase_colors = ['blue', 'cyan', 'cornflowerblue', 'green', 'purple']
# --- P-type phases ---
p_phases = ['P', 'PcP', 'PP', 'PKP', 'PKIKP']
p_dists = {
'P': np.arange(10, 100, 5),
'PcP': np.arange(10, 80, 5),
'PP': np.arange(40, 170, 10),
'PKP': np.arange(140, 180, 5),
'PKIKP': np.arange(120, 180, 5),
}
# --- S-type phases ---
s_phases = ['S', 'ScS', 'SS', 'SKS']
s_dists = {
'S': np.arange(10, 100, 5),
'ScS': np.arange(10, 80, 5),
'SS': np.arange(40, 170, 10),
'SKS': np.arange(80, 140, 5),
}
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8), subplot_kw={'polar': True})
for ax, phases, dists, title in [
(ax1, p_phases, p_dists, 'P-wave phases'),
(ax2, s_phases, s_dists, 'S-wave phases'),
]:
plot_earth_layers(ax)
for phase, color in zip(phases, phase_colors):
plot_rays_on_ax(ax, model, [phase], dists[phase], source_depth, [color])
ax.set_title(title, fontsize=14, pad=20)
ax.legend(loc='lower left', fontsize=10, framealpha=0.9)
ax.plot(0, R_earth - source_depth, 'r*', markersize=15, zorder=5)
plt.tight_layout()
plt.show()
R_earth = 6371 # km
def plot_earth_layers(ax):
"""Draw Earth's major boundaries on a polar axis."""
for r in [R_earth, R_earth - 35, R_earth - 2891, R_earth - 5150]:
theta = np.linspace(0, 2 * np.pi, 300)
ax.plot(theta, np.full_like(theta, r), 'k-', lw=0.8, alpha=0.4)
# Shade outer core and inner core
ax.fill_between(np.linspace(0, 2*np.pi, 300),
0, R_earth - 5150, alpha=0.15, color='gold', label='Inner core')
ax.fill_between(np.linspace(0, 2*np.pi, 300),
R_earth - 5150, R_earth - 2891, alpha=0.10, color='orange', label='Outer core')
ax.set_ylim(0, R_earth)
ax.set_yticks([])
ax.set_xticks([])
def plot_rays_on_ax(ax, model, phase_list, distances, source_depth, colors):
"""Plot ray paths for multiple phases on a single polar axis."""
for phase_name, color in zip(phase_list, colors):
for dist in distances:
arrivals = model.get_ray_paths(source_depth_in_km=source_depth,
distance_in_degree=dist,
phase_list=[phase_name])
for arr in arrivals:
path = arr.path
theta = np.array([p[2] for p in path]) # distance in radians
r = R_earth - np.array([p[3] for p in path]) # radius
ax.plot(theta, r, color=color, lw=0.8, alpha=0.6)
# Dummy line for legend
ax.plot([], [], color=color, lw=2, label=phase_name)
# Same color palette for both panels (separate plots, no ambiguity)
phase_colors = ['blue', 'cyan', 'cornflowerblue', 'green', 'purple']
# --- P-type phases ---
p_phases = ['P', 'PcP', 'PP', 'PKP', 'PKIKP']
p_dists = {
'P': np.arange(10, 100, 5),
'PcP': np.arange(10, 80, 5),
'PP': np.arange(40, 170, 10),
'PKP': np.arange(140, 180, 5),
'PKIKP': np.arange(120, 180, 5),
}
# --- S-type phases ---
s_phases = ['S', 'ScS', 'SS', 'SKS']
s_dists = {
'S': np.arange(10, 100, 5),
'ScS': np.arange(10, 80, 5),
'SS': np.arange(40, 170, 10),
'SKS': np.arange(80, 140, 5),
}
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8), subplot_kw={'polar': True})
for ax, phases, dists, title in [
(ax1, p_phases, p_dists, 'P-wave phases'),
(ax2, s_phases, s_dists, 'S-wave phases'),
]:
plot_earth_layers(ax)
for phase, color in zip(phases, phase_colors):
plot_rays_on_ax(ax, model, [phase], dists[phase], source_depth, [color])
ax.set_title(title, fontsize=14, pad=20)
ax.legend(loc='lower left', fontsize=10, framealpha=0.9)
ax.plot(0, R_earth - source_depth, 'r*', markersize=15, zorder=5)
plt.tight_layout()
plt.show()
Comparison with real data¶
Let's download waveforms from a real earthquake recorded by the global IU network and overlay theoretical travel time curves for multiple phases.
In [23]:
Copied!
# Find a large, deep earthquake (clear phases, well-recorded globally)
client = Client("IRIS")
catalog = client.get_events(
starttime=UTCDateTime() - 365 * 24 * 3600,
endtime=UTCDateTime(),
minmagnitude=7.0,
mindepth=50,
orderby="magnitude",
limit=5
)
print("Recent large earthquakes:")
for i, ev in enumerate(catalog):
org = ev.origins[0]
mag = ev.magnitudes[0]
print(f" [{i}] M{mag.mag:.1f} - Depth {org.depth/1000:.0f} km - {org.time.strftime('%Y-%m-%d')}")
# Find a large, deep earthquake (clear phases, well-recorded globally)
client = Client("IRIS")
catalog = client.get_events(
starttime=UTCDateTime() - 365 * 24 * 3600,
endtime=UTCDateTime(),
minmagnitude=7.0,
mindepth=50,
orderby="magnitude",
limit=5
)
print("Recent large earthquakes:")
for i, ev in enumerate(catalog):
org = ev.origins[0]
mag = ev.magnitudes[0]
print(f" [{i}] M{mag.mag:.1f} - Depth {org.depth/1000:.0f} km - {org.time.strftime('%Y-%m-%d')}")
Recent large earthquakes: [0] M7.5 - Depth 229 km - 2026-03-24 [1] M7.4 - Depth 59 km - 2025-10-10 [2] M7.3 - Depth 121 km - 2026-03-30 [3] M7.0 - Depth 629 km - 2026-02-22
In [24]:
Copied!
# Select the first earthquake
event = catalog[0]
origin = event.origins[0]
eq_lat = origin.latitude
eq_lon = origin.longitude
eq_depth = origin.depth / 1000
eq_time = origin.time
eq_mag = event.magnitudes[0].mag
print(f"Selected: M{eq_mag:.1f}, depth {eq_depth:.0f} km, {eq_time.strftime('%Y-%m-%d')}")
# Download waveforms from global stations (IU network)
print("\nDownloading waveforms from global IU network...")
try:
st = client.get_waveforms(
network="IU", station="*", location="00", channel="BHZ",
starttime=eq_time - 60,
endtime=eq_time + 30 * 60, # 30 minutes of data
)
print(f"Downloaded {len(st)} traces")
except Exception as e:
print(f"Download issue: {e}")
st = None
# Select the first earthquake
event = catalog[0]
origin = event.origins[0]
eq_lat = origin.latitude
eq_lon = origin.longitude
eq_depth = origin.depth / 1000
eq_time = origin.time
eq_mag = event.magnitudes[0].mag
print(f"Selected: M{eq_mag:.1f}, depth {eq_depth:.0f} km, {eq_time.strftime('%Y-%m-%d')}")
# Download waveforms from global stations (IU network)
print("\nDownloading waveforms from global IU network...")
try:
st = client.get_waveforms(
network="IU", station="*", location="00", channel="BHZ",
starttime=eq_time - 60,
endtime=eq_time + 30 * 60, # 30 minutes of data
)
print(f"Downloaded {len(st)} traces")
except Exception as e:
print(f"Download issue: {e}")
st = None
Selected: M7.5, depth 229 km, 2026-03-24 Downloading waveforms from global IU network...
Downloaded 62 traces
In [25]:
Copied!
if st and len(st) > 0:
# Get station locations
inv = client.get_stations(
network="IU", station="*", channel="BHZ",
starttime=eq_time, endtime=eq_time + 60,
level="station"
)
taup_model = TauPyModel(model="iasp91")
fig, ax = plt.subplots(figsize=(14, 10))
for tr in st:
try:
sta_inv = inv.select(station=tr.stats.station)
sta_lat = sta_inv[0][0].latitude
sta_lon = sta_inv[0][0].longitude
except Exception:
continue
dist_m, az, baz = gps2dist_azimuth(eq_lat, eq_lon, sta_lat, sta_lon)
dist_deg = dist_m / 1000 / 111.19
if dist_deg < 10 or dist_deg > 100:
continue
tr_proc = tr.copy()
tr_proc.detrend('demean')
tr_proc.taper(max_percentage=0.05)
tr_proc.filter('bandpass', freqmin=0.5, freqmax=2.0)
times = tr_proc.times() + (tr_proc.stats.starttime - eq_time)
data_norm = tr_proc.data / np.max(np.abs(tr_proc.data)) * 2
ax.plot(times / 60, data_norm + dist_deg, 'k-', lw=0.4)
# Overlay theoretical travel times — only first arrival per distance to avoid zigzags
phases_record = [
('P', 'blue', '-', 2.0),
('pP', 'blue', '--', 1.5),
('S', 'red', '-', 2.0),
('sS', 'red', '--', 1.5),
('PcP', 'cyan', '-', 1.5),
('ScS', 'orange', '-', 1.5),
('PP', 'blue', ':', 1.5),
('SS', 'red', ':', 1.5),
]
dist_range = np.arange(10, 101, 0.5)
for phase_name, color, ls, lw in phases_record:
d_arr, t_arr = [], []
for dist in dist_range:
arrivals = taup_model.get_travel_times(source_depth_in_km=eq_depth,
distance_in_degree=dist,
phase_list=[phase_name])
if arrivals:
d_arr.append(dist)
t_arr.append(arrivals[0].time / 60) # first arrival only
if d_arr:
ax.plot(t_arr, d_arr, color=color, ls=ls, lw=lw, label=phase_name)
ax.set_xlabel('Time since origin (minutes)', fontsize=13)
ax.set_ylabel('Epicentral distance (degrees)', fontsize=13)
ax.set_title(f'Record section: M{eq_mag:.1f}, depth {eq_depth:.0f} km\n'
f'IU network, filtered 0.5-2 Hz', fontsize=14)
ax.set_xlim(0, 30)
ax.legend(fontsize=10, loc='lower right', ncol=2)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
else:
print("Could not download waveforms.")
if st and len(st) > 0:
# Get station locations
inv = client.get_stations(
network="IU", station="*", channel="BHZ",
starttime=eq_time, endtime=eq_time + 60,
level="station"
)
taup_model = TauPyModel(model="iasp91")
fig, ax = plt.subplots(figsize=(14, 10))
for tr in st:
try:
sta_inv = inv.select(station=tr.stats.station)
sta_lat = sta_inv[0][0].latitude
sta_lon = sta_inv[0][0].longitude
except Exception:
continue
dist_m, az, baz = gps2dist_azimuth(eq_lat, eq_lon, sta_lat, sta_lon)
dist_deg = dist_m / 1000 / 111.19
if dist_deg < 10 or dist_deg > 100:
continue
tr_proc = tr.copy()
tr_proc.detrend('demean')
tr_proc.taper(max_percentage=0.05)
tr_proc.filter('bandpass', freqmin=0.5, freqmax=2.0)
times = tr_proc.times() + (tr_proc.stats.starttime - eq_time)
data_norm = tr_proc.data / np.max(np.abs(tr_proc.data)) * 2
ax.plot(times / 60, data_norm + dist_deg, 'k-', lw=0.4)
# Overlay theoretical travel times — only first arrival per distance to avoid zigzags
phases_record = [
('P', 'blue', '-', 2.0),
('pP', 'blue', '--', 1.5),
('S', 'red', '-', 2.0),
('sS', 'red', '--', 1.5),
('PcP', 'cyan', '-', 1.5),
('ScS', 'orange', '-', 1.5),
('PP', 'blue', ':', 1.5),
('SS', 'red', ':', 1.5),
]
dist_range = np.arange(10, 101, 0.5)
for phase_name, color, ls, lw in phases_record:
d_arr, t_arr = [], []
for dist in dist_range:
arrivals = taup_model.get_travel_times(source_depth_in_km=eq_depth,
distance_in_degree=dist,
phase_list=[phase_name])
if arrivals:
d_arr.append(dist)
t_arr.append(arrivals[0].time / 60) # first arrival only
if d_arr:
ax.plot(t_arr, d_arr, color=color, ls=ls, lw=lw, label=phase_name)
ax.set_xlabel('Time since origin (minutes)', fontsize=13)
ax.set_ylabel('Epicentral distance (degrees)', fontsize=13)
ax.set_title(f'Record section: M{eq_mag:.1f}, depth {eq_depth:.0f} km\n'
f'IU network, filtered 0.5-2 Hz', fontsize=14)
ax.set_xlim(0, 30)
ax.legend(fontsize=10, loc='lower right', ncol=2)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
else:
print("Could not download waveforms.")
Summary¶
| Concept | Key result |
|---|---|
| Snell's law | $\sin\theta/v = p$ (ray parameter) is conserved along the ray path |
| Travel time formulas | $\Delta x = p\Delta z/\eta$, $\Delta t = u^2\Delta z/\eta$ where $\eta = \sqrt{u^2 - p^2}$ |
| Continuous models | Sums become integrals; linear gradient → circular arc ray paths |
| Velocity structure | Sharp jumps → triplications; low velocity zones → shadow zones |
| PREM | Standard 1D Earth model with Moho, LVZ, CMB, ICB |
| Seismic phases | P, S, PcP, ScS, PP, SS, PKP, PKIKP — each probes different depths |
| P-wave shadow zone | 100°–140°: caused by the liquid outer core |
The big picture¶
Everything connects:
- Snell's law governs ray bending at interfaces
- The ray parameter is conserved and controls the geometry
- Travel time formulas let us compute T(X) for any layered model
- Continuous models extend this to realistic velocity gradients
- Velocity structure creates distinctive T(X) signatures (triplications, shadow zones)
- The real Earth exhibits all of these — and matching observed travel times to models is how we discovered the Earth's internal structure