A Novel Immersed Boundary Approach for Irregular Topography with Acoustic Wave Equations

Irregular terrain has a pronounced effect on the propagation of seismic and acoustic wavefields but is not straightforwardly reconciled with structured finite-difference (FD) methods used to model such phenomena. Accurate wavefield simulation is paramount in subsurface imaging applications such as reverse-time migration (RTM) and full-waveform inversion (FWI), requiring suitable topography handling. Methods currently detailed in the literature are generally limited in scope application-wise or non-trivial to apply to real-world geometries. A general immersed boundary treatment capable of imposing a range of boundary conditions in a relatively equation-agnostic manner has been developed, alongside a framework implementing this approach to complement emerging code-generation paradigms. The approach is distinguished by the use of N-dimensional Taylor-series extrapolants constrained by boundary conditions imposed at some suitably distributed set of surface points. The extrapolation process is encapsulated in modified derivative stencils applied in the vicinity of the boundary, utilizing hyperspherical support regions. This method ensures boundary representation is consistent with the FD discretization. Furthermore, high-dimensional and vector boundary conditions can be applied without approximation prior to discretization. A consistent methodology can thus be applied across free and rigid surfaces with first and second-order acoustic wave equation formulations. Application to both equations is demonstrated, and numerical examples based on analytic and real-world topography implementing free and rigid surfaces in 2D and 3D are presented. Numerical examples and convergence tests demonstrate the accuracy of boundary treatments devised by the prescribed approach, their suitability to practical applications, and the feasibility of automatically generating treatments to suit each case.

💡 Research Summary

The paper introduces a versatile immersed‑boundary (IB) methodology designed to handle irregular topography in acoustic wave simulations without sacrificing the simplicity and efficiency of structured finite‑difference (FD) schemes. The authors observe that conventional FD methods struggle with complex surfaces because the grid does not conform to the terrain, leading to large discretization errors near the boundary. Existing IB techniques either target very specific problems (e.g., only free surfaces) or require cumbersome preprocessing steps that are difficult to automate for realistic geological models.

To overcome these limitations, the authors develop a general, equation‑agnostic framework based on N‑dimensional Taylor‑series extrapolation. For each point in a hyperspherical support region surrounding the boundary, the acoustic field (pressure and, for second‑order formulations, velocity components) is expressed as a Taylor expansion about a set of surface points where the boundary conditions are known. The expansion coefficients are constrained directly by the prescribed boundary condition—zero traction for a free surface or zero velocity for a rigid surface—so that the extrapolated values automatically satisfy the physics at the interface. These extrapolated values replace the standard FD stencil entries that would otherwise cross the interface, yielding modified derivative stencils that are consistent with the underlying grid spacing and retain the formal order of accuracy of the original scheme.

A key strength of the approach is its flexibility: the order N of the Taylor series can be chosen independently of the spatial dimension, allowing straightforward extension from the first‑order acoustic pressure equation to the second‑order pressure‑velocity system. Moreover, vector‑valued boundary conditions are handled without any prior approximation; the constraints are imposed before discretization, eliminating the need for ad‑hoc “ghost‑cell” values or penalty terms.

The authors also integrate the IB treatment into a modern code‑generation paradigm. By treating the extrapolation and stencil modification as a reusable module, the system can automatically generate the appropriate finite‑difference operators for any given topography and boundary condition set. Users supply only the terrain description and the type of surface (free or rigid), and the code‑generator produces the complete solver, dramatically reducing the manual effort required to adapt existing FD codes to new geometries.

Numerical validation is performed on both synthetic and real‑world topographies in two and three dimensions. For a sinusoidal hill with a free surface, convergence studies show second‑order L2 error reduction as the grid is refined, confirming that the modified stencils preserve the nominal accuracy of the base FD scheme. In a second test, a rugged mountainous model with a rigid surface demonstrates correct reflection coefficients and wave‑field patterns, indicating that the rigid‑boundary condition is enforced accurately. A third, fully three‑dimensional experiment combines free and rigid surfaces on a realistic geological model and is used in a reverse‑time migration (RTM) scenario; the resulting migrated images exhibit markedly reduced boundary‑related artefacts compared with a conventional stair‑step FD implementation. Across all cases, the computational overhead introduced by the IB treatment remains modest because the hyperspherical support regions are limited to a few grid points around the interface.

In summary, the paper delivers a robust, high‑order immersed‑boundary technique that (1) embeds boundary conditions directly into Taylor‑series extrapolants, (2) works for both first‑ and second‑order acoustic equations, (3) supports free and rigid surfaces without additional approximations, (4) integrates seamlessly with automatic code generation, and (5) demonstrates accurate, convergent results on complex 2‑D and 3‑D topographies. The method’s generality and automation potential make it highly attractive for large‑scale seismic imaging workflows such as full‑waveform inversion (FWI) and reverse‑time migration, where accurate handling of irregular terrain is essential for reliable subsurface interpretation.