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

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

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 the transmitted wave in a faster medium becomes evanescent. For an incident P-wave crossing into a layer where \(v_{P2} > v_{P1}\), the P-critical angle is

\[\theta_c^{P} = \arcsin\!\left(\frac{v_{P1}}{v_{P2}}\right).\]

Similarly, when \(v_{S2} > v_{P1}\) an S-to-P critical angle exists:

\[\theta_c^{S} = \arcsin\!\left(\frac{v_{P1}}{v_{S2}}\right).\]

For an incident SV-wave the same logic applies with \(v_{S1}\) in the numerator. Beyond the critical angle the corresponding vertical slowness becomes imaginary, the coefficient becomes complex, and total reflection occurs for that mode.

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, 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.

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#

All the formulae above are derived in a 2-D vertical plane containing both the source and the receiver. 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).

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

  • Transmission coefficients — depend on ray parameter and layer impedances, not on azimuth.

Thus, the layered-media solver need only operate in the 2-D ray plane; the full 3-D solution is recovered by geometry alone.

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)

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]

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