Skip to main content

Science – Society – Technology

A fracture flow permeability and stress dependency simulation applied to multi-reservoirs, multi-production scenarios analysis

Abstract

The use of the subsurface and the exploitation of subsurface resources require prior knowledge of fluid flow through fracture networks. For nuclear waste disposal, for the enhancement of hydrocarbon recovery from a field, or the development of an enhanced geothermal system (EGS), it is fundamental to constrain the fractures and the fracture network. This study is part of the GEMex project, an international collaboration of two consortia, one from Europe and one from Mexico. The research is based on exploration, characterization and assessment of two geothermal systems located in the Trans-Mexican volcanic belt, Los Humeros and Acoculco. In Acoculco, two wells reached very high temperatures, but did not find any fluids. For that reason, the Acoculco Caldera is foreseen as an EGS development site, hoping to connect existing wells to a productive zone. This implies that the fluid flow through the geothermal reservoir would be mainly fracture dominated. This study investigates the dependency of fracture permeability, constrained by fracture lengths and apertures, with stress field conditions. Simulations are computed in 2D, using COMSOL Multiphysics\(^{\textregistered }\) Finite Elements Method Software, populated with mechanical data obtained in the rock physics laboratory and with dense discrete fracture networks generated from 1D scanline surveys measured in Las Minas analogue outcrops for Acoculco reservoir. The method offers a prediction for multiple scenarios of the reservoir flow characteristics which could be a major improvement in the development of the EGS technology.

Introduction

Background

The understanding of fractures is crucial for several subsurface activities such as oil and gas, nuclear waste disposal, enhanced geothermal system or CO2 sequestration (Fox et al. 2015). In the case of fluid flow circulation, fractures can act either as flow boundaries or fluid pathways.

The first general question here is to evaluate to what extent fractures in subsurface can be characterized. Two main methods are commonly used to predict the impact of the fractures: one method is using a discrete fracture model (DFM) and the second method is using a continuum described as the representative volume element (RVE) (Geiger and Matthäi 2012; Surrette 2006). In the first case, using DFM, all fractures are explicitly described as boundary elements. Therefore, a DFM simulates the geometry of real fractures and is a closer representation of reality (Li et al. 2013). Computationally, handling a DFM drastically increases the number of elements and the complexity in the mesh. Because of that, the calculation takes longer to converge to a solution (Garipov et al. 2016; Zhang et al. 2007). In the second case, using RVE consists of modelling the bulk impact of the fractures on the reservoir properties. These properties are averaged to represent the matrix and the fractures together in a continuum model (Flemisch et al. 2018; Singhal and Gupta 2010; Surrette 2006; Dershowitz et al. 2004). Because working with DFM is computationally very expensive, it is more common to use the RVE approaches to model fluid flow through fractured media. However, working with DFM allows a higher level of details, such as separating fractures properties based on the structural fracture sets. Additionally, a DFM allows the coupling of the mechanical properties and stress field to the flow properties of the fracture network. The characterization of the fractures is, therefore, of utmost importance to the results.

Study context

In this study, a combination of a field study fracture characterization in Acoculco, with porosity, permeability, rock strength laboratory tests, stress field and coupled thermal-hydraulic (TH) modelling is carried out. A streamlined integration of a DFM generated from field data, mechanical rock properties of field samples measured in the rock physics laboratory, stress field data and simulation of a geothermal field production has not been presented before to the best of our knowledge.

This study is in the context of analysing the feasibility of developing an enhanced geothermal system (EGS) in the Acoculco geothermal field, which is being studied within the GEMex project. Acoculco is located in the Trans-Mexican volcanic belt (TMVB), about 100 km northeast of Mexico City. Two exploration wells were drilled, EAC1 and EAC2, at about 500 m apart, in 1995 and 2008, respectively. Their lithological sequence is similar and can be summarized as described by Viggiano-Guerra et al. (2011) and López-Hernández et al. (2009), in Table 1.

Table 1 Fracture sets, colors and orientations

The carbonate section regroups from dolomitic limestones to carbonates metamorphosed into marbles or metasomatized into skarns due to the intrusion of a granodioritic pluton. Acoculco lithological reservoir section is composed of limestones, marbles and skarns. The main assumption here is that these reservoirs are not permeable because the fractures are not connected enough to allow a meaningful production of fluid through the fractures. The problem is then to quantify the role of the fractures for these three formations. This implies an understanding of the fracture properties in terms of mechanics and their impact regarding the fluid circulation (i.e. determining whether fractures behave as fluid barriers or fluid pathways). How are fractures spatially distributed in the formation? How are these fractures connected? To what extent does the stress field deform the aperture of the fracture (DAF) and, therefore, how much does it influence the fracture permeability?

