How Equivalent Are Equivalent Porous Media?

Geoenergy and geoengineering applications usually involve fluid injection into and production from fractured media. Accounting for fractures is important because of the strong poromechanical coupling that ties pore pressure changes and deformation. A possible approach to the problem uses equivalent porous media to reduce the computational cost and model complexity instead of explicitly including fractures in the models. We investigate the validity of this simplification by comparing these two approaches. Simulation results show that pore pressure distribution significantly differs between the two approaches even when both are calibrated to predict identical values at the injection and production wells. Additionally, changes in fracture stability are not well captured with the equivalent porous medium. We conclude that explicitly accounting for fractures in numerical models may be necessary under some circumstances to perform reliable coupled thermohydromechanical simulations, which could be used in conjunction with other tools for induced seismicity forecasting.

Introduction

Fractures are ubiquitous in rock masses and control their transport and deformation properties (Cornet, 2015; Jaeger et al., 2007). On the one hand, fractures represent preferential paths for fluid flow because their permeability is usually several orders of magnitude higher than that of intact rock (Rutqvist, 2015; Witherspoon et al., 1980; Zimmerman & Bodvarsson, 1996). On the other hand, fractures represent weak planes with lower strength and higher deformability than the intact rock (Bandis et al., 1983; Barton, 1976). The hydraulic and mechanical responses of fractured media are strongly coupled, meaning that the deformation induced by pore pressure changes causes permeability changes as a result of fracture opening or closure, which in turn affects the pore pressure (Barton et al., 1985; Stephanson et al., 1997; Tsang, 1999; Watanabe et al., 2012).

Many geoengineering and geoenergy applications, including enhanced geothermal systems, geologic carbon storage, energy storage, and nuclear waste disposal, take place in fractured media (Rutqvist & Stephansson, 2003) where they induce pore pressure, temperature, and stress changes that are mostly controlled by the physical properties and topology of the fractures. Moreover, the presence of fractures introduces additional nonlinearity and discontinuities in the displacement field. The combined effect of pressure and stress changes affects fracture stability and may induce seismicity (Clarke et al., 2014; Ellsworth et al., 2019; Häring et al., 2008; National Research Council, 2013). The magnitude of induced seismicity is proportional to the rupture area and its frequency distribution obeys a power law called the Gutenberg‐Richter law (Kanamori & Brodsky, 2004).

Including fractures in numerical thermohydromechanical (THM) models is possible and several approaches have been proposed (e.g., Keilegavlen et al., 2021; Pan et al., 2014; T. Wang et al., 2019; Yao et al., 2018; Zareidarmiyan, Salarirad, et al., 2020). Yet, explicitly accounting for fractures in numerical models of poromechanical problems is a challenging task and often comes at the expenses of a very high computational cost (Lei et al., 2017; McDermott & Kolditz, 2006; Pandey et al., 2017; Salimzadeh et al., 2018; Thomas et al., 2020; Yoshioka et al., 2019). In order to facilitate the calculation for large problems, dense fractured networks are often simplified as equivalent porous media, like the ubiquitous joint model, and computational resources are used for other purposes, such as reactive transport, phase changes, and dynamics or complex plasticity models (e.g., Birdsell et al., 2018; Liu et al., 2009; Parisio, Vilarrasa, et al., 2019; Rinaldi et al., 2015; Rutqvist, Wu, et al., 2002; Rutqvist et al., 2005; Taron & Elsworth, 2009; Vallier et al., 2020; Y. Wang et al., 2019).

Here, we investigate the validity of this simplification made by representing a naturally fractured reservoir with an equivalent porous medium. We compare the pore pressure, temperature, and stress changes caused by injection of cold water and production from a distant well between a model that explicitly includes fractures and one that implicitly includes them as an equivalent porous medium. In particular, we assess the induced seismicity potential resulting from both approaches. As a case study, we use the data from Hontomin Technology Development Plant (TDP) for CO2 geological storage. The storage formation at Hontomin is a fractured carbonate rock with a very low‐permeable rock matrix and two perpendicular dominant fracture sets (de Dios et al., 2017). We investigate the validity of equivalent porous media to reproduce the fluid flow and geomechanical response to injection and production into a fractured reservoir.

Materials and Methods

Geometry

