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.

@en

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.

@en

This note provides the derivation of closed-form expressions for elastic displacements, strains, and stresses inside an inclusion. Jansen et al. (2019) and Wu et al. (2021) obtained correct expressions for the stresses inside an inclusion, but their derivation of these expressions contained mistakes. In this note, the correct derivation of expressions for the stresses inside an inclusion is presented and some of the results of the aforementioned studies are clarified.@en

Producing geothermal heat from production water causes cooldown from the reservoir temperature up to 250C at fluid pressures from over 100 bar to 10 bar.During the process degassification of CO2 and methane cause reduction in pH and by that dissolution and precipitation of minerals.At depth, mineral precipitation in the reservoir restricts flow paths through the cyclic system, resulting into injectivity loss, by that higher injection pressures result in additional costs.Due the large number of timesteps,numerically modeling mineralization, accounting for the reaction kinetics, can be computationally expensive. These simulations are less expensive when assuming a local equilibrium between the reactants and reaction-products. As described in Meulenbroek et all. (2020) we present an analytical model for mineral precipitation in a low-enthalpy geothermal reservoir.The three different reaction regimes are (1) fast reactions (2) very slow reactions (3) reaction/transport intermediate zone.We focus on the near-wellbore region in the reservoir, where precipitation can behave as a ‘skin’ and has a more dramatic impact on the injectivity than precipitation further downstream. Our numerical model uses a coupling approach between PHREEQC and COMSOL utilizing the qualification of the different reaction regimes. This methodology was validated by using an analytical solution of a specific mineralization case. In addition it was compared to a field case. @en

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.

@en

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

@en

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.

@en

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.

@en

Compressing a porous material or injecting fluid into a porous material can induce changes in the pore space, leading to a change in porosity and permeability. In a continuum scale PDE model, such as Biot’s theory of linear poroelasticity, the Kozeny–Carman equation is commonly used to determine the permeability of the porous medium from the porosity. The Kozeny–Carman relation assumes that there will be flow through the porous medium at a certain location as long as the porosity is larger than zero at this location. In contrast, from discrete network models it is known that percolation thresholds larger than zero exist, indicating that the fluid will stop flowing if the average porosity becomes smaller than a certain value dictated by these thresholds. In this study, the difference between the Kozeny–Carman equation and the equation based on the percolation theory, is investigated.@en

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.

@en

In this work, we model the biofilm growth at the microscale using a rectangular pore network model in 2D and a cubic network in 3D. For the 2D network, we study the effects of bioclogging on porosity and permeability when we change parameters like the number of nodes in the network, the network size, and the concentration of nutrients at the inlet. We use a 3D cubic network to study the influence of the number of nodes in the z direction on the biofilm growth and on upscalability. We show that the biofilm can grow uniformly or heterogeneously through the network. Using these results, we determine the conditions for upscalability of bioclogging for rectangular and cubic networks. If there is uniform biofilm growth, there is a unique relation between permeability and porosity, K ∼ ϕ2, this relation does not depend on the volume of the network, therefore the system is upscalable. However, if there is preferential biofilm growth, the porosity-permeability relation is not uniquely defined, hence upscalability is not possible. The Damköhler number is used to determine when upscalability is possible. If the Damköhler number is less than 101, the biofilm grows uniformly and therefore the system is upscalable. However, if the Damköhler number is greater than 103, the biofilm growth exhibits a deviation from uniform biofilm growth and heterogeneous growth is observed, therefore upscalability is not possible. There is a transition from uniform growth to preferential growth if the Damköhler number is between 101 and 103.@en

After the primary extraction in oil reservoirs up to 60 % of the oil remains trapped in the reservoir (Sen, 2008). Therefore, different mechanisms have been developed to get the oil out to the reservoir. One of these techniques is Microbial Enhanced Oil Recovery (MEOR) which is a technique used to produce more oil in a secondary extraction by using microbes in the reservoir. The main effects caused by microbes in oil recovery is the reduction of the interfacial tension between oil and water, wettability change of the rock and bioclogging caused by the growth and development of biofilm. Among these mechanisms, interfacial tension reduction and biclogging is thought to have the greatest impact on recovery (Sen, 2008). In this work, we describe the growth of biofilm, the growth of the microbial population and the transport of nutrients using a pore network model. We follow the previous models of Thullner et al. (Thullner, 2008) and Ezeuko et al. (Ezeuko, 2011) in which the biofilm is considered as a permeable layer. We consider the biofilm and the bacteria separately. Additionally, we assume that once a tube is full with biofilm, this biofilm can spread to the neighboring tubes. Finally, we study the changes in the hydrodynamic properties of the medium caused by the plugging of the pores and we study the flow diversion of water caused by plugging of the high permeability zones. @en