Innovation

To answer this question, we propose a numerical method to calculate the fluid flow and heat exchange in a fracture controlled reservoir. The models are built using real data from the field. The discrete fracture models (DFMs) are obtained from scanline surveys measured in the field on Las Minas analogue outcrops and processed using the SkaPy script and then extrapolated using the multiple point statistic method as presented in Bruna et al. (2019). The rock properties have been measured in the rock physics laboratory on rocks sampled on the analogue outcrops of Las Minas. The numerical model runs two studies: (i) a static analysis measuring the influence of the stress field on the fracture aperture; (ii) a transient analysis which couples fluid flow in porous media and thermal exchange through the DFMs. Because Acoculco reservoir is composed of limestone, marble and skarn, a comparison is made between these three formations. Therefore, we created three models representing the three formations. Each formation is populated with its own mechanical properties and its own DFM.

The main innovations of this study are to cover a geothermal doublet reservoir scale of 600 × 600 m, with a very dense fracture system, with an amount of tens of thousands of fractures (exact amount is given in Table 2), and material properties with very low rock matrix permeability (model populated with a matrix permeability of \(1 \times 10^{-17}\,\text {m}^{2})\). The second innovation concerns the method implemented to be able to use these real field data. We solve this problem by looking at 2D (xy) models and then we solve the mechanical part in a steady state study. Then, we calculate the fluid flow and thermal exchange using the deformed aperture fracture (DAF) from the steady state study. Finally, a field development perspective is applied by simulating production over 100 years. For that purpose, the model is implemented in multiple reservoir geometries and multiple wells positioning. This method helps to evaluate which formation, and more precisely, which part of each formation, is best suited for the development of an EGS. For each of the three analysed reservoirs, four development scenarios are considered. Acoculco case study is given to illustrate our results.

This multi-scenario model is simulated using the COMSOL Multiphysics software for Finite Elements Method (FEM) analysis. All data are freely available on: https://github.com/BatLep/Sim_THM_Multi_Res-Sc.git.

Table 2 Fracture numbers, and initial apertures per fracture sets and per formation

Methods: simulation model implementation

General implementation

Assuming a 2D horizontal domain D, composed of the rock matrix defining a first sub-domain \(\mathbf {\Omega }\), and several fractures that all together constitute a second sub-domain \(\mathbf {\Gamma }\), such as \(\mathbf {D = \{ \Omega + \Gamma \} }\). The solid mechanics for 2D plane stress and plane strain, the subsurface fluid flow in porous media and the heat transfer are solved using the Finite Elements Method (FEM), implemented in the COMSOL Multiphysics\(^{\textregistered }\) software. The simulation computes two studies: the first one is steady state and solves the solid mechanics; the second one is transient and solves the fluid flow and thermal exchange. The time-dependent solver used the implicit Backward Differentiation Formula (BDF), and especially the Backward Euler, where the fraction of the initial step is set at 0.001. The rock matrix sub-domain \(\Omega \) and the linear elements for fractures in 2D model sub-domain \(\Gamma \) are discretized into triangular elements following a conform mesh which honors the geometrical characteristics of the fracture system as defined by Flemisch et al. (2018). The fractures are modelled as “single layers material” (zero thickness) to which specific material properties are applied, independently of the rock matrix \(\Omega \) material properties. Boundary conditions are applied to the outer boundaries of the domain D such that no displacement is allowed, but thermal transfer and fluid flow are free along these edges. Injection and production of fluid from the well bores are idealized as points and placed on a fracture for calculation simplicity.

Governing equations

Study 1: solid mechanics, steady state

For the solid mechanics analysis, the physical descriptions are based on the laws for the balance of forces and the constitutive relations that relate the stresses to strains (Jaeger et al. 2007; Zoback 2007; Fossen 2013). In this study, a relation between the fluid flow and the stress field is induced using the fracture aperture. The effect of fracture slippage and shearing and consequently its influence of fracture permeability is not considered in this study. The stress field acts as a load on the fracture aperture. The deformed aperture of the fracture (DAF) follows the linear elasticity behaviour of a spring as described with the constitutive relation of the Hooke’s law as

$$\begin{aligned} f_\text{s} = -K(u-u_0) \end{aligned}$$
(1)

where \(f_\text{s}\) is a force/unit area, u is the displacement deforming the spring, and K is the stiffness matrix. \(u_0\) is an optional offset, which describes the stress-free state of the spring. As in the case of analysing the fracture elasticity, the stiffness is a function of the fracture material properties and the fracture width (aperture) \(d_\text{f}\). The stiffness in the normal direction is computed based on a state of plane strain, so that