We adopt a plane strain geometry for the Hontomin TDP site. The model consists of three horizontal layers that represent the caprock, the reservoir, and the bedrock with thickness of 160, 120, and 50 m, respectively. The reservoir is a fractured carbonate with two perpendicular fracture sets, one subhorizontal and the other one subvertical (de Dios et al., 2017). We represent the two sets with 10‐cm thick planar fractures embedded into the reservoir layer dipping 15° and 75°. Water is injected into the reservoir through a horizontal well on the left‐hand side of the model and produced from another horizontal well on the right‐hand side of the model. The two wells are horizontally spaced by 500 m and a mass flow rate of 5·10−3 kg/s/m is applied (Figure S1). In addition to the fractured reservoir model, we model an equivalent porous medium that yields the same injection pressure in the long term after calibrating the equivalent permeability.

THM Model

Nonisothermal fluid flow in deformable geologic media requires solving the equations of mass, energy, and momentum conservation. Neglecting the inertial term, that is, quasistatic condition, the momentum balance results in: where is the total stress tensor and is the body forces vector. Fluid mass conservation can be expressed as (Bear, 1972) where is the porosity, t(T) is the time, is the fluid density, q(L 3 L −2 T −1) is the volumetric flux, and r(ML −3 T) is the source/sink term of fluid. The volumetric flux is calculated using Darcy's law as where k(L 2) is the intrinsic permeability, is the dynamic viscosity, p(ML −1 T −2) is the fluid pressure, g(LT −2) is the acceleration of gravity, and z[L] is the elevation from a reference plane.

Energy balance is given by (Nield & Bejan, 2006) where is the solid phase density, h (L 2 T −2) is the solid phase enthalpy, h (L 2 T −2) is the fluid phase enthalpy, is thermal conductivity of the geological medium, and is temperature. We assume thermal equilibrium of all phases at every point.

We consider the linear thermoporoelasticity theory in porous media to account for the effect of pore pressure and temperature changes on rock deformation. Elastic strain tensor depends on the total stress, overpressure, and temperature as where is the mean total stress, I (‐) is the identity matrix, is the linear thermal expansion coefficient, ν(−) and E(ML −1 T −2) are the drained Poisson's ratio and the Young's modulus of the rock, respectively, and is the Biot coefficient. The effective stress tensor is .

Fractures are modeled using a finite thickness representation with a higher permeability and lower normal and shear stiffness than the surrounding rock matrix (Zareidarmiyan, Salarirad et al., 2020). For the sake of comparing the fractured media and the equivalent porous medium, we consider that fracture permeability remains constant throughout the whole simulation. Otherwise, an equivalent porous medium cannot be obtained. We assume linear thermoelastic behavior up to failure (and subsequent induced microseismicity) described by the Mohr‐Coulomb failure criterion (Labuz & Zang, 2012) where and are, respectively, the shear stress and the effective normal stress acting on a plane of a given orientation, c(ML −1 T −2) is cohesion, and is the friction coefficient, which is expressed in terms of the friction angle as . Theoretically, the effective stress for failure analyses is the Terzaghi effective stress, that is, with Biot coefficient equal to one (Makhnenko & Labuz, 2015). Because the assumption of = 1 has been made for all rocks (Table S1), the adopted formulation is consistent and does not need to differentiate between the two possible effective stress choices. Considering cohesionless fractures, the mobilized friction coefficient can be defined as the ratio of the shear stress to the effective normal stress

The mobilized friction coefficient quantifies the shear slip tendency and how close the stress state is to shear failure conditions. Shear failure occurs when the mobilized friction coefficient is equal to the actual friction coefficient. As an alternative representation that does not imply a specific orientation in space, one can express the mobilized friction resistance in terms of invariants of the effective stress tensor . Considering the mean effective stress and the deviatoric stress , the mobilized friction angle in a normal faulting stress regime can be expressed as where .

Numerical Models

The model represents a two‐dimensioanl plane strain vertical section of the formation. Initial conditions correspond to a hydrostatic pore pressure with 15 MPa at the center of the reservoir, a geothermal gradient that yields a temperature of 45°C in the reservoir at a depth of 1.44 km, and a normal faulting stress regime, with the horizontal to vertical effective stress ratio being equal to 0.65, that is, . Tables S1 and S2 show the material and fluid properties used in the simulations of the fractured reservoir. For the equivalent porous medium, permeability has been set at 2.75·10−16 m2 and the rest of the properties are equal to those of the rock matrix.

