Traveltime module

Overview

The traveltime module provides numerical computation of P- and S-wave travel times in a horizontally layered 1D velocity model. It implements a robust ray-based solver that accounts for both direct and head-wave arrivals. The algorithm is fully analytical in layered segments and uses an iterative shooting method on the ray parameter \(p\) to match a given horizontal offset.

The method assumes a horizontally stratified medium and straight-ray propagation within each layer. Reflection phases are not included, and refraction (head-wave) branches are handled explicitly by Snell’s law at layer interfaces.

The module is a core component of ShakeLab’s Seismicity package, used to compute expected travel times or source–receiver distances for both real-time and synthetic scenarios.

Theory

In a 1D layered medium, the travel time \(T\) and horizontal offset \(X\) for a given ray parameter \(p\) are expressed as:

\[\begin{split}X(p) = \sum_i h_i \frac{p v_i}{\sqrt{1 - (p v_i)^2}} \\ T(p) = \sum_i h_i \frac{1}{v_i \sqrt{1 - (p v_i)^2}}\end{split}\]

where \(v_i\) is the wave velocity and \(h_i\) the vertical thickness of the i-th layer. The shooting method iteratively adjusts \(p\) to satisfy \(X(p) = X_{target}\) using bisection within the physically admissible range \(p \in [0, 1/v_{max}]\).

Head-wave branches are generated when Snell’s law is satisfied at an interface such that the critical angle is reached, producing a refracted phase that propagates horizontally in the lower layer and returns to the surface.

API Reference

LayeredModel1D

class shakelab.seismicity.traveltime.LayeredModel1D(thickness, vp, vs, rho, qp=None, qs=None)

Defines a horizontally layered velocity model.

Parameters:

  • thickness (list of float): Layer thicknesses in meters. The last element can be np.inf for a half-space.

  • vp, vs (list of float): P- and S-wave velocities (m/s).

  • rho (list of float): Density (kg/m³).

  • qp, qs (optional): Quality factors (dimensionless).

Attributes:

  • n_layers: number of layers.

  • total_thickness(): total model depth.

  • find_layer(z): returns index of layer at depth z.

Main solver functions

shakelab.seismicity.traveltime.solve_traveltime(model, offsets, source_depth, receiver_depths=None, phases=('P',), arrivals='all', tol=1e-6, max_iter=60, require_convergence=True, depth_eps=0.5, verbose=False)

Computes travel times for given source–receiver geometry.

Parameters:

  • model: a LayeredModel1D instance.

  • offsets: array-like, receiver horizontal distances (m).

  • source_depth: float, source depth (m).

  • receiver_depths: optional array, one per offset (m).

  • phases: list or tuple of {“P”, “S”}.

  • arrivals: “direct” or “all” (includes head-waves).

  • tol: numerical tolerance for offset matching.

  • max_iter: maximum bisection iterations.

  • require_convergence: discard non-converged solutions if True.

  • depth_eps: tolerance to detect equal source/receiver depth.

  • verbose: print iteration info.

Returns:

Dictionary keyed by phase (“P”, “S”), each containing an array of travel times (s) corresponding to offsets.

shakelab.seismicity.traveltime.solve_with_heads(model, offset, source_depth, receiver_depth, phase='P', tol=1e-6, max_iter=60)

Computes both direct and head-wave branches for a single offset and phase, returning the physically valid (minimum-time) solution.

Secondary utilities

shakelab.seismicity.traveltime.compute_branches(model, offset, source_depth, receiver_depth=0.0, phase='P', include_direct=True, include_heads=True, tol=1e-6, max_iter=60)

Returns all geometrically valid branches (direct and/or head) as a list of dictionaries with keys:

{"t", "x", "kind", "p", "converged", "dx"}

shakelab.seismicity.traveltime.tt_direct(model, offset, source_depth, receiver_depth, phase='P', tol=1e-6, max_iter=60)

Computes only the direct branch, using iterative bisection.

shakelab.seismicity.traveltime.tt_direct_surface(model, offset, source_depth, phase='P', tol=1e-6, max_iter=60)

Direct arrival from depth source_depth to a surface receiver.

shakelab.seismicity.traveltime.tt_direct_curve(model, offsets, source_depth, receiver_depth, phase='P', tol=1e-6, max_iter=60)

Returns the direct traveltime curve for multiple offsets.

shakelab.seismicity.traveltime.curve_with_heads(model, offsets, source_depth, receiver_depth=0.0, phase='P', tol=1e-6, max_iter=60)

Computes a continuous curve of traveltimes including both direct and head-wave segments.

Usage examples

Example 1 – Source at 10 km, receivers at surface

from shakelab.seismicity.traveltime import solve_traveltime, LayeredModel1D
import numpy as np

model = LayeredModel1D(
    thickness=[2000, 3000, 5000, np.inf],
    vp=[3000, 4000, 6000, 8000],
    vs=[1500, 2300, 3400, 4500],
    rho=[2100, 2300, 2500, 2700]
)

offsets = np.linspace(0, 200e3, 21)
results = solve_traveltime(model, offsets, source_depth=10000.0,
                           phases=("P", "S"), arrivals="direct")

print(results["P"])

Example 2 – Source and receivers at 10 km, including head waves

results = solve_traveltime(model, offsets, source_depth=10000.0,
                           receiver_depths=np.ones_like(offsets)*10000.0,
                           phases=("P",), arrivals="all")

import matplotlib.pyplot as plt
plt.plot(offsets/1000, results["P"], label="P-wave (direct + head)")
plt.xlabel("Offset [km]")
plt.ylabel("Travel time [s]")
plt.legend()
plt.show()

Notes and limitations

  • Valid for horizontally layered media only (1D).

  • No reflection or multiple refraction phases.

  • For high-contrast interfaces, the solver may fail to converge if the offset exceeds the critical limit for a direct branch; in such cases, use arrivals="all".

  • Numerical stability depends on tol and max_iter; typical values (1e-6, 60) are sufficient for most engineering models.

References

  • Cervený, V. (2001). Seismic Ray Theory. Cambridge University Press.

  • Aki, K., & Richards, P. G. (2002). Quantitative Seismology. University Science Books.