The effects of stress perturbations on friction are crucial for understanding earthquake triggering. Previous experimental studies have primarily been conducted at room temperature, where fault gouge materials typically exhibit velocity-strengthening and frictionally stable behaviour. In this study, we investigate how variations in effective normal stress ((Formula presented.)) influence fault (in-)stability by performing (Formula presented.) -perturbation experiments on simulated carbonate fault gouges under fluid-drained, hydrothermal conditions. Our results indicate that in the velocity-neutral or -weakening regime, perturbing (Formula presented.) can reinforce frictional instability, leading to accelerated slow slips or enhanced stick-slip events. This effect is particularly pronounced when the excitation period ((Formula presented.)) approaches or exceeds the characteristic recurrence period ((Formula presented.)) associated with pre-perturbation instabilities. Stress drops of the resulting events can have larger amplitudes than expected from quasi-steady-state. When cyclic perturbations are imposed, slip events tend to synchronize at specific phases when (Formula presented.) is close to (Formula presented.) — notably between 0.5π and 1π in radian, corresponding to maximum destressing rate and minimum (Formula presented.), respectively. Additionally, short-period ((Formula presented.)) perturbations can induce significant shear stress reduction (or fault weakening), with magnitudes comparable to the stress drops from stick slips, yet they are surprisingly associated with acoustically quiet slow slip, suggesting a stabilizing effect. These findings underscore the critical role of perturbation period in controlling fault response. In the context of induced seismicity, our results imply that cyclic or monotonic fluid injections should be carefully designed, considering both perturbation amplitude and period. Properly turned cyclic injections could potentially mitigate seismic risk by promoting quiet, slow slip over seismic fault slip.

Depletion-induced fault slip and seismicity in the Groningen natural gas field are known to be caused by compaction of reservoir rock, most likely at locations in faults where reservoir rock juxtaposes non-reservoir rock leading to severely-peaked shear stresses at the reservoir-fault corners. The resulting fault slip is probably initially aseismic until a critical nucleation length is reached. Under the assumption of slip-weakening friction, the nucleation length can be approximated with a classic stability criterion developed by Uenishi and Rice (U&R) in 2003 for a single-peaked stress distribution. Earlier work revealed that the validity of this criterion breaks down when the fault offset exceeds approximately 70% of the reservoir height because interaction effects between neighboring stress peaks can no longer be ignored. We therefore extended the U&R criterion to cope with such a double-peaked shear stress. The key mathematical innovation involves a singular functional form of the slip gradient which allows for the formulation of slip-patch end conditions that can be directly extended to multiple patches. We derived an exact double-patch eigenvalue criterion and an approximate closed-form “Extended Uenishi and Rice (EU&R) criterion” that is dependent on a parameter representing the scaled distance between the slip patches. For examples with parameter values roughly based on those of the Groningen field, we found a good agreement between our approximate EU&R criterion and an exact eigenvalue-based approach. Our results can serve as a robust code-independent termination criterion during numerical simulation of depletion-induced onset of seismicity resulting from natural gas or geothermal energy production.

A general trend in the use of the Dutch deep subsurface is a shift from hydrocarbon production to geothermal energy production and subsurface storage of CO2 and H2. A broad mistrust by the general public, and many local governments, of any deep-subsurfacerelated activity leads to an increasing tendency to block any development that might potentially cause harm to humans or the environment. A complicating factor in this respect is the uncertainty surrounding new technology and the lack of related historical data (notably seismic records). We support the opinion of Boulder and Lofstedt (2024) that to effectively address this problem, we should avoid a dichotomy of acceptable versus unacceptable risks, as stimulated by the use of the precautionary principle. Instead, we should use a Tolerability Of Risk (TOR) approach with the As Low As Reasonably Practicable (ALARP) principle as key element.

We present expressions to compute the inverse of a Cauchy-type singular integral equation representing the relation between a double-peaked Coulomb stress in a fault or fracture and the resulting slip gradient in two distinct collinear slip patches. In particular we consider a situation where the patches are close enough to account for the influence of the slip gradient in one patch on the slip-induced shear stress in the other patch and vice versa. This situation can occur during depletion-induced or injection-induced fault slip in subsurface reservoirs for, e.g., natural gas production, hydrogen or CO₂ storage, or geothermal operations. The theory for a single slip patch is well-developed but the situation is less clear for a configuration with two patches although the monographs of Muskhelishvili (1953) and Weertman (1996) provide earlier results. We show that the general inverse solution for the coupled two-patch problem requires six auxiliary conditions to ensure six physical requirements: boundedness of the slip gradient at the four end points of the slip patches and vanishing of the integrals of the slip gradient over the patches. Mathematically, the presence of two additional conditions, as compared to earlier formulations, corresponds to two undetermined coefficients in the general solution of the governing integral equation. Numerical simulation confirms that at least one of these is always non-zero in the coupled situation. For a coupled double-patch case with a symmetric pre-slip Coulomb stress pattern, the general inverse solution requires three auxiliary conditions. Moreover the conditions for the asymmetric case may be reduced to a set of four again, but these are different from the sets of four obtained earlier by Muskhelishvili (1953) and Weertman (1996). We illustrate the theory with a numerical example in which the evaluation of the Cauchy integrals is performed with a modified version of augmented Gauss–Chebyshev quadrature that relies on analytical inversion.

