Methodology#

This chapter presents the theoretical foundations of the ray tracing algorithm implemented in LayTracer.

Two-point ray tracing problem#

Given a horizontally layered velocity model with \(N\) constant-velocity layers and two points—a source S at \((x_s, z_s)\) and a receiver R at \((x_r, z_r)\)—the task is to find the ray path connecting them that satisfies Snell’s law at every interface:

\[p = \frac{\sin \theta_k}{v_k} = \text{const}\]

where \(p\) is the horizontal slowness (ray parameter) and \(\theta_k\) is the ray angle from vertical in layer \(k\).

Phase and polarization convention#

For each source-receiver pair, LayTracer uses a local ray coordinate system. The radial axis lies in the horizontal source-receiver direction, the vertical axis is normal to the horizontal interfaces, and the transverse axis is perpendicular to this radial-vertical plane.

The scalar ray-tracing equations are phase agnostic once the appropriate layer velocity is selected: \(v_P\) for P waves and \(v_S\) for both SV and SH waves. Thus SV and SH share the same ray path, ray parameter, travel time, geometrical spreading, and \(t^*\) in isotropic layered media. They differ only in polarization and interface coefficients: SV belongs to the coupled P-SV system, whereas SH is the transverse, decoupled shear mode.

Prescribed phase itineraries#

LayTracer solves a branch selected by the caller. A RayItinerary is an ordered sequence of interface interactions. reflect reverses vertical propagation direction and transmit retains it; each interaction also specifies the phase on the outgoing leg. Every monotonic segment is split into atomic physical-layer legs, so a layer crossed more than once contributes more than once to the offset and traveltime while retaining one physical model-parameter index.

The legacy reflection and refraction tuple arguments compile to the same representation. In that legacy spelling, refraction denotes a prescribed transmission or mode conversion. It does not enumerate a critically refracted head wave. The solver does not perform branch discovery, automatic phase search, or global first-arrival selection.

Dimensionless ray parameter#

Fang and Chen [2019] introduce a dimensionless parameter

\[q = \sqrt{\frac{p^2}{1/v_{\max}^2 - p^2}} = \frac{p\,v_{\max}}{\sqrt{1 - p^2\,v_{\max}^2}}\]

where \(v_{\max} = \max_k v_k\). The inverse relation is

\[p = \frac{q}{v_{\max}\,\sqrt{1 + q^2}}\]

Key advantages:

  • \(q \in [0, \infty)\) maps the full range of valid take-off angles

  • The offset function \(X(q)\) is well-behaved (monotonically increasing, smooth, avoids singularities near the critical angle)

Offset equation#

The total horizontal distance (offset) \(X\) traversed by a ray through \(N\) layers is (Eq. 3 of Fang and Chen [2019]):

\[X(q) = \sum_{k=1}^{N} \frac{q\,\lambda_k\,h_k} {\sqrt{1 + (1 - \lambda_k^2)\,q^2}}\]

where \(\lambda_k = v_k / v_{\max}\) and \(h_k\) is the thickness of the \(k\)-th traversed layer.

First derivative#

\[\frac{\mathrm{d}X}{\mathrm{d}q} = \sum_{k=1}^{N} \frac{\lambda_k\,h_k} {\bigl[1 + (1 - \lambda_k^2)\,q^2\bigr]^{3/2}}\]

Second derivative#

\[\frac{\mathrm{d}^2X}{\mathrm{d}q^2} = -3\,q\,\sum_{k=1}^{N} \frac{(1 - \lambda_k^2)\,\lambda_k\,h_k} {\bigl[1 + (1 - \lambda_k^2)\,q^2\bigr]^{5/2}}\]

Asymptotic initial estimate#

Two linear asymptotes of \(X(q)\) provide an efficient initial guess (see “Initial estimate of q” in Fang and Chen [2019]):

Near-field (\(q \to 0\)):

\[X \approx m_0\,q, \qquad m_0 = \sum_k \lambda_k\,h_k\]

Far-field (\(q \to \infty\)):

\[X \approx m_\infty\,q + b_\infty\]

where

