- Email: [email protected]
heat transport in fractured formations using discrete fracture and equivalent porous media modeling at the reservoir scale
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Physics and Chemistry of the Earth journal homepage: www.elsevier.com/locate/pce
Comparison of solute/heat transport in fractured formations using discrete fracture and equivalent porous media modeling at the reservoir scale ⁎
Miad Jarrahi , Kayla R. Moore, Hartmut M. Holländer Department of Civil Engineering, University of Manitoba, 15 Gillson Street, Winnipeg, MB, R3T 5V6, Canada
A B S T R A C T
In deep reservoirs, the formations consist of naturally or artificially fractured rocks. The fluid flow coupled with heat transfer through fractures is commonly modeled by discrete fracture approach. In this study, fluid flow and solute/heat transport phenomena were investigated using the equivalent porous media approach and compared to discrete fracture modeling. In the equivalent porous media approach, the fractures were replaced with a porous medium. Its effective porosity and permeability were adjusted to get the simulation results similar to discrete fracture modeling results. In addition, the solid deformation, in equivalent porous media approach, was coupled with heat transfer, while it was de-coupled with fluid flow to increase the simulations’ efficiency. The equivalent porous media approach was applied to two benchmark problems; geothermal doublet, and thermohaline aquifer-aquitard-aquifer model. The comparisons between the equivalent porous media and discrete fracture modeling approaches were in a good agreement and implied the feasibility and efficiency of the equivalent porous media relative to the fully coupled discrete fracture modeling approach.
1. Introduction Predictive modeling in deep reservoirs, such as oil and gas and geothermal reservoirs can be used to better understand flow rates, transport properties, their evolutions, and efficient production. These projects often have high cost and models can be used to better estimate potential issues e.g. (Clauser and Ewert, 2018). However, suitable models must be used to provide accurate results and attention must be given to the simplification of physical processes. In deep reservoirs, the formations are naturally or artificially fractured rocks. Formation water and heat exchange fluids primarily move through the fractures, however, there might be interactions with the surrounding rock matrix. Flow through fractures is represented in models by one of three methods: equivalent porous media, or discrete fracture matrix or dual continuum (Blessent et al., 2014; Dietrich et al., 2005). In the equivalent porous media approach, a zone of increased permeability is used instead of simulating individual fractures and the surrounding damage zone. The representation of fractures as equivalent porous media involves adjusting the volumetric flow rate of the porous media to the same as the fracture by adjusting the permeability e.g. (Jain et al., 2015). This is done on large scale projects, such as reservoir modeling (Jain et al., 2015). Discrete fracture modeling represents the individual fractures features. Dual continuum models combine both matrix flow through porous media and discrete fractures. Flow within a real rough fracture within a deep enhanced geothermal system, considering the detailed geometry of fracture walls and
⁎
aperture size in every location of fracture plate was numerically investigated (Jarrahi and Holländer, 2017). Identifying the features of the fracture (fracture network) in deep formations is complicated and costly. Furthermore, it requires specific tools to model the geometry of the fracture. In regional scale modeling this data is normally unavailable and would be difficult to represent. A number of benchmark models for discrete fracture zones using coupled thermal-hydraulic mechanical (THM), including a reservoir were conducted by Cacace and Jacquey (2017). In their simulations, the fully implicit approach was used to capture the nonlinear behavior of formation matrix. The aim of their modeling was to evaluate the thermo-poroelastic response induced by geothermal operations. The coupling between fluid flow and solid deformation provided robustness of the simulations (Cacace and Jacquey, 2017). In modeling a fractured formation, a discrete fracture approach requires the knowledge of the fractures aperture size distribution, the orientation of the fractures, and the fractures roughness distribution, as well as the interconnection angle between fractures and opening and closing effect of that angle on the fluid flow (Long et al., 1982). Considering the typical depth of 400–800 m below the ground level of these formations, it is hard to access this information. On the other hand, there are other parameters like the flow/heat transfer parameters such as volumetric flow rate, pressure gradient, and temperature gradient that can be measured in the field or in an experimental setup. However, these flow/heat parameters are not sufficient to conduct a discrete fracture modeling as the mentioned parameters of the fracture network
Corresponding author. E-mail address: [email protected] (M. Jarrahi).
https://doi.org/10.1016/j.pce.2019.08.001 Received 31 August 2018; Received in revised form 20 July 2019; Accepted 3 August 2019 1474-7065/ © 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: Miad Jarrahi, Kayla R. Moore and Hartmut M. Holländer, Physics and Chemistry of the Earth, https://doi.org/10.1016/j.pce.2019.08.001
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
2. Mathematical description
is unknown. In this case, the equivalent porous media approach is a proper method that can be implemented only by knowing the flow/heat transfer parameters. The use of equivalent porous media to represent fractures for contaminant transport was compared to dual continuum and discrete fracture approaches and could be calibrated to produce similar results (Jain et al., 2015). The dual continuum method overestimated concentration when groundwater velocity and degradation rates were modified to match the measured values (Daniela et al., 2014). Flow and heat transfer in enhanced geothermal systems were modeled using the equivalent porous media approach e.g. Blessent et al. (2014), and Jain et al. (2015). However, the dual continuum and equivalent porous media approaches required extensive calibration to avoid high errors in predicted mass transport in column experiment. The advantage of using calibration is to improve the accuracy of the simulations without adding any complication to the problem (Blessent et al., 2014). The equivalent porous media approach is more feasible than discrete fracture network for some deep reservoir modeling, as it does not require extensive knowledge of the fracture network. The equivalent porous media approach is a method that can be implemented only by knowing the flow/heat transfer parameters such as volumetric flow rate or pressure or temperature gradient. The main challenges in studying transport phenomena in porous media are based on pores geometry and the complex network of the pore-pore connection and pore throats (Xiong et al., 2016). A proper solution of solute exchange between the solid matrix and fracture and its effect on solute transport within the fracture is developed by (Genuchten, 1985; Harrison et al., 1992; Therrien and Sudicky, 1996). In addition, permeability and porosity of the porous media may be the function of chemical solution and the phase of the solute transport (Appelo et al., 2010). The simulation of different physical and chemical processes such as reactive transport, phase exchange, and thermodynamically consistent oil layers are carried out at reservoir scale using pore network model (Xiong et al., 2016). Accurate simulation of multiphase flow in fractured porous media was performed at the reservoir scale using a discontinuous finite element method with discrete fracture approach (Su et al., 2015). Blessent et al. (2014) used the FRAC3Dvs numerical model (Therrien and Sudicky, 1996) to compare the dualporosity (DP), equivalent porous (EP) medium, and discrete fracture matrix conceptual models to predict field-scale contaminant transport in a fractured clayey till aquitard. Their results showed that the contaminant breakthrough could be calibrated in all the methods. However, they concluded that the EP models and the DP models have limitations for predicting groundwater velocity and mass transport, due to lack of robustness with respect to variable flow and transport conditions in these models (Blessent et al., 2014). Therefore, the accuracy of the equivalent porous media approach in comparison to discrete fracture approach at the reservoir scale is yet to be improved. The objective of this study was to provide a better accuracy of using an equivalent porous media to represent solute and heat transport at the reservoir scale and to reduce the limitations (Blessent et al., 2014) of using EP model in fractured reservoir. Accordingly, the numerical model of a discrete fracture network and an equivalent porous media were compared. The hydraulic parameters of the equivalent porous media required to represent discrete fractures were found by calibrating the distribution of temperature and density of the porous media model to the fracture model. This approach was applied to a 3-D transient fully coupled thermo-hydro-mechanical model to simulate fluid flow and heat transport in a geothermal doublet, considering the deformation of the reservoir. In addition to fluid flow and heat transfer, solute transport was modeled through an equivalent porous matrix, which was replaced with the fracture in a thermohaline aquifer-aquitard-aquifer system. The equivalent porous medium provided simulation results that were in a good agreement with their benchmark results.
Geomechanics involves the complex nature of geomaterials in a heterogeneous system of flow, heat, and stress-strain conditions. The presence of a mobile fluid within a solid matrix alters the solid phase and its governing stress field. The altered stress field within the solid matrix may deform the grains, thus the void spaces. This may cause pore pressure evolutions throughout the matrix domain and may result in significant fluid-solid interactions (FSI). The governing equations of FSI in geomechanics were described in terms of deformation of solid skeleton and the dynamics of fluid flow. The Navier-Stokes equation determines the velocity and pressure field within a fluid domain. However, in porous matrix, Navier-Stokes equation extends to Darcy equation to define the Darcy velocity and pore pressure. The conservation of mass and momentum for fluid flow with the velocity of u (m/s) in fully saturated pores with pore pressure p (Pa) are given (Wang, 2000).
δnρ + ∇ . ρu = 0 δt u=
(1)
k (∇p + ρg∇D) μ
where n (−) is the porosity of the matrix and ρ (kg/m3), g (m/s2 ), are fluid density, and gravity, respectively. k (m2) is the intrinsic permeability, μ (Pa.s) is the dynamic viscosity, and D (m) is the elevation perpendicular to the flow, and t (s) is the time. In this paper, the 3D synthetic geothermal doublet has a plain strain condition. The governing equation for the solid matrix is shown in a plain stain format of the Hook's law (Wang, 2000).
σ σ σ 0 ⎤ v ⎡ xx xy xz ⎤ E ⎡1 − v σ = ⎢ σyx σyy σyz ⎥ = 1−v 0 ⎥ε ⎢ v (1 )(1 2 ) + − v v ⎢ 0 1 − 2v ⎦ ⎣ 0 ⎣ σzx σzy σzz ⎥ ⎦
(2)
− ∇ . σ = ρg − I ∇ . aB p Here, σ (Pa) is the stress tensor, E (Pa), and v (−) are the Young's modulus and Poisson's ratio of the solid matrix, respectively. The term I ∇ . aB p is poroelasticity coupling expression denoted as Biot-Willis pore pressure tensor, where I is the unit tensor. ε (−) is the volumetric strain. Here, the effective stress, σe (Pa) can be defined as σe = σ − aB p , based on Biot consolidation theory (Biot, 1973). The constitutive equation for THM processes for a fluid-solid system including a deformable solid matrix with a volumetric strain rate of ξ (1/s ), and the fluid flow with a velocity of u, in a fully saturated pores with pore pressure p are given. The formation was considered to be homogeneous (Cacace and Jacquey, 2017).
S
δp + ∇ . (−K∇p) = −ρgaBξ δt
(3)
where S (1/m ) is the poroelastic storage coefficient, K (m/s ) is the hydraulic conductivity, and aB (−) is the Biot-Willis coefficient. The energy balance combining with Fourier's law for porous medium forms the heat transfer equation coupled with flow variables (Cacace and Jacquey, 2017).
ρCp
δT − ∇ . (kT∇T) = ∇ . (ρCp uT) δt
(4)
where Cp (J/kgK ) is the heat capacity at constant pressure, T (K) is the temperature, and kT (W/mK ) is thermal conductivity of solid. The evolution of porosity due to the deformation of solid matrix was calculated based on the strain rate and the change in pore pressure (Cacace and Jacquey, 2017).
δn a −n δp = nξ + B δt KB, S δt
(5)
where KB, S (Pa) is the bulk modulus of the solid matrix. The 2
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
conservative form of the solute transport in a porous medium is the governing equation for diffusion, convection, and dispersion of a solute concentration c (mol/m3 ), coupled with flow variables (Cacace and Jacquey, 2017).
δnc + (u−∇ . nDL ) ∇ c= 0 δt
Table 1 Material properties for geothermal doublet.
(6)
where DL (m2/s ) is the fluid's diffusion coefficient. 3. Models discription
Property name
Symbol
Porosity Intrinsic permeability Fluid bulk modulus
n k KB, F
Density
ρ
Thermal conductivity
kT Cp
Heat capacity
For two benchmarks of geological coupled phenomena, synthetic geothermal doublet and thermohaline aquifer-aquitard-aquifer, discrete fracture approach and equivalent porous media approach were performed and their corresponding results were compared to show the accuracy versus simplicity of the approaches. In discrete fracture approach of the synthetic geothermal doublet, the fully coupled thermohydro-mechanical problem is solved. The results were compared to the cited references for validation purposes. Equations (1)–(6), in mathematical description section, describe the coupling parameters in between equations. However, in equivalent porous media approach of the synthetic geothermal doublet, the fluid flow was de-coupled with solid deformation to additionally simplify the model and increase the efficiency of the simulations. The solid deformation in equivalent porous media approach was only considered and coupled with the heat transfer phenomena. The computational fluid dynamics (CFD) techniques (i.e. finite element method) were used to discretize the problem domains. The governing equations were solved using finite element based software COMSOL Multiphysics 5.3a (COMSOL, 2016). The discrete fracture features were modeled as a lower dimension interface within the matrix domain. The fracture object was delimited by assigning a fracture thickness to that lower dimension feature. In COMSOL Multiphysics the fracture corresponding fluid flux can be calculated based on fracture thickness, hydraulic conductivity of the fracture and the pressure head gradient on the fracture's lateral extent.
Value in reservoir
Value in fracture
SI unit
0.1
1
–
1 × 10−15
8.3 × 10−2
1 × 108 2600
2.5 × 109 1000
m2 Pa
3
0.65
Wm−1K−1
950
4200
Wkg−1K−1
kgm−3
considered in the model (i.e. single continuum approach) (see Fig. 1b). For initial conditions, pore pressure and temperature were set to 40 MPa and 150 °C, respectively. At the boundaries, pore pressure and temperature were kept constant and equal to their initial values. The model constrained and no displacement was considered for all the outer faces. The simulation was performed for 50 years of cold water injection. All physical properties can be found in (Cacace and Jacquey, 2017) and are summarized in Table 1. The initial model containing the discrete features is a dual continuum model. For the equivalent porous media simulation, a single permeability was used to represent the entire reservoir. 3.2. Thermohaline aquifer-aquitard-aquifer For the second model, the fluid flow, heat transfer, and solute transfer were coupled and applied to a 2-D thermohaline aquiferaquitard-aquifer (Diersch, 2014). The model domain was 100 × 100 m containing three layers of aquifer-aquitard-aquifer system (Fig. 2). The aquitard, with a thickness, b of 60 m, was centered and surrounded by other two aquifers, each with a thickness of 20 m. The properties of aquifer and aquitard are summarized in Table 2 (Diersch, 2014). A fracture is located at the center of the domain within the aquitard. In the dual continuum setting, the porosity of the fracture is 0.05 with a hydraulic conductivity of 0.1 m/s (Diersch, 2014). While the porosity of the aquifer is 0.2 with a hydraulic conductivity of 10−4 m/s. Longitudinal and transverse dispersivity were 2 m and 0.2 m in the matrix and 0.1 m in the borehole. In this model the effect of a single fracture in a deep stratified reservoir was investigated, numerically. The flow was driven by saltwater and buoyant thermal gradients. The fracture was modeled with a 1-D interface where the Darcy law was solved to obtain the Darcy velocity within the fracture. Each layer was assumed to be homogenous and isotropic. Fig. 2 shows the synthetic geothermal doublet model and its mesh structure and boundary conditions. For the initial conditions, the concentration and hydraulic head H (m) were set to 0. The initial temperature T0 (°C) was a linear function of depth, y (m).
3.1. Synthetic geothermal doublet For the first model, the fluid flow, heat transfer, and poroelastic deformation processes were coupled and applied to a 3-D synthetic geothermal doublet (Cacace and Jacquey, 2017). The reservoir was at depth of 4 km below sea level and extended 500 × 500 × 200 m in x, y, and z directions. The reservoir contained two vertical hydraulic fractures, separated at a distance of 200 m by a porous matrix, at the location of injection and production wells to enhance the efficiency of the geothermal doublet. The hydraulic fractures were modeled with a 2-D interface where the Darcy law was solved to obtain the Darcy velocity within the fracture. The reservoir was homogenous and isotropic. Fig. 1 shows the synthetic geothermal doublet models and their mesh structure and model boundary conditions. The injection and production wells are shown in blue and red lines, respectively. In Fig. 1a, a platelike vertical fracture is considered at the vicinity of each well to address the dual continuum approach. In a second approach, no fracture was
T0 = −0.5 y+ 60
(7)
Dirichlet-type boundary conditions were set to bottom (H = −0.5m ) and top (H = 0m ) of the model. No flow and no heat boundary conditions were assigned to the lateral boundaries (see Fig. 2).
Fig. 1. The synthetic geothermal doublet computational domain; a) dual continuum approach, b) single continuum approach, each containing 158,000 tetrahedral domain elements. 3
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
however, it is sensitive to effective porosity. In a constant flux and inlet temperature (assuming that the flux and inlet temperature are input in the model), according to the Darcy law (Eq. (1)), increasing the permeability would decrease the pressure gradient that maintains the Darcy velocity within the model. Table 3 shows that the increase in permeability, decreases the pressure gradient between the injection and production wells. Thus, the temperature gradient does not change when the Darcy velocity remains constant. While the Darcy velocity is independent of effective porosity, the fluid flow velocity (i.e. the velocity that fluid particles experience within the pores) is inversely related to the effective porosity of the domain. Therefore, changing the effective porosity would change the fluid velocity which in turn alters the fluid density gradient, ∇ρ. According to energy balance (Eq. (4)), the change in density gradient would change the temperature gradient. At a constant injection temperature, the temperature in the production well changes. The effective Porosity of equivalent porous medium was obtained by calibration based on temperature of the production well. Fig. 5 shows that calibration of porosity and the sensitivity of the production temperature. The sensitivity analysis of effective porosity was performed based on the variations of temperature at production well, while the deformation of reservoir was applied (Cacace and Jacquey, 2017). Finally, the sensitivity analysis of porosity of the single continuum approach was tested while the permeability was constant at 7 × 10−15 m2. According to Fig. 5, the equivalent porosity of 0.5 showed the best agreement with the variations of temperature at production well, comparing to Cacace and Jacquey (2017). The difference between the temperature variations in (Cacace and Jacquey, 2017) and our single continuum approach with the equivalent porosity of 0.5 was due to the different assumptions made for deformation of the reservoir, where the solid deformations were decoupled in the our single continuum modeling. Decoupling the solid deformation decreased the simulations run-time significantly. This is due to the flexibility of single continuum approach, where the effective porosity can be adjusted in order to get the same variations in temperature at production well over 50 years (Cacace and Jacquey, 2017). The equivalent porous media simulation used one permeability value and one porosity value for the reservoir. Permeability was adjusted to calibrate the pressure distribution in the dual porosity model and porosity was adjusted to calibrate temperature distribution. As the permeability of the formation decreases, the pressure will increase for the same pumping rate. Therefore, we chose a permeability value so that the resulting volumetric flow rate is the same as for a discrete fracture combined with the surrounding matrix. An increase in porosity increases storage available within the matrix, and is a component in the conservation of mass. Decreasing the porosity decreases the fluid stored in the matrix, and causes the temperature breakthrough of cooler water at the production well to occur faster (Fig. 5). The evolution of effective porosity in equivalent porous medium model due to deformation of solid matrix coupled with heat transfer was calculated and shown in Fig. 6. In current simulation, the deformation of the solid matrix was decoupled with fluid flow to enhance the simulation efficiency and run-time. The changes of effective porosity near the injection well were 0.35% and the production well 0.25%. Whereas in Cacace and Jacquey (2017), the changes of porosity near the injection well were 0.5% and the production well −0.01%. Cacace and Jacquey (2017) modeled a fully coupled termo-hydro-mechanical geothermal doublet considering dual continuum media (discrete fracture and porous matrix). As a result they showed that the maximum changes in porosity and permability occurs around the injection and production wells. The range of changes in porosity and permeability was reported between −0.01–0.5% and −0.03–1.75%, respectively. These ranges of changes in porosity and permability were considered to be negligible. Specially, for the purpose of calibration, the porosity, for example, was changed from 0.1 to 0.5, 0.9, and 0.95, where the physical changes of porosity due to changes in
Fig. 2. Thermohaline aquifer-aquitard-aquifer model domain including 10000 square elements. Table 2 Material properties for thermohaline aquifer-aquitard-aquifer. Property name
Symbol
Value in aquifer
Value in aquitard
SI unit
Thickness Hydraulic conductivity
20
60
m
Porosity Fluid Heat conductivity
b K n kT,f
1 × 10−4 0.2 0.65
1 × 10−8 0.35 0.65
ms−1 –
Solid Heat conductivity
kT
3
3
Wm−1K−1
Wm−1K−1
Temperature and salinity concentration on top boundary was set to 10 °C and 1, while on bottom were set to 60 °C and 0, respectively. 4. Results and discussions In this section, the geothermal doublet and thermohaline aquiferaquitard-aquifer were simulated considering each model's fracture(s) as discrete fractures in a dual continuum model then as an equivalent porous media. The results of equivalent porous medium approach were calibrated with those of discrete fracture modeling. The porous media simulation results were compared to dual porosity models benchmark results available in (Cacace and Jacquey, 2017) and (Diersch, 2014). 4.1. Synthetic geothermal doublet The thermo-poroelastic response of the geothermal doublet was obtained through dual continuum approach for 50 years and compared to reference values (Cacace and Jacquey, 2017). The fractures were then considered as an equivalent porous media and calibrated to match the discrete fracture model. The cold water at 70 °C was injected at a rate of 5 L/s into the reservoir, where a vertical hydraulic fracture was located, and produced at the same rate in the production well, where a similar fracture was located. The contours of pore pressure and temperature simulations are shown in Figs. 3b and 4b and compared to the results of (Cacace and Jacquey, 2017) in Figs. 3a and 4a, respectively. In equivalent single continuum approach, the effect of flow on pressure in the hydraulic fractures was calibrated by adjusting the permeability in the single continuum reservoir to match the discrete fracture modeling approach. The pressure gradient between injection and production wells was evaluated at each permeability that was assigned to the reservoir and summarized in Table 3. The permeability of the fractures in dual continuum modeling approach was 10−15 m2. The permeability of the single continuum porous medium varied to get the same pressure gradient at both modeling approaches. Therefore, the equivalent permeability of the single continuum was obtained to be 7 × 10−15 m2. The calibration of permeability showed that the variations of temperature at production well do not change with the permeability, 4
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
Fig. 3. Contour of pore pressure distribution for a doublet geothermal well after 50 years of stimulation. a) (Cacace and Jacquey, 2017) discrete fracture modeling, b) COMSOL, discrete fracture modeling.
permeability were negligible. It is not considered that these two parameters are independent. In fact, the changes in porosity and permeability and their effect on each other were neglected comparing to the changes that was considered for the sake of calibration. It was assumed that once the permeability is calibrated, the porosity can be adjusted to get the temperature at production well in agreement with the discrete fracture modeling results. The porosity was shown to be decreased near the production well in Cacace and Jacquey (2017) while in our simulation, it was shown to be increased. This is due to the fact that in our equivalent porous media modeling, the permeability was increased that may result a lower pressure at production well comparing to Cacace and Jacquey (2017). As the pore pressure field is different in the equivalent porous media model and Cacace and Jacquey (2017), this will alter the effective porosity change in time according to Eq. (5). Another reason for the difference between the results of the current modeling and Cacace and Jacquey (2017) is that in our modeling the solid deformation was coupled with heat transfer, however, was decoupled with fluid flow, while in Cacace and Jacquey (2017) the model was fully coupled. The decoupling was to simplify the modeling and increase the efficiency of the modeling.
Table 3 Calibration of permeability of equivalent porous medium based on pressure gradient between injection and production wells. Modeling type
Permeability (m2 )
Pressure gradient (kPa/ m )
Absolute differtence (kPa/ m )
Dual continuum matrix Fracture
10−15
73.0
0
Single continuum matrix
10−15
507.5
434.5
2 × 10−15
254.0
181.0
4 × 10−15
127.0
54.0
6 × 10−15
84.5
11.5
7 × 10−15
72.5
0.5
8 × 10−15
63.5
9.5
8.33 × 10−10
saltwater begins entering on the top aquifer and moves through a borehole to the bottom aquifer. The fracture interconnects the upper and lower aquifers, making a shortcut for governing density driven flow. The density driven flow and solute transport processes were affected by thermal buoyancy. First, the salinity and heat transport were modeled through discrete fracture modeling to match benchmark conditions. The salinity and heat transport were then modeled as an equivalent porous media and compared to the discrete fracture dual continuum results to validate the simulation. The contours of concentration and temperature are shown in Fig. 7.
4.2. Thermohaline aquifer-aquitard-aquifer A domain aquifer-aquitard-aquifer containing freshwater was exposed to an increasing thermal gradient. In the simulation, a heavy
Fig. 4. Contour of temperature after 50 years of stimulation; a) discrete fracture modeling (Cacace and Jacquey, 2017), b) discrete fracture modeling COMSOL. 5
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
itself, with 2-D and 3-D simulations varying by 2.5 days and discretization of the model varying by 2 days (Diersch, 2014). However, we were able to achieve a good calibration when comparing discrete fracture and equivalent porous media under the same model conditions (Fig. 8). Especially, in early time simulations, during the initial breakthrough, the curve fit between the fractured model and equivalent porous media approach was good. In late time, the increase in salinity concentration at the exit point at lower aquifer was faster in the equivalent porous media. This is due to the higher fluid velocity in the equivalent porous zone comparing to the fracture. As the porous zone with the effective porosity of n = 0.05 is replaced with the fracture, the fluid velocity was increased by the factor of 1/ n which resulted the faster breakthrough in comparison to the fracture model. Fig. 8b shows the log-log plot of concentration-time curve in late time simulation. It revealed the agreement between the equivalent porous medium modeling at K = 3 × 10−5m / s and discrete fracture modeling. As a result the solute transport in aquifer-aquitard-aquifer may alter by changing the hydraulic conductivity of the considered region in the middle of the domain. While in geothermal doublet model, the temperature gradient was only sensitive to the effective porosity of the replaced porous region (see Fig. 5). This showed that the sensitivity of a phenomenon in a fractured porous medium model depends on the behavior of the flow and solute/heat transport within the fractures and the fractures geometrical features (Long et al., 1982). Moreover, Fig. 8 showed that solute transport within a porous region with hydraulic conductivity of 3 × 10−5 m/s was equal to the solute transport within the fracture at early time simulation. However, due to the higher fluid velocity (i.e. Darcy velocity/porosity) within the porous region in comparison to the fluid velocity within the fracture, the breakthrough phenomenon was occurred faster. Therefore, at late time simulations, the salinity curve in equivalent porous media modeling (red line in Fig. 8) is above the salinity curve in discrete fracture modeling (solid black line in Fig. 8) with the error of 2%. Both results were converged after 50 days. Thus, it can be concluded that the porous region with the effective porosity of 0.05 and hydraulic conductivity of 3 × 10−5 m/s is equivalent to the existing fracture in the thermohaline aquifer-aquitardaquifer system with the maximum error of 2%. The representation of fracture flow in COMSOL was simulated as a porous media was calibrated to match discrete flow by adjusting only the permeability of the fractured area. By limiting the number of parameters that require adjustment we simplify the transition between discrete fractures and equivalent porous media. Variations in breakthrough curves are expected between equivalent porous media and discrete fracture models (Diersch, 2014), however total mass flux can be calibrated to be equal for both model types. Other simulations e.g. Daniela et al. (2014) have adjusted dispersivity to unrealistically high value to treat the velocity difference between porous media and fractures as an apparent dispersion process. The method of adjusting permeability and porosity provides reasonable calibration results for discrete fractures represented as equivalent porous media for a thermohaline reservoir. This method would be effective in reservoirs where limited calibration data are available.
Fig. 5. Production temperature-time plot: calibration of porosity of equivalent porous medium based on production temperature.
Fig. 6. The distribution of changes in porosity in horizontal slice of the model after 50 years.
To represent the fracture as a porous media, a porous zone with higher permeability and the effective porosity of 0.05, was defined at the center of the aquifer-aquitard-aquifer system with the thickness of 0.5 m to enhance the volumetric flow rate of the fracture. Similar to the synthetic geothermal doublet model, the sensitivity analysis of effective porosity was performed based on the variations in volumetric flow rate. Moreover, the effective porosity of 0.05 was obtained in a calibration process that provided the desired volumetric flow rate. The flow and transport through the fracture represented as a porous zone was calibrated using hydraulic conductivity values based on salinity concentrations at the exit point, where the wells meets the top of the bottom aquifer. Porosity, 0.05, and dispersivity, 0.1 m, were consistent with the corresponding values in discrete fracture modeling. The salinity concentration plots over 50 days for porous media with discrete fracture modeling and with results of Diersch (2014) are depicted in Fig. 8a (magnified in 8b), illustrating the calibration of permeability in the porous zone. Fig. 9 shows the contours of salinity concentration and temperature after 20 d for discrete fracture modeling and the equivalent porous media approach. The equivalent hydraulic conductivity was 3 × 10−5 m/s, based on comparison to the discrete fracture modeling results found in this study. The jump in the curves shows the breakthrough of the salinity through the aquitard to the lower aquifer. The initial breakthrough in the discrete fracture in this study was faster than (Diersch, 2014), approximately 1 day compared to 10 days (Figs. 7 and 8). A breakthrough time of approximately 5 days was observed for 5 × 10−6 m/s (Fig. 8), which was comparable to the 3-D case for discrete fractures in Diersch (2014), where the fracture hydraulic conductivity was 0.1 m/s. The breakthrough time was sensitive to the geometry and flow within the fracture
5. Conclusion In this paper, we presented simulations comparing equivalent porous media modeling approach to fully coupled discrete fracture (or dual continuum approach) for geothermal and thermohaline fractured reservoir applications. The efficiency and simplicity of modeling at the reservoir scale was enhanced by removing the fractures and substituting a matrix with equivalent hydraulic properties (i.e. porosity and permeability). Geological heterogeneities were not taken into account, however we were able to produce similar results to other modeling approaches that did take fractures into account e.g. Cacace and Jacquey (2017). Although the reliability of predications for reservoirs can be increased by considering accurate physical phenomena into account 6
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
Fig. 7. Contours of I) concentration of salinity, a) Diersch (2014), b) COMSOL; and II) temperature, c) Diersch (2014), d) COMSOL, at different times t = 5, 20 and 100 days in upper, middle, and lower row, respectively.
Fig. 8. Salinity concentration at exit point (50 m, 20 m); a) Concentration-time curve over 50 years, b) magnified of 8a between salinity concentration of 0.8–1.
7
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
Fig. 9. Contour of salinity concentration after 20 d, a) discrete fracture modeling, b) equivalent porous media approach.
References
(Cacace and Jacquey, 2017), simple equivalent porous media approaches can be used to produce accurate results, when adequate calibration data are available. According to the results, equivalent porous media models can be calibrated to reproduce THM and thermohaline results similar to dual porosity models by varying permeability and porosity. Equivalent porous media parameters were estimated using reservoir pressure and temperature data. The validation of this approach was confirmed by making the comparison between models of two benchmark problems in Cacace and Jacquey (2017) and Diersch (2014). Modeling the discrete fracture approach requires the detailed knowledge of fracture network in a fractured formation. However, by knowing the flow/heat transfer parameters such as volumetric flow rate and pressure/temperature gradient, that can be obtained from the field measurements, the equivalent porous media approach may be used to get numerical results that are in good agreement with the available field data. Nevertheless, the application of equivalent porous media approach instead of the discrete fracture modeling is not recommended for any fractured formation, as the nature of the fractures may not allow to find any equivalent porous region for the modeling. More details of the necessary conditions of the formations that may be modeled by equivalent porous media can be found in (Long et al., 1982).
Appelo, C.A.J., Van Loon, L.R., Wersin, P., 2010. Multicomponent diffusion of a suite of tracers (HTO, Cl, Br, I, Na, Sr, Cs) in a single sample of Opalinus Clay. Geochem. Cosmochim. Acta 74 (4), 1201–1219. Biot, 1973. Nonlinear and semilinear rheology of porous solids. J. Geophys. Res. 78 (23), 4924–4937. Blessent, D., Jorgensen, P.R., Therrien, R.e., 2014. Comparing Discrete Fracture and ContinuumModels to Predict Contaminant Transportin Fractured Porous Media. Groundwater, vol. 52. Cacace, M., Jacquey, A.B., 2017. Flexible parallel implicit modelling of coupled ThermalHydraulic-Mechanical processes in fractured rocks. Solid Earth 8, 921–941. Clauser, C., Ewert, M., 2018. The renewables cost challenge: levelized cost of geothermal electric energy compared to other sources of primary energy – review and case study. Renew. Sustain. Energy Rev. 82, 3683–3693. COMSOL, 2016. Comsol Multiphysics, User's Guide. COMSOL AB, Stockholm, Sweden. Daniela, B., Jorgensen, P.R., René, T., 2014. Comparing discrete fracture and continuum models to predict contaminant transport in fractured porous media. Gr. Water 52 (1), 84–95. Diersch, H.-J., 2014. Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media. Springer-Verlag Berlin Heidelberg. Dietrich, P., Helmig, R., Sauter, M., Hötzl, H., Köngeter, J., Teutsch, G., 2005. Flow and Transport in Fractured Porous Media. Springer-Verlag Berlin Heidelberg. Genuchten, v., 1985. A general approach for modelingsolute transport in structured soils. IAH Memories 17 (2). Harrison, B., Sudicky, E.A., Cherry, J.A., 1992. Numerical analysis of solute migration through fractured clayey deposits into underlying aquifers. Water Resour. Res. 28 (2), 515–526. Jain, C., Vogt, C., Clauser, C., 2015. Maximum potential for geothermal power in Germany based on engineered geothermal systems. Geotherm. Energy 3 (1), 15. Jarrahi, M., Holländer, H., 2017. Numerical simulation of flow heterogeneities within a real rough fracture and its transmissivity. Energy Procedia 125, 353–362. Long, J.C.S., Remer, J.S., Wilson, C.R., Wotherspoon, P.A., 1982. Porous media equivalents for networks of discontinuous fractures. Water Resour. Res. 18. Su, K., Latham, J.-P., Pavlidis, D., Xiang, J., Fang, F., Mostaghimi, P., Percival, J.R., Pain, C.C., Jackson, M.D., 2015. Multiphase flow simulation through porous media with explicitly resolved fractures. Geofluids 15 (4), 592–607. Therrien, R., Sudicky, E.A., 1996. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media. J. Contam. Hydrol. 23 (1), 1–44. Wang, H.F., 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton University Press, pp. 304. Xiong, Q., Baychev, T.G., Jivkov, A.P., 2016. Review of pore network modelling of porous media: experimental characterisations, network constructions and applications to reactive transport. J. Contam. Hydrol. 192, 101–117.
Acknowledgment and Permissions The authors would like to appreciate Dr.-Ing. Antoine Jacquey, Scientist at GFZ German research center for Geosciences, for providing required data and the permission to use the figures of (Cacace and Jacquey, 2017). The agreement with Springer Nature to reuse figures of (Diersch, 2014) is issued under license number 4643910403717. Additionally, the authors would like to thank the reviewers of this manuscript for their constructive comments and suggestions. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.pce.2019.08.001.
8
Contents lists available at ScienceDirect
Physics and Chemistry of the Earth journal homepage: www.elsevier.com/locate/pce
Comparison of solute/heat transport in fractured formations using discrete fracture and equivalent porous media modeling at the reservoir scale ⁎
Miad Jarrahi , Kayla R. Moore, Hartmut M. Holländer Department of Civil Engineering, University of Manitoba, 15 Gillson Street, Winnipeg, MB, R3T 5V6, Canada
A B S T R A C T
In deep reservoirs, the formations consist of naturally or artificially fractured rocks. The fluid flow coupled with heat transfer through fractures is commonly modeled by discrete fracture approach. In this study, fluid flow and solute/heat transport phenomena were investigated using the equivalent porous media approach and compared to discrete fracture modeling. In the equivalent porous media approach, the fractures were replaced with a porous medium. Its effective porosity and permeability were adjusted to get the simulation results similar to discrete fracture modeling results. In addition, the solid deformation, in equivalent porous media approach, was coupled with heat transfer, while it was de-coupled with fluid flow to increase the simulations’ efficiency. The equivalent porous media approach was applied to two benchmark problems; geothermal doublet, and thermohaline aquifer-aquitard-aquifer model. The comparisons between the equivalent porous media and discrete fracture modeling approaches were in a good agreement and implied the feasibility and efficiency of the equivalent porous media relative to the fully coupled discrete fracture modeling approach.
1. Introduction Predictive modeling in deep reservoirs, such as oil and gas and geothermal reservoirs can be used to better understand flow rates, transport properties, their evolutions, and efficient production. These projects often have high cost and models can be used to better estimate potential issues e.g. (Clauser and Ewert, 2018). However, suitable models must be used to provide accurate results and attention must be given to the simplification of physical processes. In deep reservoirs, the formations are naturally or artificially fractured rocks. Formation water and heat exchange fluids primarily move through the fractures, however, there might be interactions with the surrounding rock matrix. Flow through fractures is represented in models by one of three methods: equivalent porous media, or discrete fracture matrix or dual continuum (Blessent et al., 2014; Dietrich et al., 2005). In the equivalent porous media approach, a zone of increased permeability is used instead of simulating individual fractures and the surrounding damage zone. The representation of fractures as equivalent porous media involves adjusting the volumetric flow rate of the porous media to the same as the fracture by adjusting the permeability e.g. (Jain et al., 2015). This is done on large scale projects, such as reservoir modeling (Jain et al., 2015). Discrete fracture modeling represents the individual fractures features. Dual continuum models combine both matrix flow through porous media and discrete fractures. Flow within a real rough fracture within a deep enhanced geothermal system, considering the detailed geometry of fracture walls and
⁎
aperture size in every location of fracture plate was numerically investigated (Jarrahi and Holländer, 2017). Identifying the features of the fracture (fracture network) in deep formations is complicated and costly. Furthermore, it requires specific tools to model the geometry of the fracture. In regional scale modeling this data is normally unavailable and would be difficult to represent. A number of benchmark models for discrete fracture zones using coupled thermal-hydraulic mechanical (THM), including a reservoir were conducted by Cacace and Jacquey (2017). In their simulations, the fully implicit approach was used to capture the nonlinear behavior of formation matrix. The aim of their modeling was to evaluate the thermo-poroelastic response induced by geothermal operations. The coupling between fluid flow and solid deformation provided robustness of the simulations (Cacace and Jacquey, 2017). In modeling a fractured formation, a discrete fracture approach requires the knowledge of the fractures aperture size distribution, the orientation of the fractures, and the fractures roughness distribution, as well as the interconnection angle between fractures and opening and closing effect of that angle on the fluid flow (Long et al., 1982). Considering the typical depth of 400–800 m below the ground level of these formations, it is hard to access this information. On the other hand, there are other parameters like the flow/heat transfer parameters such as volumetric flow rate, pressure gradient, and temperature gradient that can be measured in the field or in an experimental setup. However, these flow/heat parameters are not sufficient to conduct a discrete fracture modeling as the mentioned parameters of the fracture network
Corresponding author. E-mail address: [email protected] (M. Jarrahi).
https://doi.org/10.1016/j.pce.2019.08.001 Received 31 August 2018; Received in revised form 20 July 2019; Accepted 3 August 2019 1474-7065/ © 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: Miad Jarrahi, Kayla R. Moore and Hartmut M. Holländer, Physics and Chemistry of the Earth, https://doi.org/10.1016/j.pce.2019.08.001
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
2. Mathematical description
is unknown. In this case, the equivalent porous media approach is a proper method that can be implemented only by knowing the flow/heat transfer parameters. The use of equivalent porous media to represent fractures for contaminant transport was compared to dual continuum and discrete fracture approaches and could be calibrated to produce similar results (Jain et al., 2015). The dual continuum method overestimated concentration when groundwater velocity and degradation rates were modified to match the measured values (Daniela et al., 2014). Flow and heat transfer in enhanced geothermal systems were modeled using the equivalent porous media approach e.g. Blessent et al. (2014), and Jain et al. (2015). However, the dual continuum and equivalent porous media approaches required extensive calibration to avoid high errors in predicted mass transport in column experiment. The advantage of using calibration is to improve the accuracy of the simulations without adding any complication to the problem (Blessent et al., 2014). The equivalent porous media approach is more feasible than discrete fracture network for some deep reservoir modeling, as it does not require extensive knowledge of the fracture network. The equivalent porous media approach is a method that can be implemented only by knowing the flow/heat transfer parameters such as volumetric flow rate or pressure or temperature gradient. The main challenges in studying transport phenomena in porous media are based on pores geometry and the complex network of the pore-pore connection and pore throats (Xiong et al., 2016). A proper solution of solute exchange between the solid matrix and fracture and its effect on solute transport within the fracture is developed by (Genuchten, 1985; Harrison et al., 1992; Therrien and Sudicky, 1996). In addition, permeability and porosity of the porous media may be the function of chemical solution and the phase of the solute transport (Appelo et al., 2010). The simulation of different physical and chemical processes such as reactive transport, phase exchange, and thermodynamically consistent oil layers are carried out at reservoir scale using pore network model (Xiong et al., 2016). Accurate simulation of multiphase flow in fractured porous media was performed at the reservoir scale using a discontinuous finite element method with discrete fracture approach (Su et al., 2015). Blessent et al. (2014) used the FRAC3Dvs numerical model (Therrien and Sudicky, 1996) to compare the dualporosity (DP), equivalent porous (EP) medium, and discrete fracture matrix conceptual models to predict field-scale contaminant transport in a fractured clayey till aquitard. Their results showed that the contaminant breakthrough could be calibrated in all the methods. However, they concluded that the EP models and the DP models have limitations for predicting groundwater velocity and mass transport, due to lack of robustness with respect to variable flow and transport conditions in these models (Blessent et al., 2014). Therefore, the accuracy of the equivalent porous media approach in comparison to discrete fracture approach at the reservoir scale is yet to be improved. The objective of this study was to provide a better accuracy of using an equivalent porous media to represent solute and heat transport at the reservoir scale and to reduce the limitations (Blessent et al., 2014) of using EP model in fractured reservoir. Accordingly, the numerical model of a discrete fracture network and an equivalent porous media were compared. The hydraulic parameters of the equivalent porous media required to represent discrete fractures were found by calibrating the distribution of temperature and density of the porous media model to the fracture model. This approach was applied to a 3-D transient fully coupled thermo-hydro-mechanical model to simulate fluid flow and heat transport in a geothermal doublet, considering the deformation of the reservoir. In addition to fluid flow and heat transfer, solute transport was modeled through an equivalent porous matrix, which was replaced with the fracture in a thermohaline aquifer-aquitard-aquifer system. The equivalent porous medium provided simulation results that were in a good agreement with their benchmark results.
Geomechanics involves the complex nature of geomaterials in a heterogeneous system of flow, heat, and stress-strain conditions. The presence of a mobile fluid within a solid matrix alters the solid phase and its governing stress field. The altered stress field within the solid matrix may deform the grains, thus the void spaces. This may cause pore pressure evolutions throughout the matrix domain and may result in significant fluid-solid interactions (FSI). The governing equations of FSI in geomechanics were described in terms of deformation of solid skeleton and the dynamics of fluid flow. The Navier-Stokes equation determines the velocity and pressure field within a fluid domain. However, in porous matrix, Navier-Stokes equation extends to Darcy equation to define the Darcy velocity and pore pressure. The conservation of mass and momentum for fluid flow with the velocity of u (m/s) in fully saturated pores with pore pressure p (Pa) are given (Wang, 2000).
δnρ + ∇ . ρu = 0 δt u=
(1)
k (∇p + ρg∇D) μ
where n (−) is the porosity of the matrix and ρ (kg/m3), g (m/s2 ), are fluid density, and gravity, respectively. k (m2) is the intrinsic permeability, μ (Pa.s) is the dynamic viscosity, and D (m) is the elevation perpendicular to the flow, and t (s) is the time. In this paper, the 3D synthetic geothermal doublet has a plain strain condition. The governing equation for the solid matrix is shown in a plain stain format of the Hook's law (Wang, 2000).
σ σ σ 0 ⎤ v ⎡ xx xy xz ⎤ E ⎡1 − v σ = ⎢ σyx σyy σyz ⎥ = 1−v 0 ⎥ε ⎢ v (1 )(1 2 ) + − v v ⎢ 0 1 − 2v ⎦ ⎣ 0 ⎣ σzx σzy σzz ⎥ ⎦
(2)
− ∇ . σ = ρg − I ∇ . aB p Here, σ (Pa) is the stress tensor, E (Pa), and v (−) are the Young's modulus and Poisson's ratio of the solid matrix, respectively. The term I ∇ . aB p is poroelasticity coupling expression denoted as Biot-Willis pore pressure tensor, where I is the unit tensor. ε (−) is the volumetric strain. Here, the effective stress, σe (Pa) can be defined as σe = σ − aB p , based on Biot consolidation theory (Biot, 1973). The constitutive equation for THM processes for a fluid-solid system including a deformable solid matrix with a volumetric strain rate of ξ (1/s ), and the fluid flow with a velocity of u, in a fully saturated pores with pore pressure p are given. The formation was considered to be homogeneous (Cacace and Jacquey, 2017).
S
δp + ∇ . (−K∇p) = −ρgaBξ δt
(3)
where S (1/m ) is the poroelastic storage coefficient, K (m/s ) is the hydraulic conductivity, and aB (−) is the Biot-Willis coefficient. The energy balance combining with Fourier's law for porous medium forms the heat transfer equation coupled with flow variables (Cacace and Jacquey, 2017).
ρCp
δT − ∇ . (kT∇T) = ∇ . (ρCp uT) δt
(4)
where Cp (J/kgK ) is the heat capacity at constant pressure, T (K) is the temperature, and kT (W/mK ) is thermal conductivity of solid. The evolution of porosity due to the deformation of solid matrix was calculated based on the strain rate and the change in pore pressure (Cacace and Jacquey, 2017).
δn a −n δp = nξ + B δt KB, S δt
(5)
where KB, S (Pa) is the bulk modulus of the solid matrix. The 2
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
conservative form of the solute transport in a porous medium is the governing equation for diffusion, convection, and dispersion of a solute concentration c (mol/m3 ), coupled with flow variables (Cacace and Jacquey, 2017).
δnc + (u−∇ . nDL ) ∇ c= 0 δt
Table 1 Material properties for geothermal doublet.
(6)
where DL (m2/s ) is the fluid's diffusion coefficient. 3. Models discription
Property name
Symbol
Porosity Intrinsic permeability Fluid bulk modulus
n k KB, F
Density
ρ
Thermal conductivity
kT Cp
Heat capacity
For two benchmarks of geological coupled phenomena, synthetic geothermal doublet and thermohaline aquifer-aquitard-aquifer, discrete fracture approach and equivalent porous media approach were performed and their corresponding results were compared to show the accuracy versus simplicity of the approaches. In discrete fracture approach of the synthetic geothermal doublet, the fully coupled thermohydro-mechanical problem is solved. The results were compared to the cited references for validation purposes. Equations (1)–(6), in mathematical description section, describe the coupling parameters in between equations. However, in equivalent porous media approach of the synthetic geothermal doublet, the fluid flow was de-coupled with solid deformation to additionally simplify the model and increase the efficiency of the simulations. The solid deformation in equivalent porous media approach was only considered and coupled with the heat transfer phenomena. The computational fluid dynamics (CFD) techniques (i.e. finite element method) were used to discretize the problem domains. The governing equations were solved using finite element based software COMSOL Multiphysics 5.3a (COMSOL, 2016). The discrete fracture features were modeled as a lower dimension interface within the matrix domain. The fracture object was delimited by assigning a fracture thickness to that lower dimension feature. In COMSOL Multiphysics the fracture corresponding fluid flux can be calculated based on fracture thickness, hydraulic conductivity of the fracture and the pressure head gradient on the fracture's lateral extent.
Value in reservoir
Value in fracture
SI unit
0.1
1
–
1 × 10−15
8.3 × 10−2
1 × 108 2600
2.5 × 109 1000
m2 Pa
3
0.65
Wm−1K−1
950
4200
Wkg−1K−1
kgm−3
considered in the model (i.e. single continuum approach) (see Fig. 1b). For initial conditions, pore pressure and temperature were set to 40 MPa and 150 °C, respectively. At the boundaries, pore pressure and temperature were kept constant and equal to their initial values. The model constrained and no displacement was considered for all the outer faces. The simulation was performed for 50 years of cold water injection. All physical properties can be found in (Cacace and Jacquey, 2017) and are summarized in Table 1. The initial model containing the discrete features is a dual continuum model. For the equivalent porous media simulation, a single permeability was used to represent the entire reservoir. 3.2. Thermohaline aquifer-aquitard-aquifer For the second model, the fluid flow, heat transfer, and solute transfer were coupled and applied to a 2-D thermohaline aquiferaquitard-aquifer (Diersch, 2014). The model domain was 100 × 100 m containing three layers of aquifer-aquitard-aquifer system (Fig. 2). The aquitard, with a thickness, b of 60 m, was centered and surrounded by other two aquifers, each with a thickness of 20 m. The properties of aquifer and aquitard are summarized in Table 2 (Diersch, 2014). A fracture is located at the center of the domain within the aquitard. In the dual continuum setting, the porosity of the fracture is 0.05 with a hydraulic conductivity of 0.1 m/s (Diersch, 2014). While the porosity of the aquifer is 0.2 with a hydraulic conductivity of 10−4 m/s. Longitudinal and transverse dispersivity were 2 m and 0.2 m in the matrix and 0.1 m in the borehole. In this model the effect of a single fracture in a deep stratified reservoir was investigated, numerically. The flow was driven by saltwater and buoyant thermal gradients. The fracture was modeled with a 1-D interface where the Darcy law was solved to obtain the Darcy velocity within the fracture. Each layer was assumed to be homogenous and isotropic. Fig. 2 shows the synthetic geothermal doublet model and its mesh structure and boundary conditions. For the initial conditions, the concentration and hydraulic head H (m) were set to 0. The initial temperature T0 (°C) was a linear function of depth, y (m).
3.1. Synthetic geothermal doublet For the first model, the fluid flow, heat transfer, and poroelastic deformation processes were coupled and applied to a 3-D synthetic geothermal doublet (Cacace and Jacquey, 2017). The reservoir was at depth of 4 km below sea level and extended 500 × 500 × 200 m in x, y, and z directions. The reservoir contained two vertical hydraulic fractures, separated at a distance of 200 m by a porous matrix, at the location of injection and production wells to enhance the efficiency of the geothermal doublet. The hydraulic fractures were modeled with a 2-D interface where the Darcy law was solved to obtain the Darcy velocity within the fracture. The reservoir was homogenous and isotropic. Fig. 1 shows the synthetic geothermal doublet models and their mesh structure and model boundary conditions. The injection and production wells are shown in blue and red lines, respectively. In Fig. 1a, a platelike vertical fracture is considered at the vicinity of each well to address the dual continuum approach. In a second approach, no fracture was
T0 = −0.5 y+ 60
(7)
Dirichlet-type boundary conditions were set to bottom (H = −0.5m ) and top (H = 0m ) of the model. No flow and no heat boundary conditions were assigned to the lateral boundaries (see Fig. 2).
Fig. 1. The synthetic geothermal doublet computational domain; a) dual continuum approach, b) single continuum approach, each containing 158,000 tetrahedral domain elements. 3
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
however, it is sensitive to effective porosity. In a constant flux and inlet temperature (assuming that the flux and inlet temperature are input in the model), according to the Darcy law (Eq. (1)), increasing the permeability would decrease the pressure gradient that maintains the Darcy velocity within the model. Table 3 shows that the increase in permeability, decreases the pressure gradient between the injection and production wells. Thus, the temperature gradient does not change when the Darcy velocity remains constant. While the Darcy velocity is independent of effective porosity, the fluid flow velocity (i.e. the velocity that fluid particles experience within the pores) is inversely related to the effective porosity of the domain. Therefore, changing the effective porosity would change the fluid velocity which in turn alters the fluid density gradient, ∇ρ. According to energy balance (Eq. (4)), the change in density gradient would change the temperature gradient. At a constant injection temperature, the temperature in the production well changes. The effective Porosity of equivalent porous medium was obtained by calibration based on temperature of the production well. Fig. 5 shows that calibration of porosity and the sensitivity of the production temperature. The sensitivity analysis of effective porosity was performed based on the variations of temperature at production well, while the deformation of reservoir was applied (Cacace and Jacquey, 2017). Finally, the sensitivity analysis of porosity of the single continuum approach was tested while the permeability was constant at 7 × 10−15 m2. According to Fig. 5, the equivalent porosity of 0.5 showed the best agreement with the variations of temperature at production well, comparing to Cacace and Jacquey (2017). The difference between the temperature variations in (Cacace and Jacquey, 2017) and our single continuum approach with the equivalent porosity of 0.5 was due to the different assumptions made for deformation of the reservoir, where the solid deformations were decoupled in the our single continuum modeling. Decoupling the solid deformation decreased the simulations run-time significantly. This is due to the flexibility of single continuum approach, where the effective porosity can be adjusted in order to get the same variations in temperature at production well over 50 years (Cacace and Jacquey, 2017). The equivalent porous media simulation used one permeability value and one porosity value for the reservoir. Permeability was adjusted to calibrate the pressure distribution in the dual porosity model and porosity was adjusted to calibrate temperature distribution. As the permeability of the formation decreases, the pressure will increase for the same pumping rate. Therefore, we chose a permeability value so that the resulting volumetric flow rate is the same as for a discrete fracture combined with the surrounding matrix. An increase in porosity increases storage available within the matrix, and is a component in the conservation of mass. Decreasing the porosity decreases the fluid stored in the matrix, and causes the temperature breakthrough of cooler water at the production well to occur faster (Fig. 5). The evolution of effective porosity in equivalent porous medium model due to deformation of solid matrix coupled with heat transfer was calculated and shown in Fig. 6. In current simulation, the deformation of the solid matrix was decoupled with fluid flow to enhance the simulation efficiency and run-time. The changes of effective porosity near the injection well were 0.35% and the production well 0.25%. Whereas in Cacace and Jacquey (2017), the changes of porosity near the injection well were 0.5% and the production well −0.01%. Cacace and Jacquey (2017) modeled a fully coupled termo-hydro-mechanical geothermal doublet considering dual continuum media (discrete fracture and porous matrix). As a result they showed that the maximum changes in porosity and permability occurs around the injection and production wells. The range of changes in porosity and permeability was reported between −0.01–0.5% and −0.03–1.75%, respectively. These ranges of changes in porosity and permability were considered to be negligible. Specially, for the purpose of calibration, the porosity, for example, was changed from 0.1 to 0.5, 0.9, and 0.95, where the physical changes of porosity due to changes in
Fig. 2. Thermohaline aquifer-aquitard-aquifer model domain including 10000 square elements. Table 2 Material properties for thermohaline aquifer-aquitard-aquifer. Property name
Symbol
Value in aquifer
Value in aquitard
SI unit
Thickness Hydraulic conductivity
20
60
m
Porosity Fluid Heat conductivity
b K n kT,f
1 × 10−4 0.2 0.65
1 × 10−8 0.35 0.65
ms−1 –
Solid Heat conductivity
kT
3
3
Wm−1K−1
Wm−1K−1
Temperature and salinity concentration on top boundary was set to 10 °C and 1, while on bottom were set to 60 °C and 0, respectively. 4. Results and discussions In this section, the geothermal doublet and thermohaline aquiferaquitard-aquifer were simulated considering each model's fracture(s) as discrete fractures in a dual continuum model then as an equivalent porous media. The results of equivalent porous medium approach were calibrated with those of discrete fracture modeling. The porous media simulation results were compared to dual porosity models benchmark results available in (Cacace and Jacquey, 2017) and (Diersch, 2014). 4.1. Synthetic geothermal doublet The thermo-poroelastic response of the geothermal doublet was obtained through dual continuum approach for 50 years and compared to reference values (Cacace and Jacquey, 2017). The fractures were then considered as an equivalent porous media and calibrated to match the discrete fracture model. The cold water at 70 °C was injected at a rate of 5 L/s into the reservoir, where a vertical hydraulic fracture was located, and produced at the same rate in the production well, where a similar fracture was located. The contours of pore pressure and temperature simulations are shown in Figs. 3b and 4b and compared to the results of (Cacace and Jacquey, 2017) in Figs. 3a and 4a, respectively. In equivalent single continuum approach, the effect of flow on pressure in the hydraulic fractures was calibrated by adjusting the permeability in the single continuum reservoir to match the discrete fracture modeling approach. The pressure gradient between injection and production wells was evaluated at each permeability that was assigned to the reservoir and summarized in Table 3. The permeability of the fractures in dual continuum modeling approach was 10−15 m2. The permeability of the single continuum porous medium varied to get the same pressure gradient at both modeling approaches. Therefore, the equivalent permeability of the single continuum was obtained to be 7 × 10−15 m2. The calibration of permeability showed that the variations of temperature at production well do not change with the permeability, 4
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
Fig. 3. Contour of pore pressure distribution for a doublet geothermal well after 50 years of stimulation. a) (Cacace and Jacquey, 2017) discrete fracture modeling, b) COMSOL, discrete fracture modeling.
permeability were negligible. It is not considered that these two parameters are independent. In fact, the changes in porosity and permeability and their effect on each other were neglected comparing to the changes that was considered for the sake of calibration. It was assumed that once the permeability is calibrated, the porosity can be adjusted to get the temperature at production well in agreement with the discrete fracture modeling results. The porosity was shown to be decreased near the production well in Cacace and Jacquey (2017) while in our simulation, it was shown to be increased. This is due to the fact that in our equivalent porous media modeling, the permeability was increased that may result a lower pressure at production well comparing to Cacace and Jacquey (2017). As the pore pressure field is different in the equivalent porous media model and Cacace and Jacquey (2017), this will alter the effective porosity change in time according to Eq. (5). Another reason for the difference between the results of the current modeling and Cacace and Jacquey (2017) is that in our modeling the solid deformation was coupled with heat transfer, however, was decoupled with fluid flow, while in Cacace and Jacquey (2017) the model was fully coupled. The decoupling was to simplify the modeling and increase the efficiency of the modeling.
Table 3 Calibration of permeability of equivalent porous medium based on pressure gradient between injection and production wells. Modeling type
Permeability (m2 )
Pressure gradient (kPa/ m )
Absolute differtence (kPa/ m )
Dual continuum matrix Fracture
10−15
73.0
0
Single continuum matrix
10−15
507.5
434.5
2 × 10−15
254.0
181.0
4 × 10−15
127.0
54.0
6 × 10−15
84.5
11.5
7 × 10−15
72.5
0.5
8 × 10−15
63.5
9.5
8.33 × 10−10
saltwater begins entering on the top aquifer and moves through a borehole to the bottom aquifer. The fracture interconnects the upper and lower aquifers, making a shortcut for governing density driven flow. The density driven flow and solute transport processes were affected by thermal buoyancy. First, the salinity and heat transport were modeled through discrete fracture modeling to match benchmark conditions. The salinity and heat transport were then modeled as an equivalent porous media and compared to the discrete fracture dual continuum results to validate the simulation. The contours of concentration and temperature are shown in Fig. 7.
4.2. Thermohaline aquifer-aquitard-aquifer A domain aquifer-aquitard-aquifer containing freshwater was exposed to an increasing thermal gradient. In the simulation, a heavy
Fig. 4. Contour of temperature after 50 years of stimulation; a) discrete fracture modeling (Cacace and Jacquey, 2017), b) discrete fracture modeling COMSOL. 5
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
itself, with 2-D and 3-D simulations varying by 2.5 days and discretization of the model varying by 2 days (Diersch, 2014). However, we were able to achieve a good calibration when comparing discrete fracture and equivalent porous media under the same model conditions (Fig. 8). Especially, in early time simulations, during the initial breakthrough, the curve fit between the fractured model and equivalent porous media approach was good. In late time, the increase in salinity concentration at the exit point at lower aquifer was faster in the equivalent porous media. This is due to the higher fluid velocity in the equivalent porous zone comparing to the fracture. As the porous zone with the effective porosity of n = 0.05 is replaced with the fracture, the fluid velocity was increased by the factor of 1/ n which resulted the faster breakthrough in comparison to the fracture model. Fig. 8b shows the log-log plot of concentration-time curve in late time simulation. It revealed the agreement between the equivalent porous medium modeling at K = 3 × 10−5m / s and discrete fracture modeling. As a result the solute transport in aquifer-aquitard-aquifer may alter by changing the hydraulic conductivity of the considered region in the middle of the domain. While in geothermal doublet model, the temperature gradient was only sensitive to the effective porosity of the replaced porous region (see Fig. 5). This showed that the sensitivity of a phenomenon in a fractured porous medium model depends on the behavior of the flow and solute/heat transport within the fractures and the fractures geometrical features (Long et al., 1982). Moreover, Fig. 8 showed that solute transport within a porous region with hydraulic conductivity of 3 × 10−5 m/s was equal to the solute transport within the fracture at early time simulation. However, due to the higher fluid velocity (i.e. Darcy velocity/porosity) within the porous region in comparison to the fluid velocity within the fracture, the breakthrough phenomenon was occurred faster. Therefore, at late time simulations, the salinity curve in equivalent porous media modeling (red line in Fig. 8) is above the salinity curve in discrete fracture modeling (solid black line in Fig. 8) with the error of 2%. Both results were converged after 50 days. Thus, it can be concluded that the porous region with the effective porosity of 0.05 and hydraulic conductivity of 3 × 10−5 m/s is equivalent to the existing fracture in the thermohaline aquifer-aquitardaquifer system with the maximum error of 2%. The representation of fracture flow in COMSOL was simulated as a porous media was calibrated to match discrete flow by adjusting only the permeability of the fractured area. By limiting the number of parameters that require adjustment we simplify the transition between discrete fractures and equivalent porous media. Variations in breakthrough curves are expected between equivalent porous media and discrete fracture models (Diersch, 2014), however total mass flux can be calibrated to be equal for both model types. Other simulations e.g. Daniela et al. (2014) have adjusted dispersivity to unrealistically high value to treat the velocity difference between porous media and fractures as an apparent dispersion process. The method of adjusting permeability and porosity provides reasonable calibration results for discrete fractures represented as equivalent porous media for a thermohaline reservoir. This method would be effective in reservoirs where limited calibration data are available.
Fig. 5. Production temperature-time plot: calibration of porosity of equivalent porous medium based on production temperature.
Fig. 6. The distribution of changes in porosity in horizontal slice of the model after 50 years.
To represent the fracture as a porous media, a porous zone with higher permeability and the effective porosity of 0.05, was defined at the center of the aquifer-aquitard-aquifer system with the thickness of 0.5 m to enhance the volumetric flow rate of the fracture. Similar to the synthetic geothermal doublet model, the sensitivity analysis of effective porosity was performed based on the variations in volumetric flow rate. Moreover, the effective porosity of 0.05 was obtained in a calibration process that provided the desired volumetric flow rate. The flow and transport through the fracture represented as a porous zone was calibrated using hydraulic conductivity values based on salinity concentrations at the exit point, where the wells meets the top of the bottom aquifer. Porosity, 0.05, and dispersivity, 0.1 m, were consistent with the corresponding values in discrete fracture modeling. The salinity concentration plots over 50 days for porous media with discrete fracture modeling and with results of Diersch (2014) are depicted in Fig. 8a (magnified in 8b), illustrating the calibration of permeability in the porous zone. Fig. 9 shows the contours of salinity concentration and temperature after 20 d for discrete fracture modeling and the equivalent porous media approach. The equivalent hydraulic conductivity was 3 × 10−5 m/s, based on comparison to the discrete fracture modeling results found in this study. The jump in the curves shows the breakthrough of the salinity through the aquitard to the lower aquifer. The initial breakthrough in the discrete fracture in this study was faster than (Diersch, 2014), approximately 1 day compared to 10 days (Figs. 7 and 8). A breakthrough time of approximately 5 days was observed for 5 × 10−6 m/s (Fig. 8), which was comparable to the 3-D case for discrete fractures in Diersch (2014), where the fracture hydraulic conductivity was 0.1 m/s. The breakthrough time was sensitive to the geometry and flow within the fracture
5. Conclusion In this paper, we presented simulations comparing equivalent porous media modeling approach to fully coupled discrete fracture (or dual continuum approach) for geothermal and thermohaline fractured reservoir applications. The efficiency and simplicity of modeling at the reservoir scale was enhanced by removing the fractures and substituting a matrix with equivalent hydraulic properties (i.e. porosity and permeability). Geological heterogeneities were not taken into account, however we were able to produce similar results to other modeling approaches that did take fractures into account e.g. Cacace and Jacquey (2017). Although the reliability of predications for reservoirs can be increased by considering accurate physical phenomena into account 6
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
Fig. 7. Contours of I) concentration of salinity, a) Diersch (2014), b) COMSOL; and II) temperature, c) Diersch (2014), d) COMSOL, at different times t = 5, 20 and 100 days in upper, middle, and lower row, respectively.
Fig. 8. Salinity concentration at exit point (50 m, 20 m); a) Concentration-time curve over 50 years, b) magnified of 8a between salinity concentration of 0.8–1.
7
Physics and Chemistry of the Earth xxx (xxxx) xxx–xxx
M. Jarrahi, et al.
Fig. 9. Contour of salinity concentration after 20 d, a) discrete fracture modeling, b) equivalent porous media approach.
References
(Cacace and Jacquey, 2017), simple equivalent porous media approaches can be used to produce accurate results, when adequate calibration data are available. According to the results, equivalent porous media models can be calibrated to reproduce THM and thermohaline results similar to dual porosity models by varying permeability and porosity. Equivalent porous media parameters were estimated using reservoir pressure and temperature data. The validation of this approach was confirmed by making the comparison between models of two benchmark problems in Cacace and Jacquey (2017) and Diersch (2014). Modeling the discrete fracture approach requires the detailed knowledge of fracture network in a fractured formation. However, by knowing the flow/heat transfer parameters such as volumetric flow rate and pressure/temperature gradient, that can be obtained from the field measurements, the equivalent porous media approach may be used to get numerical results that are in good agreement with the available field data. Nevertheless, the application of equivalent porous media approach instead of the discrete fracture modeling is not recommended for any fractured formation, as the nature of the fractures may not allow to find any equivalent porous region for the modeling. More details of the necessary conditions of the formations that may be modeled by equivalent porous media can be found in (Long et al., 1982).
Appelo, C.A.J., Van Loon, L.R., Wersin, P., 2010. Multicomponent diffusion of a suite of tracers (HTO, Cl, Br, I, Na, Sr, Cs) in a single sample of Opalinus Clay. Geochem. Cosmochim. Acta 74 (4), 1201–1219. Biot, 1973. Nonlinear and semilinear rheology of porous solids. J. Geophys. Res. 78 (23), 4924–4937. Blessent, D., Jorgensen, P.R., Therrien, R.e., 2014. Comparing Discrete Fracture and ContinuumModels to Predict Contaminant Transportin Fractured Porous Media. Groundwater, vol. 52. Cacace, M., Jacquey, A.B., 2017. Flexible parallel implicit modelling of coupled ThermalHydraulic-Mechanical processes in fractured rocks. Solid Earth 8, 921–941. Clauser, C., Ewert, M., 2018. The renewables cost challenge: levelized cost of geothermal electric energy compared to other sources of primary energy – review and case study. Renew. Sustain. Energy Rev. 82, 3683–3693. COMSOL, 2016. Comsol Multiphysics, User's Guide. COMSOL AB, Stockholm, Sweden. Daniela, B., Jorgensen, P.R., René, T., 2014. Comparing discrete fracture and continuum models to predict contaminant transport in fractured porous media. Gr. Water 52 (1), 84–95. Diersch, H.-J., 2014. Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media. Springer-Verlag Berlin Heidelberg. Dietrich, P., Helmig, R., Sauter, M., Hötzl, H., Köngeter, J., Teutsch, G., 2005. Flow and Transport in Fractured Porous Media. Springer-Verlag Berlin Heidelberg. Genuchten, v., 1985. A general approach for modelingsolute transport in structured soils. IAH Memories 17 (2). Harrison, B., Sudicky, E.A., Cherry, J.A., 1992. Numerical analysis of solute migration through fractured clayey deposits into underlying aquifers. Water Resour. Res. 28 (2), 515–526. Jain, C., Vogt, C., Clauser, C., 2015. Maximum potential for geothermal power in Germany based on engineered geothermal systems. Geotherm. Energy 3 (1), 15. Jarrahi, M., Holländer, H., 2017. Numerical simulation of flow heterogeneities within a real rough fracture and its transmissivity. Energy Procedia 125, 353–362. Long, J.C.S., Remer, J.S., Wilson, C.R., Wotherspoon, P.A., 1982. Porous media equivalents for networks of discontinuous fractures. Water Resour. Res. 18. Su, K., Latham, J.-P., Pavlidis, D., Xiang, J., Fang, F., Mostaghimi, P., Percival, J.R., Pain, C.C., Jackson, M.D., 2015. Multiphase flow simulation through porous media with explicitly resolved fractures. Geofluids 15 (4), 592–607. Therrien, R., Sudicky, E.A., 1996. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media. J. Contam. Hydrol. 23 (1), 1–44. Wang, H.F., 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton University Press, pp. 304. Xiong, Q., Baychev, T.G., Jivkov, A.P., 2016. Review of pore network modelling of porous media: experimental characterisations, network constructions and applications to reactive transport. J. Contam. Hydrol. 192, 101–117.
Acknowledgment and Permissions The authors would like to appreciate Dr.-Ing. Antoine Jacquey, Scientist at GFZ German research center for Geosciences, for providing required data and the permission to use the figures of (Cacace and Jacquey, 2017). The agreement with Springer Nature to reuse figures of (Diersch, 2014) is issued under license number 4643910403717. Additionally, the authors would like to thank the reviewers of this manuscript for their constructive comments and suggestions. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.pce.2019.08.001.
8