Elastic VTI (anisotropic) wave equation

Note

Stability: experimental – the label records that this propagator has never appeared in a release, so 0.1.0 is the first time real use meets its design. It is not a lesser promise: its interface is frozen exactly as the other eight propagators’ are, calling it emits no warning about its maturity, and every other Deepwave promise applies to it unchanged. See Stability policy.

The deepwave.elastic_vti() propagator solves the elastic wave equation in a vertically transverse isotropic (VTI) medium: an anisotropic medium with a single axis of rotational symmetry, taken to be vertical (the slowest spatial dimension – \(z\) in 3D, \(y\) in 2D). VTI is the most common form of seismic anisotropy, arising from fine horizontal layering and from aligned horizontal fractures or shale fabric. It is Deepwave’s first anisotropic propagator.

It shares the entire numerical machinery of the isotropic elastic propagator – the velocity–stress formulation on the same staggered grid, the C-PML (and optional multiaxial M-PML) absorbing boundaries, the improved-vacuum (Zeng et al., 2012) free surface, the staggering of sources and receivers, and the storage/adjoint skeleton. Only the constitutive (stress–strain) relation differs. Reading the elastic page first is recommended; this page describes only what anisotropy changes.

See Material models and physical scope for the material domain, dimensions, and limits of the interior energy argument, and Units and source conventions for source and receiver normalisation.

Minimal example

deepwave.elastic_vti() takes the VTI stiffnesses and a buoyancy model, a grid spacing, a time step, and a source. This example builds the stiffnesses from Thomsen parameters with deepwave.anisotropy.thomsen_to_cij() (see the Thomsen parameters section below):

import torch
import deepwave

vp0 = 1500 * torch.ones(40, 40)
vs0 = 800 * torch.ones(40, 40)
rho = 2200 * torch.ones(40, 40)
epsilon = 0.2 * torch.ones(40, 40)
delta = 0.1 * torch.ones(40, 40)
gamma = 0.15 * torch.ones(40, 40)
c11, c13, c33, c55, c66, buoyancy = deepwave.anisotropy.thomsen_to_cij(vp0, vs0, rho, epsilon, delta, gamma)

out = deepwave.elastic_vti(
    c11=c11, c13=c13, c33=c33, c55=c55, buoyancy=buoyancy,
    grid_spacing=5, dt=0.001,
    source_amplitudes_y=deepwave.wavelets.ricker(25, 200, 0.001, 0.05).reshape(1, 1, -1),
    source_locations_y=torch.tensor([[[5, 20]]]),
    receiver_locations_y=torch.tensor([[[35, 20]]]),
    pml_width=10,
    pml_freq=25,
)
out.receiver_amplitudes_y  # shape [1, 1, 200]

Alongside the receiver data, out also carries the final wavefields plus one opaque auxiliary_state holding the absorbing boundary’s memory variables. Passing the whole output back in as initial_state continues the propagation in a later call (see Usage).

Constitutive relation

For a general linear-elastic medium the stress is related to the strain by the fourth-order stiffness tensor \(c_{ijkl}\),

(1)\[\sigma_{ij} = c_{ijkl}\epsilon_{kl}, \qquad \epsilon_{kl} = \tfrac{1}{2}\left(\partial_k u_l + \partial_l u_k\right),\]

which in Voigt notation (\(11\!\to\!1\), \(22\!\to\!2\), \(33\!\to\!3\), \(23\!\to\!4\), \(13\!\to\!5\), \(12\!\to\!6\)) becomes a symmetric \(6\times 6\) matrix \(C_{IJ}\). For a VTI medium with the symmetry axis along direction \(3\) this matrix has five independent stiffnesses \(c_{11}, c_{13}, c_{33}, c_{55}, c_{66}\):

(2)\[\begin{split}C = \begin{pmatrix} c_{11} & c_{12} & c_{13} & 0 & 0 & 0\\ c_{12} & c_{11} & c_{13} & 0 & 0 & 0\\ c_{13} & c_{13} & c_{33} & 0 & 0 & 0\\ 0 & 0 & 0 & c_{55} & 0 & 0\\ 0 & 0 & 0 & 0 & c_{55} & 0\\ 0 & 0 & 0 & 0 & 0 & c_{66} \end{pmatrix}, \qquad c_{12} = c_{11} - 2c_{66}.\end{split}\]

The isotropic elastic medium is the special case \(c_{11}=c_{33}=\lambda + 2\mu\), \(c_{13}=c_{12}=\lambda\), \(c_{55}=c_{66}=\mu\), so deepwave.elastic_vti() reproduces deepwave.elastic() exactly in that limit (verified to machine precision in the test suite, for the forward field and for every gradient, in 2D and 3D).

3D. With the symmetry axis as \(z\) (Voigt 3), \(x\)=``1``, \(y\)=``2``, the stress-rate updates are

