Efficient numerical solution of global dynamic and quasistatic problems using a spectral-element method coupled with an infinite-element approach

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.

As an example, we compute the viscoelastic response of the Earth under glacial loading. We consider an Earth model that includes surface topography, bathymetry, and ellipticity. We use ICE-5G VM2 [10] for the loading history and viscosity profile. Figure 1b shows a typical response obtained solving a set of governing equations including Poisson’s equation. 

Authors

References

[1] D. Al-Attar and J. Tromp. Sensitivity kernels for viscoelastic loading based on adjoint methods. Geophysical Journal International, 196(1):34–77, 2014. [2]  E. Chaljub and B. Valette. Spectral element modelling of three-dimensional wave propagation in a self-gravitating Earth with an arbitrarily stratified outer core. Geophysical Journal International, 158(1):131–141, 2004.[3]  F. A. Dahlen and J. Tromp. Theoretical Global Seismology. Princeton University Press, 1998.

[4]  H. N. Gharti, D. Komatitsch, V. Oye, R. Martin, and J. Tromp. Application of an elastoplastic spectral-element method to 3D slope stability analysis. International Journal for Numerical Methods in Engineering, 91:1–26, 2012.

[5]  H. N. Gharti, V. Oye, D. Komatitsch, and J. Tromp. Simulation of multistage excavation based on a 3D spectral-element method. Computers & Structures, 100–101:54–69, 2012.

[6]  D. Komatitsch and J. Tromp. Introduction to the spectral element method for three- dimensional seismic wave propagation. Geophysical Journal International, 139:806–822, 1999.
[7]  D. Komatitsch and J. Tromp. Spectral-element simulations of global seismic wave propagation – I. validation. Geophysical Journal International, 149:390–412, 2002.
[8]  D. Komatitsch and J. Tromp. Spectral-element simulations of global seismic wave propagation – II. three-dimensional models, oceans, rotation and self-gravitation. Geophysical Journal International, 150:303–318, 2002.
[9]  K. Latychev, J. X. Mitrovica, J. Tromp, M. E. Tamisiea, D. Komatitsch, and C. C. Christara. Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation. Geophysical Journal International, 161(2):421–444, 2005.
[10]  W. R. Peltier. Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE. Annual Review of Earth and Planetary Sciences, 32(1):111– 149, 2004.
[11]  J. C. Simo and T. J. R. Hughes. Computational inelasticity. Springer, 1998.

[12]  J. Tromp and J. X. Mitrovica. Surface loading of a viscoelastic Earth–II. spherical models. Geophysical Journal International, 137(3):856–872, 1999.
[13]  S. Zhong, A. Paulson, and J. Wahr. Three-dimensional finite-element modelling of Earth’s viscoelastic deformation: effects of lateral variations in lithospheric thickness. Geophysical Journal International, 155(2):679–695, 2003.
[14]  O. C. Zienkiewicz, C. Emson, and P. Bettess. A novel boundary infinite element. International Journal for Numerical Methods in Engineering, 19(3):393–404, 1983.
[15]  O. C. Zienkiewicz and R. L. Taylor. The finite element method for solid and structural mechanics. Elsevier Butterworth-Heinemann, 2005.