The spontaneously propagating shear crack on a frictional interface has proven to be a useful idealization of a natural earthquake. The corresponding boundary value problems are nonlinear and usually require computationally intensive numerical methods for their solution. Assessing the convergence and accuracy of the numerical methods is challenging, as we lack appropriate analytical solutions for comparison. As a complement to other methods of assessment, we compare solutions obtained by two independent numerical methods, a finite difference method and a boundary integral (BI) method. The finite difference implementation, called DFM, uses a traction‐at‐split‐node formulation of the fault discontinuity. The BI implementation employs spectral representation of the stress transfer functional. The three‐dimensional (3‐D) test problem involves spontaneous rupture spreading on a planar interface governed by linear slip‐weakening friction that essentially defines a cohesive law. To get a priori understanding of the spatial resolution that would be required in this and similar problems, we review and combine some simple estimates of the cohesive zone sizes which correspond quite well to the sizes observed in simulations. We have assessed agreement between the methods in terms of the RMS differences in rupture time, final slip, and peak slip rate and related these to median and minimum measures of the cohesive zone resolution observed in the numerical solutions. The BI and DFM methods give virtually indistinguishable solutions to the 3‐D spontaneous rupture test problem when their grid spacing Δxis small enough so that the solutions adequately resolve the cohesive zone, with at least three points for BI and at least five node points for DFM. Furthermore, grid‐dependent differences in the results, for each of the two methods taken separately, decay as a power law in Δx, with the same convergence rate for each method, the calculations apparently converging to a common, grid interval invariant solution. This result provides strong evidence for the accuracy of both methods. In addition, the specific solution presented here, by virtue of being demonstrably grid‐independent and consistent between two very different numerical methods, may prove useful for testing new numerical methods for spontaneous rupture problems.

Numerical simulations of earthquake rupture dynamics are now common, yet it has been difficult to test the validity of these simulations because there have been few field observations and no analytic solutions with which to compare the results. This paper describes the Southern California Earthquake Center/U.S. Geological Survey (SCEC/USGS) Dynamic Earthquake Rupture Code Verification Exercise, where codes that simulate spontaneous rupture dynamics in three dimensions are evaluated and the results produced by these codes are compared using Web-based tools. This is the first time that a broad and rigorous examination of numerous spontaneous rupture codes has been performed—a significant advance...

Large dynamic stresses near earthquake rupture fronts may induce an inelastic response of the surrounding materials, leading to increased energy absorption that may affect dynamic rupture. We systematically investigate the effects of off‐fault plastic energy dissipation in 2‐D in‐plane dynamic rupture simulations under velocity‐and‐state‐dependent friction with severe weakening at high slip velocity. We find that plasticity does not alter the nature of the transitions between different rupture styles (decaying versus growing, pulse‐like versus crack‐like, and subshear versus supershear ruptures) but increases their required background stress and nucleation size. We systematically quantify the effect of amplitude and orientation of background shear stresses on the asymptotic properties of self‐similar pulse‐like ruptures: peak slip rate, rupture speed, healing front speed, slip gradient, and the relative contribution of plastic strain to seismic moment. Peak slip velocity and rupture speed remain bounded. From fracture mechanics arguments, we derive a nonlinear relation between their limiting values, appropriate also for crack‐like and supershear ruptures. At low background stress, plasticity turns self‐similar pulses into steady state pulses, for which plastic strain contributes significantly to the seismic moment. We find that the closeness to failure of the background stress state is an adequate predictor of rupture speed for relatively slow events. Our proposed relations between state of stress and earthquake source properties in the presence of off‐fault plasticity may contribute to the improved interpretation of earthquake observations and to pseudodynamic source modeling for ground motion prediction.

The underestimation of the size of recent megathrust earthquakes illustrates our limited understanding of their spatiotemporal occurrence and governing physics. To unravel their relation to associated subduction dynamics and long‐term deformation, we developed a 2‐D continuum viscoelastoplastic model that uses an Eulerian‐Lagrangian finite difference framework with similar on‐ and off‐fault physics. We extend the validation of this numerical tool to a realistic subduction zone setting that resembles Southern Chile. The resulting quasi‐periodic pattern of quasi‐characteristic M8–M9 megathrust events compares quantitatively with observed recurrence and earthquake source parameters, albeit at very slow coseismic speeds. Without any data fitting, surface displacements agree with GPS data recorded before and during the 2010 M8.8 Maule earthquake, including the presence of a second‐order flexural bulge. These surface displacements show cycle‐to‐cycle variations of slip deficits, which overall accommodate ∼5% of permanent internal shortening. We find that thermally (and stress) driven creep governs a spontaneous conditionally stable downdip transition zone between temperatures of ∼350°C and ∼450°C. Ruptures initiate above it (and below the forearc Moho), propagate within it, interspersed by small intermittent events, and arrest below it as ductile shearing relaxes stresses. Ruptures typically propagate upward along lithological boundaries and widen as pressures drop. The main thrust is constrained to be weak due to fluid‐induced weakening required to sustain regular subduction and to generate events with natural characteristics (fluid pressures of ∼75–99% of solid pressures). The agreement with a range of seismological, geodetic, and geological observations demonstrates the validity and strength of this physically consistent seismo‐thermo‐mechanical approach.

An important goal of computational seismology is to simulate dynamic earthquake rupture and strong ground motion in realistic models that include crustal heterogeneities and complex fault geometries. To accomplish this, we incorporate dynamic rupture modelling capabilities in a spectral element solver on unstructured meshes, the 3-D open source code SPECFEM3D, and employ state-of-the-art software for the generation of unstructured meshes of hexahedral elements. These tools provide high flexibility in representing fault systems with complex geometries, including faults with branches and non-planar faults. The domain size is extended with progressive mesh coarsening to maintain an accurate resolution of the static field. Our implementation of dynamic rupture does not affect the parallel scalability of the code. We verify our implementation by comparing our results to those of two finite element codes on benchmark problems including branched faults. Finally, we present a preliminary dynamic rupture model of the 2011 M_w 9.0 Tohoku earthquake including a non-planar plate interface with heterogeneous frictional properties and initial stresses. Our simulation reproduces qualitatively the depth-dependent frequency content of the source and the large slip close to the trench observed for this earthquake.