Fault roughness controls injection-induced seismicity.

Surface roughness ubiquitously prevails in natural faults across various length scales. Despite extensive studies highlighting the important role of fault geometry in the dynamics of tectonic earthquakes, whether and how fault roughness affects fluid-induced seismicity remains elusive. Here, we investigate the effects of fault geometry and stress heterogeneity on fluid-induced fault slip and associated seismicity characteristics using laboratory experiments and numerical modeling. We perform fluid injection experiments on quartz-rich sandstone samples containing either a smooth or a rough fault. We find that geometrical roughness slows down injection-induced fault slip and reduces macroscopic slip velocities and fault slip-weakening rates. Stress heterogeneity and roughness control hypocenter distribution, frequency-magnitude characteristics, and source mechanisms of injection-induced acoustic emissions (AEs) (analogous to natural seismicity). In contrast to smooth faults where injection-induced AEs are uniformly distributed, slip on rough faults produces spatially localized AEs with pronounced non-double-couple source mechanisms. We demonstrate that these clustered AEs occur around highly stressed asperities where induced local slip rates are higher, accompanied by lower Gutenberg-Richter b-values. Our findings suggest that real-time monitoring of induced microseismicity during fluid injection may allow identifying progressive localization of seismic activity and improve forecasting of runaway events.

Effects of reservoir mechanical properties on induced seismicity during subsurface hydrogen storage.

The intermittent storage of hydrogen in subsurface porous media such as depleted gas fields could be pivotal to a successful energy transition. Numerical simulations investigate the intermittent storage of hydrogen in a porous, depleted subsurface reservoir. Various parametric studies are performed to assess the effect of mechanical properties of the reservoir (i.e. Young's modulus, Poisson's ratio, Biot coefficient and permeability) on the induced fault slip of a single through-going fault that transverses the entire reservoir. Simulations are run using a three-dimensional, finite element, fully coupled hydromechanical code with explicit representations of layers and faults. The effect of the domain mesh refinement and fault mesh refinement on the fault slip versus operation time solution is investigated. The fault is observed to slip in two distinct events, one during the second injection period and one in the third injection period. The fault is not observed to slip during the storage or withdrawal periods. It is found that in order to minimize seismic risk, a reservoir rock with high Young's modulus (>40 GPa), high Poisson's ratio (>0.30) and high Biot coefficient (>0.65) would be preferable for hydrogen storage. Reservoir rocks of low Young's modulus (10-30 GPa), intermediate Poisson's ratio (0.00-0.30) and low-to-intermediate Biot coefficient (0.25-0.65), at high injection rates, were found to have higher potential of inducing large seismic events.This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Growth and stabilization of induced seismicity rates during long-term, low-pressure fluid injection.

We examine the temporal evolution of sequences of induced seismicity caused by long-term fluid injection using a compilation of over 20 case studies where moderate magnitude (M > 3.0) induced events have been recorded. We compare rates of seismicity with injection rates via the seismogenic index and seismic efficiency parameters, computing both cumulative and time-windowed values. We find that cumulative values tend to accelerate steeply as each seismicity sequence initiates-most cases reach a value that is within 0.5 units of their maximum value within 1-3 years. Time-windowed values tend to increase to maximum values within 25%-35% of the overall sequence, before decreasing as levels of seismicity stabilize. We interpret these observations with respect to the pore pressure changes that will be generated in highly porous, high permeability reservoirs. In such situations, the rate of pore pressure change is highest during the early phases of injection and decreases with time. If induced seismicity scales with the rate of deformation, which in turn is controlled by the rate of pore pressure change, then it is to be expected that induced seismicity is highest during the early phases of injection, and then decreases with time. This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Forecasting fluid-injection induced seismicity to choose the best injection strategy for safety and efficiency.

