While the notion that injecting fluids into the subsurface can reactivate faults by reducing frictional resistance is well established, the ensuing evolution of slip is still poorly understood. What controls whether the induced slip remains stable and confined to the fluid-affected zone or accelerates into a runaway earthquake? Are there observable indicators of the propensity to earthquakes before they happen? Here, we investigate these questions by modeling a unique fluid-injection experiment on a natural fault with laboratory-derived friction laws. We show that a range of fault models with diverging stability with sustained injection reproduce the slip measured during pressurization. Upon depressurization, however, the most unstable scenario departs from the observations, suggesting that the fault is relatively stable. The models could be further distinguished with optimized depressurization tests or spatially distributed monitoring. Our findings indicate that avoiding injection near low-residual-friction faults and depressurizing upon slip acceleration could help prevent large-scale earthquakes.

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.

Fluid-fault interactions result in many two-way coupled processes across a range of length scales, from the micron scale of the shear zone to the kilometer scale of the slip patch. The scale separation and complex coupling render fluid-fault interactions challenging to simulate, yet they are key for our understanding of experimental data and induced seismicity. Here we present spectral boundary-integral solutions for in-plane interface sliding and opening in a poroelastic solid. We solve for fault slip in the presence of rate-and-state frictional properties, inelastic dilatancy, injection, and the coupling of a shear zone and a diffusive poroelastic bulk. The shear localization zone is treated as having a finite width and non-constant pore pressure, albeit with a simplified mathematical representation. The dimension of the 2D plane strain problem is reduced to a 1D problem resulting in increased computational efficiency and incorporation of small-scale shear-zone physics into the boundary conditions. We apply the method to data from a fault injection experiment that has been previously studied with modeling. We explore the influence of bulk poroelastic response, bulk diffusivity in addition to inelastic dilatancy on fault slip during injection. Dilatancy not only alters drastically the stability of fault slip but also the nature of pore pressure evolution on the fault, causing significant deviation from the standard square-root-of-time diffusion. More surprisingly, varying the bulk’s poroelastic response (by using different values of the undrained Poisson’s ratio) and bulk hydraulic diffusivity can be as critical in determining rupture stability as the inelastic dilatancy.

Physics-based numerical modeling of earthquake source processes strives to predict quantities of interest for seismic hazard, such as the probability of an earthquake rupture jumping between fault segments. How to assess the predictive power of numerical models remains a topic of ongoing debate. Here, we investigate how sensitive are the outcomes of numerical simulations of sequences of earthquakes and aseismic slip to choices in numerical discretization and treatment of inertial effects, using a simplified 2-D crustal fault model with two co-planar segments separated by a creeping barrier. Our simulations demonstrate that simplifying inertial effects and using oversized cells significantly affects the resulting earthquake sequences, including the rate of two-segment ruptures. We find that a number of fault models with different properties and modeling assumptions can produce comparable frequency-magnitude statistics and static stress drops but have rates of two-segment ruptures ranging from 0 (single-segment ruptures only) to 1 (two-segment ruptures only). For sufficiently long faults, we find that long-term sequences of events can substantially differ even among simulations that are well-resolved by standard considerations. In such simulations, some outcomes, such as static stress drops, are stable among adequately-resolved simulations, whereas others, such as the rate of two-segment ruptures, can be highly sensitive to numerical procedures and physical assumptions, and hence cannot be reliably inferred. Our results emphasize the need to examine the potential dependence of simulation outcomes on the modeling procedures and resolution, particularly when assessing their predictive value for seismic hazard assessment.

Determining conditions for earthquake slip on faults is a key goal of fault mechanics highly relevant to seismic hazard. Previous studies have demonstrated that enhanced dynamic weakening (EDW) can lead to dynamic rupture of faults with much lower shear stress than required for rupture nucleation. We study the stress conditions before earthquake ruptures of different sizes that spontaneously evolve in numerical simulations of earthquake sequences on rate-and-state faults with EDW due to thermal pressurization of pore fluids. We find that average shear stress right before dynamic rupture (aka shear prestress) systematically varies with the rupture size. The smallest ruptures have prestress comparable to the local shear stress required for nucleation. Larger ruptures weaken the fault more, propagate over increasingly under-stressed areas due to dynamic stress concentration, and result in progressively lower average prestress over the entire rupture. The effect is more significant in fault models with more efficient EDW. We find that, as a result, fault models with more efficient weakening produce fewer small events and result in systematically lower b-values of the frequency-magnitude event distributions. The findings 1) illustrate that large earthquakes can occur on faults that appear not to be critically stressed compared to stresses required for slip nucleation; 2) highlight the importance of finite-fault modeling in relating the local friction behavior determined in the lab to the field scale; and 3) suggest that paucity of small events or seismic quiescence may be the observational indication of mature faults that operate under low shear stress due to EDW.