$$\begin{aligned} k_\text{n} = \frac{\textit{E}(1-\nu )}{d_\text{f}(1+\nu )(1-\nu )} \end{aligned}$$
(2)

where \(k_\text{n}\) is the normal stiffness (or spring constant), E is the Young’s modulus and \(\nu \) is the Poisson’s ratio. The assumption of plane strain conditions is relevant when the material of the elastic layer (in this case the fracture material) is softer than its surroundings (in this case the rock matrix).

The regional stress field is then assigned by creating a second reference system defined by the principal stress vector system. Using the solid mechanics physics analysis, initial stresses and strains are assigned to the whole domain D as the far field stress S, with

$$\begin{aligned} S=\begin{bmatrix} S_{\text{Hmax}}&0&0 \\ 0&S_{\text{hmin}}&0 \\ 0&0&S_\text{v} \end{bmatrix} \end{aligned}$$
(3)

where \(S_{\text{Hmax}}\) is the maximum horizontal stress, \(S_{\text{hmin}}\) is the minimum horizontal stress, \(S_\text{v}\) is the vertical stress and S is the stress tensor.

Study 2: coupled subsurface flow and thermal exchange, transient

The second study is a transient analysis simulating over 100 years. The heat transfer is coupled to the subsurface fluid flow. In fluid flow, heat transfer, and mass transfer, the simulations are solved based on the laws for conservation of momentum, mass, and energy. The mass transfer is governed by Darcy’s law and heat transfer is governed by Fourier’s law of heat conduction. As the fluid flow in the reservoir is mainly controlled by the fractures, a distinction can be made between the fracture permeability \(k_\text{f}\) and the intergranular rock matrix permeability \(k_\text{m}.\) The fracture permeability is calculated using the Cubic law (Snow 1969; Zhang et al. 2007; Singhal and Gupta 2010), as

$$\begin{aligned} k_\text{f}= \frac{d_\text{f}^2}{12} \end{aligned}$$
(4)

where \(k_\text{f}\) is the fracture permeability and \(d_\text{f}\) is the aperture of the fracture.

Model scenarios

The model is built with a core frame that can be populated with different properties based on conditionals. Hence, this model handles multi-analyses of multiple formations with multiple well-positioning scenarios.

Application: field example

Input data, and domain delimitation

\(\Omega \), matrix material properties

Determination of the material properties was measured in the Geoscience and Engineering Laboratory of the Civil Engineering and Geosciences department of Delft University of Technology. All tests have been run on sample taken from the analogue outcrops of Las Minas.

A total of 12 unconfined compressive strength (UCS) tests have been run to determine the elastic properties of the material, such as Young’s modulus (E) and Poisson’s ratio (\(\nu \)). The apparatus is uni-axial stress/strain device, with a capacity \(500\,\text{kN}\). Axial displacement was recorded using two LVDTs and radial displacement was recorded using a radial chain with LVDT sensor around the cylindrical samples of \(75\times 30\,\text {mm}\) size. Measurements and calculations were realized following ASTM guidelines.

Porosity measurements have been performed with the Ultrapycnometer 1000 Version 2.12 (Quantachrome Corporation). Measurements were repeated 20 times per sample and each rock type was tested on 3 to 5 specimens, depending on available material.

Rock permeability measurements were performed at German Research Centre for Geosciences (GFZ) at conditions corresponding to room temperature on oven-dry sample cores using the steady state gas flow method described by Tanikawa and Shimamoto (2009). The cores had a diameter of \(30\,\text {mm}\) and variable lengths between 15 and \(30\,\text {mm}\). At least five pore pressure levels, \(P_\text{up}\), ranging from 0.5 to \(3.5\,\text {MPa}\) were applied to the upstream side of the samples, at a constant isotropic confining pressure of \(8.0\,\text {MPa}\), with Ar as both confining and pore pressure medium. The corresponding gas flow rates through the samples were measured with one of the four different commercial flowmeters (MKS Instruments) that can measure flow rates in the range of 10, 100, 1000, and \(10,000\,\text {cm}^3\,\text {min}^{-1},\) respectively. As the gas is released from the sample to atmospheric pressure, we assume a constant pore pressure at downstream side, \(P_\text{down},\) of \(0.1\,\text {MPa}.\) The apparent gas permeability was calculated after Scheidegger (1974) as

$$\begin{aligned} k_\text{gas}=\frac{\eta L Q 2 P_\text{down}}{A\left( P^2_\text{up}-P^2_\text{down}\right) } \end{aligned}$$
(5)