Induced seismicity poses a challenge to the development of Enhanced Geothermal Systems (EGS). Improving monitoring and forecasting techniques is essential to mitigate induced seismicity and thereby fostering a positive perception of EGS projects among local authorities and population. Induced seismicity is the result of complex and coupled thermo-hydro-mechanical-chemical mechanisms. Injection flux and pressure are crucial controlling parameters for both hydraulic stimulation and circulation protocols. We develop a methodology combining a hydro-mechanical model with a seismicity rate model to estimate the magnitude and frequency of mainshocks and aftershocks induced by fluid injection. We apply the methodology to the case of the Basel EGS (2006, Switzerland) to compare the effects of progressive, cyclic and constant injections on the mechanical response of discrete faults. Results from the coupled hydro-mechanical models show that the pore pressure diffusion and consequent enhancement of fault permeability are limited to the vicinity of the injection well during cyclic injection. Additionally, constant injection induces seismicity from the start of the injection but enhances the permeability of most of the faults within a shorter duration, inducing less post-injection seismicity. The methodology can be adapted to any numerical model and allows new projects to be developed by anticipating the safest injection protocol.This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Sealing faults and fluid-induced seismicity.

Sealing faults are nearly impermeable barriers that can form boundaries between subsurface pore-pressure domains. In hydrocarbon systems, sealing faults commonly form part of a structural trap; they are thus important elements for future storage of CO2 and other gases in depleted reservoirs. The Triassic Montney Formation in western Canada hosts low-permeability gas reservoirs containing sealing faults that have previously been assumed to compartmentalize pressure domains. In this study, we show that the distribution of induced seismicity associated with hydraulic fracturing (HF) exhibits a statistically significant spatial correlation with zones of high lateral gradient in pore pressure. These high-gradient zones are interpreted as sealing fault systems. The largest induced seismicity sequence, including a 4.5 ML mainshock on 30 November 2018, occurred during HF treatments in two horizontal wells, between which there is an exceptionally large contrast (~10 MPa) in measured pore pressure. Numerical simulation of a simplified model of a hydraulic fracture intersecting a nearby vertical fault, followed by fault rupture using rate-and-state friction rheology, generates results that are in good agreement with observed strike-slip faulting near one of the HF wells. Our study demonstrates that sealing faults exhibit previously unrecognized behaviour that may be important for understanding induced seismicity risk. This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Fluid-driven aseismic fault slip with permeability enhancement and dilatancy.

Injection-induced seismicity and aseismic slip often involve the reactivation of long-dormant faults, which may have extremely low permeability prior to slip. In contrast, most previous models of fluid-driven aseismic slip have assumed linear pressure diffusion in a fault zone of constant permeability and porosity. Slip occurs within a frictional shear crack whose edge can either lag or lead pressure diffusion, depending on the dimensionless stress-injection parameter that quantifies the prestress and injection conditions. Here, we extend this foundational work by accounting for permeability enhancement and dilatancy, assumed to occur instantaneously upon the onset of slip. The fault zone ahead of the crack is assumed to be impermeable, so fluid flow and pressure diffusion are confined to the interior, slipped part of the crack. The confinement of flow increases the pressurization rate and reduction of fault strength, facilitating crack growth even for severely understressed faults. Suctions from dilatancy slow crack growth, preventing propagation beyond the hydraulic diffusion length. Our new two-dimensional and three-dimensional solutions can facilitate the interpretation of induced seismicity data sets. They are especially relevant for faults in initially low permeability formations, such as shale layers serving as caprock seals for geologic carbon storage, or for hydraulic stimulation of geothermal reservoirs.This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Inferring maximum magnitudes from the ordered sequence of large earthquakes.

The largest magnitude earthquake in a sequence is often used as a proxy for hazard estimates, as consequences are often predominately from this single event (in small seismic zones). In this article, the concept of order statistics is adapted to infer the maximum magnitude ([Formula: see text]) of an earthquake catalogue. A suite tools developed here can discern [Formula: see text] influences through hypothesis testing, quantify [Formula: see text] through maximum likelihood estimation (MLE) or select the best [Formula: see text] prediction amongst several models. The efficacy of these tools is benchmarked against synthetic and real-data tests, demonstrating their utility. Ultimately, 13 cases of induced seismicity spanning wastewater disposal, hydraulic fracturing and enhanced geothermal systems are tested for volume-based [Formula: see text]. I find that there is no evidence of volume-based processes influencing any of these cases. On the contrary, all these cases are adequately explained by an unbounded magnitude distribution. This is significant because it suggests that induced earthquake hazards should also be treated as unbounded. On the other hand, if bounded cases exist, then the tools developed here will be able to discern them, potentially changing how an operator mitigates these hazards. Overall, this suite of tools will be important for better-understanding earthquakes and managing their risks. This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

