Jean‐Philippe Avouac

These data and software supplement the following publication: Sirorattanakul, K. & Avouac, J.-P. Seismic Rhythms: Earthquake Response to Tectonic, Hydrological, and Tidal Forcing in California. Any questions can be directed to Krittanon (Pond) Sirorattanakul at krittanon.pond@gmail.com.

SUMMARY Realistic models of earthquake sequences can be simulated by assuming faults governed by rate-and-state friction embedded in an elastic medium. Exploring the possibility of using such models for earthquake forecasting is challenging due to the difficulty of integrating partial differential equation models with sparse, low-resolution observational data. This paper presents a machine-learning-based reduced-order model (ROM) for earthquake sequences that addresses this limitation. The proposed ROM captures the slow/fast chaotic dynamics of earthquake sequences using a low-dimensional representation, enabling computational efficiency and robustness to high-frequency noise in observational data. The ROM’s efficiency facilitates effective data assimilation using the Ensemble Kalman Filter, even with low-resolution, noisy observations. Results demonstrate the ROM’s ability to replicate key scaling properties of the sequence—namely the magnitude–frequency, moment–duration, and moment–area relationships—and to estimate the distributions of fault slip rate and state variable, enabling predictions of large events in time and space with uncertainty quantification. These findings underscore the ROM’s potential for forecasting and for addressing challenges in inverse problems for nonlinear geophysical systems.

ABSTRACT The magnitude–frequency distribution (MFD), which quantifies the relative frequency of large versus small earthquakes, is commonly used in seismic hazard assessment studies and is thought to characterize earthquake dynamics. The classic Gutenberg–Richter (GR) model posits that earthquake frequency decays exponentially with magnitude. The tapered Gutenberg–Richter (TGR) model is a variant that assumes a further reduced frequency of larger earthquakes. Distinguishing which of these two distributions better fits observations is important not only for a better understanding of earthquake physics but also for robust forecasting of earthquake magnitudes. Therefore, we evaluate methods used to differentiate these two distributions and their statistical significance given a set of observations. We find the likelihood-ratio test to be the most effective approach. It rarely misclassifies a GR distribution as a TGR distribution, whereas a TGR distribution can be misclassified as GR when the tail of the MFD is insufficiently sampled. We demonstrate that the probability of correctly identifying a TGR model exceeds 90% when the corner magnitude is one unit smaller than the maximum magnitude predicted by the GR distribution. Furthermore, we introduce an objective framework aimed at detecting potential temporal shifts between the two distributions. We apply this framework to global seismicity and two observational cases of induced seismicity. The MFD of global seismicity shows transitions between GR and TGR distributions over time, which might be explained by either inherent temporal variation in behavior or by the random sampling of a bilinear GR model with a larger b value for M > 7.6. Regarding the induced seismicity cases, we demonstrate significant and persistent TGR distributions in seismicity induced by geothermal well stimulations at Otaniemi, Finland. Furthermore, we find that earthquakes in the Coso geothermal field (California) exhibit TGR behavior during a specific period, likely influenced by the type of magnitude scale used.

We use remote sensing observations to document surface deformation caused by the 2025 M w 7.7 Mandalay earthquake. This event is a unique case of an extremely long (~510 km) and sustained supershear rupture probably favored by the rather smooth and continuous geometry of this section of the structurally mature Sagaing Fault. The seismic rupture involved the locked portion of the fault over its entire depth extent (0 to 13 km) with a remarkably uniform slip distribution that averages 3.3 m, and an average stress drop of 4.7 MPa. No shallow-slip deficit is observed. The rupture extent challenges usual scaling laws relating earthquake magnitude, fault length, and slip. The fault ruptured along a known seismic gap that last ruptured in 1839 and tailed off into sections that ruptured during large earthquakes in 1930 and 1946. The amplitude and spatial distribution of fault slip in the 2025 event conform only approximatively to the slip-predictable model and the segmentation inferred from the fault geometry and past ruptures. Plausible sequences of earthquakes with variable magnitude, segmentation, and return periods, including events similar to the 2025 earthquake are produced in quasidynamic simulations using a simplified but nonplanar fault geometry. Based on this simulation, M w >7.5 events return irregularly with an interevent time of ~141 y on average and a SD of ~40 y. The simulation is consistent with the historical seismicity and with the maximum magnitude ~M w 7.9 and return period (~250 y) derived from moment conservation. Data assimilation into such simulations could provide a way for time-dependent hazard assessment in the future.