where \(\eta \) is the viscosity of the pore fluid at given pressure and temperature, L is the sample length, Q is the flow rate and A is the cross-sectional area of the sample. The water permeability was derived from the gas permeabilities that were Klinkenberg corrected for the gas slip flow effect at the gas–mineral interface (Klinkenberg 1941):

$$\begin{aligned} k_\text{gas}=k_\text{w} \left( 1+\frac{2b}{P_\text{up}+P_\text{down}}\right) \end{aligned}$$
(6)

where \(k_\text{w}\) is the (intrinsic) water permeability and b is the Klinkenberg factor. In the low pressure range, gas permeabilities were significantly enhanced due to ‘slip flow’. Thus, gas permeabilities measured at differential pore pressures lower than \(0.75\,\text {MPa}\) were predominantly ignored for the conversion and water permeabilities were obtained from linear plots of gas permeabilities vs. the inverse of at least four different mean pore pressures. The measurement accuracy was better than 1.5%.

All parameters used as input in the model are given in Table 3.

Table 3 Input parameters for simulation

\(\Gamma \), the DFMs

The fractures are implemented in the model using existing DFMs. The fractures were measured in the field using the linear scanline reporting method on several outcrops. The trace of the scanline survey follows the direction of the outcrop walls, offering to capture different angles and, hence, a better prediction of the fracture distribution. The scanline surveys are then processed using the SkaPyFootnote 1 Python\(^{\mathrm{TM}}\)Footnote 2 script, together with the multiple point statistic method (Bruna et al. 2019), to extrapolate the fracture distributions into a larger spatial domain. As mentioned in "Study context", the reservoir formation of Acoculco geothermal site is composed of limestone, marble and skarn. Hence, each of these formations has its own DFM. These DFMs are separated in fractures sets based on their strike and dip values, related to the regional structural trends described in the literature (Campos-Enriquez and Garduño-Monroy 1987; López-Hernández et al. 2009; Carrasco Núñez et al. 2017; Norini et al. 2015). For each of the three formations treated here, we calculated the statistical distributions of the fracture apertures from our field measurements. As a result, we obtain for each formation and, for each fracture set, one specific aperture value. All these DFMs represent a 2D horizontal plane covering \(600\times 600\,\text {m}.\) These DFMs represent the spectrum of fracture distributions in between the two outcrops measured with the scanline surveys for each rock type. Therefore, the domain of (600 × 600) \(\text {m}^{2}\) is divided into three equivalent sub-domains of (200 × 600) \(\text {m}^{2}.\) The first sub-domain is populated using scanline survey dataset, the second is with the other dataset and the third in the middle is populated with the calculated interpolation. This explains why west and east parts of the domain present different fracture distributions.

More detail is given in Table 2 regarding the quantity of fractures per set and the value given for their aperture.

These values of aperture correspond to the fracture aperture at the surface, measured at the outcrop. To correct the aperture of the fractures at the reservoir depth, we apply stress dependency to the fracture aperture. This correction first requires an estimation of the mechanical properties of these fractures. The aperture \(d_\text{f}\) is a function of the spring constant \(k_\text{n},\) itself a function of the Young’s modulus E and the Poisson’s ratio \(\nu, \) According to Jeanne et al. (2017) and Rutqvist et al. (1998), the following can be safely used:

$$\begin{aligned} E_{\text{fracture}} &= 0.1 E_{\text{matrix}} \nonumber \\ \nu _{\text{fracture}} &= 0.4 \nu _{\text{matrix}} \end{aligned}$$
(7)

Stress field

The stress field conditions are constrained using the World Stress Map (Heidbach et al. 2016), with two data points, at stations (19,79; − 98,47) and (19,97; − 98,63) from “drilling induced fractures”, about 40 km away from the wells. Both stations indicate a normal fault stress regime, giving an average maximum horizontal stress \(S_{\text{Hmax}}\) of azimuth N55°. Vertical stress is estimated to be hydrostatic according to \(S_\text{v}=\rho g z\). Being in a normal stress regime, \(S_\text{v}> S_\text{Hmax} >S_{\text{hmin}}.\)

Problem definition and simulation implementation

In this study, we analyse three potential reservoirs, here named the limestone, the marble and the skarn formations. Each of these formations is modelled using a rock matrix domain \(\Omega \) populated using its own mechanical parameters. For each of these formations, a specific DFM is implemented and discretized. As these DFMs are varying spatially, we want to analyse their influence on production performances. Therefore, we simulate four scenarios of well doublets located at different parts of the domain D. The model implementation is summarized in Fig. 1.

