Historically, the acoustic or elastic wave equation was solved on a discrete grid using first or second-order differencing in time and space (Kelly et al., 1976). The accuracy of the finite-difference method improves by using longer, more accurate spatial differencing operators, commonly 5-33 points long (Holberg, 1987). Taking this approach to its limit, the pseudospectral method (Kosloff et al., 1985) uses multiplication by ik in the wavenumber domain to compute spatial derivatives that are accurate to spatial Nyquist. In addition to the above improvements in the accuracy of spatial discretizations, Edwards et al. (1987) describe a new time-update method that eliminates errors due to time discretization.
These accurate modeling methods are efficient, especially in 3-D, because they allow the wave field and velocity model to be coarsely sampled. Unfortunately, there is also a drawback: if we sample the velocity model too coarsely it becomes difficult to represent curved or dipping boundaries accurately.
One solution is to use finite-element methods, which attack this problem by discretizing the model with elements that follow the boundaries. Fornberg (1988) introduced a mapping method that can be applied to pseudospectral finite-differences to convert curved-grid models to vertically layered models that can be more accurately represented on rectangular grids.
The main attraction of explicit finite-difference methods on rectangular grids is their efficiency and simplicity. Finite-element methods or other distorted-grid methods require complicated setup steps while the simpler rectangular-grid methods just need 2-D or 3-D arrays of values. So from a computational viewpoint, the real importance of the interpolation method we present in this paper is the ability to handle dipping interfaces without relinquishing regularly sampled grids.