Pseudo-acoustic VTI 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.acoustic_vti() propagator solves the pseudo-acoustic VTI wave equation of Duveneck et al. (2008): the V_S = 0 limit (Alkhalifah, 1998) of the elastic VTI equation, in a vertically transverse isotropic medium with variable density. Setting the shear-wave speed to zero collapses the elastic VTI system onto a first-order velocity–stress system with only two stress-like scalars, giving the industry-standard kinematic (“acoustic”) engine for anisotropic reverse-time migration and full-waveform inversion at a fraction of the cost of full elastic VTI.

It shares the numerical machinery of the isotropic elastic and elastic VTI propagators – the velocity–stress formulation on the same staggered grid, the C-PML (and optional multiaxial M-PML) absorbing boundaries, and the staggering of sources and receivers – and is the exact c55 = c66 = 0 reduction of deepwave.elastic_vti(). Reading the elastic VTI page first is recommended; this page describes only what the pseudo-acoustic reduction changes.

Note

The compiled C/CUDA backend, which backend="auto" (the default for this propagator, matching every other compiled propagator in Deepwave) resolves to, implements both forward propagation and its hand-derived adjoint: forward matches the eager backend to near machine precision, and gradients (every stiffness and the buoyancy) match the eager backend to near machine precision when the effective model-gradient stride is one and storage is lossless. Larger strides use the same sampled-gradient approximation as deepwave.elastic_vti(); the stride is out.info.step_ratio * model_gradient_sampling_interval, so interval 1 alone does not imply exact accumulation when CFL substepping is needed. This is the approximation shared with the elastic-VTI compiled backend under the same setting, not specific to this propagator. Both backends have been checked against the full physics/gradient test battery (the elastic_vti oracle, the isotropic-limit and analytic Christoffel oracles, reciprocity, PML absorption, torch.autograd.gradcheck, and the discrete adjoint dot-product identity), not merely against each other. The eager Python backend (backend="eager") instead provides gradients via native autograd through the time loop. backend="torch_compile" runs that same eager time loop with the two step functions compiled.

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.acoustic_vti() takes the three VTI stiffnesses, 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() (vs0 = 0, as the Parameterisation and units section below explains):

import torch
import deepwave

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

out = deepwave.acoustic_vti(
    c11=c11, c13=c13, c33=c33, buoyancy=buoyancy,
    grid_spacing=5, dt=0.001,
    source_amplitudes_sigma_h=deepwave.wavelets.ricker(25, 200, 0.001, 0.05).reshape(1, 1, -1),
    source_amplitudes_sigma_v=deepwave.wavelets.ricker(25, 200, 0.001, 0.05).reshape(1, 1, -1),
    source_locations_sigma_h=torch.tensor([[[5, 20]]]),
    source_locations_sigma_v=torch.tensor([[[5, 20]]]),
    receiver_locations_sigma_h=torch.tensor([[[35, 20]]]),
    pml_width=10,
    pml_freq=25,
)
out.receiver_amplitudes_sigma_h  # 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).

Fields and governing equations

When the shear moduli vanish the elastic VTI normal stresses collapse onto two scalars: the normal stress across the symmetry axis \(\sigma_H\) (equal in every direction perpendicular to it: \(\sigma_H = \sigma_{xx} = \sigma_{yy}\) in 3D, \(\sigma_H = \sigma_{xx}\) in 2D) and the normal stress along the symmetry axis \(\sigma_V\) (\(\sigma_V = \sigma_{zz}\) in 3D, \(\sigma_V = \sigma_{yy}\) in 2D). The v and h suffixes read as vertical and horizontal only because the symmetry axis is conventionally drawn vertical; Deepwave attaches no orientation to any dimension, so read them as along and across that axis. With \(b = 1/\rho\) the buoyancy and the symmetry axis along the slowest spatial dimension \(z\) (drop \(y\) in 2D), the system is

(1)\[\begin{split}\begin{aligned} \dot\sigma_H &= c_{11}\left(\partial_x v_x + \partial_y v_y\right) + c_{13}\,\partial_z v_z + s_H,\\ \dot\sigma_V &= c_{13}\left(\partial_x v_x + \partial_y v_y\right) + c_{33}\,\partial_z v_z + s_V,\\ \dot v_x &= b\,\partial_x \sigma_H,\qquad \dot v_y = b\,\partial_y \sigma_H,\qquad \dot v_z = b\,\partial_z \sigma_V, \end{aligned}\end{split}\]

where \(s_H\) and \(s_V\) are optional stress-rate source terms. Compared with running deepwave.elastic_vti() with zero shear moduli, this stores roughly half the fields (two normal stresses instead of the full stress tensor, and only \(2\,n_\mathrm{dim}\) C-PML memory variables) and evaluates none of the shear stencils.

There is no 1D form: with a single spatial dimension only \(c_{33}\) survives, which is the 1D scalar/acoustic equation, so 1D inputs raise an error (use deepwave.acoustic() or deepwave.scalar()).

Parameterisation and units

The three stiffnesses are the true (density-normalised-free) VTI stiffnesses, identical to those of deepwave.elastic_vti(), supplied at the normal-stress (integer) grid nodes with units of pressure (Pa in SI). In terms of the axial P velocity \(v_{P0}\) and Thomsen (1986) \(\epsilon,\delta\),

