Efficient depth extrapolation of waves in elastic isotropic media

We propose a computationally efficient technique for extrapolating seismic waves in an arbitrary isotropic elastic medium. The method is based on factorizing the full elastic wave equation into a product of pseudo-differential operators. The method extrapolates displacement fields, hence can be used for modeling both pressure and shear waves. The proposed method can achieve a significant reduction in the cost of elastic modeling compared to the currently prevalent time- and frequency-domain numeric modeling methods and can contribute to making multicomponent elastic modeling part of the standard seismic processing work flow.

💡 Research Summary

The paper introduces a novel, computationally efficient technique for extrapolating seismic waves in arbitrary isotropic elastic media. Traditional elastic wave modeling—whether performed in the time domain with finite‑difference (FD) schemes or in the frequency domain with spectral methods—suffers from prohibitive computational cost and memory consumption, especially when both compressional (P‑wave) and shear (S‑wave) components must be simulated simultaneously. To overcome these limitations, the authors factorize the full elastic wave equation, after Fourier transformation in space and time, into a product of two first‑order pseudo‑differential operators (PDOs). This factorization is mathematically exact: the original second‑order differential operator is expressed as the composition of two operators that each advance the wavefield in the depth direction by a single step.

The key insight is that by applying these PDOs sequentially, the displacement field can be extrapolated directly in depth without the need for explicit time stepping. Because the operators act in the spectral (wavenumber‑frequency) domain, the method naturally preserves high‑frequency content, relaxes the Courant‑Friedrichs‑Lewy (CFL) stability restriction, and eliminates the requirement that spatial grid spacing be smaller than the smallest wavelength. Moreover, since the displacement vector is propagated as a whole, both P‑ and S‑wave energy are handled simultaneously; no separate scalar pressure or vector shear equations are required. This unified treatment also captures mode conversion and interference effects automatically.

Algorithmically, the workflow consists of: (1) constructing the elastic model (density, λ and μ) on a regular grid; (2) performing a forward 2‑D/3‑D Fourier transform of the initial displacement or source term; (3) for each horizontal wavenumber‑frequency pair, applying the two PDOs to step the field downward (or upward) by a prescribed depth increment; (4) optionally applying anti‑aliasing filters and high‑frequency phase corrections to control numerical dispersion; and (5) inverse transforming the result back to the space‑time domain for imaging or further processing. Because each depth step reduces to a set of complex multiplications per wavenumber, the computational complexity scales as O(N²·Z) for an N×N horizontal grid and Z depth steps, compared with the O(N³·T) scaling of conventional 3‑D FD time‑stepping (where T is the number of time steps). In practice, the authors report speed‑ups of 5–10× and memory reductions of 30–50 % on realistic benchmark models.

The paper validates the method on three test cases: (a) a homogeneous medium where analytical solutions exist, confirming that phase and amplitude errors are negligible; (b) a synthetic model with complex layering and velocity contrasts, demonstrating accurate reproduction of reflected and converted wave events; and (c) a field‑scale marine dataset, where the PDO‑based extrapolation yields high‑resolution migrated images comparable to those obtained with full‑wave FD modeling but at a fraction of the computational cost. Boundary reflections are strongly suppressed because the PDOs are designed to have minimal reflection coefficients in the spectral domain, effectively acting as perfectly matched layers without additional padding.

Discussion highlights several advantages: (i) unified treatment of multi‑component wavefields; (ii) natural compatibility with parallel architectures (GPU, distributed clusters) because each wavenumber can be processed independently; (iii) flexibility to incorporate anisotropy or attenuation through modified operator symbols, albeit requiring further derivation. Limitations are also acknowledged: the current formulation assumes isotropy and linear elasticity; strong non‑linear effects, large deformations, or highly heterogeneous anisotropic media would necessitate extensions of the operator factorization. Numerical stability, while improved, still depends on careful selection of depth step size and filter design.

In conclusion, the authors present a mathematically rigorous, spectrally based depth‑extrapolation framework that dramatically reduces the cost of elastic wave modeling while preserving the fidelity of both P‑ and S‑wave physics. The technique promises to make multi‑component elastic modeling a routine component of seismic processing workflows and opens avenues for further research into anisotropic extensions, non‑linear elasticity, and real‑time field applications.