Travel Times — Exercises¶

Complete the following exercises to practice the concepts from the travel times lecture.

In [ ]:

Copied!

!pip install obspy -q
!pip install obspy -q

In [ ]:

Copied!





import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from obspy.taup import TauPyModel


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
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from obspy.taup import TauPyModel


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

Exercise 1¶

Recap: Snell's law says $\sin\theta / v = p$ (ray parameter, constant along a ray). The critical angle $\theta_c = \arcsin(v_1/v_2)$ is where the refracted ray travels horizontally along the interface.: Snell's Law & Critical Angle

Consider two layers with $v_1 = 5$ km/s and $v_2 = 8$ km/s.

Recall: $\frac{\sin\theta_1}{v_1} = \frac{\sin\theta_2}{v_2} = p$

(a) Compute the critical angle $\theta_c$ and the critical ray parameter $p_c$.

(b) For an incident angle of $\theta_1 = 30°$, compute the ray parameter $p$ and the refracted angle $\theta_2$.

In [ ]:

Copied!





v1 = 5.0  # km/s
v2 = 8.0  # km/s

# (a) Critical angle and critical ray parameter
# Recall: θ_c = arcsin(v1/v2), p_c = 1/v2
theta_c = ???  # arcsin(v1/v2), convert to degrees
p_c = ???  # 1/v2 (or sin(theta_c)/v1)

print(f"Critical angle: θ_c = {theta_c:.2f}°")
print(f"Critical ray parameter: p_c = {p_c:.4f} s/km")
v1 = 5.0  # km/s
v2 = 8.0  # km/s

# (a) Critical angle and critical ray parameter
# Recall: θ_c = arcsin(v1/v2), p_c = 1/v2
theta_c = ???  # arcsin(v1/v2), convert to degrees
p_c = ???  # 1/v2 (or sin(theta_c)/v1)

print(f"Critical angle: θ_c = {theta_c:.2f}°")
print(f"Critical ray parameter: p_c = {p_c:.4f} s/km")

In [ ]:

Copied!





# (b) For incident angle θ₁ = 30°
theta1 = 30.0  # degrees

# Compute ray parameter: p = sin(θ₁) / v₁
p = ???  # sin(theta1)/v1

# Compute refracted angle: θ₂ = arcsin(p × v₂)
theta2 = ???  # arcsin(p * v2), convert to degrees

print(f"Ray parameter: p = {p:.4f} s/km")
print(f"Refracted angle: θ₂ = {theta2:.2f}°")
print(f"\nVerify Snell's law:")
print(f"  sin(θ₁)/v₁ = {np.sin(np.radians(theta1))/v1:.4f}")
print(f"  sin(θ₂)/v₂ = {np.sin(np.radians(theta2))/v2:.4f}")
# (b) For incident angle θ₁ = 30°
theta1 = 30.0  # degrees

# Compute ray parameter: p = sin(θ₁) / v₁
p = ???  # sin(theta1)/v1

# Compute refracted angle: θ₂ = arcsin(p × v₂)
theta2 = ???  # arcsin(p * v2), convert to degrees

print(f"Ray parameter: p = {p:.4f} s/km")
print(f"Refracted angle: θ₂ = {theta2:.2f}°")
print(f"\nVerify Snell's law:")
print(f"  sin(θ₁)/v₁ = {np.sin(np.radians(theta1))/v1:.4f}")
print(f"  sin(θ₂)/v₂ = {np.sin(np.radians(theta2))/v2:.4f}")

Question: If $v_2 < v_1$ (velocity decreases across the interface), is there a critical angle? What happens to the refracted ray?

Your answer here:

Exercise 2¶

Recap: For a stack of horizontal layers, the total travel time and distance are sums over layers: $X = \sum 2p\,\Delta z / \eta_k$ and $T = \sum 2u_k^2\,\Delta z / \eta_k$, where $\eta_k = \sqrt{u_k^2 - p^2}$ and $u_k = 1/v_k$.: Implement the Travel Time Formulas

Implement my_compute_XT(velocities, thicknesses, p) from scratch using the formulas:

$$\Delta x = \frac{p \, \Delta z}{\sqrt{u^2 - p^2}}, \qquad \Delta t = \frac{u^2 \, \Delta z}{\sqrt{u^2 - p^2}}$$

