An Open-Source Tool for Accurate Topography Implementation within Finite Difference Wave Solvers

The presence of irregular topography can be highly problematic in seismic modelling if not implemented with care. Seismic waves interact with the surface, reflecting, focusing, and scattering: all behaviors which must be accounted for if their propagation is to be accurately simulated. Capturing topography effects is of great importance to the emerging field of teleseismic full-waveform inversion (FWI), and as such numerical methods using unstructured meshes are favored for wavefield modelling in these applications, in the place of the structured meshes typically used for traditional FWI. Topography introduces non-grid-aligned discontinuities in material properties, making accurate representation with finite difference (FD) methods challenging. By constructing appropriate extrapolations of the subsurface wavefield, it is possible to enforce the free-surface boundary condition at non-node locations, ensuring that topography is well resolved and preventing spurious scattering. We present an open-source Python module leveraging symbolic computation to provide a relatively general, high-level abstraction for immersed boundary implementations. For a given surface and discretization, a scheme based on spatially varying stencil coefficients is automatically devised, with weights provided in a format which can be directly utilized within FD models in Devito (an open-source domain-specific language for high-performance stencil computations). To demonstrate its functionality, scattering from uneven terrain, obtained from a digital elevation map of a mountainous area, is modelled using the constant-density acoustic wave equation. Preliminary test cases of a more complex scheme for application to the variable-density acoustic wave equation are also presented.