Fig. 1
figure 1

Model implementation (from top to bottom): \(\Omega \), the material properties; \(\Gamma \), the DFMs; assignment of the stress field; well production scenarios

Results

Aperture dependency to the stress field

When assigning the stress field to the material, we expect to see the fracture width (aperture) decreasing depending on the orientation of the fracture compared to the direction of the stress field. Figure 2 gives the orientation of the seven fracture sets considered in this study. In the Acoculco area, as mentioned in "Application: field example", the stress field is oriented N55°. Therefore, the fracture sets F2 and F5 are expected to be the most affected by a decrease in aperture, while the fracture sets F1 and F7 should experience less deformation.

Fig. 2
figure 2

Stress field orientation and fracture sets: a stress field where \(S_{\text{Hmax}}\) is N55°; b stress field where \(S_{\text{Hmax}}\) is N90°

In Fig. 3, we compare the fracture sets deformation (from initial aperture measured at the outcrop, corresponding to 100%) with using two different stress field orientations: in the first case, we are using the Acoculco N55°; in the second case, we are using an East–West direction (N90°). The use of this East–West stress field is only meant to emphasize the aperture dependency to stress and observe the influence of stress field direction on the calculated DAF. The plots represent, for each formation, the percentage of aperture deformation. The results confirm the above-mentioned expectation: in all the formations, when assigning an N55° stress field, the fracture sets F2 and F5 are the most deformed and F1 and F7 the least. Similarly, when using the N90° stress field orientation, the sets F3, F5 are the most deformed and F1, F6, F7 the least.

Fig. 3
figure 3

Fracture aperture deformation, for each fracture set, in: the limestone formation; the marble formation; the skarn formation (The vertical scale shows the deformation as a percentage. Fractures sets F4 and F7 are not present in the marble formation)

Production scenarios

In this part, the field production performance based on spatial well positioning is analysed. Due to the fact that DFMs are not homogeneous (Fig. 1), the positioning of the wells can have a significant impact on the production. For instance, the fracture density is much higher on the western part of the DFMs, especially for the limestone and marble formations. Because of that, one would expect a higher fracture connectivity and, therefore, a higher bulk permeability. Here, the performance is evaluated as the estimated fluid flow between the two wells, injection and production. As the in-flow and out-flow at the wells are constrained by the boundary conditions of the model, we look at the lowest \(\Delta P\) as a representation of the best flow. As a consequence, the scenario 1 would be expected to have the best performance. Because the scenario 4 implies a fluid circulation across the fracture density variation, we would expect a medium performance compared to the others. The first result shown in Fig. 4 is that the well-positioning scenario 4 results in a medium performance compared to the other three scenarios. The second information is that the scenario 1 in the limestone formation gives a good performance, as expected. However, this in not the case in the marble formation even though the fracture density is higher in the western part. Referring back to the marble DFM (Fig. 1), even though the fracture density is lower on the eastern part, the fracture length is higher and, therefore, the degree to which these fractures are connected is much higher. This explains why the scenario 3, in both the marble and the skarn, has the best performance in terms of required \(\Delta P.\)

Fig. 4
figure 4

Pressure (\(\Delta P=W_{\text{prod}}-W_{\text{inj}}\)) and temperature profiles probed near injection and prodcution wells, for different injection scenarios (Sc1–Sc4), using Acoculco local stress field (N.55)

Heat transfer

Looking at the heat transfer (Fig. 5), a few observations can be made: (i) despite the high fracture density and the permeability trends created by the fractures, the thermal front is circular and, hence, homogeneous. This can be explained by the poor continuity and connectivity of the fractures. (ii) However, even though the fractures (\(\Gamma \)) are not creating a continuous fluid pathway, they are always showing lower temperatures than the matrix (\(\Omega \)). This confirms that the fluid flow propagates quicker in the fractures and, hence, corroborates higher permeability in the fractures. (iii) Considering the implemented conditions of pumping flow rate and the resulting bulk permeability of the reservoir: for this example of the limestone, the injected fluid, pure water at 50 °C, needs a 100 m to reach the initial reservoir temperature of 300 °C.

Fig. 5
figure 5

Map view model of heat transfer over time in the limestone formation, with well scenario 4; a \(\text {time}=10\ \text {years}\); b \(\text {time}=50\ \text {years}\); c \(\text {time}=100\ \text {years}\)

Discussion