(3)\[\begin{split}\begin{aligned} \dot\sigma_{xx} &= c_{11}\partial_x v_x + c_{12}\partial_y v_y + c_{13}\partial_z v_z,\\ \dot\sigma_{yy} &= c_{12}\partial_x v_x + c_{11}\partial_y v_y + c_{13}\partial_z v_z,\\ \dot\sigma_{zz} &= c_{13}\partial_x v_x + c_{13}\partial_y v_y + c_{33}\partial_z v_z,\\ \dot\sigma_{yz} &= c_{55}\left(\partial_y v_z + \partial_z v_y\right),\\ \dot\sigma_{xz} &= c_{55}\left(\partial_x v_z + \partial_z v_x\right),\\ \dot\sigma_{xy} &= c_{66}\left(\partial_x v_y + \partial_y v_x\right). \end{aligned}\end{split}\]

2D (P–SV in the \((y,x)\) plane, axis = \(y\)). Here \(y\) is Voigt 3 and \(x\) is 1; \(c_{66}\) and \(\sigma_{xy}\)-plane terms drop out and \(c_{66}\) must be None:

(4)\[\begin{split}\begin{aligned} \dot\sigma_{xx} &= c_{11}\partial_x v_x + c_{13}\partial_y v_y,\\ \dot\sigma_{yy} &= c_{13}\partial_x v_x + c_{33}\partial_y v_y,\\ \dot\sigma_{xy} &= c_{55}\left(\partial_x v_y + \partial_y v_x\right). \end{aligned}\end{split}\]

There is no 1D VTI propagator: with a single spatial dimension only \(c_{33}\) survives, which is the isotropic 1D elastic equation, so 1D inputs raise an error (use deepwave.elastic()).

The velocity update is unchanged from the isotropic case, \(\dot v_i = b\,\partial_j\sigma_{ij}\) with \(b=1/\rho\) the buoyancy; anisotropy enters only through the stiffnesses above.

Units and locations

All stiffnesses \(c_{IJ}\) have units of pressure (Pa in SI), the same as the Lamé parameters they generalise, and are supplied at the normal-stress (integer) grid nodes. The buoyancy is \(1/\rho\) (units \(\mathrm{m^3\,kg^{-1}}\)). As in the isotropic propagator, the shear stiffnesses are internally averaged onto the staggered shear-stress nodes and the buoyancy onto the staggered velocity nodes.

Free-surface / improved-vacuum averaging

The improved-vacuum free surface treats a cell with buoyancy = 0 (and zero stiffnesses) as a vacuum, with the model parameters averaged onto the staggered nodes exactly as in deepwave.elastic(). Each VTI stiffness is averaged in the way appropriate to the elastic modulus it generalises:

  • \(c_{55}\) generalises the shear modulus \(\mu\) for the stresses whose planes contain the symmetry axis (\(\sigma_{xz}\), \(\sigma_{yz}\) in 3D, and the single shear stress in 2D). It is harmonically averaged (with vacuum masking) onto those shear-stress nodes exactly as \(\mu\) is.

  • \(c_{66}\) generalises \(\mu\) for the horizontal-plane shear stress \(\sigma_{xy}\) (3D only) and is harmonically averaged onto that node in the same way.

  • The normal-node stiffnesses \(c_{11}, c_{12}, c_{13}, c_{33}\) are used unaveraged at the nodes where the user supplies them, exactly as \(\lambda\) and \(\lambda+2\mu\) are in the isotropic propagator.

Thomsen parameters

Anisotropy is often described by Thomsen (1986) parameters \((v_{P0}, v_{S0}, \rho, \epsilon, \delta, \gamma)\), where \(v_{P0}\) and \(v_{S0}\) are the P- and S-wave speeds along the symmetry axis, \(\epsilon\) and \(\gamma\) measure the fractional difference between horizontal and vertical P and SH velocities, and \(\delta\) controls the near-vertical curvature of the qP wavefront (and hence the short-spread NMO velocity). The helper deepwave.anisotropy.thomsen_to_cij() converts these to stiffnesses using the exact (not the weak-anisotropy linearised) relations

(5)\[\begin{split}\begin{aligned} c_{33} &= \rho v_{P0}^2, & c_{55} &= \rho v_{S0}^2,\\ c_{11} &= c_{33}(1 + 2\epsilon), & c_{66} &= c_{55}(1 + 2\gamma),\\ c_{13} &= \sqrt{2\delta\, c_{33}(c_{33}-c_{55}) + (c_{33}-c_{55})^2} - c_{55}. \end{aligned}\end{split}\]

The \(c_{13}\) expression inverts Thomsen’s exact definition \(\delta = \left[(c_{13}+c_{55})^2 - (c_{33}-c_{55})^2\right] / \left[2 c_{33}(c_{33}-c_{55})\right]\), taking the positive root, which corresponds to \(c_{13}+c_{55} > 0\); this holds for essentially all real materials. The square-root argument is non-negative iff \(\delta \ge -(c_{33}-c_{55})/(2 c_{33})\); more negative \(\delta\) cannot be represented and is rejected with a clear error. The conversion is differentiable, so inversion may be performed directly in Thomsen-parameter space by chaining autograd through it (verified by finite-difference gradient checks in the test suite).

