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.

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.

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 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.

Recently, there is an increased interest in reactive flow in porous media, in groundwater, agricultural and fuel recovery applications. Reactive flow modeling involves vastly different reaction rates, i.e., differing by many orders of magnitude. Solving the ensuing model equations can be computationally intensive. Categorizing reactions according to their speeds makes it possible to greatly simplify the relevant model equations. Indeed some reactions proceed so slow that they can be disregarded. Other reactions occur so fast that they are well described by thermodynamic equilibrium in the time and spatial region of interest. At intermediate rates kinetics needs to be taken into account. In this paper, we categorize selected reactions as slow, fast or intermediate. We model 2D radially symmetric reactive flow with a reaction-convection-diffusion equation. We show that we can subdivide the PeDa_II phasespace in three regions. Region I (slow reaction); reaction can be ignored, region II (intermediate reaction); initially kinetics need to be taken into account, region III (fast reaction); all reaction takes places in a very narrow region around the injection point. We investigate these aspects for a few specific examples. We compute the location in phase space of a few selected minerals depending on salinity and temperature. We note that the conditions, e.g., salinity and temperature may be essential for assigning the reaction to the correct region in phase space. The methodology described can be applied to any mineral precipitation/decomposition problem and consequently greatly simplifies reactive flow modeling in porous media.

Water injection in the aquifer induces deformations in the soil. These mechanical deformations give rise to a change in porosity and permeability, which results in non-linearity of the mathematical problem. Assuming that the deformations are very small, the model provided by Biot’s theory of linear poroelasticity is used to determine the local displacement of the skeleton of a porous medium, as well as the fluid flow through the pores. In this continuum scale model, the Kozeny–Carman equation is commonly used to determine the permeability of the porous medium from the porosity. The Kozeny–Carman relation states that flow through the pores is possible at a certain location as long as the porosity is larger than zero at this location in the aquifer. However, from network models it is known that percolation thresholds exist, indicating that the permeability will be equal to zero if the porosity becomes smaller than these thresholds. In this paper, the relationship between permeability and porosity is investigated. A new permeability-porosity relation, based on the percolation theory, is derived and compared with the Kozeny–Carman relation. The strongest feature of the new approach is related to its capability to give a good description of the permeability in case of low porosities. However, with this network-inspired approach small values of the permeability are more likely to occur. Since we show that the solution of Biot’s model converges to the solution of a saddle point problem for small time steps and low permeability, we need stabilisation in the finite element approximation.

1D water oil displacement in porous media is usually described by the Buckley-Leverett equation or the Rapoport-Leas equation when capillary diffusion is included. The rectilinear geometry is not representative for near well oil displacement problems. It is therefore of interest to describe the radially symmetric Buckley-Leverett or Rapoport-Leas equation in cylindrical geometry (radial Buckley-Leverett problem). We can show that under appropriate conditions, one can apply a similarity transformation (r, t) → η= r²/ (2 t) that reduces the PDE in radial geometry to an ODE, even when capillary diffusion is included (as opposed to the situation in the rectilinear geometry (Yortsos, Y.C. (Phys. Fluids 30(10),2928–2935 1987)). We consider two cases (1) where the capillary diffusion is independent of the saturation and (2) where the capillary diffusion is dependent on the saturation. It turns out that the solution with a constant capillary diffusion coefficient is fundamentally different from the solution with saturation-dependent capillary diffusion. Our analytical approach allows us to observe the following conspicuous difference in the behavior of the dispersed front, where we obtain a smoothly dispersed front in the constant diffusion case and a power-law behavior around the front for a saturation-dependent capillary diffusion. We compare the numerical solution of the initial value problem for the case of saturation-dependent capillary diffusion obtained with a finite element software package to a partially analytical solution of the problem in terms of the similarity variable η.

Successful microbial enhanced oil recovery depends on several factors like reservoir characteristics and microbial activity. In this work, a pore network is used to study the hydrodynamic evolution over time as a result of the development of a biofilm in the pores. A new microscopic model is proposed for biofilm growth which takes into account that nutrients might not fully penetrate the biofilm. An important novelty in this model is that acknowledges the continuous spreading of the biofilm over the network. The results from the current study can be used to obtain a new relation between the porosity and permeability which might be used as an alternative to the Kozeny Carman relation.