In this study, we relate the fluid flow and the stress field by calculating the deformation on the fracture aperture. The DAF follows the linear elasticity behaviour of a spring as described with the Hooke’s law. Therefore, the DAF is directly related to the spring constant \(k_\text{n}\), itself depending on the mechanical material properties, such as Young’s modulus E and the Poisson’s ration \(\nu \). Figure 6 shows the strong relation between \(k_\text{n}\) and the DAF and, hence, the direct relation between the elastic material properties (\(E, \nu \)) and the deformation. In this case, we see a threshold where \(k_\text{n} = 1\times 10^{12}\,\text {Pa}\,\text {m}^{-1}\). When \(k_n < 1 \times 10^{12}\,\text {Pa}\,\text {m}^{-1}\), the deformation happens. If \(k_n> 1 \times 10^{12}\,\text {Pa} \,\text {m}^{-1}\), almost no deformation (i.e. closure) occurs.

Fig. 6
figure 6

The spring constant, and how mechanical properties influence fracture deformation (DAF) under stresses

The geothermal site studied here is foreseen as a potential EGS. For that reason, the permeability is very low even when including the fractures. As a consequence, we do not see any thermal breakthrough at the production well. Increasing the flow rate would show a thermal breakthrough. However, we assume that having a pressure difference (\(\Delta P\)) higher than 10 MPa would be unrealistically too high as it would eventually induce some hydraulic fracturing. Unfortunately the mechanics of hydraulic fracturing is not implemented here. For that reason, we maintain the \(\Delta P\) below or around 10 MPa.

Increasing the fracture trace length would improve the fracture network connectivity and the fracture aperture would necessarily increase the global permeability (Fig. 7).

Fig. 7
figure 7

Example of simulation with connected fractures and high fracture aperture

In further works, we would investigate the hydraulic fracture stimulation prediction. Thanks to that we can implement an extra dimension to this analysis by looking at fracture connectivity and fluid flow performance before and after well stimulation.

Another characteristic which influences the model results is the thermal conductivity of the material. The matrix and the fracture materials are not always made of the same mineral composition. Therefore, the thermal conductivity in the fracture could be much higher than that in the rock matrix (Fig. 8).

Fig. 8
figure 8

Thermal conductivity contrasts in a test model

From the current results, without stimulation, the scenario 1 has the best performance in the limestone formation, while for the marble and the skarn formation, the scenario 4 would be the best.

Finally, Tables 4 and 5 give the computation information per reservoir and scenario, such as number of meshed elements, degrees of freedom, and computation memory or computation time.

Table 4 Number of elements in mesh and degrees of freedom
Table 5 Computational physical and virtual memories and computation time

The computation time so far is relatively short. The next step to be fully challenging the RVE approach would be to perform the same simulations on 3D models. The fracture prediction would be even more accurate and would certainly influence the fluid flow circulation.

Conclusion

When developing subsurface activities such as oil and gas, nuclear waste disposal, CO2 sequestration or, as in this case, enhanced geothermal system, it is fundamental to quantify the role of the fracture system present in the subsurface. Very often, fractures are up-scaled to represent the matrix and the fractures together in a continuum model. In this article, we present a stress-dependent fracture aperture model and, hence, a stress field-dependent fluid flow and heat transfer model, using field data to populate the material properties and the fracture networks. These DFMs are composed of tens of thousands of fractures. The DFMs are separated into fracture sets related to the regional structural context. Each fracture set is characterized with its own aperture value from field measurements, corrected to depth stress conditions. Thanks to these simulations, we can identify the threshold, on the mechanical properties of the fractures, below which deformation takes place. The purpose of this multiple scenarios analysis is to evaluate the risks of the project. Combining multiple formations and multiple well positioning scenarios gives a thorough evaluation of the reservoir performance. This is fundamental in the context of field production risks analysis. Based on our results, the safer scenario for this project would be developed according to the scenario 4, which simulates a well doublet across the different fracture distribution patterns instead of targeting one zone in particular. This is even more important in the case of developing the Marble and the Skarn formations. This method provides, to our knowledge, a more realistic model of the existing and expected fracture network at depth. This would certainly be a major improvement in the development of the EGS technology.

Notes

  1. https://github.com/BatLep/SkaPy.git.

  2. Python Software Foundation (2019).

Abbreviations

EGS:

enhanced geothermal system

FEM:

Finite Elements Method

DFM:

discrete fracture model

RVE:

representative volume element

TH:

thermo-hydraulic (modelling)

GEMex:

Geothermal Mexico (http://www.gemex-h2020.eu/)

TMVB:

Trans-Mexican volcanic belt

MD:

meter drilled

TD:

total depth

DAF:

deformed aperture of the fracture

D :

whole simulated domain

\(\Omega \) :

matrix domain

\(\Gamma \) :

fractures domain

\(f_\text{s}\) :

force/unit area

u :

displacement

K :

stiffness matrix

\(d_\text{f}\) :

fracture aperture

\(k_\text{n}\) :

spring constant, stiffness in the normal direction

E :

Young’s modulus

\(\nu \) :

Poisson’s ratio

S :

far field stress tensor

\(S_{\text{V}}\) :

vertical stress

\(S_{\text{Hmax}}\) :

maximum horizontal stress

\(S_{\text{hmin}}\) :

minimum horizontal stress

\(k_\text{f}\) :

fracture permeability

\(k_\text{m}\) :

matrix permeability

UCS:

unconfined compressive strength

LVDT:

Linear Variable Differential Transformer

L :

sample length for rock physics laboratory tests

A :

the cross-sectional area of the sample

Q :

flow rate for rock physics laboratory tests

\(\mu \) :

fluid viscosity

References

  • Bruna PO, Straubhaar J, Prabhakaran R, Bertotti G, Bisdom K, Mariethoz G, Meda M. A new methodology to train fracture network simulation using multiple-point statistics. Solid Earth. 2019;10(2):537–59.

    Article  Google Scholar 

  • Campos-Enriquez JO, Garduño-Monroy VH. The shallow structure of Los Humeros and Las Derrumbadas geothermal fields, Mexico. Geothermics. 1987;16(5–6):539–54.

    Article  Google Scholar 

  • Carrasco Núñez G, Hernández J, De León L, Dávila P, Norini G, Bernal JP, Jicha B, Navarro M, López-Quiroz P, Digitalis T. Geologic Map of Los Humeros volcanic complex and geothermal field, eastern Trans-Mexican Volcanic Belt Mapa geológico del complejo volcánico Los Humeros y campo geotérmico, sector oriental del Cinturón Volcánico Trans-Mexicano. Terra Digitalis. 2017;1(2):1–11.

    Article  Google Scholar 

  • Dershowitz WS, La Pointe PR, Doe TW. Advances in discrete fracture network modeling. In: Proceedings of the US EPA/NGWA fractured rock conference 2004, 2004. p. 882–94.

  • Flemisch B, Berre I, Boon W, Fumagalli A, Schwenck N, Scotti A, Stefansson I, Tatomir A. Benchmarks for single-phase flow in fractured porous media. Adv Water Resour. 2018;111:239–58.

    Article  Google Scholar 

  • Fossen H. Structural geology. No. c. Cambridge University Press. 2013. https://www.cambridge.org/core/books/structural-geology/4F99FE5530610FB31EAA89FB0E9865EF.

  • Fox DB, Koch DL, Tester JW. The effect of spatial aperture variations on the thermal performance of discretely fractured geothermal reservoirs. Geotherm Energy. 2015;3(1):1–29.

    Article  Google Scholar 

  • Garipov TT, Karimi-Fard M, Tchelepi HA. Discrete fracture model for coupled flow and geomechanics. Comput Geosci. 2016;20(1):149–60.

    Article  Google Scholar 

  • Geiger S, Matthäi S. What can we learn from high-resolution numerical simulations of single- and multi-phase fluid flow in fractured outcrop analogues? Geol Soc Lond Spec Publ. 2012;374(1):125–44.

    Article  Google Scholar 

  • Heidbach O, Rajabi M, Reiter K, Ziegler M, Team WSM. World Stress Map Database Release 2016. GFZ Data Services. 2016. http://dataservices.gfz-potsdam.de/wsm/showshort.php?id=escidoc:1680899.

  • Jaeger J, Cook N, Zimmerman R. Fundamentals of rock mechanics, vol. 38. 4th ed. New York: Wiley; 2007.

    Google Scholar 

  • Jeanne P, Guglielmi Y, Rutqvist J, Nussbaum C, Birkholzer J. Field characterization of elastic properties across a fault zone reactivated by fluid injection. J Geophys Res Solid Earth. 2017;122(8):6583–98.

    Article  Google Scholar 

  • Klinkenberg L. The permeability of porous media to liquids and gases. In: Drilling production practice. 1941. p. 200–13.

  • Li W, Yost K, Sousa R. Heat transfer between fluid flow and fractured rocks. In: GRC 2013 annual meeting. 2013;37:165–71.

  • López-Hernández A, García-Estrada G, Aguirre-Díaz G, González-Partida E, Palma-Guzmán H, Quijano-León JL. Hydrothermal activity in the Tulancingo–Acoculco Caldera Complex, central Mexico: exploratory studies. Geothermics. 2009;38(3):279–93.

    Article  Google Scholar 

  • Norini G, Groppelli G, Sulpizio R, Carrasco-Núñez G, Dávila-Harris P, Pellicioli C, Zucca F, De Franco R. Structural analysis and thermal remote sensing of the Los Humeros Volcanic Complex: implications for volcano structure and geothermal exploration. J Volcanol Geotherm Res. 2015;301:221–37. https://doi.org/10.1016/j.jvolgeores.2015.05.014.

    Article  Google Scholar 

  • Python Software Foundation. Python. Copyright \(\copyright \)2001–2018. Python Software Foundation (1). https://www.python.org/. 2019.

  • Rutqvist J, Noorishad J, Tsang C-F, Stephansson O. Determination of fracture storativity in hard rocks using high-pressure injection testing. Water Resour Res. 1998;34(10):2551–60.

    Article  Google Scholar 

  • Scheidegger A. The physics of flow through porous media. 3rd ed. Toronto: Univ. of Toronto Press; 1974.

    Google Scholar 

  • Singhal B, Gupta R. Applied hydrogeology of fractured rocks, vol. 2. Dordrecht: Springer; 2010.

    Book  Google Scholar 

  • Snow D. Anisotropic permeability of fractured media. Water Resour Res. 1969;5(6):1273–89. https://doi.org/10.1029/WR005i006p01273.

    Article  Google Scholar 

  • Surrette MJ. Quantifying heterogeneity in variably fractured rock using a hydrostructural domain approach. Simon Fraser University. 2006.

  • Tanikawa W, Shimamoto T. Frictional and transport properties of the Chelungpu fault from shallow borehole data and their correlation with seismic behavior during the 1999 Chi-Chi earthquake. J Geophys Res Solid Earth. 2009;114(1):1–15.

    Google Scholar 

  • Viggiano-Guerra JC, Armenta MF, Silva GRR. Evolucion del sistema geotérmico de acoculco, Pue., México: Un estudio con base en estudios petrográficos del pozo EAC-2 y en otras consideraciones. Geotermia. 2011;24(1):14–24.

    Google Scholar 

  • Zhang J, Standifird WB, Roegiers JC, Zhang Y. Stress-dependent fluid flow and permeability in fractured media: from lab experiments to engineering applications. Rock Mech Rock Eng. 2007;40(1):3–21.

    Article  Google Scholar 

  • Zoback MD. Reservoir geomechanics. Cambridge University Press. 2007. https://www.cambridge.org/core/books/reservoir-geomechanics/71717D4521278CE16D1DED0476ACEDBE.

Download references

Authors' contributions

BL designed the concept of the study, collected the samples in the field, analysed the samples in the rock physics laboratory, implemented the modelling part in COMSOL Multiphysics and wrote the manuscript. AD provided input to the concept, assisted in the implementation in COMSOL Multiphysics and provided input to the manuscript. NDG helped in implementing the COMSOL Multiphysics model and provided input to the manuscript. JK realized all permeability measurements and provided input to the manuscript. POB contributed to the computation of the DFMs and provided input to the manuscript. DB provided input to the manuscript. All authors of this article contributed to review of the manuscript. All authors read and approved the final manuscript.

Acknowledgements

This work is part of the GEMex consortium, a European–Mexican collaboration, which aims at investigating non-conventional geothermal systems: (i) an enhanced geothermal system, as studied in the case of Acoculco site; (ii) and a supercritical system, as studied in the geothermal field of Los Humeros. We would like to thank all the GEMex team, and especially our Mexican partners who guided us in the field. GEMex has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant agreement No. 727550 and the Mexican Energy Sustainability Fund CONACYT-SENER, project 2015-04-268074. More information is available on: http://www.gemex-h2020.eu. We would like to extend this acknowledgment with thanking COMSOL Multiphysics in The Netherlands for providing the technical support; more information is available on: http://www.comsol.com.

Availability of data and materials

The datasets generated during and/or analysed during the current study are available in the repository: https://github.com/BatLep/Sim_THM_Multi_Res-Sc.githttps://github.com/BatLep/SkaPy.git

Competing interests

The authors declare that they have no competing interests.

Funding

GEMex has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant agreement No. 727550 and the Mexican Energy Sustainability Fund CONACYT-SENER, project 2015-04-268074. More information is available on: http://www.gemex-h2020.eu.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Baptiste Lepillier.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lepillier, B., Daniilidis, A., Doonechaly Gholizadeh, N. et al. A fracture flow permeability and stress dependency simulation applied to multi-reservoirs, multi-production scenarios analysis. Geotherm Energy 7, 24 (2019). https://doi.org/10.1186/s40517-019-0141-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40517-019-0141-8

Keywords