A laboratory study on fault slip caused by fluid injection directly versus indirectly into a fault: implications for induced seismicity in EGSs.

Enhanced geothermal systems (EGSs) developed by hydraulic stimulation are promising for exploiting petrothermal heat by improving fluid pathways in low-permeable geothermal reservoir rocks. However, fluid injection into the subsurface can potentially cause large seismic events by reactivating pre-existing faults, which is a significant barrier to EGSs. The management of injection-induced seismicity is, therefore, essential for the success of EGSs. During the hydraulic stimulation of an EGS, fluid can be injected into a fault zone or into the rock matrix containing pre-existing faults adjacent to the injection well. The differences in hydromechanical responses between fluid injection into and adjacent to a fault have not been investigated in detail. Here, we performed triaxial fluid injection experiments involving injecting fluid directly and indirectly into a fault in granite rock samples to analyse the distinct hydromechanical responses and estimate the injection-induced seismicity in both cases. Our results suggest that in addition to directly injecting fluid into a critically stressed fault, injecting into nearly intact granite adjacent to the fault could also cause injection-induced seismic hazards owing to the high fluid pressure required to create new fractures in the granite matrix. It is, therefore, important to carefully identify pre-existing faults within tight reservoirs to avoid injecting fluid adjacent to them. Additionally, once prior unknown faults are delineated during hydraulic stimulation, appropriate shut-in strategies should be implemented immediately to mitigate seismic risks. This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Injection-induced basement seismicity beneath the Raton Basin: constraints from refined fault architectures and basin structure.

Intraplate earthquakes induced by anthropogenic fluid injection present unexpected seismic risk to previously quiescent or low seismicity-rate regions. Despite many studies of induced seismicity, there are relatively few with detailed openly accessible constraints on the interaction between seismic sources and subsurface structures. In this study of the Raton Basin, we refine source observations from a dense nodal array and constrain basin structure using teleseismic receiver functions. The cross-correlation-based relocated hypocentres and a new set of focal mechanisms light up active fault segments and show clear spatiotemporal patterns. The geometric complexity of reactivated fault clusters appears greatest near higher rate injection wells. Simpler normal fault structure is found farther from injection wells and near abrupt structural transitions suggested by receiver functions. While less induced seismicity in the crystalline basement is expected when injection is >1 km from the top of the basement (like Raton), our receiver function analysis identified a basin thickness ~3 km beneath the nodal array and lateral variations in sedimentary structures. Our results explain potential fluid connectivity between the injection depths focused at ~1-1.5 km below the surface and basement fault activity that begins at ~3 km and reaches peak activity at ~4-8 km depths. This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Widespread deep seismicity in the Delaware Basin, Texas, is mainly driven by shallow wastewater injection.

Industrial activity away from plate boundaries can induce earthquakes and has evolved into a global issue. Much of the induced seismicity in the United States' midcontinent is attributed to a direct pressure increase from deep wastewater disposal. This mechanism is not applicable where deep basement faults are hydraulically isolated from shallow injection aquifers, leading to a debate about the mechanisms for induced seismicity. Here, we compile industrial, seismic, geodetic, and geological data within the Delaware Basin, western Texas, and calculate stress and pressure changes at seismogenic depth using a coupled poroelastic model. We show that the widespread deep seismicity is mainly driven by shallow wastewater injection through the transmission of poroelastic stresses assuming that unfractured shales are hydraulic barriers over decadal time scales. A zone of seismic quiescence to the north, where injection-induced stress changes would promote seismicity, suggests a regional tectonic control on the occurrence of induced earthquakes. Comparing the poroelastic responses from injection and extraction operations, we find that the basement stress is most sensitive to shallow reservoir hydrogeological parameters, particularly hydraulic diffusivity. These results demonstrate that intraplate seismicity can be caused by shallow human activities that poroelastically perturb stresses at hydraulically isolated seismogenic depths, with impacts on seismicity that are preconditioned by regional tectonics.

Dynamic rupture initiation and propagation in a fluid-injection laboratory setup with diagnostics across multiple temporal scales.

