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. Key Points:• We develop a new optimal strong stability preserving Runge-Kutta (SSPRK) method for solving the anelastic wave equation • We perform theoretical and numerical analyses of the SSPRK method • We demonstrate the accuracy and efficiency of the SSPRK method in anelastic wavefield modeling in the presence of surface topographySupporting Information:• Supporting Information S1 ). Modeling three-dimensional wave propagation in anelastic models with surface topography by the optimal strong stability preserving Runge-Kutta method. Withers et al., 2015) for a weak frequency dependence of Q is sometimes of interest. We follow a conventional way to calculate the anelastic coefficients determined by the Q(ω) law (Emmerich & Korn, 1987;Käser et al., 2007;Kristek & Moczo, 2003), but we constrain the coefficients to be positive (Blanc et al., 2016;Yang et al., 2016). This positivity ensures the decay of total energy over time as proved in the supporting information (SI) section A. Assuming the attenuation factors at a reference frequency f r for P and S waves is Q P and Q S , the frequency range of interest in the classical approach to calculate the anelastic coefficients is defined as [f min , f max ] = [f r /f amp , f r × f amp ],

At seismic discontinuities in the crust and mantle, part of the compressional wave energy converts to shear wave, and vice versa. These converted waves have been widely used in receiver function (RF) studies to image discontinuity structures in the Earth. While generally successful, the conventional RF method has its limitations and is suited mostly to flat or gently dipping structures. Among the efforts to overcome the limitations of the conventional RF method is the development of the wave‐theory‐based, passive‐source reverse‐time migration (PS‐RTM) for imaging complex seismic discontinuities and scatters. To date, PS‐RTM has been implemented only in 2D in the Cartesian coordinate for local problems and thus has limited applicability. In this paper, we introduce a 3D PS‐RTM approach in the spherical coordinate, which is better suited for regional and global problems. New computational procedures are developed to reduce artifacts and enhance migrated images, including back‐propagating the main arrival and the coda containing the converted waves separately, using a modified Helmholtz decomposition operator to separate the P and S modes in the back‐propagated wavefields, and applying an imaging condition that maintains a consistent polarity for a given velocity contrast. Our new approach allows us to use migration velocity models with realistic velocity discontinuities, improving accuracy of the migrated images. We present several synthetic experiments to demonstrate the method, using regional and teleseismic sources. The results show that both regional and teleseismic sources can illuminate complex structures and this method is well suited for imaging dipping interfaces and sharp lateral changes in discontinuity structures.

To date developments of seismic attenuation models of the Earth have lagged those of velocity models. This is partly due to difficulties in isolating waveform perturbations caused by attenuation and velocity heterogeneities and partly due to different theories (e.g., ray theory, finite frequency theory) used in most previous and current studies to approximate sensitivity kernels and invert for the attenuation structure. We present in this paper a new method for computing the 3D waveform Fréchet kernels that account for full physical‐dispersion and dissipation attenuation. For solving the 3D isotropic anelastic wave equation described by the generalized Maxwell body model, we extend our previously proposed full waveform modeling method in cartesian coordinates to that in spherical coordinates, which can provide stable numerical solutions even in the presence of strong attenuation. Then, we apply the scattering‐integral method for calculating 3D travel time and amplitude sensitivity kernels with respect to velocity and attenuation structures. We demonstrate the accuracy of our forward method and the effectiveness of the implementation of absorbing and free surface boundary conditions through numerical tests. Moreover, by choosing the Northwestern United States region as a realistic example, we verify the accuracy of the computed 3D sensitivity kernels through comparing the waveform measurements with predictions from the kernels. Finally, we discuss the importance of calculating full anelastic sensitivity kernels including both effects of physical dispersion and dissipation, where we specially explore the effect of scattering due to random velocity and attenuation heterogeneities on waveform measurements.