We present an efficient numerical method for earthquake cycle simulations that employs a finite difference discretization of the off‐fault material to accommodate spatially variable elastic properties. The method is developed for the two‐dimensional antiplane shear problem of a vertical strike‐slip fault with rate‐and‐state friction. We compare earthquake cycles in a homogeneous half‐space with those in which the upper portion of the fault cuts through a sedimentary basin. In both cases, we assume velocity‐weakening behavior over the full seismogenic depth, even in the basin, to isolate the influence of elastic heterogeneity. In a homogeneous half‐space, events rupturing the full seismogenic depth occur periodically. Event sequences are more complex in basin models, with one or several subbasin events confined to the lower section of the fault followed by a much larger, surface‐rupturing event that breaks through the basin. This phenomenology emerges only for sufficiently compliant and deep basins. Predicted surface velocities are essentially identical before subbasin events and surface‐rupturing events, suggesting that geodetic observations would not be useful in predicting the rupture mode. The alternating sequence of subbasin and surface‐rupturing events would complicate interpretation of paleoseismic data. Our results also offer one potential explanation for the shallow slip deficit observed in many recent earthquakes, namely, that these events, which lack appreciable surface slip, are simply one style of rupture. Subsequent events on these faults might be larger, with slip extending all the way to the surface. The 1940 Mw 7.0 and 1979 Mw 6.5 Imperial Valley events might be considered as examples of these two rupture styles.

Abstract.Conditions under which a single oscillator model coupled with Dieterich-Ruina's rate and state dependent friction exhibits chaotic dynamics is studied. Properties of spring-block models are discussed. The parameter values of the system are explored and the corresponding numerical solutions presented. Bifurcation analysis is performed to determine the bifurcations and stability of stationary solutions and we find that the system undergoes a Hopf bifurcation to a periodic orbit. This periodic orbit then undergoes a period doubling cascade into a strange attractor, recognized as broadband noise in the power spectrum. The implications for earthquakes are discussed.

S U M M A R YWe investigate the emergent dynamics when the slip law formulation of the non-linear rateand-state friction law is attached to a Burridge-Knopoff spring-block model. We derive both the discrete equations and the continuum formulation governing the system in this framework. The discrete system (ODEs) exhibits both periodic and chaotic motion, where the system's transition to chaos is size-dependent, that is, how many blocks are considered. From the discrete model we derive the non-linear elastic wave equation by taking the continuum limit. This results in a non-linear partial differential equation (PDE) and we find that chaos ensues when the same parameter is increased. This critical parameter value needed for the onset of chaos in the continuous model is much smaller than the value needed in the case of a single block and we discuss the implications this has on dynamic modelling of earthquake rupture with this specific friction law. Most importantly, these results suggest that the friction law is scale-dependent, thus caution should be taken when attaching a friction law derived at laboratory scales to full-scale earthquake rupture models. Furthermore, we find solutions where the initial slip pulse propagates like a travelling wave, or remains localized in space, suggesting the presence of soliton and breather solutions. We discuss the significance of these pulse-like solutions and how they can be understood as a proxy for the propagation of the rupture front across the fault surface during an earthquake. We compute analytically the conditions for soliton solutions and by exploring the resulting parameter space, we introduce a possible method for determining a range of suitable parameter values to be used in future dynamic earthquake modelling.

We have developed an efficient computational framework for simulating multiple earthquake cycles with off-fault plasticity. The method is developed for the classical antiplane problem of a vertical strike-slip fault governed by rate-and-state friction, with inertial effects captured through the radiationdamping approximation. Both rate-independent plasticity and viscoplasticity are considered, where stresses are constrained by a Drucker-Prager yield condition. The off-fault volume is discretized using finite differences and tectonic loading is imposed by displacing the remote side boundaries at a constant rate. Time-stepping combines an adaptive Runge-Kutta method with an incremental solution process which makes use of an elastoplastic tangent stiffness tensor and the return-mapping algorithm. Solutions are verified by convergence tests and comparison to a finite element solution. We quantify how viscosity, isotropic hardening, and cohesion affect the magnitude and off-fault extent of plastic strain that develops over many ruptures. If hardening is included, plastic strain saturates after the first event and the response during subsequent ruptures is effectively elastic. For viscoplasticity without hardening, however, successive ruptures continue to generate additional plastic strain. In all cases, coseismic slip in the shallow sub-surface is diminished compared to slip accumulated at depth during interseismic loading. The evolution of this slip deficit with each subsequent event, however, is dictated by the plasticity model. Integration of the off-fault plastic strain from the viscoplastic model reveals that a significant amount of tectonic off

Numerical simulations of Sequences of Earthquakes and Aseismic Slip (SEAS) have made great progress over the past decades to address important questions in earthquake physics and fault mechanics. However, significant challenges in SEAS modeling remain in resolving multiscale interactions between aseismic fault slip, earthquake nucleation, and dynamic rupture; and understanding physical factors controlling observables such as seismicity and ground deformation. The increasing capability and complexity of SEAS modeling calls for extensive efforts to verify codes and advance these simulations with rigor, reproducibility, and broadened impact. In 2018, we initiated a community code-verification exercise for SEAS simulations, supported by the Southern California Earthquake Center (SCEC). Here we report the findings from our first two benchmark problems (BP1 and BP2), designed to test the capabilities of different computational methods in correctly solving a mathematically well-defined, basic problem in crustal faulting. These benchmarks are for a 2D antiplane problem, with a 1D planar vertical strike-slip fault obeying rate-and-state friction, embedded in a 2D homogeneous, linear elastic half-space. Sequences of quasi-dynamic earthquakes with periodic occurrences (BP1) or bimodal sizes (BP2) and their interactions with aseismic slip are simulated. The comparison of >70 simulation results from 11 groups using different numerical methods, uploaded to our online platform, show excellent agreements in long-term and coseismic evolution of fault properties. In BP1, we found that the truncated domain boundaries influence interseismic fault stressing, earthquake recurrence, and coseismic rupture process, and that agreement between models is only achieved with sufficiently large domain sizes. In BP2, we found that complexity of long-term fault behavior depends on how well important physical length scales related to spontaneous nucleation and rupture propagation are resolved. Poor numerical resolution can result in the generation of artificial complexity, impacting simulation results that are of potential interest for characterizing seismic hazard, such as earthquake size distributions, moment release, and earthquake recurrence times. These results inform the development of more advanced SEAS models, contributing to our further understanding of earthquake system dynamics.