Remember: multiply by 2 for the two-way path, and the ray turns back when $p \geq u_i$.

Hint: Loop over layers. For each layer, check if $p < u_i$. If yes, compute $\eta = \sqrt{u^2 - p^2}$ and accumulate $\Delta x$ and $\Delta t$. If no, the ray has already turned — stop.

In [ ]:

Copied!





def my_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):
        # YOUR CODE HERE
        # 1. Check if p >= u → ray can't enter, break
        # 2. Compute eta = sqrt(u^2 - p^2)
        # 3. Accumulate X += 2 * p * dz / eta
        # 4. Accumulate T += 2 * u^2 * dz / eta
        pass
    return X, T
def my_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):
        # YOUR CODE HERE
        # 1. Check if p >= u → ray can't enter, break
        # 2. Compute eta = sqrt(u^2 - p^2)
        # 3. Accumulate X += 2 * p * dz / eta
        # 4. Accumulate T += 2 * u^2 * dz / eta
        pass
    return X, T

In [ ]:

Copied!





# Test your implementation against the provided compute_XT
velocities = [4.0, 6.0, 8.0]
thicknesses = [3.0, 3.0, 3.0]

for p_test in [0.05, 0.10, 0.15, 0.20]:
    X_ref, T_ref = compute_XT(velocities, thicknesses, p_test)
    X_my, T_my = my_compute_XT(velocities, thicknesses, p_test)
    match = "✓" if abs(X_my - X_ref) < 0.01 and abs(T_my - T_ref) < 0.01 else "✗"
    print(f"p = {p_test}: X = {X_my:.2f} (ref: {X_ref:.2f}), T = {T_my:.2f} (ref: {T_ref:.2f}) {match}")
# Test your implementation against the provided compute_XT
velocities = [4.0, 6.0, 8.0]
thicknesses = [3.0, 3.0, 3.0]

for p_test in [0.05, 0.10, 0.15, 0.20]:
    X_ref, T_ref = compute_XT(velocities, thicknesses, p_test)
    X_my, T_my = my_compute_XT(velocities, thicknesses, p_test)
    match = "✓" if abs(X_my - X_ref) < 0.01 and abs(T_my - T_ref) < 0.01 else "✗"
    print(f"p = {p_test}: X = {X_my:.2f} (ref: {X_ref:.2f}), T = {T_my:.2f} (ref: {T_ref:.2f}) {match}")

Question: For the 3-layer model $v = [4, 6, 8]$ km/s, what is the maximum ray parameter $p$ for which the ray still enters the second layer? What is the physical meaning of this value?

Your answer here:

Exercise 3¶

Recap: The T(X) curve shows travel time vs. epicentral distance. Different ray parameters $p$ produce different (X, T) pairs. Triplications appear when velocity decreases with depth.: Build a T(X) Curve

Approximate a linear velocity gradient ($v = 4$ to $8$ km/s over 30 km depth) by discretizing it into 10 thin layers. Sweep through many ray parameters and build the complete $T(X)$ curve.

Use the provided compute_XT function.

In [ ]:

Copied!





# Discretize a linear velocity gradient into thin layers
n_layers = 10
velocities = np.linspace(4.0, 8.0, n_layers)  # km/s
thicknesses = np.full(n_layers, 3.0)           # 3 km each

# Range of ray parameters
slownesses = 1.0 / np.array(velocities)
p_values = np.linspace(slownesses.min() + 1e-4, slownesses.max() - 1e-4, 300)

# TODO: Compute X and T for each p
X_arr, T_arr = [], []
for p in p_values:
    # YOUR CODE HERE: use compute_XT, append results
    pass

# Plot T(X)
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(X_arr, T_arr, c='b', s=4)
ax.set_xlabel('Distance X (km)')
ax.set_ylabel('Travel time T (s)')
ax.set_title('T(X) curve — discretized linear gradient')
ax.set_xlim(0, 100)
ax.set_ylim(0, 20)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Discretize a linear velocity gradient into thin layers
n_layers = 10
velocities = np.linspace(4.0, 8.0, n_layers)  # km/s
thicknesses = np.full(n_layers, 3.0)           # 3 km each

# Range of ray parameters
slownesses = 1.0 / np.array(velocities)
p_values = np.linspace(slownesses.min() + 1e-4, slownesses.max() - 1e-4, 300)

