Frequency domain finite-difference wave propagation modeling in 2D acoustic media

Seismic wave propagation modeling is helpful in understanding waveform behavior involving velocity anomalies due to complicated geological structures. On the other hand, forward modeling is the kernel of inversion and tomography procedures, so developing forward modeling algorithms is an inevitable need. In order to model wave propagation, it is necessary to solve complete wave equation for 3D media. Because of limitations in computational resources, solving complete wave equation with all ideal considerations such as heterogeneity, anisotropy, and absorption, for 3D media is difficult. Thus by considering easier conditions we can solve wave equation with normal computers. Solving wave equation with computers needs to discretization techniques such as boundary integral, finite-element or finite-difference methods. Boundary integral methods are suitable for simple models while finite-element or finite-difference methods are predominantly used for heterogeneous models. Depending on the domain for which wave equation is going to be solved, we can categorize methods to time-space, frequency-space, Laplace, slowness-space and etc. Recently, the frequency domain finite-difference (FDFD) method has found extensive application in multi-source experiment modeling, especially in waveform tomography. Because of the special form of wave equation in the frequency domain, the modeling of multi-source experiments is a straightforward job. On the other hand, considering absorption mechanisms is easy. This study deals with wave propagation modeling in a 2D acoustic heterogeneous media, using the second-order approximation of acoustic wave equation in the frequency domain. The acoustic wave equation is formulated as the first-order hyperbolic system involving fields of pressure and particle velocities. This system is discretized using the second-order staggered grid stencil. To avoid spurious reflections from the model boundaries, a sponge-like perfectly matched layer (PML) is implemented. Solving wave equation in the frequency domain leads to a large matrix equation. The key step in the frequency domain finite-difference modeling that controls computational efficiency is the numerical inversion of the massive matrix. The matrix structure depends on the spatial derivatives approximation. In order to solve this system, different direct solvers can be used. The UMFPACK (unsymmetric multifrontal sparse LU factorization package) solver, which is embedded in MATLAB and has acceptable performance in solving the general system of equations, was selected for this study. For each frequency component, it is necessary to solve a large system of equations to obtain a single frequency component of the pressure wavefield. In this paper we review criteria to avoid numerical dispersion and errors during the finite-difference approximation, and time aliasing during frequency sampling. Choosing appropriate spatial and frequency discretization intervals is very important in FDFD modeling. Choosing large Δ (spatial discretization interval) values will cause the pressure field to be inadequately sampled in space and numerical dispersion. For the scheme presented here, we should have more than 10 grid points per minimum wavelength to keep dispersion errors small. On the other hand according to the sampling theorem, if, where tmax is maximum time and df is frequency components sampling interval, we encounter time aliasing. The present FDFD package is written in MATLAB programming language. MATLAB supports sparse algebra and makes possible the solution of large system of equations with minimum usage of memory, which is the main concern in FDFD algorithms. Examples of waveform modeling using FDFD are shown. In the first example, the headwave originating from a high velocity layer is modeled; in the second example, the wave behavior in a trap model is shown; and in the last example, a salt dome model is studied.