Accurate viscoelastic modeling by frequency‐domain finite differences using rotated operators

Abstract

The viscoelastic wave equation is an integro‐differential equation that requires special methods when using time‐domain numerical finite‐difference methods. In the frequency domain, the integral terms are easily represented by complex valued elastic media properties. There are further significant advantages to using the frequency domain if the forward or the inverse problem requires modeling or inverting a large number of prestack source gathers. Numerical modeling is expensive for seismic data because of the large number of wavelenghths typically separating sources from receivers, which results in a need for a large number of grid points. A major obstacle to using frequency‐domain methods is the consequent storage requirements. To reduce these, we maximize the accuracy and simultaneously minimize the spatial extent of the numerical operators. We achieve this by extending earlier published methods introduced for the viscoacoustic case to the viscoelastic case. This requires the formulation of two new numerical operators: a differencing operator in a rotated coordinate frame and a lumped mass term. The new operators are combined with ordinary second‐ order, finite‐difference operators in an optimal manner to minimize numerical errors without increasing the size of the numerical operator. For a fixed number of grid points, the resulting second‐ order differencing scheme is no more expensive than an ordinary second‐order differencing scheme, but a numerical dispersion analysis shows that the number of grid points required per smallest wavelength is reduced from approximately 15 to approximately 4. The new scheme is also capable of handling embedded fluid layers without instability. We demonstrate that no further improvement in performance can be achieved using higher order spatial operators because of the associated computational overheads associated with the larger differencing operators. The new viscoelastic modeling scheme is used to study a crosshole data set in which the exact nature of the seismic coda is unclear. The results of the modeling study indicate this coda is likely related to the generation of mode‐converted shear waves within the complicated, finely layered sediments at the site.