The mechanical boundary conditions are the lithostatic stress on the upper boundary, no displacement perpendicular to the lateral boundaries, and fixed displacement on the lower boundary. The hydraulic boundary conditions are no flow on all boundaries except from the wells where water is injected or produced at a prescribed flow rate of 5.10−3 kg/s/m for 30 years. The wells maintain their location in the equivalent porous medium model. We solve the problem both in isothermal conditions and considering an injection temperature being 20°C colder than that of the reservoir.

An unstructured mesh with 17,775 quadrilateral elements is used. The mesh is refined around fractures, especially at fracture intersections. For the equivalent porous medium, we have used the same mesh as for the fractured reservoir and a more homogeneous mesh with refinement around the wells, obtaining the same results in both cases. Simulations are performed with CODE_BRIGHT (Olivella et al., 1996), an extensively validated finite element numerical code that solves coupled THM problems in porous media. The models are available at Zareidarmiyan, Parisio, et al. (2020).

Material Properties

The reservoir rock properties are taken from the measurements performed on limestone Sopeña Formation from Hontomin injection well (de Dios et al., 2017). The physical properties of Sopeña Formation are highly heterogeneous and the fluctuations are likely connected to the presence of fractures in some of the tested specimens, that is, porosity measurements oscillate from 0.002 to 0.16. In the numerical analyses, we assumed mean values for the elastic and strength properties (Table S1). The permeability of the rock matrix is reported to be ∼10−18 m2, while for the fractures it can go up to ∼10−13 m2. Marly Lias constitutes the caprock (permeability <10−18 m2) at Hontomin and given the lack of information, we have assumed Opalinus Clay as an analogue for the physical properties of the ductile caprock and base rock. Opalinus Clay is a low‐permeable Jurassic shale from northern Switzerland with high CO2 entry pressure values (Makhnenko et al., 2017). Undrained poroelastic parameters and permeability perpendicular to the bedding planes are measured for the fully saturated shale (Makhnenko & Podladchikov, 2018) and are included in the model. Failure properties of shale are adopted from the literature (Gräsle, 2011).

Results

The pressure change is positive or negative at the injection or production wellbore, respectively (Figure 1). The steady‐state conditions in isothermal injection yield the same pressure buildup at the injection well for both the fractured and the equivalent porous medium, although a slight difference is observed at the production highlighting that flow through fractures differs from that of an equivalent porous medium (Figure 1a). The rate of pressure change is more rapid for the equivalent porous medium, which reaches steady‐state conditions within approximately 6 months, while it takes around 1 yr for the fractured medium because of the higher diffusion time across the low‐permeable matrix blocks. The injection overpressure difference between equivalent and fractured medium is greater when cold fluid is injected (Figure 1b) because the effect of the increase in fluid viscosity (Equation S2) on pressure buildup is smaller in the fractured medium as a result of the localized preferential flow paths in comparison to the distributed head losses in the equivalent porous medium. At the production well, the pressure change has not reached a steady‐state condition after 30 yrs of injection (Figure 1b) due to the continuously propagating cooling front. The thermal diffusion process in fractured rock, which includes conduction and advection, takes longer time than the pore pressure diffusion.

Figure 1

Pressure change evolution at the injection and production points considering a fractured reservoir and an equivalent porous medium for (a) isothermal conditions and (b) injection of cold water.

In all cases, the difference between equivalent and fractured medium is small (by design) at the injection and production wells and the equivalent approach does a consistent job in approximating fractures. However, the spatial distribution of pore pressure within the reservoir significantly differs between the fractured medium and the equivalent porous medium (Figures S2 and S3). The difference between the two approaches is as high as 50% when normalized with respect to the pressure change occurring at the injection well, which is the largest pressure change within the reservoir (Figure 2). The difference in pore pressure distribution between the fractured and equivalent approach is concentrated around the fractures during early time of the simulation because pressure has not completely diffused through the matrix blocks (Figures 2a and S2). In contrast, pore pressure diffusion has a smearing effect in the long term, as the difference is more continuously distributed in space and does not conform with the fracture topology, but caused by the fact that the pressure gradient is lower along fractures and is distributed over a longer distance than through the equivalent porous medium (Figures 2b and S2).

