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

Abstract

Physics-based simulations provide a path to overcome the lack of observational data which is 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 more than one fault, complex geometries, and elastic heterogeneities are challenging. We present a symmetric interior penalty discontinuous Galerkin (SIPG) method to perform SEAS simulations accounting for the complex geometries and heterogeneity of the subsurface. Due to the discontinuous nature of the approximation, the spatial discretisation natively provides a mean to impose boundary and interface conditions associated with geometrically complex domains and embedded faults. 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 and interior boundaries and interfaces including intersecting 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. Further flexibility and efficiency is provided by optionally defining the displacement evaluation via a discrete Green’s function, which is evaluated once in an embarrassingly parallel pre-computation step using algorithmically optimal and scalable sparse parallel solvers and pre-conditioners.We illustrate the characteristics of the SIPG formulation via an extensive suite of verification problems (analytic, manufactured, and code comparison) for elasto-static and seismic cycle 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 geometry.Lastly, 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 multi-fault scenario of elasto-static instantaneous displacement due to 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. Our demonstrator models exploit advantages of both the boundary integral and volumetric methods and open 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).…”

Fused GEMMs towards an efficient GPU implementation of the ADER‐DG method in SeisSol

“…YATeTo therefore plays two important roles in SeisSol (and in similar applications 26 ): (1) it provides a convenient way for expressing DG-like numerical schemes for various wave propagation models 4,21,22 and (2) it works as an abstract performance portability layer between application and hardware.…”

Modeling Sequences of Earthquakes and Aseismic Slip (SEAS) in Elasto-plastic Fault Zones with A Hybrid Finite Element Spectral Boundary Integral scheme

“…The flexibility of domainbased methods allows for handling small-scale heterogeneities, material nonlinearities, as well as complexities of fault geometry (Kuna, 2013;Taborda & Bielak, 2011;Aagaard et al, 2013;Kaneko et al, 2011;Allison & Dunham, 2018;Erickson & Dunham, 2014;Thakur et al, 2020;Barbot, 2019;Mia et al, 2022;Erickson et al, 2017). However, modeling sequences of earthquakes and seismic slip with domain-based methods requires substantial computational effort due to the different spatial, and temporal scales (Tong & Lavier, 2018;Biemiller & Lavier, 2017;Kaneko et al, 2008;Allison & Dunham, 2018;Van Dinther et al, 2013;Mckay et al, 2019;Uphoff et al, 2022). Alternatively, boundary integral techniques reduce the dimensions of the problem by one, limiting the computations to the fault surface; thus, alleviating the computational cost (Aliabadi, 1997;Lapusta et al, 2000;Lapusta & Liu, 2009;Barbot, 2019;Cattania & Segall, 2021;Nie & Barbot, 2022;Heimisson & Rinaldi, 2022).…”