Summary
Although observation of gravity perturbations induced by earthquakes is possible, simulation of seismic wave propagation in a self-gravitating, rotating Earth model with 3D heterogeneity is challenging due to the numerical complexities associated with the unbounded Poisson/Laplace equation that governs gravity perturbations. Therefore, gravity perturbations are generally omitted, and only the background gravity is taken into account using the so-called Cowling approximation. However, gravity perturbations may be significant for large earthquakes (Mw ≥ 6.0) and long-period responses. In this study, we develop a time-domain solver based on the spectral-infinite-element approach, which combines the spectral-element method inside the Earth domain with a mapped-infinite-element method in the infinite space outside. This combination allows us to solve the complete, coupled momentum-gravitational equations in a fully discretized domain while accommodating complex 3D Earth models. We compute displacement and gravity perturbations considering various Earth models, including PREM and S40RTS and conduct comprehensive benchmarks of our method against the spherical harmonics normal-mode approach and the direct radial integration method. Our 3D simulations accommodate topography, bathymetry, rotation, ellipticity, and oceans. Results show that our technique is accurate and stable for long simulations. Our method provides a new scope for incorporating earthquake-induced gravity perturbations into source and adjoint tomographic inversions.