Portal de Periódicos da CAPES

Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics

2011; American Geophysical Union; Volume: 116; Issue: B10 Linguagem: Inglês

10.1029/2011jb008395

2156-2202

Yoshihiro Kaneko, Jean‐Paul Ampuero, N. Lapusta,

High-pressure geophysics and materials

Journal of Geophysical Research: Solid EarthVolume 116, Issue B10 SeismologyFree Access Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics Y. Kaneko, Y. Kaneko [email protected] Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California, San Diego, La Jolla, California, USASearch for more papers by this authorJ.-P. Ampuero, J.-P. Ampuero Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California, USASearch for more papers by this authorN. Lapusta, N. Lapusta Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California, USA Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California, USASearch for more papers by this author Y. Kaneko, Y. Kaneko [email protected] Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California, San Diego, La Jolla, California, USASearch for more papers by this authorJ.-P. Ampuero, J.-P. Ampuero Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California, USASearch for more papers by this authorN. Lapusta, N. Lapusta Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California, USA Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California, USASearch for more papers by this author First published: 29 October 2011 https://doi.org/10.1029/2011JB008395Citations: 54AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinkedInRedditWechat Abstract [1] We develop a spectral element method for the simulation of long-term histories of spontaneous seismic and aseismic slip on faults subjected to tectonic loading. Our approach reproduces all stages of earthquake cycles: nucleation and propagation of earthquake rupture, postseismic slip and interseismic creep. We apply the developed methodology to study the effects of low-rigidity layers on the dynamics of the earthquake cycle in 2-D. We consider two cases: small (M ∼ 1) earthquakes on a fault surrounded by a damaged fault zone and large (M ∼ 7) earthquakes on a vertical strike-slip fault that cuts through shallow low-rigidity layers. Our results indicate how the source properties of repeating earthquakes are affected by the presence of a damaged fault zone with low rigidity. Compared to faults in homogeneous media, we find (1) reduction in the earthquake nucleation size, (2) amplification of slip rates during dynamic rupture propagation, (3) larger recurrence interval, and (4) smaller amount of aseismic slip. Based on linear stability analysis, we derive a theoretical estimate of the nucleation size as a function of the width and rigidity reduction of the fault zone layer, which is in good agreement with simulated nucleation sizes. We further examine the effects of vertically-stratified layers (e.g., sedimentary basins) on the nature of shallow coseismic slip deficit. Our results suggest that low-rigidity shallow layers alone do not lead to coseismic slip deficit. While the low-rigidity layers result in lower interseismic stress accumulation, they also cause dynamic amplification of slip rates, with the net effect on slip being nearly zero. Key Points Spectral element modeling of spontaneous earthquake cycles is presented Source properties of repeating earthquakes are affected by damaged fault zones Near-surface low-rigidity layers do not lead to the reduction of coseismic slip 1. Introduction [2] Earthquake cycle simulations are important for understanding earthquake mechanics and physics-based hazard analysis. Modeling long-term slip histories of faults remains, however, quite challenging due to a wide range of spatial and temporal scales involved. To simulate a spontaneous earthquake sequence on a fault, a model needs to incorporate and resolve slow tectonic loading during the interseismic periods, nucleation and propagation of rupture during earthquakes that involve rapid changes in stress and slip rate at the propagating dynamic rupture tips, and the subsequent postseismic deformation and aseismic afterslip. In addition, destructive large earthquakes occur on faults that extend tens to hundreds of kilometers while variations in stress changes and slip rate at the rupture tip occur over distances of the order of meters. [3] Many approaches to modeling long-term histories of fault slip have been proposed [e.g., Shibazaki and Matsu'ura, 1992; Rice, 1993; Cochard and Madariaga, 1996; Tullis, 1996; Ward, 1997; Rundle et al., 1999; Kato, 2004; Duan and Oglesby, 2005; Liu and Rice, 2005; Hillers et al., 2006; Dieterich and Richards-Dinger, 2010] but all of them used simplified treatments of either aseismic slip processes (e.g., nucleation, fault creep, and afterslip) or inertial effects during dynamic rupture. Lapusta et al. [2000] and Lapusta and Liu [2009] developed 2-D and 3-D boundary integral methods (BIMs) capable of capturing both seismic and aseismic slip and the gradual process of earthquake nucleation. Those studies are restricted to planar faults embedded in a uniform elastic space or a half-space. At the same time, observations indicate complex crustal structures with variable bulk properties, fault damage zones, and non-planar fault geometries. Hence it is important to include those factors into long-term earthquake cycle models, combining them with laboratory-derived fault constitutive relations. [4] In this work, we develop a spectral element method (SEM) that can enable us to simulate long-term history of spontaneous seismic and aseismic slip on a fault embedded into heterogeneous bulk and subjected to laboratory-derived rate and state friction and slow tectonic loading. In particular, we present a quasi-static SEM with a time updating scheme that can be used to model long-term deformation histories and that is suitable for a fault boundary governed by a rate and state friction formulation. Our model merges the explicit scheme presented by Kaneko et al. [2008] for simulating dynamic rupture propagation with the quasi-static SEM developed in this work for modeling slow tectonic loading and the associated crustal deformation and fault slip. The SEM with an explicit time scheme was built upon prior studies by Komatitsch and Vilotte [1998], Komatitsch and Tromp [1999], and Ampuero [2002]. The combined algorithm is able to resolve all stages of an earthquake cycle, including gradual nucleation processes, dynamic rupture propagation, postseismic slip, and aseismic processes throughout the loading period. Since the presented approach is based on finite element methodology, it can be adapted to simulations of non-planar faults and fault systems. [5] Our methodology is described in sections 2 and 3. Section 3 also illustrates its potential by presenting simulations of long-term slip on a fault segment with relatively simple distributions of fault friction properties. To verify the developed SEM approach, we compare SEM and BIM simulation results in a 2-D model of small (M ∼ 1) repeating earthquakes. We then consider two application examples to explore the potential effects of variable bulk properties on repeating earthquakes (sections 4 and 5). [6] In section 4, we consider small repeating earthquakes on a fault that bisects a meter-scale fault-parallel low-rigidity zone embedded in undamaged or damaged but higher-rigidity host rock. Such a model configuration is motivated by localized damaged zones surrounding fault cores often found on exhumed faults [e.g., Chester et al., 1993]. Recently, meter-thick foliated fault gouge of extremely low P- and S-wave speeds embedded within a 200-m wide damaged zone at 2.6–2.8 km depths was found in the EarthScope's San Andreas Fault Observatory at Depth (SAFOD) project [Zoback et al., 2010]. Simulations of single earthquake rupture suggest that such a mechanical configuration leads to perturbed rupture speeds and slip velocity of propagating rupture, resulting in high-frequency oscillations in the slip function near the rupture front [Harris and Day, 1997]. Here we examine additional effects of a fault-parallel low-rigidity zone on earthquake source properties of simulated repeating earthquakes. We compare the model response in a layered bulk with that in a homogeneous bulk, and investigate how earthquake source properties, such as stress drop, recurrence intervals, and nucleation sizes, depend on the width of a low-rigidity layer. Given the recent successes of rate and state models in explaining several properties of repeating earthquake sequences [Chen and Lapusta, 2009; Chen et al., 2010], it is important to consider these additional effects for proper interpretation of observations. [7] In section 5, we consider the effects of stratified bulk properties (e.g., a sedimentary basin) on simple repeating earthquakes that rupture the entire seismogenic depth. Fialko et al. [2005] pointed out that, for several large (M∼7) strike-slip earthquakes, coseismic slip in the uppermost crust is systematically lower than that at seismogenic depth. A reduction of coseismic slip at shallow depths (<3–4 km), referred to as 'shallow slip deficit', has been inferred for several large strike-slip crustal earthquakes [e.g., Simons et al., 2002; Fialko et al., 2005; Bilham, 2010], including the 1992 M7.3 Landers earthquake, the 1999 M7.1 Hector Mine earthquake, the 2005 M6.5 Bam earthquake, and the 2010 M7.0 Haiti earthquake. Several mechanisms have been proposed to explain the shallow slip deficit, including the presence of velocity-strengthening friction at shallow depths that releases accumulated strain by fault creep [e.g., Marone et al., 1991; Marone, 1998; Kaneko et al., 2008] and low pre-stress in low-rigidity shallow bulk materials resulting from uniform tectonic strain [e.g., Rybicki, 1992; Rybicki and Yamashita, 1998]. Using the developed SEM, we investigate whether a vertical strike-slip fault embedded in a vertically stratified bulk structure can cause shallow coseismic slip deficit without the presence of a shallow velocity-strengthening region. 2. A Quasi-Static SEM Algorithm for Simulations of Long-Term Deformation Histories Discretized Elastodynamic Relations [8] The SEM dynamic model presented by Kaneko et al. [2008] relies on an explicit time updating scheme, the approach commonly used in SEMs for wave propagation [e.g., Komatitsch and Vilotte, 1998; Komatitsch and Tromp, 1999]. However, the explicit time scheme limits the maximum length of each time step Δt by the Courant stability condition. For dynamic rupture simulations, the Courant condition can be rewritten in terms of the cohesive zone size Λ divided by the wave speeds of the medium: where D is the dimension of the problem, C is a stability parameter that depends only on the time scheme and is of order one, and N is the number of fault plane node points within the cohesive zone. (In SEM, the critical time step is actually smaller due to the non-uniform distribution of the Gauss-Lobatto-Legendre nodes inside each spectral element.) N should be at least 3–5 for well-resolved simulations of dynamic rupture in the cases with slip-weakening or weakly rate-dependent friction laws [Day et al., 2005; Kaneko et al., 2008; Kaneko and Lapusta, 2010]. For dynamic rupture simulations with cohesive zone sizes of 1–100 meters, Δt ∼ 10−4 − 10−2 s, and hence simulating tens to thousands of years of deformation histories is not computationally feasible. To take a longer time step, one needs to use an implicit time updating scheme. SEMs with implicit schemes have been used to solve elastic and acoustic wave equations [e.g., Ampuero, 2002; Zampieri and Pavarino, 2006; Dupros et al., 2010]. Here we develop a quasi-static SEM with an adaptive time stepping and merge it with the fully dynamic SEM to simulate long-term slip histories on a rate-and-state fault. [9] As in the work of Kaneko et al. [2008], we start from the discretized weak form of the equation of motion in its matrix form: where M and K are the mass and stiffness matrix respectively, B is the fault-boundary matrix (described in Appendix A), τ = τtot − τo is the relative traction vector on the fault, τtot is the total traction and τo is the traction in the reference static-equilibrium state. Vectors u, , and collect the values of displacements, particle velocities and accelerations, respectively, of all the computational nodes of the bulk mesh. [10] In the case of quasi-static problems, equation (2) becomes Let us decompose the displacement vector u into the values on fault nodes, denoted by uf, and the values on nodes within the medium, denoted by um. Then where K11 and K12 are the parts of the stiffness matrix corresponding to uf, and K21 and K22 are the parts corresponding to um. From equation (5), we have Given the displacement on the fault, uf, this equation yields the corresponding displacement field in the medium, um. [11] In equation (4), let us introduce f ≡ K11uf + K12um. We now write equation (4) for the fault nodes with the ± signs indicating the values of field variables on the two sides of the fault (Figure 1): Subtracting the minus side from the plus side, and using the sign convention τ = −τ+ = τ−, where τ± are defined with respect to the outward normal from the fault boundary Γ± (Figure 1), we obtain Equation (6) allows to eliminate the off-fault degrees of freedom, um, to obtain a formulation (8) involving only the degrees of freedom on the fault (τ and uf). This procedure is known as static condensation or sub-structuring in computational mechanics. Because the B± matrices are diagonal, the expression (8) is a local relation which can be computed node by node on the fault once the f± terms are known or predictor values are assumed. It is convenient to rewrite equation (8) in terms of total traction, τtot = τo + τ: Note that for the cases we consider in this study, the fault-normal components of traction τ remain unchanged, and hence the fault-normal components of τ are neglected in the calculations. Figure 1Open in figure viewerPowerPoint The fault divided into two non-overlapping surfaces Γ±. [12] The developed quasi-static time stepping algorithm is summarized in Appendix B. The algorithm is written in general terms, independent of specific forms of fault constitutive relations. The quasi-static algorithm is similar to Heun's method, which can be seen as an extension of the Euler method into a two-stage second-order Runge-Kutta method. While Heun's method is usually qualified as an explicit method, we still need to solve a large linear system during static condensation. Hence we prefer to qualify our algorithm as implicit. In section 2.3, we discuss the quasi-static time updating scheme coupled with a specific form of friction laws in more detail. Fault Constitutive Response: Rate and State Friction Laws [13] The fault resistance to sliding is described by laboratory-derived rate and state friction laws, which were developed to incorporate observations of rock friction experiments at low slip rate [Dieterich, 1978, 1979; Ruina, 1983; Blanpied et al., 1995, 1998; Marone, 1998]. For time-independent effective normal stress , the shear strength on the fault is expressed as where a and b are rate and state constitutive parameters with magnitudes of the order of 0.01, is the magnitude of slip velocity, f0 is a reference friction coefficient corresponding to a reference slip velocity o, θ is a state variable which is typically interpreted as the average age of the population of contacts between two surfaces, and L is the characteristic slip for state evolution [Dieterich, 1978, 1979; Rice and Ruina, 1983; Ruina, 1983; Dieterich and Kilgore, 1994]. Two types of state-variable evolution laws are commonly used in modeling: [14] The parameter combination a − b < 0 corresponds to steady state velocity-weakening friction and can lead to unstable slip, whereas a − b > 0 corresponds to steady state velocity-strengthening and leads to stable sliding [Rice and Ruina, 1983; Ruina, 1983]. Throughout this article, we omit the words "steady state" and simply refer to velocity weakening/strengthening. [15] In expression (10), shear frictional strength is undefined for slip velocities = 0, which is unphysical. To regularize (10) near = 0, we follow the approach of Rice and Ben-Zion [1996] and Lapusta et al. [2000] in using a thermally activated creep model of the direct effect term a ln (/o) to obtain Solving (13) for gives This regularization produces a negligible change from equation (10) in the range of slip velocities explored by laboratory experiments; the difference in at ∼ o is of the order of exp(−2f0/a) or less, and the typical value of f0/a in this study is 40. [16] Under slow tectonic loading, frictional instability (i.e., an earthquake) is able to develop only if the velocity-weakening region of the fault exceeds the nucleation size h* [Rice and Ruina, 1983; Rice, 1993; Rubin and Ampuero, 2005]. Two theoretical estimates of the earthquake nucleation size for 2-D problems are given by where μ* = μ for mode III and μ* = μ/(1 − ν) for mode II. The estimate h*RR was derived from the linear stability analysis of steady sliding by Rice and Ruina [1983], while h*RA was obtained for the parameter regime a/b > 0.5 by Rubin and Ampuero [2005] on the basis of energy balance for a quasi-statically expanding crack. Note that Rubin and Ampuero [2005] gave formulae for half of the nucleation size but we use full nucleation sizes here. Updating Scheme: Advancing One Evolution Time Step [17] We have developed an updating scheme appropriate for the rate and state fault boundary condition. Here, we discuss how values of field variables are updated over one evolution time step. We adopt a multistage predictor-corrector strategy to solve the statically condensed problem. Suppose that the discretized values of displacement u and particle velocity are known at the nth time step. To find the values of the field variables at the (n + 1)th time step, we perform the following steps. [18] 1. Predict the values of displacements on the fault uf, based on the known values at the nth time step: [19] 2. Solve for the displacement field in the medium u*n+1m using equation (6): This is solved by a preconditioned conjugate gradient method, an iterative method for solving symmetric-positive-definite systems of equations. The algorithm we use is based on Hestenes and Stiefel [1952] and is summarized by Trefethen and Bau [1997]. Because the stiffness matrix K22 is large (∼105 by 105), a direct method such as Gaussian elimination cannot be used. Fortunately, the matrix K22 is sparse, and the product Ku is always computed at a local elemental level as in the case of the dynamic SEM [Kaneko et al., 2008]. This is why we use an iterative method. [20] 3. Compute f* = K11u*n+1f + K12u*n+1m and τ*n+1tot in equation (9): [21] 4. Determine the first prediction of the state variable, *n+1. By integrating the evolution law (11) or (12) with the constant magnitude n of slip velocity n = un+ − un− during the time step, we obtain for the aging law, and for the slip law. [22] 5. Find the first prediction of slip velocity *n+1 by equating the magnitude of shear stress in equation (19) and strength in equation (13). The directions of shear traction vector τn+1tot and slip velocity vector n+1 have to coincide. From equation (19), the traction τ*n+1tot and f*n+1 have the same direction. By projecting (19) onto that direction and using equation (14), we obtain Using the directional cosines constructed from the components of τ*n+1tot, we obtain the components of *n+1. [23] 6. Calculate the final prediction of displacement and slip on the fault, u**n+1f = at the (n + 1)th time step by [24] 7. Make the corresponding prediction u**n+1m of the displacement in the medium using the u**n+1f as in step 2. This involves another conjugate gradient solution, with u*n+1m as an initial guess. [25] 8. Make the corresponding prediction τn+1tot** and **n+1 by repeating steps 3 and 4 and by replacing n in equation (20) or (21) with (n + *n+1)/2. [26] 9. Find the final prediction **n+1 and the components of **n+1 by repeating step 5 with τn+1tot** and θ**n+1 instead of τn+1tot* and θ*n+1. [27] 10. Declare the values of n+1, θn+1, and τn+1tot on the fault, and the values of displacement of the entire medium un+1, to be equal to the predictions with the superscript double asterisks. [28] In steps 2 and 7, we terminate the conjugate gradient iteration when ∥LHS − RHS∥/∥RHS∥ < ɛCG = 10−5, where LHS and RHS are the left-hand side and right-hand side of equation (18), respectively. The convergence rate in the conjugate gradient iteration greatly depends on the type of preconditioner. Currently, we use the Jacobi preconditioner, one of the simplest forms of preconditioning, for the conjugate gradient method. For SEMs, a special preconditioner called 'the Schwarz method' has been shown to significantly increase the convergence rate for elastic problems [Lottes and Fischer, 2004; Zampieri and Pavarino, 2006]. Implementing the Schwarz preconditioner remains a subject of future work. 3. Implementation Example Formulation of a 2-D Model [29] The response of faults to tectonic loading is characterized by long periods of quasi-static deformation combined with short periods of fast dynamic slip. To simulate such response, we adopt the variable time stepping procedure of Lapusta et al. [2000], in which the time step is set to be inversely proportional to slip velocity on the fault as described in Appendix C. As a result, relatively large time steps, a fraction of a year, are used in the interseismic periods, while small time steps, a fraction of a second or smaller, are used to simulate fast seismic slip. Note that the efficiency of the time stepping procedure (Appendix C) depends on the degree of the positive direct effect in the rate and state formulation [Lapusta et al., 2000], a feature that has ample laboratory confirmation. [30] The updating scheme introduced in the previous section can be merged with the explicit time stepping scheme of the fully dynamic SEM problem (Appendix B). The main challenge is to find proper criteria for switching from the quasi-static implicit scheme to the dynamic explicit scheme and vice versa. At the onset of earthquakes (or at the end of nucleation processes), slip velocities abruptly increase from values much smaller than typical plate loading rates (∼10−10–10−9 m/s) to coseismic values (∼1–10 m/s), and the time step progressively becomes smaller. Hence we switch from one scheme to the other based on the values of the maximum slip velocity. For the problems discussed below, we switch from the quasi-static to dynamic scheme at maxQD = 0.5 mm/s and from the dynamic to quasi-static scheme at maxDQ = 0.2 mm/s. Ideally, one could formulate a switching criterion based on the relative importance of the inertial term in the governing equation (2). However, we find that such quantity is numerically more unstable than the criteria based on the values of the maximum slip velocity (Appendix D). We confirm that the criteria we adopt here ensures that, at the times of the switch, the inertial term in the governing equation is much smaller (∼10−6) relative to the other terms (Appendix D) and that the results compare well with BIM methods as discussed below. [31] To demonstrate how the ideas outlined so far are combined to produce long-term deformation histories, let us consider the response of a 2-D model of a vertical strike-slip fault embedded in an elastic medium (Figure 2). On the fault, a potentially seismogenic patch borders regions steadily moving with the prescribed slip rate Vpl = 2 mm/yr, as illustrated in Figure 2. That steady motion provides loading. The fault motion is in the along-strike direction y, but only variations with depth z are considered, so that the fault behavior is described by strike-parallel slip δ(z, t), slip velocity (or slip rate) (z, t) = ∂δ(z, t)/∂t, and the relevant component of shear stress τtot(z, t). The symmetries of the problem allow us to restrict the computational domain to the medium on one side of the fault (x ≥ 0). Figure 2Open in figure viewerPowerPoint 2-D models of a vertical strike-slip fault. Small repeating earthquakes at seismogenic depths in a region indicated by a black rectangle are simulated using these models. The fault is 90-m long. Slip evolution is computed, based on the assumed friction law, on the 45-m long central portion, and the prescribed slip rate Vpl = 2 mm/yr is imposed on the two 22.5-m long outer portions. The fault is divided into three segments: a 19-m long central velocity-weakening (VW) patch surrounded by two 13-m long velocity-strengthening (VS) regions. By symmetry consideration, the medium across the fault boundary has equal and opposite motion. Unless otherwise noted, the medium is assumed to be homogeneous. In section 4, we investigate the effect of a fault-parallel low-rigidity layer of a width H and rigidity μD on seismic and aseismic slip. [32] It is convenient to express the formulae in terms of variables (u(x, z, t) − Vplt/2) and ((x, z, t) − Vpl/2), in which case τo (x, z, t) becomes independent of time and equal to the initial stress τo (x, z). This approach was used for the BIM model of Lapusta et al. [2000]. For the 2-D problems we consider here, the medium across the fault boundary has equal and opposite motion by symmetry consideration. Then the relation (9) on the fault becomes Note that our mesh is conformal and hence B = B+ = B−. [33] The SEM model consists of a 90 m by 60 m rectangular domain (Figure 2). To allow comparison to the BIM, the domain is replicated using periodic boundary conditions on both sides of the domain (Figure 2). The fault boundary obeys rate and state friction with the aging law (11). The model contains variations in steady state friction properties that create rheological transitions (Figure 3). The parameters used in the simulations are listed in Table 1. The effective normal stress and characteristic slip L are uniform along the fault. Figure 3Open in figure viewerPowerPoint (a) Depth-variable distribution of friction parameters a and (a − b). (b) Distribution of the ratio h*/Δx. Two theoretical estimates h* of the nucleation size by Rubin-Ampuero (RA) and Rice-Ruina (RR) are shown. Table 1. Parameters Used in the 2-D SEM and 2-D BIM Models of Small Repeating Earthquakes Parameter Symbol Value Shear modulus μ 32.0 GPa Shear wave speed Vs 3.464 km/s Reference slip rate 0 10−6 m/s Reference friction coefficient f0 0.6 Characteristic slip distance L 84.0 micron Effective normal stress 120 MPa Rate and state parameter a a 0.0144 Rate and state parameter b b 0.0191a a The indicated value of b is for the velocity-weakening region in Figure 2. [34] We use the criteria for spatial discretizations developed in the work by Perfettini and Ampuero [2008] and Lapusta and Liu [2009], which showed that resolving a cohesive zone size is a more stringent requirement than resolving the nucleation size, for the aging formulation of rate and state friction and typical rate and state parameters. The spatial discretization required to resolve dynamic rupture is largely sufficient for properly resolving the interseismic deformation. An optimal mesh for the quasi-static problem would involve mesh coarsening as a function of distance from the fault. Here we do not consider such optimization and, to avoid mapping between different meshes, we adopt the same mesh for both the dynamic and quasi-static problems. We discretize the domain into quadrilateral spectral elements with an average node spacing Δx = 0.25 m in each element. Note that for SEMs, the nodes are generally nonuniformly distributed, and hence we report the node spacing in terms of their average value. The average node spacing Δx = 0.25 m results in Λo/Δx ≈ μL/(bΔx) ≈ 5 where Λo is the cohesive zone size at the rupture speed Vr → 0+. Such resolution has shown to be adequate in the work of Day et al. [2005] and Lapusta and Liu [2009], and it leads to stable results in our simulations that do not change due to finer discretizations. The selected spatial discretization corresponds to h*RA/Δx ≈ 50 (Figure 3b), where h*RA is the estimate of the nucleation size obtained by Rubin and Ampuero [2005] for a/b ⪆0.5, see equation (16). Comparison of Simulation Results Obtained With 2-D SEM and 2-D BIM [35] To assess the accuracy of numerical results, we conduct comparison of simulation results obtained using the developed SEM model with those of the BIM spectral formulation of Lapusta et al. [2000], which resolves all stages of each earthquake episode under a single computational scheme. Figure 2 illustrates the geometry of the antipl