Fluids are known to trigger a broad range of slip events, from slow, creeping transients to dynamic earthquake ruptures. Yet, the detailed mechanics underlying these processes and the conditions leading to different rupture behaviors are not well understood. Here, we use a laboratory earthquake setup, capable of injecting pressurized fluids, to compare the rupture behavior for different rates of fluid injection, slow (megapascals per hour) versus fast (megapascals per second). We find that for the fast injection rates, dynamic ruptures are triggered at lower pressure levels and over spatial scales much smaller than the quasistatic theoretical estimates of nucleation sizes, suggesting that such fast injection rates constitute dynamic loading. In contrast, the relatively slow injection rates result in gradual nucleation processes, with the fluid spreading along the interface and causing stress changes consistent with gradually accelerating slow slip. The resulting dynamic ruptures propagating over wetted interfaces exhibit dynamic stress drops almost twice as large as those over the dry interfaces. These results suggest the need to take into account the rate of the pore-pressure increase when considering nucleation processes and motivate further investigation on how friction properties depend on the presence of fluids.

Applying the extended parallel process model to understand households' responses to tornado and earthquake risks in Oklahoma.

Oklahoma is a multihazard environment where both natural (e.g., tornadoes) and technological hazards (e.g., induced seismicity) are significant, making Oklahoma a unique setting to better understand how to manage and prepare for multiple hazards. While studies have attempted to understand drivers of hazard adjustments, few have focused on the overall number of adjustments undertaken instead of individual adjustments or adjustments in a multihazard environment. To address these gaps, we employ a survey sample of 866 households in Oklahoma to understand households' danger control responses (protective hazard adjustments) for tornado and earthquake risks in Oklahoma. We apply the extended parallel processing model (EPPM) to categorize respondents according to their relative level of perceived threat and efficacy of protective actions in predicting the number of hazard adjustments they intend to or have adopted in response to tornadoes and induced earthquakes. In line with the EPPM, we found that households have the highest number of danger control responses when their perceived threat and efficacy are both high. Counter to the EPPM literature, we found low threat coupled with high efficacy moved some individuals toward the adoption of danger control responses in response to both tornadoes and earthquakes. When households have high efficacy, threat appraisals matter in tornado danger control responses but not in earthquake danger control responses. This EPPM categorization opens new research approaches for studies of natural and technological hazards. This study also provides information for local officials and emergency managers making mitigation and preparedness investments and policies.

Methane Emissions from Non-producing Oil and Gas Wells and the Potential Role of Seismic Activity: A Case Study in Northeast British Columbia, Canada.

Increasing seismic activity due to fluid injections for oil and gas production may be contributing to leakage along non-producing oil and gas wells and emitting methane, a potent greenhouse gas. However, the extent to which nearby seismicity may drive or exacerbate methane emissions and cause well integrity issues is unknown. Therefore, we analyze field evaluations at 448 non-producing oil and gas wells in Northeast British Columbia (NEBC) and geospatially analyze oil and gas well and fluid injection data alongside locations of 3515 earthquakes from 2001 to 2021 and 130 faults. Through analysis of ground and helicopter-based field evaluations of non-producing wells in NEBC, we show that methane emission rates of non-producing wells average at 8301 mg/h/well but vary by 10 orders of magnitude. We find that higher methane emission rates (milligrams of methane/h/well) are observed at wells with larger flowing pressures at the wellhead during completion (kPa) and with shorter distances (m) to earthquakes, particularly at plugged wells. These results imply that seismicity may increase the likelihood of non-producing well integrity issues and methane leakage, thereby also exacerbating groundwater contamination and environmental degradation risks.

Active strike-slip faults and an outer frontal thrust in the Himalayan foreland basin.