We critically review the derivation of closed-form analytical expressions for elastic displacements, strains, and stresses inside a subsurface reservoir undergoing pore pressure changes using inclusion theory. Although developed decades ago, inclusion theory has been used recently by various authors to obtain fast estimates of depletion-induced and injection-induced fault stresses in relation to induced seismicity. We therefore briefly address the current geomechanical relevance of this method, and provide a numerical example to demonstrate its use to compute induced fault stresses. However, the main goal of our paper is to correct some erroneous assumptions that were made in earlier publications. While the final expressions for the poroelastic stresses in these publications were correct, their derivation contained conceptual mistakes due to the mathematical subtleties that arise because of singularities in the Green's functions. The aim of our paper is therefore to present the correct derivation of expressions for the strains and stresses inside an inclusion and to clarify some of the results of the aforementioned studies. Furthermore, we present two conditions that the strain field must satisfy, which can be used to verify the analytical expressions.

Quantification of the poromechanical response of subsurface formations due to human-induced pore pressure fluctuations is critical for the performance and stability assessment of many geo-energy systems. In particular, natural faults in the subsurface introduce the hazard of induced seismicity. Numerical modeling of fault reactivation is challenging, while the specific details of induced stresses and fault slip in reservoirs with displaced (i.e. non-zero offset) faults may cause additional challenges depending on the type of numerical formulation employed. To facilitate the systematic development and testing of numerical tools for the simulation of induced seismicity in faulted reservoirs we developed a set of semi-analytical test problems of increasing complexity, based on inclusion theory and Cauchy singular integral equations. With these we investigate the accuracy of two recently developed Finite Volume (FV) schemes with collocated and staggered arrangements of unknowns. One of them employs a conformal discrete fault model (DFM) which can guarantee sufficient accuracy at the cost of adaptive mesh refinement but may suffer from modelling and computational challenges when addressing large-scale realistic geological configurations. The second one employs an embedded (or non-conformal) discrete fault model (EDFM) which avoids the need for excessive mesh refinement, but of which the accuracy and the range of applicability are still to be investigated. We found that both numerical schemes accurately represent the pre-slip Coulomb stresses, but show different degrees of accuracy in representing the resulting depletion-induced fault slip. The semi-analytical benchmark data are available via DOI 10.4121/22240309.

We consider steady-state single-phase confined flow through a subsurface porous layer containing a displaced, fully conductive fault causing a sudden jump in the flow path, and we employ (semi-)analytical techniques to compute the corresponding pressures and fault stresses. In particular, we obtain a new solution for the pressure field with the aid of conformal mapping and a Schwarz–Christoffel transformation. Moreover, we use an existing technique to compute the poro-elastic stress field with the aid of inclusion theory. The additional resistance to fluid flow provided by a displaced fault, relative to the resistance in a layer without a fault, is a function of dip angle, fault throw divided by reservoir height, and reservoir width divided by reservoir height. Fluid flow has a larger effect on fault stresses in case of injection than in case of depletion, where injection with up-dip flow results in increased zones of fault slip near the bottom of the reservoir. Opposedly, injection with down-dip flow results in increased slip near the top of the reservoir. An order-of-magnitude estimate of the effect of steady-state flow across displaced faults in the Groningen natural gas reservoir shows that the effect on fault stresses is probably negligible. A similar estimate of the effect in low-enthalpy geothermal doublets indicates that steady-state flow may possibly play a small role, in particular close to the injector, but site-specific assessments will be necessary to quantify the effect.