\[m_\infty = \sum_{k:\,\lambda_k=1} h_k, \qquad b_\infty = \sum_{k:\,\lambda_k<1} \frac{\lambda_k\,h_k}{\sqrt{1 - \lambda_k^2}}\]

The initial estimate \(q_0\) is chosen from the appropriate asymptote based on whether the target offset \(X_R\) falls in the near-field or far-field regime.

Quadratic Newton iteration#

The two-point problem \(X(q) = X_R\) is solved by second-order Newton iteration (see Fang and Chen [2019]). At each step, \(X(q)\) is expanded to second order about the current iterate \(q_i\):

\[\tfrac{1}{2}\,X''(q_i)\,\Delta q^2 + X'(q_i)\,\Delta q + \bigl[X(q_i) - X_R\bigr] = 0\]

This quadratic equation in \(\Delta q\) yields two roots. The root minimising \(|X(q_i + \Delta q) - X_R|\) is selected. The iteration is usually rapid, but convergence is not inferred from a fixed iteration count; the independently recomputed endpoint residual below is the acceptance test.

Safeguarding and numerical diagnostics#

The production solve accepts a quadratic-iteration result only after independently recomputing \(X(q)-X_R\). If the requested offset tolerance is not met, monotonicity of \(X(q)\) is used to construct a bracket on \(q\in[0,\infty)\) and Brent’s method supplies a safeguarded root. The returned SolveDiagnostics record identifies the accepted route (q_newton or q_brentq), primary and fallback iteration counts, initial and final \(q\) and \(p\), signed and absolute endpoint residuals, the criticality margin \(1-pv_{\max}\), and the dimensionless logarithmic slope

\[\kappa_q = \left|\frac{q}{X}\frac{\mathrm{d}X}{\mathrm{d}q}\right|.\]

Diagnostics are opt-in and do not allocate per-ray diagnostic objects for ordinary calls. Invalid itineraries and missing properties required by a requested output raise an explicit exception; they are not reported as a successful trace.

Travel time#

Once the ray parameter \(p\) is determined, the travel time through each layer is

\[\Delta t_k = \frac{h_k}{v_k^2\,\eta_k}\]

where \(\eta_k = \sqrt{1/v_k^2 - p^2}\) is the vertical slowness in layer \(k\). The total travel time is \(t = \sum_k \Delta t_k\).

Fixed-topology analytic sensitivities#

Ray tracing predicts an arrival time; many practical tasks also need to know why that prediction changes. For example, an event locator needs the Jacobian of travel time with respect to source position for a Gauss–Newton update, and a network-design or uncertainty study needs the Jacobian with respect to layer velocities and interface depths. Those Jacobians turn small model perturbations into predicted arrival-time perturbations,

\[\delta\mathbf{t} \approx \mathbf{J}_m\,\delta\mathbf{m} + \mathbf{J}_x\,\delta\mathbf{x}.\]

They also permit first-order propagation of velocity-model uncertainty into an arrival-time covariance, for example \(\mathbf{C}_t = \mathbf{J}_m\mathbf{C}_m\mathbf{J}_m^\mathsf{T} + \mathbf{C}_{\mathrm{pick}}\). Without analytic derivatives, each column of these matrices normally requires several perturbed ray traces; the resulting finite differences are slower and can be noisy when the nonlinear ray solve is poorly conditioned.

LayTracer therefore differentiates the already solved ray. "Fixed topology" means that, over the small perturbation being considered, the ray retains the same ordered sequence of layer legs, phase conversions, and reflections. It is the appropriate local approximation for a specified arrival (for example, a PP reflection from a selected interface) rather than an automatic comparison or switch between competing branches. LayTracer returns the derivatives; selecting phases, defining an objective function, and performing the inversion remain the caller’s responsibility.

The ray parameter must change when a velocity, interface, source, or receiver is perturbed so that the ray still joins the same endpoints. For atomic leg \(j\), define

\[c_j=\sqrt{1-p^2v_j^2},\qquad X_p=\sum_j\frac{h_jv_j}{c_j^3}.\]

Here \(p\) is the physical horizontal slowness, \(h_j\) and \(v_j\) are the leg thickness and phase velocity, and \(X_p\) measures the local change in offset with ray parameter. If \(X_m\) and \(T_m\) denote the partial offset and travel-time changes at fixed \(p\), implicit differentiation of the endpoint-offset constraint and Fermat stationarity give

\[\frac{\mathrm{d}p}{\mathrm{d}m}=-\frac{X_m}{X_p},\qquad \left.\frac{\mathrm{d}T}{\mathrm{d}m}\right|_X=T_m-pX_m.\]

The first expression is the small refraction-angle adjustment required to keep the endpoints fixed. The second is the corresponding travel-time derivative: the direct effect of changing the parameter plus the effect of re-bending the ray, without tracing a separate perturbed ray.

After this fixed-offset correction, the per-leg velocity and vertical-thickness building blocks reduce to

\[\frac{\mathrm{d}T}{\mathrm{d}v_j} =-\frac{h_j}{v_j^2c_j},\qquad \frac{\mathrm{d}T}{\mathrm{d}h_j}=\frac{c_j}{v_j}.\]

The implementation aggregates repeated leg contributions into sparse derivatives for \(V_P\) and \(V_S\) in each physical layer. Thus a multiple that crosses one layer several times contributes to that layer’s single model parameter rather than creating duplicate columns. Signed incidence coefficients map thickness derivatives to interface-depth derivatives. Endpoint derivatives are returned for source and receiver x, y, and z coordinates, together with derivatives of the physical ray parameter \(p\). Calculating the derivatives in \(p\) avoids introducing an artificial discontinuity when the identity of the fastest traversed leg changes.

These are deliberately local derivatives. They are valid only while the endpoint layer membership, prescribed phase itinerary, and propagating branch remain unchanged. A perturbation that moves an endpoint across an interface, creates or removes a leg, changes the selected arrival, or reaches a critical or evanescent path is outside this linearization; such sensitivity records are marked invalid rather than being silently extrapolated.

Topology-safe traveltime approximation#

The endpoint derivatives can replace repeated exact tracing in a dense local source-receiver table. For an exact anchor ray joining \(\mathbf{s}_a\) and \(\mathbf{r}_a\), LayTracer evaluates

\[T(\mathbf{s},\mathbf{r}) \approx T_a + \nabla_{\mathbf{s}}T_a\mathbin{\cdot} (\mathbf{s}-\mathbf{s}_a) + \nabla_{\mathbf{r}}T_a\mathbin{\cdot} (\mathbf{r}-\mathbf{r}_a).\]

This is the default first-order interpolation of one selected arrival, not a new ray solve. The velocity model and any explicit RayItinerary are fixed for the lifetime of a TravelTimeApproximator. Changing either requires refitting.

Second-order endpoint prediction is available with predict(order=2). Let \(R\) be the anchor’s horizontal offset, \(p\) its physical ray parameter, \(X_p=\partial X/\partial p\), and \(p'\) the directional change in ray parameter obtained from the endpoint derivatives. If \(\delta\mathbf{r}_\perp\) is the component of the relative horizontal endpoint displacement perpendicular to the anchor azimuth, the directional curvature is

\[T'' = X_p(p')^2 + \frac{p}{R}\|\delta\mathbf{r}_\perp\|^2.\]

The zero-offset limit replaces the second term by \(\|\delta\mathbf{r}\|^2/X_p\). LayTracer adds \(T''/2\) to the linear prediction. Direct endpoints in the same physical layer use the exact homogeneous straight-path traveltime instead, including horizontal and coincident limits. Second order remains a local, fixed-topology approximation; it does not make extrapolation across a layer boundary or ray branch valid.

Anchor coverage#

select_anchors uses deterministic farthest-point coverage of an arbitrary three-dimensional point cloud. It starts from the lexicographically smallest point in each group and adds the farthest point until every target lies within the requested Euclidean distance of an anchor. During approximator fitting, source and receiver anchors are initially selected independently within their physical endpoint layers. The sets are then refined until every reachable fitting pair has a topology-compatible anchor pair inside both distance limits. A limit of None keeps every supplied endpoint of that type as an exact anchor.

The final Cartesian product of source and receiver anchors is traced once with sensitivities and diagnostics. When one endpoint type is retained exactly, coverage is refined independently for each exact endpoint and topology; this avoids repeatedly scanning the full anchor-pair product. For each prediction, candidate anchor pairs are ranked using the source and receiver distances normalized by their respective limits. The nearest compatible valid sensitivity supplies the first- or second-order prediction above.

Topology and validity#

Compatibility is defined by the ordered physical path legs. Each leg records its model-layer index, phase, and vertical direction. Consequently, the approximator does not use an anchor that would move an endpoint across an interface, reverse the vertical direction of a direct ray, change a prescribed itinerary, or switch phase sequence.

Unsupported predictions are not silently retraced. They are returned as NaN with a validity mask and one of four reason codes:

  • unreachable_topology: the requested fixed itinerary cannot connect the endpoints;

  • no_topology_match: no fitted anchor ray represents the requested branch;

  • outside_anchor_distance: compatible anchors exist but are beyond a distance limit;

  • invalid_anchor_sensitivity: compatible nearby anchor rays have unusable derivatives.

Anchor distance controls locality but is not an intrinsic error estimate. Approximation error depends on model contrasts, geometry, and distance from the anchor, and can be largest near rapidly varying configurations such as very short source-receiver offsets. The second-order mode reduces this local curvature error but does not provide an error bound. Applications should select distance limits using representative validation retraces. LayTracer deliberately does not perform those retraces, adapt anchors to an error tolerance, or fall back to exact tracing automatically. The visual accuracy–cost experiments in 06. Visualizing and using sensitivities illustrate this validation workflow.

Attenuation operator \(t^*\)#

The attenuation operator (Aki and Richards [2002], Ch. 5) measures the cumulative dissipative loss of wave amplitude along the ray path. Since the spatial path length in layer \(k\) is \(\Delta s_k = v_k \Delta t_k\), the spatial integral representing intrinsic absorption corresponds exactly to:

\[t^* = \sum_{k=1}^{N} \frac{\Delta t_k}{Q_k}\]

where \(Q_k\) is the quality factor in layer \(k\). This gives the spectral decay \(\exp(-\pi f\,t^*)\) at frequency \(f\).

For a vertical ray: \(t^* = \sum h_k / (v_k\,Q_k)\).

For uniform \(Q\): \(t^* = t / Q\).

Geometrical spreading#

In a 1-D layered medium with cylindrical symmetry (3-D point source), the classical geometrical spreading factor \(L\) relates the solid angle of the ray tube at the source to its cross-sectional area at the receiver (Červený [2001], Aki and Richards [2002]):

\[L = \sqrt{\frac{X \cdot \cos\theta_s \cdot \cos\theta_r}{p} \left| \frac{\partial X}{\partial p} \right|}\]

where \(\theta_s, \theta_r\) are the ray angles at source and receiver. Equation above defines the relative geometrical spreading (see Červený [2001], Eq. 4.10.22), which measures the ray-tube geometrical divergence strictly from the ray curvature, without the source-point velocity multiplier \(1/v_r\).

The derivative \(\partial X / \partial p\) is computed analytically via the chain rule:

\[\frac{\partial X}{\partial p} = \frac{\mathrm{d}X}{\mathrm{d}q} \cdot \frac{\mathrm{d}q}{\mathrm{d}p}, \qquad \frac{\mathrm{d}q}{\mathrm{d}p} = \frac{v_{\max}}{(1 - p^2\,v_{\max}^2)^{3/2}}\]

Reflection and transmission coefficients#

In layered media, wave amplitudes are modified at every crossed interface. LayTracer computes interface coefficients using the full angle-dependent P-SV Zoeppritz formulation for P and SV, and the decoupled SH-SH formulation for SH. Two variants are available via the transcoef_method parameter:

  • "standard" — displacement-amplitude coefficients (Zoeppritz); this is the default.

  • "normalized" — energy-flux-normalized coefficients following Červený [2001], Eq. 5.3.10.

Angle-dependent P-SV formulation (welded solid-solid interface)#

For horizontal slowness (ray parameter) \(p\), the vertical slownesses are

\[\eta_{\alpha i} = \sqrt{\frac{1}{v_{Pi}^2} - p^2}, \qquad \eta_{\beta i} = \sqrt{\frac{1}{v_{Si}^2} - p^2},\]

and the auxiliary quantities

\[a = \rho_2\!\left(1-2v_{S2}^2p^2\right)-\rho_1\!\left(1-2v_{S1}^2p^2\right),\]
\[b = \rho_2\!\left(1-2v_{S2}^2p^2\right)+2\rho_1v_{S1}^2p^2,\]
\[c = \rho_1\!\left(1-2v_{S1}^2p^2\right)+2\rho_2v_{S2}^2p^2,\]
\[d = 2\left(\rho_2v_{S2}^2-\rho_1v_{S1}^2\right).\]

Define the cosine-dependent intermediate terms

\[E=b\,\eta_{\alpha1}+c\,\eta_{\alpha2}, \quad F=b\,\eta_{\beta1}+c\,\eta_{\beta2}, \quad G=a-d\,\eta_{\alpha1}\eta_{\beta2}, \quad H=a-d\,\eta_{\alpha2}\eta_{\beta1},\]

and the system determinant

\[D = EF + GH\,p^2.\]

The complete \(4\times 4\) scattering matrix (Aki and Richards [2002], Eqs. 5.38-5.40) is computed by LayTracer. The eight independent P-SV coefficients are listed below.

Incident P-wave — reflection and transmission:

\[R_{PP} = \frac{\left(b\,\eta_{\alpha1}-c\,\eta_{\alpha2}\right)F - \left(a+d\,\eta_{\alpha1}\eta_{\beta2}\right)H\,p^2}{D},\]
\[R_{PS} = -\frac{2\,\eta_{\alpha1}\left(ab+cd\,\eta_{\alpha2}\eta_{\beta2}\right) p\,(v_{P1}/v_{S1})}{D},\]
\[T_{PP} = \frac{2\rho_1\,\eta_{\alpha1}\,F\,(v_{P1}/v_{P2})}{D},\]
\[T_{PS} = \frac{2\rho_1\,\eta_{\alpha1}\,H\,p\,(v_{P1}/v_{S2})}{D}.\]

Incident SV-wave — reflection and transmission:

\[R_{SP} = -\frac{2\,\eta_{\beta1}\left(ab+cd\,\eta_{\alpha2}\eta_{\beta2}\right) p\,(v_{S1}/v_{P1})}{D},\]
\[R_{SS} = -\frac{\left(b\,\eta_{\beta1}-c\,\eta_{\beta2}\right)E - \left(a+d\,\eta_{\alpha2}\eta_{\beta1}\right)G\,p^2}{D},\]
\[T_{SP} = -\frac{2\rho_1\,\eta_{\beta1}\,G\,p\,(v_{S1}/v_{P2})}{D},\]
\[T_{SS} = \frac{2\rho_1\,\eta_{\beta1}\,E\,(v_{S1}/v_{S2})}{D}.\]

For references and details on the derivation of these formulas, see Lay and Wallace [1995] (Table 3.1, note the sign error in the second term of \(b\)) and Aki and Richards [2002] (Equations 5.38-5.40).

For post-critical incidence the coefficients may become complex; for amplitude modelling the software uses \(|T_l|\).

Angle-dependent SH formulation (welded solid-solid interface)#

In isotropic layered media, SH is decoupled from the P-SV system (Aki and Richards [2002], Section 5.2.3; Lay and Wallace [1995], Eq. 3.109). With vertical slowness

\[\eta_i = \sqrt{\frac{1}{v_{Si}^2} - p^2},\]

the SH displacement coefficients at a welded interface are

\[R_{SHSH} = \frac{\rho_1 v_{S1}^2 \eta_1 - \rho_2 v_{S2}^2 \eta_2} {\rho_1 v_{S1}^2 \eta_1 + \rho_2 v_{S2}^2 \eta_2},\]
\[T_{SHSH} = \frac{2\rho_1 v_{S1}^2 \eta_1} {\rho_1 v_{S1}^2 \eta_1 + \rho_2 v_{S2}^2 \eta_2}.\]

At normal incidence these reduce to the shear-impedance forms with \(Z_S = \rho v_S\). Since SH does not couple to P or SV in this model class, LayTracer rejects P/SV-to-SH mode conversions.

Energy-flux-normalized coefficients#

The standard (displacement-amplitude) Zoeppritz coefficients \(\bar{R}_{mn}\) do not conserve energy flux across an interface. For applications where energy conservation is essential (e.g. seismic-moment estimation in layered media), Červený [2001] (Eq. 5.3.10) defines normalized reflection/transmission coefficients:

\[\mathcal{R}_{mn} = \bar{R}_{mn}\, \sqrt{\frac{v_{\mathrm{out}}\,\rho_{\mathrm{out}}\, \cos\theta_{\mathrm{out}}} {v_{\mathrm{in}}\,\rho_{\mathrm{in}}\, \cos\theta_{\mathrm{in}}}},\]

where \(\cos\theta = \sqrt{1 - v^2 p^2}\) and subscripts “in”/”out” refer to the incident and scattered wave, respectively (for reflections, “out” uses properties of the same side of the interface). For subcritical incidence at a lossless interface, the normalized coefficients ensure conservation of energy flux, with reflected and transmitted energy fractions summing to unity.

Critical angles#

A critical angle occurs when a reflected or transmitted wave mode in a faster medium becomes evanescent, meaning it no longer propagates away from the interface but decays normal to it. For any incident mode \(m\) and outgoing mode \(n\), the critical angle is

\[\theta_c^{m \to n} = \arcsin\!\left(\frac{v_{m,\mathrm{in}}}{v_{n,\mathrm{out}}}\right),\]

provided \(v_{n,\mathrm{out}} > v_{m,\mathrm{in}}\). For P-SV scattering, \(v_m\) and \(v_n\) are selected from the P or SV wave speeds on the incident and outgoing sides. For SH-SH transmission, the same rule applies with \(v_S\) on both sides:

\[\theta_c^{SH \to SH} = \arcsin\!\left(\frac{v_{S,\mathrm{in}}}{v_{S,\mathrm{out}}}\right).\]

Beyond the critical angle the corresponding vertical slowness becomes imaginary, the coefficient becomes complex, and the outgoing wave mode carries no propagating normal energy flux.

Brewster angles#

By analogy with optics, a Brewster angle is an incidence angle at which a reflection or transmission coefficient passes through zero or a deep minimum. In electromagnetic theory, Brewster’s angle is the incidence angle at which \(R_p = 0\) for p-polarised light at a dielectric interface. In elastodynamics, the same phenomenon occurs: certain combinations of elastic parameters produce incidence angles where one of the P-SV scattering coefficients vanishes.

Unlike critical angles, which depend only on velocity ratios, P-SV Brewster angles depend on all six elastic parameters (\(v_{P1}\), \(v_{S1}\), \(\rho_1\), \(v_{P2}\), \(v_{S2}\), \(\rho_2\)). Physically, they arise from destructive interference between the P and SV displacement potentials at the welded interface: the two potential contributions to a particular scattered mode cancel exactly, driving that coefficient to zero.

For example, the reflected P coefficient \(R_{PP}\) may vanish at an angle well below the critical angle. At this Brewster angle the incident energy is partitioned entirely into the transmitted P-wave and the mode-converted waves, with no same-mode reflection.

For SH, the analogous zero has a simpler interpretation because SH is decoupled from the P-SV potentials. The reflected SH coefficient vanishes when the oblique shear impedances match across the interface:

\[\rho_1 v_{S1}^2 \eta_1 = \rho_2 v_{S2}^2 \eta_2,\]

or equivalently

\[\rho_1 v_{S1}\cos\theta_1 = \rho_2 v_{S2}\cos\theta_2.\]

This SH null depends only on shear velocities, densities, and incidence angle.

LayTracer provides the function find_brewster_angles which numerically detects these minima in the computed coefficient curves by searching for local minima of \(|C(\theta)|\) whose value falls below a user-specified threshold.

Extension to 3-D layered media#

The formulae above are written in a local incidence coordinate system. The ray geometry is two-dimensional because the velocity model varies only with depth, so the horizontal slowness direction is fixed along the path. P and SV particle motions are in the radial-vertical incidence plane; SH motion is polarized along the transverse axis perpendicular to that plane.

The 2-D reduction is therefore a reduction of geometry and material dependence, not a statement that every displacement vector lies in the incidence plane. Because horizontally layered media possess cylindrical symmetry about the vertical axis through the source, any source-receiver pair in three-dimensional space can be reduced to an equivalent 2-D problem without loss of generality (Červený [2001], Ch. 3; Aki and Richards [2002], Ch. 4). The local incidence plane is obtained by rotating the 2-D problem into the source-receiver azimuth.

Coordinate projection#

Let the source position be \(\mathbf{s} = (x_s, y_s, z_s)\) and the receiver position \(\mathbf{r} = (x_r, y_r, z_r)\). The epicentral distance (horizontal distance) is

\[\Delta = \sqrt{(x_r - x_s)^2 + (y_r - y_s)^2}\]

and the unit direction vector from source to receiver in the horizontal plane is

\[\begin{split}\hat{\mathbf{u}} = \frac{1}{\Delta} \begin{pmatrix} x_r - x_s \\ y_r - y_s \end{pmatrix}, \qquad \Delta > 0\end{split}\]

(for vanishing \(\Delta\) the azimuth is arbitrary, and we default to \(\hat{\mathbf{u}} = (1,\,0)\)).

The 2-D ray tracing problem is then solved exactly as described in the preceding sections, with the epicentral distance \(\Delta\) playing the role of the target offset \(X_R\) and the depth coordinates \(z_s, z_r\) determining the traversed layers.

Back-projection to 3-D#

Once the 2-D ray path \(\{(x_k^{(2\mathrm{D})},\, z_k)\}_{k=0}^M\) has been computed, it is mapped back to 3-D Cartesian coordinates via

\[\begin{split}\begin{aligned} X_k &= x_s + x_k^{(2\mathrm{D})}\;\hat{u}_x, \\ Y_k &= y_s + x_k^{(2\mathrm{D})}\;\hat{u}_y, \\ Z_k &= z_k. \end{aligned}\end{split}\]

This simply sweeps the 2-D ray along the source-receiver azimuth.

Validity of amplitude attributes#

Because the medium properties depend only on depth, every quantity computed from the 2-D ray remains valid in 3-D:

  • Ray parameter \(p\) — determined solely by Snell’s law across horizontal interfaces.

  • Travel time — sum of vertical-slowness contributions \(\Delta t_k\), unchanged by azimuth.

  • Attenuation operator \(t^*\) — depends only on \(\Delta t_k\) and \(Q_k\).

  • Geometrical spreading — the formula incorporates the epicentral distance \(X\) to capture the 3-D cylindrical divergence of the ray-tube out of the incidence plane (Červený [2001], §4.10).

  • Interface coefficients - scalar reflection/transmission coefficients depend on ray parameter and layer elastic properties, not on absolute azimuth. P and SV use the coupled P-SV coefficients; SH uses the decoupled SH-SH coefficients. The SH polarization direction itself rotates with the source-receiver azimuth.

Thus, the layered-media solver operates in the 2-D incidence plane; the 3-D ray geometry and scalar amplitude attributes are recovered by rotation. Full vector receiver components or source radiation patterns require an additional polarization/projection step. The legacy trans_product is a real, nonnegative product of coefficient magnitudes. Request complex_coefficient_product separately when polarity and post-critical phase are needed; for the same compatible path, trans_product == abs(complex_coefficient_product). Neither coefficient product nor the spreading factor is a complete synthetic amplitude because source radiation, free-surface response, polarization projection, and finite-frequency effects are outside this operator. Density and quality factor columns are validated only when an output that depends on them is requested.

References#

[1] (1,2,3,4)

X. Fang and X. Chen. A fast and robust two-point ray tracing method in layered media with constant or linearly varying layer velocity. Geophysical Prospecting, 67(7):1811–1824, 2019. doi:10.1111/1365-2478.12799.

[2] (1,2,3,4,5,6)

K. Aki and P. G. Richards. Quantitative Seismology. University Science Books, 2nd edition, 2002.

[3] (1,2,3,4,5,6)

V. Červený. Seismic Ray Theory. Cambridge University Press, 2001. doi:10.1017/CBO9780511529399.

[4] (1,2)

T. Lay and T. C. Wallace. Modern Global Seismology. Academic Press, 1995.