(2)\[c_{33} = \rho v_{P0}^2, \qquad c_{11} = c_{33}(1 + 2\epsilon), \qquad c_{13} = c_{33}\sqrt{1 + 2\delta}.\]

These are exactly what deepwave.anisotropy.thomsen_to_cij() produces with \(v_{S0} = 0\) (its general c13 relation reduces to \(c_{33}\sqrt{1+2\delta}\), and c55 = c66 = 0), so inversion may be performed differentiably in Thomsen space by chaining autograd through it; the returned c55 (identically zero) is simply discarded, while the buoyancy it also returns is passed straight to the propagator. That buoyancy \(b = 1/\rho\) is internally averaged onto the staggered velocity nodes exactly as in deepwave.elastic() (a buoyancy = 0 cell, with zero stiffnesses, is a vacuum, giving a free surface); the stiffnesses are used unaveraged.

Sources and receivers

Here \(h=\texttt{out.info.internal_dt}\) is the internal step. Sources and receivers may be placed on any velocity component (z – 3D only –, y, x; force-density sources, scaled internally by the buoyancy and \(h\) as in deepwave.elastic()) and on either normal stress (sigma_h, sigma_v; stress-rate sources \(s_H\), \(s_V\), scaled internally by \(h\)). An isotropic explosive/pressure source is the same wavelet applied to both the sigma_h and sigma_v channels.

Stability

The pseudo-acoustic VTI system is well-posed only where \(\epsilon \ge \delta\) everywhere, equivalently

(3)\[c_{11} c_{33} \ge c_{13}^2 \quad\text{(elementwise)}.\]

This is the Grechka–Zhang–Rector (2004) / Duveneck (2008) condition; deepwave.acoustic_vti() validates it and raises ValueError otherwise. The elliptical case \(\epsilon = \delta\) (\(c_{13} = \sqrt{c_{11}c_{33}}\)) is permitted and is a useful, artifact-free test point. For the CFL time step and PML profiles the maximum qP phase speed is bounded (when max_vel is not supplied) by \(\sqrt{\max(c_{11}, c_{33})\,b}\) pointwise, which is attained on a symmetry axis. Strong anisotropy is a known trigger of C-PML corner instability; the multiaxial-PML option (mpml_fraction, typically 0.15–0.2) cures it at the cost of slightly larger boundary reflections (see Boundary conditions).

Pseudo-shear source artifact

Where \(\epsilon \ne \delta\), a point source in the pseudo-acoustic VTI system emits, in addition to the qP wave, a spurious diamond/rhombic pseudo-shear wave – the kinematic remnant of the shear wave that the V_S = 0 approximation is supposed to remove. This is intrinsic to the formulation, not a numerical defect. The optional helper deepwave.anisotropy.source_taper() mitigates it, following Duveneck et al. (2008), by blending \(c_{13}\) towards the locally elliptical value \(\sqrt{c_{11}c_{33}}\) (\(\delta \to \epsilon\), for which no pseudo-shear is emitted) within a small raised-cosine radius of each source.

What it returns is a modified \(c_{13}\) model, not a mask: the propagator you then hand it to is running on a medium that differs from the one you built, within radius cells of every source, where the anisotropy has been made elliptical. That is what suppresses the artefact, but it also means the near-source wavefield is not your model’s wavefield, and a gradient taken with respect to the tapered \(c_{13}\) is not the gradient with respect to yours there. Keep the radius to a few cells, and leave the helper alone if the near-source medium is what you are trying to recover.

Correctness

Because deepwave.acoustic_vti() is the exact c55 = c66 = 0 reduction of deepwave.elastic_vti(), its forward field is bit-identical to elastic_vti with zero shear moduli (\(\sigma_V = \sigma_{zz}\), \(\sigma_H = \sigma_{xx} = \sigma_{yy}\), velocities equal) in 2D and 3D – the primary correctness check in the test suite. Independently, the measured qP phase velocity across angle (from band-limited plane waves in a homogeneous strongly-anisotropic medium) matches the exact V_S = 0 Christoffel solution within the grid-dispersion budget, and every model and source gradient matches finite differences.

References

Duveneck, E., Milcik, P., Bakker, P. M., and Perkins, C. (2008), “Acoustic VTI wave equations and their application for anisotropic reverse-time migration”, SEG Technical Program Expanded Abstracts, 2186–2190, doi:10.1190/1.3059320.

Alkhalifah, T. (1998), “Acoustic approximations for processing in transversely isotropic media”, Geophysics 63, 623–631.

Grechka, V., Zhang, L., and Rector, J. W. (2004), “Shear waves in acoustic anisotropic media”, Geophysics 69, 576–582.

Specifying sources and receivers

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

out = deepwave.acoustic_vti(
    ...,
    sources={"y": (amplitudes, source_locations)},
    receivers={"y": receiver_locations},
)

This propagator’s source channels are pressure_rate, z, y, x, sigma_v, sigma_h and its receiver channels are p, z, y, x, sigma_v, sigma_h (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.

Both scalar channels are synthesised, because this propagator has no pressure variable. The pressure_rate source drives both stress channels, each receiving -source_amplitudes_pressure_rate. The p receiver records minus the mean normal stress, -(sigma_v + (ndim - 1) * sigma_h) / ndim – the pressure a hydrophone in this medium measures, and the same quantity every other propagator’s p receiver records.

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