Alongside the five stiffnesses, the helper returns the buoyancy \(b = 1/\rho\), because that – not \(\rho\) – is what the propagators take, so its whole return value can be unpacked straight into deepwave.elastic_vti(). A vacuum cell (\(\rho = 0\)) is mapped to \(b = 0\) exactly, which is the vacuum sentinel the propagators expect, rather than to the infinity a plain 1 / rho would produce. deepwave.anisotropy.cij_to_thomsen() inverts the whole conversion, taking the buoyancy in the same place.

P-wave surface data alone cannot resolve \(\epsilon\) and \(\delta\) separately – time processing (NMO, DMO, time migration) is controlled only by the P-wave NMO velocity \(v_{nmo}\) and the anellipticity parameter \(\eta\) (Alkhalifah and Tsvankin, 1995). The helper deepwave.anisotropy.alkhalifah_to_thomsen() converts \((v_{P0}, v_{nmo}, \eta)\) to \((\epsilon, \delta)\) via the exact relations \(\delta = (v_{nmo}^2/v_{P0}^2 - 1)/2\) and \(\epsilon = \eta(1 + 2\delta) + \delta\), so a model estimated from NMO-velocity analysis can be chained straight into deepwave.anisotropy.thomsen_to_cij() (choosing an \(v_{S0}\), since shear information is not recoverable from P-wave-only data) to build the Cij stiffness model that deepwave.elastic_vti() (or deepwave.acoustic_vti(), with \(v_{S0}=0\)) consumes – also a differentiable chain, so FWI may be performed directly in \((v_{P0}, v_{nmo}, \eta)\) parameter space. deepwave.anisotropy.thomsen_to_alkhalifah() is its inverse, converting \((v_{P0}, \epsilon, \delta)\) back to \((v_{nmo}, \eta)\) via \(v_{nmo} = v_{P0}\sqrt{1 + 2\delta}\) and \(\eta = (\epsilon - \delta)/(1 + 2\delta)\), for comparing a Thomsen model with, or initialising it from, a conventional NMO-velocity model.

Phase velocities

Independent tests also propagate 3D SH, qP and qSV plane-wave initial states with unequal c55 and c66 and an oblique horizontal azimuth. Refining space and time converges to the continuum Christoffel solution. Separate normal-incidence P and SV tests cross a welded interface between two VTI materials and converge to the analytic reflection/transmission waveforms. These checks cover those modes and geometries; they do not establish an error bound for arbitrary oblique interfaces or staircase topography.

For a plane wave with unit propagation direction \(\mathbf{n}\) the phase velocities follow from the Christoffel matrix \(\Gamma_{ik} = c_{ijkl} n_j n_l\); the three eigenvalues \(\rho v^2\) give the qP, qSV and (in 3D) SH phase speeds, which depend on the angle to the symmetry axis. The propagator reproduces the analytic Christoffel qP and qSV phase velocities across angle (measured from band-limited plane waves in a homogeneous strongly-anisotropic medium) to within the grid-dispersion budget; this is the test that detects axis-convention or constitutive-coupling errors, which the isotropic limit (\(c_{11}=c_{33}\)) cannot.

Stability

Readable backgrounds are checked for finite values and positive represented elastic energy. Active cells require positive normal stiffnesses and nonnegative shear stiffnesses. Positive shear stiffnesses require a positive normal block; supported zero-shear fluid limits permit a semidefinite block. Vacuum requires zero buoyancy and zero stiffness. See Material models and physical scope for the exact 2D/3D conditions and tensor-transform preconditions.

When max_vel is omitted, setup derives a conservative speed from the constitutive tensor and staggered coefficient averaging, including heterogeneous contrast. A pointwise phase-speed maximum alone need not bound the discrete operator. An explicit speed must cover the resolved bound; it cannot bypass the material or CFL checks. See Discretisation: the exact scheme.

Strong anisotropy is a known trigger of C-PML corner instability. The multiaxial-PML option (mpml_fraction, typically 0.15–0.2) adds a small amount of cross-axis damping that cures it at the cost of slightly larger boundary reflections; see Boundary conditions.

Specifying sources and receivers

The channel-keyed mapping form is the canonical way to specify a survey:

out = deepwave.elastic_vti(
    ...,
    sources={"pressure_rate": (amplitudes, source_locations)},
    receivers={"p": receiver_locations},
)

This propagator’s source channels are pressure_rate, z, y, x and its receiver channels are p, z, y, x (the z channel is absent in 2-D). "pressure_rate" names the quantity the scalar source injects; see Units and source conventions. "p" is still accepted there, with a deepwave.DeepwaveFutureWarning, and is removed in 0.4.0. It is unchanged as a receiver channel name.

The equivalent flat source_amplitudes_*/source_locations_*/ receiver_locations_p keywords are equally permanent and produce bit-identical results; giving one channel in both forms raises. See Usage for the full vocabulary.