In recent laboratory experiments, varying nucleation locations of accelerating slip with changing nucleation lengths were observed. Spatial variations in effective normal stress, due to the controlling influence on fault strength and fracture energy, play an important role. We quantitatively explain how spatially heterogeneous effective normal stresses affect earthquake nucleation and slip behavior. We simulate a meter-scale laboratory experiment in a numerical earthquake sequence model with stochastically variable normal stresses. We identify five regimes of earthquake nucleation and slip behaviors, controlled by the ratio of the heterogeneity wavelength (\(\lambda\)) to the nucleation length (\(L_c\)). When \(\lambda\) is much smaller than \(L_c\), full ruptures are observed. Slip rates and recurrence intervals are similar to those on homogeneous faults with comparable averaged normal stress. When \(\lambda\) is much larger than \(L_c\), slow slip events and partial ruptures occur frequently and the nucleation length of each earthquake depends on the local stress level. We find locations of nucleation and arrest in both low and high normal stress regions (LSR and HSR, respectively) when \(\lambda\) and \(L_c\) are of the same magnitude. When \(\lambda\) is larger than \(L_c\), earthquakes nucleate in LSRs, and arrest in HSRs. However, HSRs and LSRs switch these roles when \(\lambda\) is smaller than \(L_c\). Interestingly, we observe that nucleation migrates from an LSR to its neighboring HSR in one earthquake, when \(\lambda\) is between the minimum and maximum local nucleation lengths. We observe a large amount of aseismic slip and associated stress drop in the initial LSR, which might be linked to the migration of foreshocks as documented in natural and laboratory observations. This improved understanding of earthquake nucleation is important in estimating the seismic potential of different fault patches for natural and induced seismicity.

Numerical simulations of Sequences of Earthquakes and Aseismic Slip (SEAS) have rapidly progressed to address fundamental problems in fault mechanics and provide self-consistent, physics-based frameworks to interpret and predict geophysical observations across spatial and temporal scales. To advance SEAS simulations with rigor and reproducibility, we pursue community efforts to verify numerical codes in an expanding suite of benchmarks. Here we present code comparison results from a new set of benchmark problems BP6-QD-A/S/C that consider a single aseismic slip transient induced by changes in pore fluid pressure consistent with fluid injection and diffusion in fault models with different treatments of fault friction. Ten modeling groups participated in problems BP6-QD-A and BP6-QD-S considering rate-and-state fault models using the aging and slip law formulations for frictional state evolution, respectively, allowing us to explore these ingredients across multiple codes and better understand how various computational factors affect the simulated evolution of pore pressure and aseismic slip. Comparisons of problems using the aging versus slip law illustrate how models of aseismic slip can differ in the timing and amount of slip achieved with different treatments of fault friction given the same perturbations in pore fluid pressure. We achieve excellent quantitative agreement across participating codes, with further agreement being found by ensuring sufficiently fine time-stepping and consistent treatment of remote boundary conditions. Our benchmark efforts offer a community-based example to reveal sensitivities of numerical modeling results, which is essential for advancing multi-physics SEAS models to better understand and construct reliable predictive models of fault dynamics.

High-resolution computer simulations of earthquake sequences in three or even two dimensions pose great demands on time and energy, making lower-cost simplifications a competitive alternative. We systematically study the advantages and limitations of simplifications that eliminate spatial dimensions, from 3D down to 0/1D in quasi-dynamic earthquake sequence models. We demonstrate that, when 2D or 3D models produce quasi-periodic characteristic earthquakes, their behavior is qualitatively similar to lower-dimension models. Certain coseismic characteristics like stress drop and fracture energy are largely controlled by frictional parameters and are thus largely comparable. However, other observations are quantitatively clearly affected by dimension reduction. We find corresponding increases in recurrence interval, coseismic slip, peak slip velocity, and rupture speed. These changes are to a large extend explained by the elimination of velocity-strengthening patches that transmit tectonic loading onto the velocity-weakening fault patch, thereby reducing the interseismic stress rate and enhancing the slip deficit. This explanation is supported by a concise theoretical framework, which explains some of these findings quantitatively and effectively estimates recurrence interval and slip. Through accounting for an equivalent stressing rate at the nucleation size h* into 2/3D models, 0/1D models can also effectively estimate these earthquake cycle parameters. Given the computational efficiency of lower-dimensional models that run more than a million times faster, this paper aims to provide qualitative and quantitative guidance on economical model design and interpretation of modeling studies.

Dynamic modeling of sequences of earthquakes and aseismic slip (SEAS) provides a self-consistent, physics-based framework to connect, interpret, and predict diverse geophysical observations across spatial and temporal scales. Amid growing applications of SEAS models, numerical code verification is essential to ensure reliable simulation results but is often infeasible due to the lack of analytical solutions. Here, we develop two benchmarks for three-dimensional (3D) SEAS problems to compare and verify numerical codes based on boundary-element, finite-element, and finite-difference methods, in a community initiative. Our benchmarks consider a planar vertical strike-slip fault obeying a rate- and state-dependent friction law, in a 3D homogeneous, linear elastic whole-space or half-space, where spontaneous earthquakes and slow slip arise due to tectonic-like loading. We use a suite of quasi-dynamic simulations from 10 modeling groups to assess the agreement during all phases of multiple seismic cycles. We find excellent quantitative agreement among simulated outputs for sufficiently large model domains and small grid spacings. However, discrepancies in rupture fronts of the initial event are influenced by the free surface and various computational factors. The recurrence intervals and nucleation phase of later earthquakes are particularly sensitive to numerical resolution and domain-size-dependent loading. Despite such variability, key properties of individual earthquakes, including rupture style, duration, total slip, peak slip rate, and stress drop, are comparable among even marginally resolved simulations. Our benchmark efforts offer a community-based example to improve numerical simulations and reveal sensitivities of model observables, which are important for advancing SEAS models to better understand earthquake system dynamics.

Geodetic data provide an opportunity to improve our understanding of the processes and parameters controlling the dynamics of deformation during the earthquake cycle at subduction zones. However, the observations contain noise and are temporally and spatially sparse, whereas dynamical models are unequivocally imperfect. Also, the relative contributions from various drivers of surface deformation are poorly constrained by independent observations. Some drivers may be static or vary slowly in time (e.g., plate motion), whilst others vary significantly during the earthquake cycle (e.g., viscoelastic relaxation). Data assimilation combines prior estimates of dynamical models, with the likelihood of observations into posterior estimates of the state evolution and time-independent parameters of a physical process. We explore the usefulness of data assimilation by using a particle filter to estimate the (spatially variable) elastic thickness of the overriding plate and the extent of the locked zone. We assimilate vertical interseismic surface displacements into a 2D elastic flexural model. The particle filter uses a Monte Carlo approach to represent the state probability distribution by a finite number of realizations (“particles”). We use sequential importance resampling to preserve particles statistically close to the observations and duplicate and perturb them. Synthetic experiments demonstrate that the particle filter effectively estimates 1D elastic thickness from synthetic observations. However, elastic thickness estimates for models with a landward increase in plate thickness show larger uncertainty near the coast as the sensitivity of surface displacements reduces with increasing plate thickness. Interestingly, the effectiveness of the elastic thickness estimation is highly sensitive to network aperture, including GPS/A. Assimilation of interseismic vertical velocities prior to the 2011 Tohoku-Oki earthquake yields estimates of upper plate thickness that agree with previous studies. However, results of the locked zone extent are not as expected, which could be due to missing physics in the relatively conceptual model. These results demonstrate the potential of the particle filter to better understand the geodynamic process parameters of the earthquake cycle at subduction zones.