Figure 2

Pressure change difference between the fractured reservoir, Δp frac, and the equivalent porous medium, Δp eq, normalized with respect to the pressure difference at the injection well, Δp inj, in (a) early times (30 days of operation) and (b) late times (30 years of operation).

The late‐time (30 years of operation) mobilized friction angle has a similar space distribution for the equivalent and fractured media approaches (Figures 3c and 3d). For cold water injection, assuming a friction angle of 30°, resistance is exceeded at the injection wellbore because of pressure buildup and cooling. The region undergoing failure is entirely localized within the injection formation, although a portion of the interface between the aquifer and the caprock is exceeding frictional resistance in both cases. The smeared nature of the equivalent porous medium approach does not distinguish between critically and noncritically oriented fractures (Figures 3b and 3d). On the contrary, the enhanced information of the fractured medium approach shows that the mobilized friction angle along subvertical fractures is exceeding resistance, whereas the subhorizontal fracture set remains stable (inset Figures 3a and 3c).

Figure 3

Mobilized friction angle development due to cold water injection at early times (30 days) for (a) The fractured reservoir and (b) The equivalent porous medium, and at late times (30 years) for (c) The fractured reservoir and (d) The equivalent porous medium.

The differences between the two approaches in terms of fracture stability can be highlighted by plotting the mobilized friction coefficient in all the considered cases and along fractures A, B, and C (Figure 4). Fracture A is connected to the injection wellbore and, at isothermal conditions, the mobilized friction coefficient is the highest closer to the wellbore in the equivalent porous medium, while it becomes smaller than one of the fractured mediums at greater distances (Figure 4a). Cold fluid injection causes the increase in the mobilized friction coefficient in both media, with values in the equivalent medium potentially exceeding failure, while they never exceed the threshold of μ = 0.6 in the fractured medium (Figure 4d). Note that fracture A is a subhorizontal fracture, thus it is not critically oriented in a normal faulting stress regime like the one considered in this study. Contrary to fracture A, fracture B is close to the critical orientation (Figures 4b and 4e) and the mobilized friction coefficient after injection is higher for the fractured medium. In this case, the mobilized friction coefficient remains below failure threshold for the isothermal injection (Figure 4b), while the cooling effect reduces the mean effective stress and mobilized friction coefficient has a sharp increment over large patches. Fluid production has a beneficial effect because the pressure reduction with time stabilizes fracture C (Figure 4c), with a more pronounced effect in the fractured medium. An inversion is observed between short‐ and long‐time behavior during cold injection and mobilized friction coefficient slightly increases close to the production point (Figure 4f). Nevertheless, no failure occurs and the mobilized friction coefficient in the noncritically oriented fracture C never exceeds failure conditions.

Figure 4

Mobilized friction coefficient, that is, deviatoric stress to effective normal stress ratio, along fractures A, B, and C for early times (30 days) and late times (30 yrs) for isothermal conditions (a)–(c) and injection of cold water (d)–(f) when considering a fractured reservoir and an equivalent porous medium.

Discussion

Equivalent porous media can misleadingly reproduce the pore pressure evolution measured at lower‐dimensionality observation points of a fractured medium as a result of injection or production (Figures 1 and 2). Even though an equivalent porous medium can be used to reproduce fluid flow through fractured media under certain conditions, like for densely fractured media (Long et al., 1982; Scanlon et al., 2003), our simulations show that the pressure distribution within the fractured medium can differ as much as 50% of the pressure change undergone at the wells. Such difference occurs both in the short term, when pore pressure has not completely diffused through the rock matrix, and in the long term, when pore pressure diffusion has already occurred. As a result, heat transport also differs between an equivalent porous medium and the fractured medium.

Fractures are usually more permeable and consequently heat advection dominates along fractures, whereas heat conduction prevails within the matrix: The temperature change is larger along fractures in the fractured medium compared to the one in the matrix in the equivalent porous medium (Figure S4 and S5). The temperature changes are restricted to a small rock volume around the injection well because of the lower thermal diffusivity compared to the hydraulic diffusivity (De Simone et al., 2017). The differences in the pore pressure and temperature distributions highlight the limitations of the equivalent porous media to accurately reproduce fluid flow and heat transport within a fractured rock.

