Modeling Three-Dimensional Wave Propagation in Anelastic Models With Surface Topography by the Optimal Strong Stability Preserving Runge-Kutta Method

Accurate and efficient forward modeling methods are important for simulation of seismic wave propagation in 3D realistic Earth models and crucial for high-resolution full waveform inversion. In the presence of attenuation, wavefield simulation could be inaccurate or unstable over time if not well treated, indicating the importance of the implementation of a strong stability preserving time discretization scheme. In this study, to solve the anelastic wave equation, we choose the optimal strong stability preserving Runge-Kutta (SSPRK) method for the temporal discretization and apply the fourth-order MacCormack scheme for the spatial discretization. We approximate the rheological behaviors of the Earth by using the generalized Maxwell body model and use an optimization procedure to calculate the anelastic coefficients determined by the Q() law. This optimization constrains positivity of the anelastic coefficients and ensures the decay of total energy with time, resulting in a stable viscoelastic system even in the presence of strong attenuation. Moreover, we perform theoretical and numerical analyses of the SSPRK method, including the stability criteria and the numerical dispersion. Compared with the traditional fourth-order Runge-Kutta method, the SSPRK has a larger stability condition number and can better suppress numerical dispersion. We use the complex-frequency-shifted perfectly matched layer for the absorbing boundary conditions based on the auxiliary difference equation and employ the traction image method for the free-surface boundary condition on curvilinear grids representing the surface topography. Finally, we perform several numerical experiments to demonstrate the accuracy of our anelastic modeling in the presence of surface topography.