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

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.