Realistic models of earthquake sequences can be simulated by assuming faults governed by rate-and-state friction embedded in an elastic medium. Exploring the possibility of using such models for earthquake forecasting is challenging due to the difficulty of integrating Partial Differential Equation (PDE) models with sparse, low-resolution observational data. This paper presents a machine-learning-based reduced-order model (ROM) for earthquake sequences that addresses this limitation. The proposed ROM captures the slow/fast chaotic dynamics of earthquake sequences using a low-dimensional representation, enabling computational efficiency and robustness to high-frequency noise in observational data. The ROM's efficiency facilitates effective data assimilation using the Ensemble Kalman Filter (EnKF), even with low-resolution, noisy observations. Results demonstrate the ROM's ability to replicate key scaling properties of the sequence —namely the magnitude–frequency, moment–duration, and moment–area relationships— and to estimate the distributions of fault slip rate and state variable, enabling predictions of large events in time and space with uncertainty quantification. These findings underscore the ROM's potential for forecasting and for addressing challenges in inverse problems for nonlinear geophysical systems.

Mobile barchan dunes are well-developed in a crater in the Hellespontus Montes region on Mars. Previous studies of their temporal evolution show that the barchans maintain their shape and migrate in a uniform pattern. Whereas barchans are typically associated with unidirectional wind regimes, the crater experiences seasonal changes in wind regime, driven by large-scale circulation patterns. Using a multi-scale modelling approach we demonstrate that the effect of upwind mesas are minimal to steering regional wind conditions, beyond the extent of the mesas themselves which limits the effect on the development and maintenance of barchan dunes further downwind. The results of high resolution CFD modelling showed individual barchan dunes had the capability to locally steer oblique wind flows along the orientation of the barchan dunes. We hypothesise that this ability of barchan dunes to ‘steer’ near surface local wind flows, combined with the uni-directional source of sediment at the site allows barchan morphology to persist in Hellespontus Montes, despite being subject to a range of incident wind directions.

