A three-dimensional hybrid finite element -- spectral boundary integral method for modeling earthquakes in complex unbounded domains

We present a 3D hybrid method which combines the Finite Element Method (FEM) and the Spectral Boundary Integral method (SBIM) to model nonlinear problems in unbounded domains. The flexibility of FEM is used to model the complex, heterogeneous, and nonlinear part -- such as the dynamic rupture along a fault with near fault plasticity -- and the high accuracy and computational efficiency of SBIM is used to simulate the exterior half spaces perfectly truncating all incident waves. The exact truncation allows us to greatly reduce the domain of spatial discretization compared to a traditional FEM approach, leading to considerable savings in computational cost and memory requirements. The coupling of FEM and SBIM is achieved by the exchange of traction and displacement boundary conditions at the computationally defined boundary. The method is suited to implementation on massively parallel computers. We validate the developed method by means of a benchmark problem. Three more complex examples with a low velocity fault zone, low velocity off-fault inclusion, and interaction of multiple faults, respectively, demonstrate the capability of the hybrid scheme in solving problems of very large sizes. Finally, we discuss potential applications of the hybrid method for problems in geophysics and engineering.

zone physics 30 as well as extension to the quasi-dynamic limit and cycle simulations. 31 The current extension to the full three 66 dimensional case represents the culmination of these efforts.

67
The remainder of this paper is organized as follows. In Section 2, we describe the physical model (Section 2.1), and the 68 numerical methods to solve it, which includes the finite-element method (Section 2.2), the spectral boundary integral method 69 (Section 2.3), and their coupling -the hybrid method (Section 2.4). In Section 3, we validate the hybrid method using the 70 benchmark problem TPV3 32 of the Southern California Earthquake Center. Next, we demonstrate the capabilities of the new 71 hybrid method on more complex problems. We consider a low velocity fault zone in Section 4, a low velocity inclusion at a 72 distance from the fault in Section 5, and interacting faults in Section 6. Finally, we discuss the advantages of the hybrid method 73 in terms of computational cost in Section 7 and draw conclusions in Section 8.
where is the material density, the displacement vector, with the "dot" being the derivative with respect to time , the 79 Cauchy stress tensor, and the coordinate vector. Body forces are neglected. Dirichlet boundary conditions are applied on 80 and Neumann boundary conditions are applied on where is the normal vector to the surface . Initially, the domain is assumed to be in equilibrium, and, hence, the initial 82 conditions are given by (0) = 0 anḋ (0) = 0. We assume linear elastic material behavior: with the infinitesimal strain tensor = ( ∕ + ∕ )∕2 and the Lamé parameters , describing the elastic properties 84 of the material. 85 The fault transmits stresses from one half-space to the other through interface tractions. We focus on tangential (friction) 86 interaction and impose non-penetration/non-opening conditions to the normal component of the fault surfaces ± . The local 87 slip vector, which corresponds to the tangential fault opening vector, is given by where is the global-to-local rotation matrix and + and − indicate the upper and lower fault sides, respectively. Local slip 89 is the amplitude of . The fault is governed by a stick-slip behavior that is described by two states. A sticking section of the 90 fault is described by: wherė is slip rate, the amplitude of the fault shear traction vector and the fault strength. A sliding fault section is described 92 by: A constitutive law is applied to model the fault strength evolution. For simplicity, we apply a linear slip-weakening friction 94 law, 33 which is given by:  wherë denotes the second time derivative of the displacement vector, and are the mass and stiffness matrix, respectively, 108 is a fault rotation-area matrix, is the fault traction vector, and is the force vector from Neumann boundary conditions. 109 We apply an explicit central-difference time integration formulation with a predictor-corrector formulation. The step-by-step 110 procedure follows: where the subscript indicates the time step and Δ is the current incremental time step, which is required to satisfy the Courant- where is the fault impedance matrix given by −1 = Δ where t+1 and̃ t+1 are individual entries in t+1 and̃ t+1 , respectively, and t+1 is the fault strength at each split-node (node 125 indicator is omitted for simplicity) and is governed by Eq. (8).
where is a diagonal matrix with 11 = 33 = 1 and 22 = s ∕ p . p and s are the longitudinal and shear wave speeds of 136 the material, respectively. Eq. (19) states that the traction at the surface of the half space, SBI , equals the far field traction, ∞ , 137 plus a "radiation damping" term, ṡ ± , and a spatiotemporal integral term . In this formulation the elastodynamic response 138 is separated between local and nonlocal contributions. represents the nonlocal elastodynamic long-range interaction between 139 different parts of the surface, and the local effect ṡ represents wave radiation from the surface. 140 We use the spectral approach 49 for computing which involves a Fourier transform in space and a convolutions in time,

