De Novellis et al. (2020, hereafter DN20) studied the effect of mass extraction from a quarry on the occurrence of the Mw 4.9 Le Teil, France, earthquake of November 11 2019. This topic was also the focus of the report of the French working group mandated by CNRS INSU (“Groupe de Travail Teil” of Institut National des Sciences de l’Univers du Centre National de la Recherche Scientifique; Delouis et al, 2019). Despite using similar data and methods, these two independent research efforts reached contrasting conclusions. While both concluded the earthquake was possibly a triggered event (i.e. its initiation was possibly promoted by the quarry activity but its further rupture growth was primarily enabled by natural pre-existing stresses), DN20 deemed realistic the hypothesis that the earthquake was an induced event (i.e. both the earthquake initiation and its further growth, up to its final size, were caused by the quarry activity). This distinction is critical for our understanding of future anthropogenic hazards in the region and in similar settings elsewhere, and may have significant social, economical and legal repercussions. Here, we show that a severe error in the calculations carried by DN20 undermines their conclusion.

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.

The subduction zone of the Cocos Plate beneath Southern Mexico has major variations in the megathrust geometry and behavior. The subduction segment beneath the Oaxaca state of Mexico has relatively frequent large earthquakes on the shallow part of the megathrust and within the subducting slab, and it also has large aseismic slow-slip events. The slab geometry under Oaxaca includes part of the subhorizontal “flat-slab” zone extending far from the trench beneath southern Mexico and the beginning of its transition to more regular subduction geometry to the southeast. We study the rupture of the 16 February 2018 Mw 7.2 Pinotepa earthquake near Pinotepa Nacional in Oaxaca that was a thrust event on the subduction interface. The Pinotepa earthquake was about 350 km away from the 8 September 2017 Mw 8.2 Tehuantepec earthquake in the subducting slab offshore Oaxaca and Chiapas; it was in an area of Coulomb stress decrease from the M8.2 quake, so it seems unlikely to be a regular aftershock and was not triggered by the static stress change. Geodetic measurements from interferometric analysis of synthetic aperture radar (InSAR) and time-series analysis of GPS station data constrain finite-fault slip models for the M7.2 Pinotepa earthquake. We analyzed InSAR data from Copernicus Sentinel-1A and -1B satellites and JAXA ALOS-2 satellite. Our Bayesian (AlTar) static slip model for the Pinotepa earthquake shows all of the slip confined to a very small (10-20 km diameter) rupture, similar to some early seismic waveform fits. The Pinotepa earthquake ruptured a portion of the Cocos megathrust that has been previously mapped as partially coupled and shows that at least small asperities in that zone of the subduction interface are fully coupled and fail in high-stress drop earthquakes. The previous 2012 Mw 7.4 Ometepec earthquake is another example of asperity in the partially coupled zone but was not imaged by InSAR so the rupture extent is not so well constrained. The preliminary NEIC epicenter for the Pinotepa earthquake was about 40 km away (NE) from the rupture imaged by InSAR, but the NEIC updated epicenter and Mexican SSN location are closer. Preliminary analysis of GPS data after the Pinotepa earthquake indicates rapid afterslip on the megathrust in the region of coseismic slip. Atmospheric noise masks the postseismic signal on early InSAR data.

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.

Anthropogenic fluid injections at depth induce seismicity which is generally organized as swarms, clustered in time and space, with moderate magnitudes. While some similarities between swarms have already been observed, whether they are driven by the same mechanism is still an open question. Pore fluid pressure or aseismic processes are often proposed to explain observations, while recent models suggest that seismicity is triggered by fluid-induced aseismic slip. Using 22 natural and anthropogenic swarms, we observe that duration, migration velocity and total moment scale similarly for all swarms. This confirms a common driving process for natural and induced swarms and highlights the ubiquity of aseismic slip. We propose a method to estimate the seismic-to-total moment ratio, which is then compared to a theoretical estimation that depends on the migration velocity, the effective stress drop and the slip velocity. Our findings lead to a generic explanation of swarms driving process.

Near-field ground motion is the major blind spot of seismic hazard studies, mainly because of the challenges in accounting for source effects. Initial stress heterogeneity is an important component of physics-based approaches to ground motion prediction that represent source effects through dynamic earthquake rupture modeling. We hypothesize that stress heterogeneity on a fault primarily originates from past background seismicity. We develop a new method to generate stochastic stress distributions as a superposition of residual stresses left by previous ruptures that are consistent with regional distributions of earthquake size and hypocentral depth. We validate our method on Mw 7 earthquake models suitable for California, by obtaining a satisfactory agreement with empirical earthquake scaling laws and ground motion prediction equations. To avoid the excessive seismic radiation produced by dynamic models with abrupt arrest at preset rupture borders, we achieve spontaneous rupture arrest by incorporating a scale-dependent fracture energy adjusted with fracture mechanics theory. Our analyses of rupture and ground motion reveal particular signatures of the initial stress heterogeneity: rupture can locally propagate at supershear speed near the highly-stressed areas; the position of high-stress and low-stress areas due to initial stress heterogeneity determines how the peak ground motion amplitudes and polarization spatially vary along the fault, as low-stress areas slows down the rupture, decrease stress drop, and change the radiation distribution before the rupture arrest. We also find that the medium stratification amplifies the moment rate spectrum at frequencies above 2 Hz, which requires understanding the interaction between site effects and rupture dynamics; therefore, we highlight the need to consider a realistic fault medium on future studies of rupture dynamics. Our approach advances our understanding of the relations between dynamic features of earthquake ruptures and the statistics of regional seismicity, and our capability to model source effects for near-field ground motion prediction studies.

Some observations of repeating earthquakes show an unusual, non self-similar scaling between seismic moment and corner frequency, a source property related to rupture size. These observations have been mostly reported in regions at the transition from stable to unstable slip, in geothermal reservoirs and subduction zones. What controls the non self-similarity of these ruptures and how this is linked to the frictional stability of the interface are still open questions. Here we develop seismic cycle simulations of a single unstable slipping patch to investigate the mechanics underlying this behavior. We show that temporal changes of normal stress on a fault can produce ruptures that exhibit the observed anomalous scaling. Our results highlight the role of fault zone fluid pressure in modulating the effective normal stress and contributing to the sliding stability of the fault.

Supershear earthquakes are rare compared to their subshear counterparts, but the cause for their paucity remains to be understood. We investigate for the first time the prevalence of supershear ruptures across multiple earthquake cycles on long faults using rate-and-state friction and a 2.5D approximation that accounts for the finite seismogenic width W. We find supershear events occur only in a narrow range of friction parameters that is not commonly observed in laboratory experiments, which may explain its rarity in nature. Particularly, the ratio between direct and evolution effects of rate-and-state friction needs to be low (a/b<0.4) and the nucleation length has to be sufficiently large compared to W, but not too large that it causes a transition from seismic to aseismic slip. Our numerical and analytical developments contribute fundamentally to understanding the state of stress on long faults over multiple earthquake cycles and their potential for hosting supershear earthquakes.