An accurate and efficient recursive convolution method to simulate viscoacoustic and viscoelastic waves

Seismic wavefield forward modeling in anelastic (attenuating) media is a fundamental tool for both data processing and interpretation in modern seismic exploration. We propose a generalized recursive convolution formula to calculate the temporal convolutions directly, rather than solving many auxiliary partial differential equations of the memory variables when dealing with a viscoelastic medium. The new formula is obtained in terms of the Taylor series expansion and offers approximations of the convolutions to arbitrary order by the number of terms retained. We conduct theoretical and numerical comparisons of the new method with the commonly used memory variable method and other traditional recursive convolution methods. The comparisons show that the new method has the highest accuracy of all these recursive convolution methods and yields better performance with various stress relaxation times and time steps than the common leapfrog time-stepping scheme to solve the auxiliary partial differential equations of the memory variables. Our numerical examples verify the versatility and feasibility of the new method for viscoacoustic and viscoelastic wave modeling.