141
where the displacement history is convolved with the elastodynamic kernels. Please refer to independent formulation in 49 for 142 the derivation of the kernels and details for computing the nonlocal term, .

144
The hybrid method consists in coupling the FEM and the SBIM at the boundaries ± , where the FE-domain is truncated (see 145 Fig. 1c). At ± , which we refer to as the virtual boundary, we apply an exact elastodynamic transparent boundary condition 146 using the SBIM, which accounts for wave propagation in the infinite half-space beyond the FE truncation. Alternative coupling methods, e.g., Langrange multiplier, could also be applied. However, as we will show in Sec. 3, the simple 163 staggered approach proposed here provides excellent accuracy and is optimal in terms of computational efficiency. Further, in 164 the finite-element domain, we apply 8-node linear hexagonal elements in a regular mesh in all presented problems.

187
We present the results of the hybrid method with a virtual strip width 2 = 4Δ (see Fig. 2b) and compare it with the reference higher precision spectral boundary method. Therefore, the convergence rate of the hybrid method corresponds to the rate of the 197 least accurate of the methods it combines, i.e., the linear finite element method. Hence, the hybrid method does not loose any 198 accuracy compared to a fully FEM model.

199
Since we are dealing with a dynamic problem, we also show the temporal evolution of the error at station C (see Fig. 5b). The  Overall, we find that the error remains mostly constant over the duration of the simulation and decreases with mesh refinement.  The previous example was a benchmark problem and could have been solved by a boundary-element approach without any 215 discretization of the bulk. Hence, the hybrid method was not required. In the following, we will consider more complex problems, 216 which require volumetric discretization. First, we consider a slip-weakening fault with a low velocity fault zone (LVFZ). LVFZ 217 are found in most mature faults, where the near fault rock is considerably damaged and, as a consequence, has a reduced wave 218 speed ranging from 20% to 60% with respect to the host rock. 50,51,52,53,54 In 2D setups, when the reduction is high enough, the 219 rupture behaves like a pulse. The results presented here will confirm this behavior on a 3D setup. 220 We consider a velocity reduction of 20% with respect to the surrounding host rock, which has the same elastic properties as

226
As a result of the nucleation procedure, the rupture front quickly propagates radially and eventually spans the entire fault (see 227 Fig. 7). When a dynamic rupture propagates, it radiates elastic waves, which are then reflected at the boundary of the LVFZ. 54

228
Depending on the incident angle, the reflected wave can have an inverted polarization and cause unloading of the fault and 229 generate a slip pulse. This effect is also observed in our 3D simulations and is shown in Fig. 7b, station C, and in Fig. 8. The 230 reflected wave causes the rupture to split into a pulse-like rupture, followed by a crack-like rupture.

231
The rapid acceleration and deceleration of a slip pulse are a source of high frequencies (see Fig. 9b Station C) and cause 232 oscillations in slip velocity, trailing the rupture front (see Fig. 8). These oscillations do not affect the rupture propagation and 233 disappear with further mesh refinement and regularized friction laws. 55 Additionally, numerical damping, which is not used 234 here, is often applied to minimize such high frequencies. Since this problem cannot be solved with the SBI method, we validate while beyond a distance ∕2 from the fault plane the material has a 20% velocity reduction.