# TODO: Compute X and T for each p
X_arr, T_arr = [], []
for p in p_values:
    # YOUR CODE HERE: use compute_XT, append results
    pass

# Plot T(X)
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(X_arr, T_arr, c='b', s=4)
ax.set_xlabel('Distance X (km)')
ax.set_ylabel('Travel time T (s)')
ax.set_title('T(X) curve — discretized linear gradient')
ax.set_xlim(0, 100)
ax.set_ylim(0, 20)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Question: The slope of the $T(X)$ curve equals the ray parameter: $dT/dX = p$. At large distances, is the slope steeper or shallower — and why?

Your answer here:

Exercise 4: Ray Tracing in a Continuous Model¶

For a linear velocity gradient $v(z) = v_0 + kz$ with $v_0 = 6$ km/s and $k = 0.01$ s$^{-1}$:

(a) Compute the turning depth $z_p$ for $p = 0.14$ s/km analytically.

(b) Trace a ray numerically by discretizing the model into thin layers.

Recall: The turning depth is where $u(z_p) = p$, i.e., $v(z_p) = 1/p$.

In [ ]:

Copied!





v0 = 6.0    # km/s
k = 0.01    # 1/s
p = 0.14    # s/km

# (a) Turning depth: v(z_p) = 1/p → v0 + k*z_p = 1/p → z_p = ?
z_turn = ???  # analytical turning depth for linear v(z)

print(f"Turning depth for p = {p} s/km: z_p = {z_turn:.1f} km")
print(f"Check: v(z_p) = {v0 + k * z_turn:.2f} km/s = 1/p = {1/p:.2f} km/s")
v0 = 6.0    # km/s
k = 0.01    # 1/s
p = 0.14    # s/km

# (a) Turning depth: v(z_p) = 1/p → v0 + k*z_p = 1/p → z_p = ?
z_turn = ???  # analytical turning depth for linear v(z)

print(f"Turning depth for p = {p} s/km: z_p = {z_turn:.1f} km")
print(f"Check: v(z_p) = {v0 + k * z_turn:.2f} km/s = 1/p = {1/p:.2f} km/s")

In [ ]:

Copied!





# (b) Trace a ray numerically
z_max = 200
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)

# Trace the downgoing ray
x, z = [0.0], [0.0]
for u in u_layers:
    if p >= u:
        break
    # YOUR CODE HERE: compute horizontal offset dx for this layer
    dx = ???  # p * dz / sqrt(u**2 - p**2)
    x.append(x[-1] + dx)
    z.append(z[-1] + dz)

x, z = np.array(x), np.array(z)

# Mirror for the upgoing path
X_total = 2 * x[-1]
x_full = np.concatenate([x, X_total - x[::-1]])
z_full = np.concatenate([z, z[::-1]])

# Plot
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(x_full, z_full, 'b-', lw=2)
ax.axhline(z_turn, color='r', ls='--', label=f'Analytical turning depth = {z_turn:.0f} km')
ax.set_ylim(z_max, -1)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Depth (km)')
ax.set_title(f'Ray path: v(z) = {v0} + {k}z, p = {p} s/km')
ax.legend()
plt.tight_layout()
plt.show()

print(f"Total distance: X = {X_total:.1f} km")
print(f"Maximum depth reached: {z[-1]:.1f} km")
# (b) Trace a ray numerically
z_max = 200
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)

# Trace the downgoing ray
x, z = [0.0], [0.0]
for u in u_layers:
    if p >= u:
        break
    # YOUR CODE HERE: compute horizontal offset dx for this layer
    dx = ???  # p * dz / sqrt(u**2 - p**2)
    x.append(x[-1] + dx)
    z.append(z[-1] + dz)

x, z = np.array(x), np.array(z)

# Mirror for the upgoing path
X_total = 2 * x[-1]
x_full = np.concatenate([x, X_total - x[::-1]])
z_full = np.concatenate([z, z[::-1]])

# Plot
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(x_full, z_full, 'b-', lw=2)
ax.axhline(z_turn, color='r', ls='--', label=f'Analytical turning depth = {z_turn:.0f} km')
ax.set_ylim(z_max, -1)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Depth (km)')
ax.set_title(f'Ray path: v(z) = {v0} + {k}z, p = {p} s/km')
ax.legend()
plt.tight_layout()
plt.show()

