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.

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. Convergence is typically achieved within 2–3 iterations.

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\).

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.

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):1648–1661, 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.