246
Similar as in the LVFZ setup of Sec. 4, elastic waves are reflected at the boundary of the low velocity inclusion and affect the 247 shear stress at the interface. However, the reflected waves have the same polarity as the incident ones and, hence, increase the 248 shear stress in front of the propagating rupture front (see Fig. 9b). This increasing shear stress peak eventually causes the rupture 249 to transition from subRayleigh to supershear velocities (see Fig. 9a and Fig. 10). SubRayleigh propagation occurs when the 250 rupture speed is lower than the Rayleigh wave speed, ≈ 0.9 , and can be observed at stations A and B in Fig. 9. Supershear 251 propagation, however, refers to ruptures propagating faster than and their speed can approach the limiting speed, . 56,57 In 252 our simulations, supershear rupture occurs within the domain surrounding station C in Fig. 9a.

253
In this 3D simulation, we observe a supershear transition through the Burridge-Andrews mechanism, 58,6 where a shear stress 254 peak in front of the existing crack nucleates the supershear rupture (see Fig. 9b, station B and Fig. 10). In contrast to 2D 255 setups, 50,54 the extent of the supershear rupture is confined to a triangular shaped region, which surrounds station C. Additionally, 256 the transition occurs progressively: first at 3 = 0 km at ≈ 7 s, then it expands towards the ± 3 direction, and finally, at ≈ 12 s 257 it spans the entire seismogenic depth. Finally, we present an example of interaction between nearby faults, i.e. two fracture planes side-by-side. We consider a dilational 265 step-over geometry with a system of two faults that overlap each other (see Fig. 11). The dilational step-over implies that

279
Using the hybrid method, we can successfully reproduce the results of Bai and Ampuero 59 (see Fig. 12): after nucleation, the  sizes, which confirms the linear relationship between computational cost and the number of degrees of freedom (see Fig. 13).

299
The computational saving of the hybrid method compared to a standard FEM lies in the reduction of 2 due to the truncation of 300 the FE domain. Moreover, the added overhead cost of the SBI as wave absorption algorithm is in the same order of magnitude 301 of only one layer of FE elements (see Fig. 13b). Therefore, it is practically negligible. For example, for a full FE simulation, 2 must be in the order of 1 to prevent artificial reflections at the domain boundary.

303
However, using the hybrid method the domain size can be truncated up to the extent of the nonlinear region or the extent 304 of the elastic heterogeneity, which are usually one to two orders of magnitude smaller than 1 . 50 Assuming a regular spatial 305 discretization, the domain truncation results in a reduction of 2 by one order of magnitude and so will the computational cost.

306
The savings may be even higher in other applications.

307
All our simulations were performed using distributed memory parallel computing with 48 threads and for the largest sim-308 ulations 96 threads. Therefore, in Fig. 13, we report the computational time multiplied by the number of parallel processes.

309
However, the SBIM library that we are using also supports shared memory parallelism and it is designed to be easily coupled 310 to any FEM library written in C++. The only requirement is that the FEM mesh at the virtual boundary ± is a regular grid, 311 due to the spectral representation of the boundary integral equations.

312
These computational savings represent an important step towards feasible modelling of complex temporal and spatial multi-313 scale 3D problems such as earthquake cycle simulations with near field heterogeneities, nonlinear material behavior and 314 plasticity, as well as a networks of interacting faults, including fault branches and non-planar fault geometry. The major chal-315 lenge of earthquake cycle simulations is that they involve very long interseismic loading time (years) while the dynamic rupture 316 happens extremely rapidly (seconds). An advantage of the hybrid method is that the SBIM is already capable of absorbing elas-317 tic waves in the dynamic as well as in the quasi-dynamic limit and these approaches can be combined in a variable time stepping

326
We developed a three dimensional hybrid method combining the finite-element method with the spectral boundary-integral 327 method. We validated the hybrid method using a benchmark problem and illustrated its potential for solving complex earthquake 328 propagation on various example problems including near systems with field heterogeneity and multiple interacting faults. The 329 hybrid method is suitable for cases where the spatial extent of near field nonlinearity and heterogeneity is too large to be lumped 330 into an effective fault constitutive law, but is still considerably smaller than the domain of interest for the wave propagation.

331
In these cases, the hybrid method allows for a reduction of computational cost by at least one order of magnitude with respect 332 to a full finite-element implementation, while maintaining the same level of accuracy. The high accuracy and computational 333 efficiency of the hybrid method enable the investigation of complex failure problems such as multi-physics fault zone problems.

335
The authors thank Dr. Nicolas Richart for technical help.