It should be noted that we have neglected the effect of changes in fracture aperture on their permeability. Including this characteristic, that is, permeability enhancement due to fracture opening when inertial losses become negligible (Rutqvist, 2015; Zhou et al., 2019), would have significantly increased the differences in pore pressure (Birdsell et al., 2018) and temperature distribution between the fractured and the equivalent porous medium. Furthermore, due to the nonlinearity of the problem, an equivalent permeability that reflects the pore pressure evolution at the injection and production wells would not be possible to reproduce.

Induced (micro) seismicity assessment is significantly affected by the differences in the pore pressure and temperature distributions, even when considering constant fracture permeability (Figure 4). In addition to the potential local stress variations around fractures previous to injection/production (Gao et al., 2019; Lei & Gao, 2018), pore pressure and temperature changes induce local rotation of the principal stress tensor along fractures that can be reproduced only when fractures are explicitly included in the model (De Simone et al., 2013; Zareidarmiyan et al., 2018, Zareidarmiyan, Salarirad, et al., 2020). The lower stiffness of fractures compared to the rock matrix (Pyrak‐Nolte & Morris, 2000) generates a lower cooling‐induced stress reduction and affects both the induced poro and thermomechanical stresses and thus, fracture stability.

Measurements of the physical properties of intact and fractured rock in the lab at representative pressure and temperature provide fundamental insights on multiphysical processes but are difficult to perform (Braun et al., 2020; Faoro et al., 2009; Samuelson et al., 2009). Additionally, fractures introduce a length‐scale dependency and lab‐sized results must be properly upscaled to reservoir conditions (Parisio, Tarokh, et al., 2019). In this regard, recent experiments at intermediate scale with comprehensive monitoring systems can provide insight on behavior of decameter‐long fractures and bridge various scales (Amann et al., 2018; Guglielmi et al., 2020). Fluid‐solid interactions and chemical reactions also play an important role and in the case of CO2 injection can significantly change fracture properties (Pluymakers et al., 2014), and the additional effects on induced seismicity are to be further studied in the lab (Elsworth & Yasuhara, 2006; Li et al., 2007; Samuelson & Spiers, 2012).

Numerical modeling of fractured media also possesses some challenges. Computationally, at the characteristic length scale of many reservoir applications, the fracture density can become excessively high and including all the discontinuities in numerical models may become unfeasible. Homogenization and other upscaling techniques are often the only solution for overcoming this issue (Bonnet et al., 2001; Gan & Elsworth, 2016; Lei et al., 2015). However, our numerical simulation results show that the information loss from neglecting the explicit representation of fractures at a spacing that is in the order of the reservoir's size can lead to unacceptable approximation on stability.

The choice of using either the fractured or equivalent media approach should be carefully pondered on a case‐by‐case basis and should heavily rely on available structural geology information. At very high fracture density, the equivalent approach and other homogenization‐based schemes would be more convenient and, assuming the rock mass representative elementary volume to contain enough fractures, potentially almost no information would be lost. Modeling bedding planes in shale, which are tightly spaced discontinuities that result in overall anisotropy, is an example where continuous approaches (Parisio et al., 2018) yield similar predictions to high‐resolution discrete ones (Lisjak et al., 2015; Kim et al., 2020). Multiscale methods have the potential to bridge the gap between the explicit fracture modeling and equivalent continuum approaches and, along with machine learning and data‐driven computational schemes, could resolve the large spatial scale separation typical of fractured reservoirs (Wang & Sun, 2018). Nevertheless, the high computational costs required from modeling explicit fractured media must be justified by a proper knowledge about the characteristics of the subsurface fracture network. Otherwise, a high uncertainty on the fracture network geometry and characteristics hinders the effort of highly realistic models.

Conclusions

We consider a problem of fluid injection and production from a fractured reservoir and compare pore pressure and temperature distributions between a fractured medium and its equivalent porous medium and how they affect fracture stability. An equivalent porous medium that reproduces the pressure evolution at the injection and production wells can be obtained, provided that fracture permeability remains constant. However, the pore pressure distribution significantly differs between the two approaches as a result of the higher permeability of fractures. Differences are also evident when evaluating fracture stability. Therefore, achieving an accurate simulation of the THM processes occurring in fractured media requires explicitly including fractures in the reservoir models whenever they introduce inhomogeneities at the reservoir scale.

Supporting Information S1

Click here for additional data file.