Abstract The seismicity rate model of J. Dieterich (1994, https://doi.org/10.1029/93jb02581 ) has been used extensively in recent years to model induced seismicity in fluid injection settings. In this study, we highlight its major assumptions and examine how they may bias the interpretations of inferred model parameters when applied to induced seismicity observed in real reservoirs. We do so by comparing the model to numerical simulations of induced earthquakes in a discrete fault network (DFN). The seismicity patterns show significant differences between the DFN and Dietrich model. In particular, DFN simulations show the development of a backfront during the injection due to the exhaustion of available nucleation sources and a shut‐off of near‐well seismicity due to aseismic slip at high pore pressure. When matched to the DFN catalogs, both the background seismicity rate parameter, R b , and the direct effect rate‐and‐state friction parameter, a , of the Dieterich model both significantly underestimate the same parameters of the DFN. We recall the seismicity induced by the 1993 stimulation of the GPK1 injection well in Soultz‐Sous‐Forts for evidence of co‐injection backfronts and aseismic slip among small faults. We do not discount the use of the Dieterich model—they successfully reproduce seismicity simulated from a wide range of initial conditions and fault network conditions—but emphasize that interpretations of inferred parameters must take into account finite size effects that are neglected by the model's original assumptions.

In this study, we apply control theory to mitigate earthquake hazards to a stress-based model of enhanced geothermal stimulation. The model considers pore pressure diffusion as the main stressing mechanism and rate-and-state friction as the shear failure mechanism. The controller is designed to follow a given average pressure and the probability of exceedance of a red-light earthquake (the magnitude at which the stimulation would have to stop by regulation) within chosen volumes surrounding the injection source and within a target time. We rigorously prove that the proposed controller can effectively force two output types within the system to given references, despite the presence of model uncertainties, and with minimal system information, using a continuous control signal. This framework is applied to a validated model of the 2018 Otaniemi geothermal stimulation. We use a suite of simulations to identify injection scenarios that outperform the 2018 Otaniemi stimulation. The optimal stimulation achieves higher average pressure in a shorter time with lower seismic hazard. The controller can help determine whether a combination of safety thresholds and optimization targets is feasible and economical. The control framework could be used to design stimulation schedules for enhanced geothermal systems .

Coseismic ruptures release stored elastic strain by a combination of shear displacement along localized, principal faults and distributed bulk inelastic failure of the surrounding material. How<br/>inelastic strain localizes as fault systems mature and structurally develop is less well understood due to the difficulty of measuring the complex, near-field and high-strain regions of coseismic surface ruptures. Here, we use radar and optical images to measure the near-field surface displacement field and magnitude of off-fault inelastic strain from 16 historic strike-slip earthquakes that occurred on faults spanning almost three orders of magnitude in cumulative displacement and fault slip rates. Our results show that inelastic shear deformation does localize<br/>as fault systems mature, where the magnitude of off-fault inelastic strain is largest (34-67%) for fault systems with the lowest cumulative displacement (&lt; 3 km) that then rapidly decays to values that saturate around 13-19% for the most ‘mature’ fault systems with cumulative displacement exceeding ~20 km. We find that more localized coseismic ruptures host faster ruptures, generate fewer aftershocks and occur along geometrically simpler fault networks.

Efforts to secure and decarbonize the energy sector are driving various subsurface reservoir operations. These operations carry a risk of inducing surface deformation and earthquakes. To assess these risks, modeling tools integrating fluid flow, geomechanical and seismicity modeling are needed. Here, we demonstrate the use of an efficient Vertical Flow Equilibrium (VFE) multiphase fluid flow model in an integrated framework for deformation and seismicity modeling both under fluid extraction or injection configurations. The VFE-computed spatio-temporal pressure evolution is fed to a geomechanical module to compute surface deformation and stress changes in and around the reservoir. Stress changes feed a seismicity module to calculate earthquake probabilities. First, we apply the benchmarked model to gas extraction from Groningen. There, we can reduce the variance of pressure measurements by ∼ 38% with respect to a pre-existing single phase flow model while remaining computationally efficient. The surface deformation and seismicity simulations show remarkable agreement with observed data. Second, we study induced seismicity due to CO 2 sequestration in the Decatur phase 1 project. We find that, for the Decatur phase 1 project, poroelastic stress changes can account for most of the non-clustered observed seismicity within modeling uncertainties. Finally we simulate scenarios for CO 2 sequestration using the Quest field. The sloping reservoir topography significantly impacts the predicted position of the CO 2 plume but the effects on geomechanical deformation (and seismicity) are minimal. Incorporating VFE models with geomechanical and seismicity forecasts with real-world case applications can allow real-time hazard assessment and mitigation procedures.

Abstract Slow slip events (SSEs) have been observed in many subduction zones and are understood to result from frictional unstable slip on the plate interface. The diversity of their characteristics and the fact that interplate slip can also be seismic suggest that frictional properties are heterogeneous. We are however lacking methods to determine spatial variations of frictional properties. In this paper, we employ a Physics‐Informed Neural Network (PINN) to achieve this goal using a synthetic model inspired by the long‐term SSEs observed in the Bungo channel. PINN is a deep learning technique that can be used to solve the differential equations representing the physics of the problem and determine the model parameters from observations. We start with an idealized case where it is assumed that fault slip is directly observed. We next move to a more realistic case where the observations consist of synthetic surface displacement velocity data measured by virtual GNSS stations. We find that the geometry and friction properties of the velocity weakening region, where the slip instability develops, are well estimated, especially if surface displacement velocity above the velocity weakening region is observed. Our PINN‐based method can be seen as an inversion technique with the regularization constraint that fault slip obeys a particular friction law. This approach remediates the issue that standard regularization techniques are based on non‐physical constraints. Our results show that the PINN‐based method is a promising approach for estimating the spatial distribution of friction parameters from GNSS observations.

ABSTRACT We present a method for earthquake causal attribution, which allows us to quantify the probability that an event is due to tectonic loading, a previous earthquake, or a fluid injection. The method is an extension of the stochastic declustering algorithm of Marsan and Lengliné (2008). Earthquake triggering is represented by nonparametric, mean-field kernels, which scale linearly with the seismic moment or hydraulic energy of the trigger. The kernels are estimated based on a linear intensity model via expectation–maximization, with uncertainties derived from Gaussian approximation of the incomplete-data likelihood. Some general implications of the resulting probabilistic causal structure, including an explicit algorithm to quantify the cascading effects, are illustrated. The estimators are validated using synthetic catalogs generated with an extended epidemic-type aftershock sequence model, which accounts for injection-induced earthquakes. Application to southern California seismicity and comparisons with the nearest-neighbor distance declustering method support the linearity assumption in the seismic moment. Application to seismicity related to CO2 injection in the Illinois Basin-Decatur Project (for the period 2011–2014) reveals that 11% of the earthquakes were directly triggered by injection, 89% were due to previous earthquakes, whereas the contribution from tectonic loading was negligible (&amp;lt;1%). The earthquake interaction kernels in both cases show ∼1/t decay in time and indicate triggering by elastic static stress transfer; the injection kernels in the Decatur case suggest pore-pressure diffusion as a more likely mechanism than poroelasticity. The Gutenberg–Richter b-value is estimated to be larger for anthropogenic events (∼1.4) than natural ones (∼1.0). Deviations from the model suggest spatial anisotropy of earthquake interaction in both natural and induced settings.

Dunes and ripples are markers of eolian activity. Dunes arise and evolve from the action of wind blowing on sand grains and can thus provide information on past and current wind regime. They constantly adjust and adapt their shape through feedback between the bed topography and the near-surface air flow. This interaction modulates erosion of the stoss side and deposition on the lee side, and eventually results in the dune migration. Here, we present a workflow that quantitatively relates the rate of barchan dunes migration, which can be measured from remote sensing, to the wind velocity, either measured at a meteorological station or extracted from reanalysis data. We validate this workflow using data from Earth and apply it on Mars. The workflow requires the selection of a sand transport law and the use of computational fluid dynamic (CFD) modeling. This modeling is used to estimate the effect of the local topography on the near surface airflow, namely the speed-up effect. We compare the dune migration rate predicted through the workflow to remote sensing observations, at two barchan dune fields located along the southern rim of the Arabia Gulf. After validating this workflow on Earth, we apply it to a barchan dune field on Mars. The dune migration is used to derive a wind speed distribution, averaged over one Martian year. Finally, we use ripple migration, that is much faster than dune migration, to derive the sub-annual variation of the wind speed.

Abstract A variety of geo‐energy operations involve extraction or injections of fluids, including hydrocarbon production or storage, hydrogen storage, CO 2 sequestration, and geothermal energy production. The surface deformation resulting from such operations can be a source of information on reservoir geomechanical properties as we show in this study. We analyze the time‐dependent surface deformation in the Groningen region in northeastern Netherlands using a comprehensive geodetic data set, which includes InSAR (Radarsat2, TerraSAR‐X, Sentinel‐1), GNSS, and optical leveling spanning several decades. We resort to an Independent Component Analysis (ICA) to isolate deformation signals of various origins. The signals related to gas production from the Groningen gas field and from seasonal storage at Norg Underground Gas Storage are clearly revealed. Surface deformation associated to the Groningen reservoir show decadal subsidence, with spatially variable subsidence rates dictated by local compressibility. The ICA reveals distinct seasonal fluctuations at Norg, closely mirroring the variations of gas storage. By comparing the observed long‐term subsidence within the Groningen reservoir and seasonal oscillations at Norg from a linear poroelastic compaction model, we quantify the fraction of inelastic deformation of the reservoir in space and time and constrain the reservoir compressibility. In Groningen, increased compressibility indicates inelastic compaction that has built over time and might account for as much as 20% of the total compaction cumulated until 2021, while Norg shows no signs of inelastic deformation and a constant compressibility. This study provides a methodology to monitor and calibrate models of the subsurface deformation induced by geo‐energy operations or aquifer management.

The magnitude-frequency distribution (MFD) quantifies the relative frequency of large versus small earthquakes, serving as a critical tool for investigating earthquake dynamics and assessing seismic hazard. The classic Gutenberg-Richter (GR) model posits that earthquake occurrences follow an exponential decay. The tapered Gutenberg-Richter (TGR) distribution, a variant of the GR model, assumes a reduced frequency of larger earthquakes. It might be a more appropriate model, as there should be a physical upper limit to the magnitude of both natural and induced seismicity. Distinguishing these two distributions is important not only for a better understanding of earthquake physics but also for both precise estimation of the MFD parameters and robust forecasting of earthquake magnitudes. We therefore evaluate various methods which can be used to differentiate these two distributions and their statistical significance given a set of observations. We find the likelihood-ratio test to be the most effective approach and demonstrate that the probability of correctly identifying the TGR catalog exceeds 90% when the corner magnitude is one unit smaller than the maximum magnitude predicted by the GR distribution. Furthermore, we introduce an objective framework aimed at detecting potential temporal shifts between the two distributions. We apply this framework to revisit several observational cases, including global earthquakes and induced seismicity. We observe diverse behaviors of MFD and significant temporal variations. Some of these signals may imply special earthquake physics, while others can be the artifacts of catalog generating.

Abstract We analyze the evolution of the rupture radius and maximum magnitude (Mmax) of injection-induced earthquakes on faults obeying rates and state friction. We define the radii of two different slip modes, aseismic (Ra) and seismic slip (Rs), and derive an expression for maximum magnitude evolution. If the flow rate is sufficiently high, the seismic moment grows with the scaled injection volume, Qt/wS, as M∼Cf(Qt/wS)3/2, in which Cf depends on the initial stress level, S is storage coefficient, and w is the thickness of the reservoir. These findings are confirmed using numerical simulations conducted with varied initial states. The simulations show that Rs behaves as a rupture arrest radius and Ra behaves as the minimum possible radius of aseismic creep at a given injection volume. The Mmax evolution curve can be steeper if the fault is slightly critically stressed. A high-flow rate results in frequent seismic events, starting at relatively low-injected volume, which helps track the evolution of Mmax, providing a way to anticipate the risk of a large event. Conversely, a low-flow rate allows for a larger volume injection without seismic events but may lead to sudden large events without precursory events.

SUMMARY Seismic and aseismic slip events result from episodic slips on faults and are often chaotic due to stress heterogeneity. Their predictability in nature is a widely open question. In this study, we forecast extreme events in a numerical model. The model, which consists of a single fault governed by rate-and-state friction, produces realistic sequences of slow events with a wide range of magnitudes and interevent times. The complex dynamics of this system arise from partial ruptures. As the system self-organizes, the state of the system is confined to a chaotic attractor of a relatively small dimension. We identify the instability regions within this attractor where large events initiate. These regions correspond to the particular stress distributions that are favourable for near complete ruptures of the fault. We show that large events can be forecasted in time and space based on the determination of these instability regions in a low-dimensional space and the knowledge of the current slip rate on the fault.

This study presents and validates a workflow that quantitatively links the rate of barchan dunes migration, which can be measured from remote sensing, to the wind velocity, either measured at a meteorological station or extracted from reanalysis data. The workflow requires the selection of a sand transport law and a procedure to estimate the effect of the local topography on the near surface airflow, namely the speed-up effect, that results from the compression of streamlines as the wind climbs up the dune topography. Additionally, the estimate of sand flux under natural conditions needs to account for short duration wind gusts which are usually not fully accounted for or sampled in climatic models. Those spatial and temporal variations of wind speed have a strong influence on the local sand flux due to the non-linearity of the sand transport models. We investigate these effects by using computational fluid dynamic (CFD) modeling to estimate the speed-up effect on airflow and sand transport. We next include that effect to compare the predicted dune migration rate with remote sensing observations, at two desert barchan dune fields located along the southern rim of the Arabia Gulf. We find that, at the two sites, the speed-up effect increases the predicted sand flux by a factor of ∼3 and that the measured and predicted dune migration rates agree well if the sand transport law of Kok et al. (2012) is used, combined with the cessation threshold from Pähtz and Durán (2023) along with reanalysis data ERA5-Land with an hourly sampling. The proposed workflow is applicable to any barchan dune field on Earth or Mars. • We present a workflow that relates dune migration to wind data, both in-situ and from climatic models. • We introduce a new relationship to estimate the local speed-up effect on barchan dunes. • An hourly sampling of wind data is sufficient to predict sand flux through an aeolian transport law. • ERA5-Land reanalysis data can be used to infer sand flux and dune migration rate.

Abstract Gas extraction from the Groningen gas field resulted in significant induced seismicity. We analyze the magnitude‐frequency distribution of these earthquakes in space, time and in view of stress changes calculated based on gas production and reservoir properties. Previous studies suggested variations related to reservoir geometry and stress. While we confirm the spatial variations, we do not detect a clear sensitivity of b‐value to Coulomb stress changes. However, we find that b‐value correlates positively with the rate of Coulomb stress changes. This correlation is statistically significant and robust to uncertainties related to stress change calculation. This study thus points to a possible influence of stress change rate on the probability of the magnitude of induced earthquakes.

Observations of coseismic surface ruptures are widely used for examining faulting mechanics and kinematics, and in fault displacement hazard analysis. Here, first, we provide new remote sensing observations of 20 historic surface rupturing continental strike-slip earthquakes (6.1 ≤ M w ≤ 7.9). To measure the near-field surface deformation pattern, we have processed a range of raw radar and optical images from satellite and aerial platforms with pixel tracking techniques, and in two cases using standard interferometric synthetic aperture radar (InSAR). We provide a range of observations, including surface displacement maps that constrain at least the horizontal component of surface motion (2D), and for three events the full 3D components. Second, we provide strain invariants for each event, which includes the finite dilatational strain, strain magnitudes, the vorticity, and maximum shear strain. Third, we provide a total of 134 surface fault rupture traces that are mapped manually from the displacement maps. Finally, we provide 2648 measurements of coseismic horizontal fault-parallel surface fault slip that are measured objectively using swath profiles with a function fitting method that we have developed.

Slow slip events (SSEs) have been observed in many subduction zones and are understood to result from frictional unstable slip on the plate interface. The diversity of their characteristics and the fact that interplate slip can also be seismic suggest that frictional properties are heterogeneous. We are however lacking methods to constrain spatial distribution of frictional properties. In this paper, we employ Physics-Informed Neural Networks (PINNs) to achieve this goal using a synthetic model inspired by the long-term SSEs observed in the Bungo channel. PINN is a deep learning technique which can be used to solve the physics-based differential equations and determine the model parameters from observations. To examine the potential of our proposed method, we execute a series of numerical experiments. We start with an idealized case where it is assumed that fault slip is directly observed. We next move to a more realistic case where the synthetic surface displacement velocity data are observed by virtual GNSS stations. The geometry and friction properties of the velocity weakening region, where the slip instability develops, are well estimated, especially if surface displacement velocities above the velocity weakening region are observed. Our PINN-based method can be seen as an inversion technique with the regularization constraint that fault slip obeys a particular friction law. This approach remediates the issue that standard regularization techniques are based on non-physical constraints. These results of numerical experiments reveal that the PINN-based method is a promising approach for estimating the spatial distribution of friction parameters from GNSS observation.

Abstract Gas extraction from the Groningen gas reservoir, located in the northeastern Netherlands, has led to a drop in pressure and drove compaction and induced seismicity. Stress-based models have shown success in forecasting induced seismicity in this particular context and elsewhere, but they generally assume that earthquake clustering is negligible. To assess earthquake clustering at Groningen, we generate an enhanced seismicity catalog using a deep-learning-based workflow. We identify and locate 1369 events between 2015 and 2022, including 660 newly detected events not previously identified by the standard catalog from the Royal Netherlands Meteorological Institute. Using the nearest-neighbor distance approach, we find that 72% of events are background independent events, whereas the remaining 28% belong to clusters. The 55% of the clustered events are swarm-like, whereas the rest are aftershock-like. Among the swarms include five newly identified sequences propagating at high velocities between 3 and 50 km/day along directions that do not follow mapped faults or existing structures and frequently exhibit a sharp turn in the middle of the sequence. The swarms occurred around the time of the maximum compaction rate between November 2016 and May 2017 in the Zechstein layer, above the anhydrite caprock, and well-above the directly induced earthquakes that occur within the reservoir and caprock. We suggest that these swarms are related to the aseismic deformation within the salt formation rather than fluids. This study suggests that the propagating swarms do not always signify fluid migration.

The 2022 Mw 7.0 Chihshang (Taiwan) earthquake, captured by almost a dozen near-fault strong-motion seismometers, high-rate GPS and satellite data, offers a rare opportunity to examine dynamic fault rupture in detail. Using dynamic rupture simulations, we investigate the particle motions recorded at near-fault strong-motion and 1 Hz GPS stations surrounding the main asperity. Some of these stations were as close as 250 m from the fault trace as determined by sub-pixel correlation of Sentinel-2 images. Our model reproduces the observed strong asymmetry in the ground motions on either side of the fault rupture, which results from along-dip spatial variability in rake angle on the steeply dipping fault (70°) at shallow depth (2 km). Observed near-fault, pulse-like fault-parallel ground velocity larger than fault-normal velocity can be explained by a model with a sub-shear rupture speed, which may be due to shallow rupture propagation within low-velocity material and to free surface reflections. In addition, we estimate a slip-weakening distance Dc of ~0.7-0.9m from strong-motion seismogram recorded at Station F073, which is located ~250 m from the fault rupture, and the results of dynamic rupture modeling. The inferred Dc is similar to other empirically derived estimates found for crustal earthquakes. These results have important implications for near-fault ground-motion hazard.

The 2023 M7.8 and M7.5 earthquake doublet near Kahramanmaraş, Turkey, provides insight regarding how large earthquakes rupture complex faults. Here we determine the faults geometry using surface ruptures and Synthetic Aperture Radar measurements, and the rupture kinematics from the joint inversion of high-rate Global Navigation Satellite System (GNSS), strong-motion waveforms, and GNSS static displacement. The M7.8 event initiated on a splay fault and subsequently propagated along the main East Anatolian Fault with an average rupture velocity between 3.0 and 4.0 km/s. In contrast, the M7.5 event demonstrated a bilateral supershear rupture of about 5.0-6.0 km/s over an 80 km length. Despite varying strike and dip angles, the sub-faults involved in the mainshock are nearly optimally oriented relative to the local stress tensor. The second event ruptured a fault misaligned with respect to the regional stress, also hinting at the effect of local stress heterogeneity in addition to a possible free surface effect.

ABSTRACT: Much progress has recently been made in the development of stress-based models for forecasting induced earthquakes. Models based on a point-source representation of earthquake nucleation can already be used to estimate seismicity rates. Forecasting magnitudes with stress-based models remains a challenge that requires taking fault finite size and network geometry into account. Quake-DFN, an open-source earthquake simulator, was developed to address this challenge. It allows simulating sequences of earthquakes in a 3-D Discrete Fault Network governed by rate and state friction, a phenomenological law established based on laboratory observations. Our simulation method aligns with the widely used quasi-dynamic earthquake simulators, but it also has the unique capability to simulate realistic discrete fault geometry and inertial overshoot effect. Quake-DFN was benchmarked against three publicly available simulation results: (1) the rupture of a planar fault with uniform prestress, (2) the propagation of a rupture across a stepover separating two parallel planar faults, and (3) a branch fault system with a secondary fault splaying from a main fault. Next, we explored the factors that determine the magnitudes of injection-induced earthquakes for various fault geometries and loading mechanisms. Firstly, we investigated a single planar fault system. Depending on the initial state, the simulations demonstrate both self-arrested ruptures with a log-linear evolution of maximum magnitude with injection volume and runaway ruptures where the entire fault ruptured early stage. We showed that these behaviors can be theoretically explained with fracture mechanics. Using the same geometry, we additionally conduct simulations with frictional heterogeneity and find that simple heterogeneity patterns can result in a Gutenberg-Richter-like magnitude distribution, making the earthquake sequence more realistic and potentially allowing us to investigate the origin of the magnitude distribution. We then test a slightly more complex fault system – a uniformly distributed discrete faults. In this system, weakly interacting small faults are distributed with different initial states corresponding to Dietrich's model (1994), a widely used earthquake model in induced earthquake studies. We find some inherent assumptions in the Dieterich model (e.g., all faults are critically stressed) may lead to a biased interpretation of induced earthquakes. Lastly, we conducted simulations with varied initial states using a fault network and stress field similar to the one that was activated during the 2011 Prague, Oklahoma, earthquake sequence. The simulations produce realistic earthquake sequences, and a few simulations successfully reproduce the foreshock-mainshock pattern observed in the actual earthquake sequence. We note that Quake-DFN can easily be coupled (one-way) with existing geomechanical models and, hence, can further accommodate inhomogeneous permeability structures. Quake-DFN uses laboratory-measured friction parameters and further allows exploring the uncertainty of laboratory measurements by simulating a wide range of parameter space due to its low computational cost. These examples show that Quake-DFN is a useful tool to forecast a large variety of earthquake sequences and, most importantly, magnitudes induced by a fluid injection near a known fault system.

Gas extraction from the Groningen gas field resulted in significant induced seismicity. We analyze the magnitude-frequency distribution of these earthquakes in space, time and in view of stress changes calculated based on gas production and reservoir properties. Previous studies suggested variations related to reservoir geometry and decreasing b-value with increasing Coulomb stress. While we confirm the spatial variations, we do not detect a clear sensitivity of b-value to Coulomb stress. However, we find that b-value correlates positively with Coulomb stress rate. This correlation is statistically significant and robust to uncertainties related to stress calculation. This study thus points to a possible influence of stress rate on the magnitude probability of induced earthquakes.

ABSTRACT We present an earthquake simulator, Quake-DFN, which allows simulating sequences of earthquakes in a 3D discrete fault network governed by rate and state friction. The simulator is quasi-dynamic, with inertial effects being approximated by radiation damping and a lumped mass. The lumped mass term allows for accounting for inertial overshoot and, in addition, makes the computation more effective. Quake-DFN is compared against three publicly available simulation results: (1) the rupture of a planar fault with uniform prestress (SEAS BP5-QD), (2) the propagation of a rupture across a stepover separating two parallel planar faults (RSQSim and FaultMod), and (3) a branch fault system with a secondary fault splaying from a main fault (FaultMod). Examples of injection-induced earthquake simulations are shown for three different fault geometries: (1) a planar fault with a wide range of initial stresses, (2) a branching fault system with varying fault angles and principal stress orientations, and (3) a fault network similar to the one that was activated during the 2011 Prague, Oklahoma, earthquake sequence. The simulations produce realistic earthquake sequences. The time and magnitude of the induced earthquakes observed in these simulations depend on the difference between the initial friction and the residual friction μi−μf, the value of which quantifies the potential for runaway ruptures (ruptures that can extend beyond the zone of stress perturbation due to the injection). The discrete fault simulations show that our simulator correctly accounts for the effect of fault geometry and regional stress tensor orientation and shape. These examples show that Quake-DFN can be used to simulate earthquake sequences and, most importantly, magnitudes, possibly induced or triggered by a fluid injection near a known fault system.

Abstract Fault maturity has been proposed to exert a first order control on earthquake rupture, yet direct observations linking individual rupture to long‐term fault growth are rare. The 2021 Mw 7.4 Maduo earthquake ruptured the east‐growing end of the slow‐moving (∼1 mm/yr) Jiangcuo fault in north Tibet, providing an opportunity to examine the relation between rupture characteristics and fault structure. Here we combine field and multiple remote sensing techniques to map the surface rupture at cm‐resolution and document comprehensively on‐fault offsets and off‐fault deformation. The 158 km‐long surface rupture consists of misoriented structurally inherited N110°‐striking segments and younger optimally oriented N093°‐striking segments, relative to the regional stress field. Despite being comparatively newly formed, the ∼N093°‐striking fault segments accommodate more localized strain, with up to 3 m on‐fault left‐lateral slip and 25%–50% off‐fault deformation, and possibly faster rupture speed. These results are in contrast with previous findings showing more localized strain and faster rupture speed on more mature fault segments; instead, our observations suggest that fault orientation with respect to the regional stress can exert a more important control than fault maturity on coseismic rupture behavior when both factors are at play.

Seismic (‘earthquakes’) and aseismic (‘slow earthquakes’) slip events result from episodic slips on faults and are often chaotic due to stress heterogeneity. Their predictability in nature is a widely open question. Here, we forecast extreme events in a numerical model of a single fault governed by rate-and-state friction, which produces realistic sequences of slow events with a wide range of magnitudes and inter-event times. The complex dynamics of this system arise from partial ruptures. As the system self-organizes, prestress is confined to a chaotic attractor of a relatively small dimension. We identify the instability regions (corresponding to particular stress distributions) within this attractor which are precursors of large events. We show that large events can be forecasted in time and space based on the determination of these instability regions in a low-dimensional space and the knowledge of the current slip rate on the fault.

Fault maturity has been proposed to exert a first-order control on earthquake rupture, yet direct&amp;#160;observations linking individual rupture to long-term fault growth are rare. The 2021 Mw 7.4 Maduo earthquake ruptured the east-growing end of the slow-moving (~1 mm/yr) Jiangcuo fault&amp;#160;in north Tibet, providing an opportunity to examine the relation between rupture characteristics&amp;#160;and fault structure. Here, we combine field and multiple remote sensing techniques to map the&amp;#160;surface rupture at cm-resolution and document comprehensively on-fault offsets and off-fault&amp;#160;deformation. The 158 km-long surface rupture consists of misoriented structurally inherited&amp;#160;N110&amp;#176;-striking segments and younger optimally oriented N093&amp;#176;-striking segments, relative to&amp;#160;the regional stress field. Despite being comparatively newly formed, the ~N093&amp;#176;-striking fault&amp;#160;segments accommodate more localized strain, with up to 3 m on-fault left-lateral slip and 25-50% off-fault deformation, and possibly faster rupture speed. These results are in contrast with&amp;#160;previous findings showing more localized strain and faster rupture speed on more mature fault&amp;#160;segments; instead, our observations suggest that fault orientation with respect to the regional&amp;#160;stress can exert a more important control than fault maturity on coseismic rupture behaviors&amp;#160;when both factors are at play.

As faults accumulate displacement, they are thought to mature from disorganized and distributed fracture networks to more simplified throughgoing fault structures with a more localized zone of inelastic strain. Understanding the degree of inelastic strain localization holds importance for seismic hazard, as smoother faults are thought to host faster rupture velocities and have different seismic shaking intensities from ruptures along rougher, less mature faults. However, quantifying this evolutionary process of strain localization along major fault systems has been difficult due to a lack of near-field coseismic measurements. Here we test if such an evolutionary process exists by measuring the near-field surface deformation pattern of 17 large (6.0 &lt; Mw &lt; 7.9) continental strike-slip surface ruptures. To do this we use a range of geodetic imaging techniques including, a new 3D optical pixel tracking method, and pixel tracking of radar amplitude data acquired by satellite and UAVSAR platforms. With these geodetic imaging data we measure the total coseismic offset across the surface rupture and difference them from the displacements recorded by field surveys, which we assume captures the on-fault, discrete component of deformation. This differencing allows us to obtain an average magnitude of off-fault deformation for each surface rupturing event, which we compare to a number of known source parameters to test the notion of progressive fault localization. Our results show that progressively smaller amounts of off-fault strain occur along fault systems with higher cumulative displacements, supporting the notion that faults systems localize as they mature. We also find strong correlations of off-fault deformation with the long-term fault slip-rate and the geometrical complexity of the mapped surface rupture, and a moderate correlation with rupture velocity. However, we find a weak-no correlation of off-fault deformation with the fault initiation age and the moment-scaled radiated energy. We also present comparisons of off-fault strain with other known seismic source parameters.

Shallow aquifers are the largest freshwater bodies in the North African Sahara and the Arabian Peninsula. Their groundwater dynamics and response to climatic variability and anthropogenic discharge remain largely unquantified due to the absence of large-scale monitoring methods. Currently, the assessment of groundwater dynamics in these aquifer systems is made primarily from sporadic well logs that barely cover a few percent of the geographical extent of these water bodies. To address this deficiency, we develop the use of an ultra-wideband (UWB) very high frequency (VHF) interferometric airborne sounding radar, under a collaboration between NASA and the Qatar Foundation, to characterize the depth and geometry of the shallowest water table in large hyperarid hydrological basins in North Africa and the Arabian Peninsula. Herein, we describe the science objectives, measurement requirements, instrument design, expected performance, flight implementation scenarios, primary targets for investigation, and the first technology demonstration of the concept. Our performance analyses suggest that an airborne, nadir-looking sounding radar system operating at a 70-MHz center frequency with a linearly polarized folded-dipole antenna array—enabling a bandwidth (BW) of 50 MHz—and a surface signal-to-noise ratio (SNR) of 85 dB flying at an altitude of 500–2,000 m can map the uppermost water table depths of aquifer systems spanning tens of kilometers at a vertical resolution of 3 m in desiccated terrains to an average penetration depth of 50 m, with a spatial resolution of 200 m. For the first time, this airborne concept will allow time-coherent high-resolution mapping of the uppermost water tables of major aquifer systems in hyperarid areas, providing unique insights into their dynamics and responses to increasing climatic and anthropogenic stressors, which remain largely uncharacterized. The aforementioned significantly surpasses the existing capabilities for mapping shallow aquifers in these harsh and remote environments, which relies today on data collected on different timescales from sparse well logs that do not cover their geographic extents. A list of key abbreviations for this article can be found in “The Key Abbreviations Used in This Article.”