An accurate and efficient scheme for wave propagation in linear viscoelastic media

The problem of wave propagation in a linear viscoelastic medium can be described mathematically as an exponential evolution operator of the form e (super M) (sub t) acting on a vector representing the initial conditions, where M is a spatial operator matrix and t is the time variable. Techniques like finite difference, for instance, are based on a Taylor expansion of this evolution operator. We propose an optimal polynomial approximation of e (super M) (sub t) based on the powerful method of interpolation in the complex plane, in a domain which includes the eigenvalues of the matrix M.The new time-integration technique is implemented to solve the isotropic viscoacoustic equation of motion. The algorithm is tested for the problem of wave propagation in a homogeneous medium and compared with second-order temporal differencing and the spectral Chebychev method. The algorithm solves efficiently, and with machine accuracy, the problem of seismic wave propagation with a dissipation mechanism in the sonic band, a characteristic of sedimentary rocks. The computational effort in forward modeling based on the new technique is half compared to that of temporal differencing when two-digit precision is required. Accuracy is very important, particularly for propagating distances of several wavelengths, since the anelastic effects should not be confused with nonphysical phenomena, such as numerical dispersion in time-stepping methods. Finally, we compute the seismic response to a single shot in a realistic geologic model which includes a gas-cap zone in an anticlinal fold. The results clearly show the importance of the attenuation and dispersion effects for an appropriate interpretation of the seismic data. In particular, the gas-cap response (a bright spot) suffers significant variations in amplitude, phase, and arrival time compared to the purely acoustic case.