The Himalayan foreland basin formed by flexure of the Indian Plate below the advancing orogen. Motion on major thrusts within the orogen has resulted in damaging historical seismicity, whereas south of the Main Frontal Thrust (MFT), the foreland basin is typically portrayed as undeformed. Using two-dimensional seismic reflection data from eastern Nepal, we present evidence of recent deformation propagating >37 km south of the MFT. A system of tear faults at a high angle to the orogen is spatially localized above the Munger-Saharsa basement ridge. A blind thrust fault is interpreted in the subsurface, above the sub-Cenozoic unconformity, bounded by two tear faults. Deformation zones beneath the Bhadrapur topographic high record an incipient tectonic wedge or triangle zone. The faults record the subsurface propagation of the Main Himalayan Thrust (MHT) into the foreland basin as an outer frontal thrust, and provide a modern snapshot of the development of tectonic wedges and lateral discontinuities preserved in higher thrust sheets of the Himalaya, and in ancient orogens elsewhere. We estimate a cumulative slip of ∼100 m, accumulated in <0.5 Ma, over a minimum slipped area of ∼780 km2 These observations demonstrate that Himalayan ruptures may pass under the present-day trace of the MFT as blind faults inaccessible to trenching, and that paleoseismic studies may underestimate Holocene convergence.

Cluster Analysis of Seismicity in the Eastern Gulf of Corinth Based on a Waveform Template Matching Catalog.

The Corinth Rift, in Central Greece, is one of the most seismically active areas in Europe. In the eastern part of the Gulf of Corinth, which has been the site of numerous large and destructive earthquakes in both historic and modern times, a pronounced earthquake swarm occurred in 2020-2021 at the Perachora peninsula. Herein, we present an in-depth analysis of this sequence, employing a high-resolution relocated earthquake catalog, further enhanced by the application of a multi-channel template matching technique, producing additional detections of over 7600 events between January 2020 and June 2021. Single-station template matching enriches the original catalog thirty-fold, providing origin times and magnitudes for over 24,000 events. We explore the variable levels of spatial and temporal resolution in the catalogs of different completeness magnitudes and also of variable location uncertainties. We characterize the frequency-magnitude distributions using the Gutenberg-Richter scaling relation and discuss possible b-value temporal variations that appear during the swarm and their implications for the stress levels in the area. The evolution of the swarm is further analyzed through spatiotemporal clustering methods, while the temporal properties of multiplet families indicate that short-lived seismic bursts, associated with the swarm, dominate the catalogs. Multiplet families present clustering effects at all time scales, suggesting triggering by aseismic factors, such as fluid diffusion, rather than constant stress loading, in accordance with the spatiotemporal migration patterns of seismicity.

Multifractal Properties of Time Series of Synthetic Earthquakes Obtained from a Spring-Block Model.

With the spring-block model proposed by Olami, Feder, and Christensen (OFC), we obtained a time series of synthetic earthquakes with different values of the conservation level (ß), which measures the fraction of the energy that a relaxing block passes to its neighbors. The time series have multifractal characteristics, and we analyzed them with the Chhabra and Jensen method. We calculated the width, symmetry, and curvature parameters for each spectrum. As the value of conservation level increases, the spectra widen, the symmetric parameter increases, and the curvature around the maximum of the spectra decreases. In a long series of synthetic seismicity, we located earthquakes of the greatest magnitude and built overlapping windows before and after them. For the time series in each window, we performed multifractal analysis to obtain multifractal spectra. We also calculated the width, symmetry, and curvature around the maximum of the multifractal spectrum. We followed the evolution of these parameters before and after large earthquakes. We found that the multifractal spectra had greater widths, were less skewed to the left, and were very pointed around the maximum before rather than after large earthquakes. We studied and calculated the same parameters and found the same results in the analysis of the Southern California seismicity catalog. This suggests that there seems to be a process of preparation for a great earthquake and that its dynamics are different from the one that occurs after this mainshock based on the behavior of the parameters mentioned before.

Tsallis Entropy and Mutability to Characterize Seismic Sequences: The Case of 2007-2014 Northern Chile Earthquakes.

Seismic data have improved in quality and quantity over the past few decades, enabling better statistical analysis. Statistical physics has proposed new ways to deal with these data to focus the attention on specific matters. The present paper combines these two progressions to find indicators that can help in the definition of areas where seismic risk is developing. Our data comes from the IPOC catalog for 2007 to 2014. It covers the intense seismic activity near Iquique in Northern Chile during March/April 2014. Centered in these hypocenters we concentrate on the rectangle Lat-22-18 and Lon-68-72 and deepness between 5 and 70 km, where the major earthquakes originate. The analysis was performed using two complementary techniques: Tsallis entropy and mutability (dynamical entropy). Two possible forecasting indicators emerge: (1) Tsallis entropy (mutability) increases (decreases) broadly about two years before the main MW8.1 earthquake. (2) Tsallis entropy (mutability) sharply decreases (increases) a few weeks before the MW8.1 earthquake. The first one is about energy accumulation, and the second one is because of energy relaxation in the parallelepiped of interest. We discuss the implications of these behaviors and project them for possible future studies.

