Notes on finite-difference wave propagation

2 minute read

Finite-difference method

“If one wants to solve an ODE or PDE to high accuracy on a simple domain, and if the data defining the problem are smooth, then spectral methods are usually the best tool. They can often achieve ten digits of accuracy where a finite-difference or finite element method would get two or three. At lower accuracies, they demand less computer memory than the alternatives.” Trefethen (2000)

O(T,S) Medium PPW PD Reference
2-2 Homo 10 10 (Kelly et al., 1976)
2-4 Homo 5 30 (Levander, 1989)
2-2 Heter 15-20 10  
2-4 Heter 10-15 30  

where O(T,S) is order of accuracy in time and space, PPW is number of points per wavelength and PD is propagation distance in terms of wavelength.

If the gridpoint spacing is too coarse, dipping interfaces appear as stair-steps, where the edge of each step acts as a strong diffractor. If the waves propagate more than 30 wavelengths, empirical evidence suggests higher-order schemes of 2-8 or higher are needed. Stork (2013) claims that numerical dispersion errors should be no more than π/20 after propagating more than 50 wavelengths, which is equivalent to a 0.05% phase error. Otherwise there can be a fatal buildup of errors when iterative methods such as FWI or least squares migration are used. He suggests that longer 25-point stencils using optimized coefficients (Zhang and Yao, 2013) are much better at controlling numerical dispersion.

Today’s best practice is to use higher-order finite-difference stencils, with recommendations of 8th-order in space and 2nd-order in time. For media with strong contrasts in physical properties (water next to elastic sediments or air), the staggered grid scheme is recommended over the conventional FD solutions to the second-order wave equation.

Pseudospectral method

Spatial derivative in the Fourier domain. (Kreiss and Oliger, 1972; Kosloff and Baysal, 1982; Kosloff et al., 1982). In theory, only two points per wavelength are needed for an accurate solution in a smooth medium.

ince the FFT represents the derivative in the wavenumber domain by the exact multiplication of ikx, then its spatial representation is equivalent to a discrete derivative operator that is as wide as the computational domain. It requires a periodic boundary condition, but if the absorbing boundary conditions are effective (see Appendix 8.9 and Liu and Sen (2010, 2012)) then this should not be a problem.

Spectral Element VS Finite Element

The SE method differs from the standard FE element in that its nodes are at the unevenly spaced quadrature points (GLL), and therefore removes the need to make the lumped mass approximation.

Sponge boundary

The simplest absorbing boundary condition is that of a damping zone (Cerjan et al., 1985), also known as a sponge zone, with a thickness of about 5-10 wavelengths located next to the sides of the model. Inside the leftside sponge region the wavefield is multiplied by an exponential damping function f(i) at each time step, where i is the gridpoint index for the x-coordinate. A modified damping function is proposed by Wang and Qin (1997)

For staggered grid, damping just the particle-velocity component after each time step is less effective than damping both the stress and particle-velocity components after each time step. However, an instability has been observed when both stress and velocity components are damped for the elastic wave equation.


For a viscoelastic fluid, the Maxwell model predicts that stress decays exponentially to zero, which is accurate for a viscoelastic fluid, not a solid.