Recent laboratory and field studies suggest that temporal variations in injection patterns (e.g., cyclic injection) might trigger less seismicity than constant monotonic injection. This study presents results from uniaxial compressive experiments performed on Red Felser sandstone samples providing new information on the effect of stress pattern and rate on seismicity evolution. Red Felser sandstone samples were subjected to three stress patterns: cyclic recursive, cyclic progressive (CP), and monotonic stress. Three different stress rates (displacement controlled) were also applied: low, medium, and high rates of 10⁻⁴ mm/s, 5 × 10⁻⁴ mm/s, and 5 × 10⁻³ mm/s, respectively. Acoustic emission (AE) waveforms were recorded throughout the experiments using 11 AE transducers placed around the sample. Microseismicity analysis shows that (i) Cyclic stress patterns and especially cyclic progressive ones are characterized by a high number of AE events and lower maximum AE amplitude, (ii) among the three different stress patterns, the largest b-value (slope of the log frequency-magnitude distribution) resulted from the cyclic progressive (CP) stress pattern, (iii) by reducing the stress rate, the maximum AE energy and final mechanical strength both decrease significantly. In addition, stress rate remarkably affects the detailed AE signature of the events classified by the distribution of events in the average frequency (AF)—rise angle (RA) space. High stress rates increase the number of events with low AF and high RA signatures. Considering all elements of the AE analysis, it can be concluded that applying cyclic stress patterns in combination with low-stress rates may potentially lead to a more favourable induced seismicity effect in subsurface-related injection operations.

We address aseismic fault slip and the onset of seismicity resulting from depletion-induced or injection-induced stresses in reservoirs with pre-existing vertical or inclined faults. Building on classic results, we discuss semi-analytical modelling techniques for fault slip including dislocation theory, Cauchy-type singular integral equations and the use of Chebyshev polynomials for their solution and an eigenvalue-based stability analysis. We consider slip patch development during depletion for faults with zero, constant static and slip-weakening friction, and our results confirm earlier findings based on numerical simulation, in particular the aseismic growth of two slip patches that may subsequently merge and/or become unstable resulting in nucleation of seismic slip. New findings include improved approximate expressions for the induced seismic moment per unit strike length and a description of the effect of coupling between the slip patches which affects both forward simulation and eigenvalue computation for high values of the ratio of fault throw to reservoir height. Our implementation based on analytical inversion and semi-analytical integration with Chebyshev polynomials is more efficient and more robust than our numerical integration approach. It is not yet well suited for Monte Carlo simulation, which typically requires sub-second simulation times, but with some further development that option seems to be within reach. Moreover, our results offer a possibility for embedded fault modelling in large-scale numerical simulation tools.

A smoothed embedded finite-volume modeling (sEFVM) method is presented for faulted and fractured heterogeneous poroelastic media. The method casts a fully coupled strategy to treat the coupling between fault slip mechanics, deformation mechanics, and fluid flow equations. This ensures the stability and consistency of the simulation results, especially, as the fault slip is implicitly found through an iterative prediction-correction procedure. The computational grid is generated independently for embedded faults and rock matrix. The efficiency is further enhanced by extending the finite-volume discrete space by introducing only one degree of freedom per fault element. The embedded approach can lead to an oscillatory stress field at the fault, which damages the robustness of the implicit slip detection strategy. To resolve this challenge, a smoothed embedded strategy is devised, in which the stress and slip profiles are post processed within the iterative loops by fitting the best curve based on a least-square error criterion. The sEFVM provides locally conservative mass flux and stress fields, on staggered grid. Its performance is further investigated for several proof-of-the-concept test cases, including a multiple fault system in a heterogeneous domain. Results indicate that the method develops a promising approach for field-scale relevant simulation of induced seismicity.

We present a scalable collocated Finite Volume Method (FVM) to simulate induced seismicity as a result of pore pressure changes. A discrete system is obtained based on a fully-implicit fully-coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a flexible unstructured mesh. The cell-centered collocated scheme leads to a convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a generic multi-point flux approximation is formulated to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable processing algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate and efficient; and provides a promising framework for accurate and efficient simulation of induced seismicity in various geoscientific applications.

We develop a collocated Finite Volume Method (FVM) to study induced seismicity as a result of pore pressure fluctuations. A discrete system is obtained based on a fully-implicit coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a fully unstructured mesh. The cell-centered collocated scheme leads to convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a multi-point flux approximation is formulated in a general procedure to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable parallel solution algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate, stable and efficient; and as such promising for accurate and efficient simulation of induced seismicity.

We explore and develop a Proper Orthogonal Decomposition (POD)-based deflation method for the solution of ill-conditioned linear systems, appearing in simulations of two-phase flow through highly heterogeneous porous media. We accelerate the convergence of a Preconditioned Conjugate Gradient (PCG) method achieving speed-ups of factors up to five. The up-front extra computational cost of the proposed method depends on the number of deflation vectors. The POD-based deflation method is tested for a particular problem and linear solver; nevertheless, it can be applied to various transient problems, and combined with multiple solvers, e.g., Krylov subspace and multigrid methods.