A discontinuous Galerkin method for sequences of earthquakes and aseismic slip on multiple faults using unstructured curvilinear grids

Abstract

Summary Physics-based simulations provide a path to overcome the lack of observational data hampering a holistic understanding of earthquake faulting and crustal deformation across the vastly varying space-time scales governing the seismic cycle. However, simulations of sequences of earthquakes and aseismic slip (SEAS) including the complex geometries and heterogeneities of the subsurface are challenging. We present a symmetric interior penalty discontinuous Galerkin (SIPG) method to perform SEAS simulations accounting for the aforementioned challenges. Due to the discontinuous nature of the approximation, the spatial discretisation natively provides a means to impose boundary and interface conditions. The method accommodates two- and three-dimensional domains, is of arbitrary order, handles sub-element variations in material properties and supports isoparametric elements, i.e., high-order representations of the exterior boundaries, interior material interfaces and embedded faults. We provide an open-source reference implementation, Tandem, that utilises highly efficient kernels for evaluating the SIPG linear and bilinear forms, is inherently parallel and well suited to perform high resolution simulations on large scale distributed memory architectures. Additional flexibility and efficiency is provided by optionally defining the displacement evaluation via a discrete Green’s function approach, exploiting advantages of both the boundary integral and volumetric methods. The optional discrete Green’s functions are evaluated once in a pre-computation stage using algorithmically optimal and scalable sparse parallel solvers and preconditioners. We illustrate the characteristics of the SIPG formulation via an extensive suite of verification problems (analytic, manufactured, and code comparison) for elastostatic and quasi-dynamic problems. Our verification suite demonstrates that high-order convergence of the discrete solution can be achieved in space and time and highlights the benefits of using a high-order representation of the displacement, material properties, and geometries. We apply Tandem to realistic demonstration models consisting of a 2D SEAS multi-fault scenario on a shallowly dipping normal fault with four curved splay faults, and a 3D intersecting multi-fault scenario of elastostatic instantaneous displacement of the 2019 Ridgecrest, CA, earthquake sequence. We exploit the curvilinear geometry representation in both application examples and elucidate the importance of accurate stress (or displacement gradient) representation on-fault. This study entails several methodological novelties. We derive a sharp bound on the smallest value of the SIPG penalty ensuring stability for isotropic, elastic materials; define a new flux to incorporate embedded faults in a standard SIPG scheme; employ a hybrid multi-level preconditioner for the discrete elasticity problem; and demonstrate that curvilinear elements are specifically beneficial for volumetric SEAS simulations. We show that our method can be applied for solving interesting geophysical problems using massively parallel computing. Finally, this is the first time a discontinuous Galerkin method is published for the numerical simulations of sequences of earthquakes and aseismic slip, opening new avenues to pursue extreme scale 3D SEAS simulations in the future.

Geometric Control on Seismic Rupture and Earthquake Sequence Along the Yingxiu‐Beichuan Fault With Implications for the 2008 Wenchuan Earthquake

“…For example, a recent dynamic rupture simulation of the 2008 Wenchuan earthquake (Tang et al., 2021) did not include the Xiaoyudong fault in their multi‐fault system because of numerical stability issues. This problem may be addressed using the discontinuous galerkin (DG) method that better handles fault geometry complexities in dynamic rupture simulations (e.g., Ulrich, Gabriel, et al., 2019; Ulrich, Vater, et al., 2019; Wollherr et al., 2019), and also possibly for quasi‐dynamic earthquake sequence models (e.g., Kozdon et al., 2019; Uphoff et al., 2022). Besides the fault geometrical complexities, previous numerical studies have also indicated that a larger ratio between the width of seismogenic zone and the critical nucleation size ( h *) along with the variation of ( b − a )/ a will introduce more complexities in earthquake cycles and rupture pattern (Barbot, 2019; Cattania, 2019).…”

Modeling Multi‐Material Structural Patterns in Tectonic Flow With a Discontinuous Galerkin Level Set Method

“…Besides its wide application in fluid dynamics (e.g., Rhebergen et al, 2013), the DG method has also been applied in the geophysical community for modeling seismic wave propagation (Chung et al, 2015;Y. Zhang et al, 2019) and fault aseismic slips (Uphoff et al, 2022). It has also been utilized to discretize the level set formulation for multiphase fluid dynamics modeling (e.g., Jibben & Herrmann, 2017;Mousavi, 2014;Mousavi et al, 2016;Owkes & Desjardins, 2013;J.…”

“…The DG discretization is originally formulated for the neutron transport problem (Reed & Hill, 1973) and has since been widely applied to the modeling of hyperbolic conservation laws (Cockburn & Shu, 1991; Shu, 2014). Besides its wide application in fluid dynamics (e.g., Rhebergen et al., 2013), the DG method has also been applied in the geophysical community for modeling seismic wave propagation (Chung et al., 2015; Y. Zhang et al., 2019) and fault aseismic slips (Uphoff et al., 2022). It has also been utilized to discretize the level set formulation for multiphase fluid dynamics modeling (e.g., Jibben & Herrmann, 2017; Mousavi, 2014; Mousavi et al., 2016; Owkes & Desjardins, 2013; J. Zhang & Yue, 2019).…”

The role of three-dimensional fault interactions in creating complex seismic sequences

“…Segmentation may well be the most important factor that modulate the large earthquakes. Two recent implications of SEAS models hinted the possibility of more insight into the emergence of power-law distribution in large-scale and non-planar fault systems (Uphoff et al, 2022;Ozawa et al, 2022), where we could test the contribution of various factors to the frequency magnitude distribution.…”