print(f"Total distance: X = {X_total:.1f} km")
print(f"Maximum depth reached: {z[-1]:.1f} km")

Question: If you increase the gradient $k$ (velocity increases faster with depth), how would the ray paths change? Would they curve more or less?

Your answer here:

Exercise 5: Exploring Seismic Phases with TauP¶

Use ObsPy's TauP to explore how seismic phases travel through the real Earth.

(a) Compute P and S travel times at epicentral distances from 10° to 100°.

(b) Plot the S-P time as a function of distance — this is how seismologists estimate the distance to an earthquake.

In [ ]:

Copied!





model = TauPyModel(model="iasp91")
source_depth = 10  # km (shallow earthquake)

# (a) Compute P and S travel times at many distances
distances = np.arange(10, 101, 2)  # degrees

P_times = []
S_times = []
for dist in distances:
    # TODO: Use model.get_travel_times() to get P and S arrival times
    # Hint: phase_list=['P'] for P, phase_list=['S'] for S
    # arrivals[0].time gives the first arrival time in seconds
    p_arr = model.get_travel_times(source_depth_in_km=source_depth,
                                    distance_in_degree=dist,
                                    phase_list=['P'])
    s_arr = model.get_travel_times(source_depth_in_km=source_depth,
                                    distance_in_degree=dist,
                                    phase_list=['S'])
    P_times.append(p_arr[0].time if p_arr else np.nan)  # arrival time of first P phase
    S_times.append(s_arr[0].time if s_arr else np.nan)  # arrival time of first S phase

P_times = np.array(P_times)
S_times = np.array(S_times)

# Plot P and S travel time curves
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(distances, P_times / 60, 'b-', lw=2, label='P')
ax.plot(distances, S_times / 60, 'r-', lw=2, label='S')
ax.set_xlabel('Epicentral distance (degrees)')
ax.set_ylabel('Travel time (minutes)')
ax.set_title('P and S travel times (IASP91 model)')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
model = TauPyModel(model="iasp91")
source_depth = 10  # km (shallow earthquake)

# (a) Compute P and S travel times at many distances
distances = np.arange(10, 101, 2)  # degrees

P_times = []
S_times = []
for dist in distances:
    # TODO: Use model.get_travel_times() to get P and S arrival times
    # Hint: phase_list=['P'] for P, phase_list=['S'] for S
    # arrivals[0].time gives the first arrival time in seconds
    p_arr = model.get_travel_times(source_depth_in_km=source_depth,
                                    distance_in_degree=dist,
                                    phase_list=['P'])
    s_arr = model.get_travel_times(source_depth_in_km=source_depth,
                                    distance_in_degree=dist,
                                    phase_list=['S'])
    P_times.append(p_arr[0].time if p_arr else np.nan)  # arrival time of first P phase
    S_times.append(s_arr[0].time if s_arr else np.nan)  # arrival time of first S phase

P_times = np.array(P_times)
S_times = np.array(S_times)

# Plot P and S travel time curves
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(distances, P_times / 60, 'b-', lw=2, label='P')
ax.plot(distances, S_times / 60, 'r-', lw=2, label='S')
ax.set_xlabel('Epicentral distance (degrees)')
ax.set_ylabel('Travel time (minutes)')
ax.set_title('P and S travel times (IASP91 model)')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

In [ ]:

Copied!





# (b) Plot S-P time vs distance
# TODO: Compute the S-P time difference at each distance
sp_time = ???  # S_time - P_time at the chosen distance

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(distances, sp_time, 'g-', lw=2)
ax.set_xlabel('Epicentral distance (degrees)')
ax.set_ylabel('S-P time (seconds)')
ax.set_title('S-P time vs. epicentral distance')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# (b) Plot S-P time vs distance
# TODO: Compute the S-P time difference at each distance
sp_time = ???  # S_time - P_time at the chosen distance

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(distances, sp_time, 'g-', lw=2)
ax.set_xlabel('Epicentral distance (degrees)')
ax.set_ylabel('S-P time (seconds)')
ax.set_title('S-P time vs. epicentral distance')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Question: The S-P time increases roughly linearly with distance. A seismologist measures an S-P time of 300 seconds. Using your plot, approximately how far away is the earthquake (in degrees)?

Your answer here: