On the modelling of self-gravitation for full 3D global seismic wave propagation

The last decade has witnessed an unprecedented increase in high quality long period seismic data. This is because of the occurrence of several very large Earthquakes that were recorded on an exponentially growing number of broadband seismic stations that are installed in very dense networks such as the USarray. While for surface and body waves, tomography based on full waveform observables and 3D numerical simulations has become a standard tool in seismology, this is not the case for the normal modes. The main reason for this discrepancy is the fact that there is no established method to model the full physics of long period seismology in 3D and compute gradients with respect to material properties at reasonable computational cost.Here, we present a new approach to the inclusion of the full gravitational response to spectral element wave propagation solvers. We leverage the Salvus meshing software to include the external domain using adaptive mesh refinement and high order shape mapping. Together with Neumann boundary conditions based on a multipole expansion of the right hand side this minimizes the number of additional elements needed. Initial conditions for the iterative solution of the Poisson equation based on temporal extrapolation from previous time steps together with a polynomial multigrid method reduce the number of iterations needed for convergence. In summary this approach reduces the extra cost for simulating full gravity to a similar order as the elastic forces.We demonstrate the efficacy of the proposed method using the displacement from an elastic global wave propagation simulation at 200s period as the right hand side in the Poisson equation.