The complete set of governing equations for global dynamic and quasistatic problems –such as post-seismic and post-glacial rebound, tidal loading, and long-period seismology– involves a coupling between the conservation laws of continuum mechanics and Poisson/Laplace’s equation [3]. For dynamic problems, such as seismic wave propagation and the free oscillations of the Earth, it is possible to decouple Poisson’s equation using an explicit time marching scheme so that it can be solved independently [2]. For quasistatic problems, such as glacial isostatic adjustment and tidal loading, inertia is neglected, requiring an implicit time marching scheme [12]. In the latter case, Poisson’s equation cannot be decoupled. Although an explicit time scheme with an independent Poisson’s solver is generally fast, such an approach is limited by conditional stability, such that a very large number of time steps is often necessary. On the other hand, an implicit time scheme coupled with Poisson’s equation is generally slow but unconditionally stable. In both cases, the unbounded and large-scale nature of the problem poses numerical challenges, particularly for 3D Earth models. Most of the existing methods use spherical harmonics to solve the unbounded Poisson/Laplace’s equation. Such methods are often limited to spherically-symmetric models or have to rely on iterative procedures [13, 9, 1]. In view of these challenges, we develop a parallel software package based on the spectral-element method [6, 8, 7, 4, 5] combined with a mapped infinite-element approach [14]. While the spectral-element method is used within the Earth model, the infinite-element approach is employed in the outer region. In the infinite element approach, a so-called infinite-element layer is used to mimic all of space. The outermost edges of an element in the infinite-element layer are mapped to infinity in order to reproduce the behavior of gravitational potential outside the domain of interest, such that the potential decays to zero at infinity. Gauss- Lobatto-Legendre (GLL) quadrature is used for numerical integration in spectral elements. Since GLL quadrature cannot be used in infinite elements due to a singularity, we use Gauss- Radau quadrature instead. Spectral and infinite elements share identical quadrature points on infinite-element boundaries, thereby providing a natural coupling of the infinite-element method with the spectral-element method (Figure 1a). We use a generalized Maxwell rheology [11, 15] for viscoelastic deformation and accommodate topography and ellipticity. Both explicit and implicit time schemes are implemented in order to address a range of problems, including long-period seismology, glacial rebound (Figure 1b), tidal loading, etc.

References