Defining the Scale to Build Complex Networks with a 40-Year Norwegian Intraplate Seismicity Dataset.

We present a new complex network-based study focused on intraplate earthquakes recorded in southern Norway during the period 1980-2020. One of the most recognized limitations of spatial complex network procedures and analyses concerns the definition of adequate cell size, which is the focus of this approach. In the present study, we analyze the influence of observational errors of hypocentral and epicentral locations of seismic events in the construction of a complex network, looking for the best cell size to build it and to develop a basis for interpreting the results in terms of the structure of the complex network in this seismic region. We focus the analysis on the degree distribution of the complex networks. We observed a strong result of the cell size for the slope of the degree distribution of the nodes, called the critical exponent γ. Based on the Abe-Suzuki method, the slope (γ) showed a negligible variation between the construction of 3- and 2-dimensional complex networks. The results were also very similar for a complex network built with subsets of seismic events. These results suggest a weak influence of observational errors measured for the coordinates latitude, longitude, and depth in the outcomes obtained with this particular methodology and for this high-quality dataset. These results imply stable behavior of the complex network, which shows a structure of hubs for small values of the cell size and a more homogeneous degree distribution when the cell size increases. In all the analyses, the γ parameter showed smaller values of the error bars for greater values of the cell size. To keep the structure of hubs and small error bars, a better range of the side sizes was determined to be between 8 to 16 km. From now on, these values can be used as the most stable cell sizes to perform any kind of study concerning complex network studies in southern Norway.

Optimized Traffic Light System with AIC and Application to the 2021 M6.7 Yangbi Earthquake Sequence.

One important question in earthquake prediction is whether a moderate or large earthquake will be followed by an even bigger one. Through temporal b-value evolution analysis, the traffic light system can be used to estimate if an earthquake is a foreshock. However, the traffic light system does not take into account the uncertainty of b-values when they constitute a criterion. In this study, we propose an optimization of the traffic light system with the Akaike Information Criterion (AIC) and bootstrap. The traffic light signals are controlled by the significance level of the difference in b-value between the sample and the background rather than an arbitrary constant. We applied the optimized traffic light system to the 2021 Yangbi earthquake sequence, which could be explicitly recognized as foreshock-mainshock-aftershock using the temporal and spatial variations in b-values. In addition, we used a new statistical parameter related to the distance between earthquakes to track earthquake nucleation features. We also confirmed that the optimized traffic light system works on a high-resolution catalog that includes small-magnitude earthquakes. The comprehensive consideration of b-value, significance probability, and seismic clustering might improve the reliability of earthquake risk judgment.

Forecasting Strong Subsequent Earthquakes in Greece with the Machine Learning Algorithm NESTORE.

Aftershocks of earthquakes can destroy many urban infrastructures and exacerbate the damage already inflicted upon weak structures. Therefore, it is important to have a method to forecast the probability of occurrence of stronger earthquakes in order to mitigate their effects. In this work, we applied the NESTORE machine learning approach to Greek seismicity from 1995 to 2022 to forecast the probability of a strong aftershock. Depending on the magnitude difference between the mainshock and the strongest aftershock, NESTORE classifies clusters into two types, Type A and Type B. Type A clusters are the most dangerous clusters, characterized by a smaller difference. The algorithm requires region-dependent training as input and evaluates performance on an independent test set. In our tests, we obtained the best results 6 h after the mainshock, as we correctly forecasted 92% of clusters corresponding to 100% of Type A clusters and more than 90% of Type B clusters. These results were also obtained thanks to an accurate analysis of cluster detection in a large part of Greece. The successful overall results show that the algorithm can be applied in this area. The approach is particularly attractive for seismic risk mitigation due to the short time required for forecasting.