The distributional finite difference method: an efficient method for modeling seismic wave propagation through 3D heterogeneous geological media.

We present a novel numerical method called the Distributional Finite Difference method (DFD) for modeling elastodynamic wave propagation in complex heterogeneous media. Efficient wave propagation modeling is crucial for solving the inverse problem where thousands of simulations are needed to infer the Earth's internal structure. The proposed method elegantly combines advantages of the finite difference method (FD) and of the spectral element method (SEM). In the past decades, the Spectral Element Method has become a popular alternative to the Finite Difference method for modeling wave propagation through the Earth. Though this can be debated, SEM is often considered to be more accurate and flexible than FD. This is because SEM has exponential convergence, it allows to accurately model material discontinuities, and complex structures can be meshed using multiple elements. In the mean time, FD is often thought to be simpler and more computationally efficient, in particular because it relies on structured meshed that are well adapted to computational architectures.  The DFD method divides the computational domain in multiple blocks or elements that can be arbitrary large. Within each block, the computational operations needed to model wave propagation are very similar that of FD which leads to high efficiency. When using smaller elements, the DFD approach allows to mesh certain regions of space having complex geometry, thus ensuring high flexibility. The DFD method permits simple specification of boundary conditions and accurately account for free surfaces and solid/fluid interfaces. Further, depending on the chosen basis functions, the DFD method can achieve "spectral like accuracy",  only a few (say 3) points per wavelength are need to model wave propagation accurately, i.e., with acceptable numerical dispersion. This results in significant reduction in memory usage. We present numerical examples demonstrating the advantages of the DFD method in various situations. We show that the DFD method si well adapted for modeling 3D wave propagation through the Earth.