Comparison of numerical codes for coupled thermo-hydro-mechanical simulations of fractured media

Ahma Zariarmiyan, Hossin Salarira, Victor Vilarrasa, Kwang-Il Kim,Jawon L, Ki-Bok Min

a Department of Mining Engineering, Amirkabir University of Technology -Tehran Polytechnic (AUT), Tehran, Iran

b Institute of Environmental Assessment and Water Research (IDAEA), Spanish National Research Council (CSIC), Barcelona, Spain

c Associated Unit: Hydrogeology Group UPC-CSIC, Barcelona, Spain

d Department of Energy Resources Engineering, Seoul National University, Seoul, Republic of Korea

e Radioactive Waste Disposal Research Division, Korea Atomic Energy Research Institute (KAERI), Daejeon, Republic of Korea

ABSTRACT Geo-energy and geo-engineering applications, such as improved oil recovery (IOR), geologic carbon storage, and enhanced geothermal systems (EGSs), involve coupled thermo-hydro-mechanical (THM)processes that result from fluid injection and production. In some cases, reservoirs are highly fractured and the geomechanical response is controlled by fractures. Therefore, fractures should explicitly be included into numerical models to realistically simulate the THM responses of the subsurface. In this study, we perform coupled THM numerical simulations of water injection into naturally fractured reservoirs (NFRs) using CODE_BRIGHT and TOUGH-UDEC codes. CODE_BRIGHT is a finite element method(FEM)code that performs fully coupled THM analysis in geological media and TOUGH-UDEC sequentially solves coupled THM processes by combining a finite volume method (FVM) code that solves nonisothermal multiphase flow (TOUGH2) with a distinct element method (DEM) code that solves the mechanical problem(UDEC).First,we validate the two codes against a semi-analytical solution for water injection into a single deformable fracture considering variable permeability based on the cubic law.Then,we compare simulation results of the two codes in an idealized conceptual model that includes one horizontal fracture and in a more realistic model with multiple fractures. Each code models fractures differently. UDEC calculates fracture deformation from the fracture normal and shear stiffnesses, while CODE_BRIGHT treats fractures as equivalent porous media and uses the equivalent Young’s modulus and Poisson’s ratio of the fracture.Finally,we obtain comparable results of pressure,temperature, stress and displacement distributions and evolutions for the single horizontal fracture model. Despite some similarities, the two codes provide increasingly different results as model complexity increases. These differences highlight the challenging task of accurately modeling coupled THM processes in fractured media given their high nonlinearity.

Keywords:Coupled thermo-hydro-mechanical (THM)analysis Improved oil recovery (IOR)Naturally fractured reservoir (NFR)CODE_BRIGHT TOUGH-UDEC

1. Introduction

Thermo-hydro-mechanical (THM) processes take place in a number of scientific and engineering geo-applications, including improved oil recovery (IOR) (Wang et al., 2017), geologic carbon storage (Rutqvist et al., 2016), and enhanced geothermal systems(EGSs)(Majer et al.,2007).These geo-applications are performed in deep geological formations,which are usually fractured rock masses(Rutqvist and Stephansson, 2003). Coupled THM processes in fractured media involve deformation, heat and fluid diffusion, permeability enhancement and,eventually,plasticity and/or shear slip that may induce (micro-)seismicity (Min et al., 2004; Rutqvist et al.,2016; Vilarrasa, 2016; Norbeck and Horne, 2018; Vilarrasa et al.,2019a). The main cause of induced (micro-)seismicity is pressure buildup (Parotidis et al., 2004; Peacock et al., 2017; Rinaldi and Rutqvist, 2019). Additionally, since the injected fluids generally reach the storage formation at a temperature lower than that of the reservoir,effective stresses are reduced even further because of the temperature difference between the injected fluid and the reservoir,which produces a significant contraction of both the fracture and the surrounding rock matrix (Ghassemi, 2012; De Simone et al., 2013;Pandey et al., 2017; Paluszny et al., 2017; Yao et al., 2018; Rutqvist et al., 2019). Hence, geomechanical stability issues of naturally fractured reservoirs (NFRs) may arise due to a combination of pressure buildup and cooling (Kim and Hosseini, 2014; Vilarrasa et al., 2019b). Consequently, to accurately assess the geomechanical response of NFRs and investigate potential instabilities,hydro-mechanical (HM) and thermo-mechanical (TM) effects should be examined(Segall and Fitzgerald,1998;Li et al.,2019).

The behavior of fractures within NFRs is a critical issue for many geo-engineering applications (Rutqvist, 2015; Tomac and Sauter,2018). The role of geomechanics in fractured reservoirs is complex as a result of strong nonlinear interactions between the fluids(pressure and temperature changes) and the geomechanics of blocky systems (deformation and stress changes). The resulting block movements within NFRs can cause major alterations of the flow properties (Gale,1982; Salimzadeh and Nick, 2019). Pressure buildup induces shear displacement,which opens up fractures due to dilatancy,enhancing permeability(Barton et al.,1985;Yeo et al.,1998). In turn, permeability increase reduces pore pressure and propagates pressure perturbation along fractures more rapidly.Such nonlinear response of NFRs is not straightforward to predict and requires to be solved numerically.

Numerical models usually consider NFRs as equivalent continuous porous media in order to simplify an already complex coupled problem (Long, 1983; Nelson, 2001). If the fractures are not continuous and not interconnected or if they are highly interconnected, the reservoir can be considered as equivalent porous media.But if the characteristics of fractures are between these two extreme cases or the reservoir matrix is of very low-permeability,fractures entail a significant effect on the THM response of NFRs.Therefore, certain reservoirs, like carbonate oil reservoirs or crystalline rock, should not be considered as an equivalent continuum(Long,1983;Long et al.,1991;Narr et al.,2006).Simulation of THM processes in fractured media requires of a numerical code capable of modeling fluid and heat flow, as well as mechanical deformations, in both fractures and the surrounding porous media(Tsang,1999;Rutqvist and Tsang,2012;Zareidarmiyan et al.,2018).

Simulators based on continuum mechanics cannot easily consider the THM processes occurring in fractures because of their small thickness, among others. Nevertheless, fractures can still be modeled as an equivalent continuum by including thin elements that represent fractures with codes like CODE_BRIGHT(CB)(Olivella and Alonso, 2008) and OpenGeoSys (Wang and Kolditz, 2007;Watanabe et al., 2012). On the other hand, codes based on the distinct element method(DEM),e.g.universal distinct element code(UDEC)(Itasca,2014a),3DEC(Itasca,2014b)and PFC(Itasca,2014c),are capable of modeling the geomechanical behavior of fractures but present limitations to simulate interactions between fluid flow and fracture mechanical behavior. To overcome this limitation, a geomechanical code and a multiphase fluid flow code can be coupled to simulate THM processes(e.g.TOUGH-UDEC(TU)(Lee et al.,2019)).

Two different strategies, i.e. fully coupled and sequentially coupled, can be used to solve coupled THM problems (Kim et al.,2011). On the one hand, the fully coupled method solves the coupled equations simultaneously, providing unconditional stability and convergence for well-posed problems.On the other hand,the sequential coupling partitions the coupled problem and solves sub-problems sequentially. This sub-division implies smaller systems of equations to be solved than the fully coupled method.While the fully coupled method requires an integrated THM simulator, sequential methods typically provide flexible and efficient code development and the use of already available and developed codes (Kim et al., 2011).

The goals of this study are to investigate THM effects of fluid injection into NFRs and to compare two different coupling schemes and their simulation results, i.e. a fully coupled code, CB (Olivella et al., 1996), which treats fractures as equivalent porous media,and a sequentially coupled code,TU(Lee et al.,2019),which models fractures as discontinuities. TU has been selected among the very few available options for sequential coupling as an alternative to develop a single fully coupled code in continuum media including discontinuities. The ways in which each code considers and treats the fractures are not the same; however, both codes are based on the same balance equations for their calculations. Thus, this study improves our understanding of THM responses of NFRs by explicitly simulating fractures in our models and investigating their geomechanical response. Moreover, this study enables us to have a better understanding and evaluation of different coupling approaches, i.e. fully coupled and sequentially coupled. We first validate the two codes by comparing their results with those of a semi-analytical solution. Then, to address the challenging task of modeling coupled THM processes induced by fluid injection/production into fractured media, we initially model an idealized conceptual model that includes one horizontal fracture and subsequently,we increase model complexity by including two sets of fractures to represent a realistic NFR scenario. Finally, we compare the results obtained with the two codes and discuss the potential sources of differences.

2. Simulators

2.1. CODE_BRIGHT

CB is a finite element method (FEM) code that performs fully coupled THM analysis in geological media (Vaunat and Olivella,2018). This code solves five governing equations: stress equilibrium, liquid mass balance, gas mass balance, energy balance and balance of a conservative solute.Models can account for nonlinear problems, large displacements, dissolution process, advective and non-advective fluxes of species, and energy and fluid phase changes. As a result, this code includes a significant number of constitutive laws for defining both uncoupled and coupled behaviors for all the variables (Olivella et al., 1994, 1996; Vaunat and Olivella, 2018).

In CB, fractures are modeled as single fractures embedded in a continuous porous medium (Fig.1). Fractures are characterized by their aperture,b, spacing between fractures,a, and the intrinsic permeability of the fracture rock matrix,kfm,which is considered to be the unfractured porous material surrounding the fracture(Fig.1).The intrinsic permeability perpendicular to the fractures is equal to the fracture rock matrix permeability, and the equivalent intrinsic permeability of the embedded fracture in the direction parallel to the fractures is (Olivella and Alonso, 2008):

The aperture of the fracture can be calculated as a function of deformation as (Olivella and Alonso, 2008):

Fig.1. Schematic of embedded fracture(s) in a continuous porous medium in CB for (a) a single fracture and (b) multi fractures. The blue arrows depict the fluid flow within the fracture.

where ε is the volumetric strain and ε0is the threshold value for volumetric strain above which fractures open up.The permeability of fracture elements depends on the fracture spacing (a) as a characteristic size of the material, but it is independent of the element thickness (s)(see Eqs.(1)and(2)).In our simulations,we consider one fracture in each element, i.e.a=s(Fig.1a).

Assuming a linear elastic behavior of fractures, CB requires the Young’s modulus,the Poisson’s ratio,and the reference porosity of the fractures as input data (Vaunat and Olivella, 2018). Failure of fractures can also be modeled by using available elastoplastic and viscoplastic constitutive laws(Vilarrasa et al., 2017).

2.2. TOUGH-UDEC

2.2.1. General description

TU is a simulator that sequentially couples TOUGH2 and UDEC(Lee et al., 2019). TOUGH2 is a finite volume method (FVM) simulator for multicomponent and multiphase fluid and heat transport in fractured and porous media (Pruess and Narasimhan, 1985;Pruess et al., 2012). UDEC is a two-dimensional (2D) numerical software that simulates the response of loading on materials containing multiple, intersecting fractures (Itasca, 2014a). Both codes are well tested and widely applied in their respective fields.

UDEC represents NFRs as an assemblage of discrete blocks and treats the fractures as boundary conditions between blocks. Using the contact area model, UDEC provides a linear representation of discontinuity stiffness and yield limit based on normal stiffness(kn), shear stiffness (ks), friction, cohesion and tensile strength properties, and dilation of rock discontinuities (Itasca, 2014a).UDEC can model large block displacement along fractures, block rotations, nonlinear and complex behaviors, giving rise to fracture aperture changes (Fig. 2) (Israelsson, 1996; Itasca, 2014a). These fracture aperture changes calculated by UDEC are used by TOUGH2,which represents fractures as equivalent porous media, to update permeability (using the cubic law) and porosity as follows:

Fig. 2. Schematic of fracture modeling in UDEC. The blue arrow depicts the fluid flow within the fracture. Note that blocks can rotate, leading to variable fracture aperture along the fracture.

wherekeis the equivalent permeability and φeis the equivalent porosity. In TOUGH2, the spacing coincides with the fracture thickness.

2.2.2. Mesh definition

For linking TOUGH2 and UDEC, the model geometry and element numbering along fractures and within the matrix blocks should be identical and consistent in both codes. To this end, a MATLAB routine is used to generate the mesh in TOUGH2 with consistent element numbering from UDEC mesh. UDEC mesh uses zones and domains to represent matrix blocks and fractures,respectively. The elements, which are triangular in the zones(matrix blocks) and quadrilateral in the domains (fractures), are defined by their corner nodes (Fig. 3). In contrast, the numerical grid in TOUGH2 is defined with finite volume elements,which have one node located in the center of the element, being connected to the center nodes of the neighboring elements (Fig. 3).

The fracture domains in UDEC have small sizes, i.e. small apertures,compared with the zones in the matrix blocks.This small size could cause numerical problems in temperature and pressure calculation by TOUGH2.To overcome these problems,the thickness of the fracture elements (a) in TOUGH2 is assumed to be 1 m.

2.2.3. Coupling logic

The way TOUGH2 and UDEC are coupled is similar to TOUGHFLAC coupling logic (Rutqvist, 2011; Lee et al., 2019). Balance equations are solved sequentially and the primary and coupling variables are passed between the two codes at individual intervals(Fig. 4). First, TOUGH2 solves the TH equations and provides pressures,p, and temperatures,T, to UDEC zones and domains.Pressures are directly stored at the center of UDEC zones and domains while temperatures are interpolated from the TOUGH2 element centers to the UDEC corner nodes (Fig. 4) (see Lee et al.(2019) for more details). Then, UDEC solves the mechanical equation by considering the updated pressures and temperatures,and the updated apertures of each fracture domain are stored when equilibrium with minimum unbalanced forces is reached.After that, if variable fracture permeability is considered, the stored apertures are converted to the equivalent permeabilities and porosities (Eqs. (3) and (4), respectively) and transferred to the TOUGH2 for the next step calculation (Fig. 4).

Fig.3. Schematic of UDEC and TOUGH2 mesh definition.Note that while UDEC considers the real fracture aperture,the thickness of the fracture elements in TOUGH2 is assumed to be 1 m to overcome numerical problems.

Fig. 4. Schematic of TOUGH-UDEC coupling logic and data interchange.

3. Methods

3.1. Mathematical model

To solve the THM coupling, mass conservation of each phase(Bear, 1972), energy balance (Nield and Bejan, 2017), and momentum balance (Biot, 1956; Mctigue, 1986) have to be solved.Assuming that there are no chemical reactions, fluid conservation can be expressed as (Bear,1972):

wheretis the time,φ is the porosity,ρ is the fluid density,qis the volumetric flux of the injected fluid (volume of fluid per unit area per unit time), andris the source/sink term of fluid. Momentum conservation of the fluid phase in a porous medium is expressed using Darcy’s law:

wherekis the intrinsic permeability;μ is the dynamic viscosity of the injected fluid,which is a function of pressure and temperature;gis the gravitational acceleration; andzis the elevation.

The flow rate per unit width of a fracture can be calculated from Darcy’s law. Assuming that the permeability of fractures is controlled by fracture aperture and that fracture spacing(a)equals fracture aperture(Eq.(3)),the flow rate per unit width follows the cubic law as (Witherspoon et al.,1980):

Non-isothermal effects are governed by energy conservation as(Nield and Bejan, 2017):

wherecsandcare the specific heat capacities of the solid and the fluid, respectively; ρsis the solid density; λ is the equivalent thermal conductivity;andfQis the sink/source energy term.

To solve the mechanical problem,the momentum equation for a porous medium has to be satisfied. Assuming quasi-static conditions, i.e. neglecting inertial terms, momentum balance results in the equilibrium of stresses:

where σ is the stress tensor andbis the body forces vector.

Based on the theory of linear thermo-poroelasticity,stresses are a function of changes in strain, pressure, and temperature (Biot,1956; Mctigue,1986):

where Δσ is the total stress tensor;εVis the volumetric strain;Iis the identity matrix;Δε is the strain tensor;K=E/[3(1-2ν)]is the bulk modulus;Eis the Young’s modulus; ν is the Poisson’s ratio;G=E/[2(1+ν)] is the shear modulus; β is the Biot coefficient,which has been assumed to be 1; and αTis the linear thermal expansion coefficient. It is worth mentioning that cooling implies an isotropic reduction of stresses equal to 3KαTΔT.

Both fluid pressure and temperature changes cause stress and strain variations, inducing porosity and intrinsic permeability changes. The temperature distribution and evolution depend on how fluid flows through the fractures and matrix blocks, which in turn is controlled by temperature because of the dependency of fluid density and viscosity on pressure, but specially on temperature.These coupled non-isothermal effects are considered in Eqs. (6), (8) and (10).

Fig. 5. Schematic of fluid injection into a discrete fracture-rock matrix system, indicating the geometry and boundary conditions. Note that since the model thickness is 1 m, the fracture has been modeled with its real aperture.

3.2. Equivalency in the treatment of discontinuities

For the purpose of having comparable simulations in both TU and CB, use of equivalent properties for discontinuities (e.g. fractures,and faults) is a requirement. While UDEC considers discontinuities with the actual aperture, CB assumes thicker equivalent porous media.When dealing with discontinuities with the actual aperture,the mechanical behavior of discontinuities can be typified by normal and shear stiffnesses which can be written as(Yu et al.,2015):

whereknandksare the normal and shear stiffnesses of the fracture,respectively; σnis the total normal stress to the fracture; τ is the shear stress;andunandusare the normal and shear displacements,respectively.

If, instead of modeling the actual fracture aperture, an equivalent porous media with a thickness ofais used to represent a discontinuity, the required equivalent elastic parameters to reproduce fracture stiffness can be written as (Yu et al., 2015; Puig Damians, 2016):

whereEnis the oedometric modulus. Assuming that deformation normal to the fracture occurs without lateral deformation of the material, the modulusEnis defined as

Fromknandks, the equivalent Young’s modulus and Poisson’s ratio of the fracture can be determined, and vice versa, by combining Eqs. (13)—(15).

3.3. Fractures stability

To assess fracture stability, we consider the Mohr-Coulomb failure criterion (Jaeger et al., 2007). The numerical simulations calculate deformation and stress changes due to low-temperature fluid injection and isothermal fluid production. Linear thermoelasticity is assumed for the whole model and in order to evaluate the potential fracture instability and subsequent induced microseismicity, we perform a slip tendency analysis of preexisting fractures (Morris et al.,1996; Streit and Hillis, 2004). Failure conditions are given by

where τcris the critical shear stress;cis the cohesion; μ is the friction coefficient,which is often expressed in terms of the friction angle,φ,i.e. μ = tanφ;andis the effective normal stress.Shear failure and thus, sliding,occurs when the shear stress τ equals the critical shear stress τcr. For cohesionless fractures (c= 0), shear failure occurs when the ratio of the shear stress to effective normal stress equals the friction coefficient:

We use Eq.(17)to estimate the mobilized friction angle φmobon pre-existing fracture planes to quantify the shear slip tendency along these planes. The mobilized friction angle represents how close the stress state is to failure conditions. In particular, failure occurs when the mobilized friction angle equals the actual friction angle.

3.4. Validation

TOUGH2, UDEC, and CB have been extensively validated, i.e.TOUGH2 for multiphase fluid flow and heat transport, UDEC for geomechanical problems including mechanical behaviors of rock matrix and fractures,and CB for performing coupled THM analysis in geological media. Here, we aim at validating TU and CB for simulating water injection into a single fracture with variable permeability by comparing their results with the semi-analytical solution of Wijesinghe (1985). This semi-analytical solution has been used to verify other numerical codes, such as ROCMAS II(Noorishad et al., 1992), Geocrack (Swenson et al., 1995), FEHM(Bower and Zyvoloski, 1997), and OpenGeoSys (Watanabe et al.,2012).

Wijesinghe(1985)considered the cubic law to calculate fracture permeability and analyzed a one-dimensional transient fluid flow coupled with deformation in a single, uniform horizontal fracture confined by a very low-permeable rock in 2D space(Fig.5).Initially,fracture aperture is uniform and equal to 1× 10-5m, and fluid pressure is 11 MPa all along the fracture.The fracture is subjected to a uniform isotropic stress equal to 50 MPa. Fluid is injected at the leftmost edge of the fracture with 0.9 MPa overpressure. Thus, a constant pressure boundary condition is fixed to 11.9 MPa at the fracture on the left boundary and a constant pressure equal to the initial pressure, i.e. 11 MPa, is prescribed on the right boundary.Material parameters are listed in Table 1.

Figs.6 and 7 depict the simulation results compared with semianalytical solution for pressure and fracture aperture along thefracture, respectively. Both TU and CB match well with the semianalytical solution. The pressure increase caused by the constant pressure boundary condition gradually propagates toward the lowpressure edge of the fracture (right side) (Fig. 6), which causes fracture opening(Fig.7)and permeability enhancement,which,in turn, affect pressure evolution (Fig. 6). This nonlinearity is well captured by both codes under the conditions of this semi-analytical solution.

Table 1Material properties used in the code validation.

3.5. Benchmarking exercises

3.5.1. Single fracture model

We consider a 200 m×100 m 2D plane strain model with three horizontal layers that represent the caprock,reservoir,and bedrock,with thicknesses of 20 m,60 m,and 20 m,respectively(Fig.8).We model one horizontal fracture embedded in the middle of the reservoir matrix. All layers are saturated with water. The initial conditions correspond to hydrostatic pressure, a geothermal gradient of 25°C/km with a surface temperature of 20°C and a normal faulting stress regime.The considered stress regime is such that the vertical stress,which equals the lithostatic stress,is greater than the horizontal stresses,following σh= 0.8σv.Table 2 lists the THM properties of the materials considered in the model. Recall Section 3.2 for the equivalency of the fracture properties between the two models. We assume that fracture permeability remains constant and equal to 10-14m2to minimize complexity and nonlinearity and make this simple case as comparable as possible.Note that the permeability of the reservoir rock matrix is four orders of magnitude higher than that assumed in the example for validation(Section 3.4)and thus,pressure will diffuse into the rock matrix surrounding the fracture.

Fig. 6. Pressure profiles along the fracture for the semi-analytical solution of Wijesinghe (1985), TU and CB after 500 s and 2000 s of water injection.

The mechanical boundary conditions are the lithostatic stress on the upper boundary, no displacement normal to the lateral boundaries and fixed displacement on the lower boundary. All boundaries are considered as no-flow boundaries.Water is injected for two months into the fracture on the left boundary at a constant overpressure of 6 MPa and a temperature of 30°C, which means that it is 21.25°C colder than the reservoir (Fig. 8).

3.5.2. Multi-fracture model

Model complexity is increased by modeling a scenario representative of a real NFR(Fig.9).We consider a model with two sets of orthogonal fractures, with dip angels of -30°and 60°degrees,embedded into the reservoir matrix.The initial conditions and the material properties of the model are the same as in the single fracture model (Fig. 9 and Table 2). The main difference is that fracture permeability depends on fracture aperture in this model.The initial permeability is 10-14m2and it can increase up to 10-13m2as a result of fracture opening. Water is injected into a fracture on the left boundary for 1 year at a constant pressure buildup of 1.5 MPa and a temperature of 30°C.Additionally,water is produced simultaneously at a constant pressure drop of-1.5 MPa from a fracture on the right boundary (Fig. 9).

Each code uses its own mesher.Since the geometry is complex,it is difficult to reproduce the same meshes in both codes. The fractures are discretized with 1 m-thick element both in CB and TOUGH2. The mesh in CB is made of unstructured quadrilateral elements,with refinement around fracture intersection.The size of the elements in the fracture is 1 m and it progressively increases further away, reaching a size of 5 m in the reservoir, caprock, and bedrock.On the other hand,the mesh is structured in TU,with the blocks being discretized into triangular finite difference zones with the maximum edge length of 5 m.

4. Results

4.1. Single fracture model

Fluid injection causes fluid pressure increase (Fig. 10) and cooling around the injection well (Fig. 11), which reduce the effective stresses, causing deformation (Fig. 12). While pressure increase produces an expansion of the fracture and reservoir,cooling induces contraction.Simulation results of CB and TU have a reasonable fitting (Figs. 10—13). The pressure propagation front along the fracture is identical after 60 d of cold water injection in both codes. However, it should be noted that for previous times,pressure diffuses more rapidly in CB than in TU.This faster pressure diffusion within the fracture generates a higher pore pressure in the matrix (reservoir) in CB than in TU (Fig. 10b). The pressure differences are 0.8 MPa at the top and bottom of the reservoir near to the injection point and 0.2 MPa at 175 m away from the injection point. These differences in pressure distribution may result from dissimilar storage coefficient calculations between the fully coupled approach, which reproduce the non-local nature of storage, and the sequentially coupled approach, which fails to include this effect (De Simone and Carrera, 2017) (see the Discussion for more details).

Fig. 7. Fracture aperture along the fracture for the semi-analytical solution of Wijesinghe (1985), TU and CB after 500 s and 2000 s of water injection.

Fig.8. Schematic of the single fracture model geometry and boundary conditions.The dark red line indicates the position of the fracture,which is embedded in the middle of the reservoir rock.

The higher pore pressure buildup in the matrix in CB than in TU causes a slightly higher cooling in CB than in TU (Fig. 11) and a higher vertical displacement at the top of caprock in CB than in TU because of higher expansion of the reservoir matrix(Fig.12).Total normal stress profiles along the fracture are very similar (Fig.13).The decrease in the total normal stress on the left side of the profiles, i.e. near to the injection point, is the result of the thermal stress reduction caused by cooling. Further away, outside of the cooled region,the total stress increase is due to the poromechanical response of the rock.

4.2. Multi-fracture model

4.2.1. Fluid pressure evolution

Changes in fluid pressure as a result of injection/production propagate much faster in the fractures than in the reservoir matrix due to the large permeability contrast between them. During the first days of operation, pressure changes mainly occur in the fractures (Figs.14 and 15). In the long term, pore pressure within the matrix blocks equilibrates with that in the fractures(Fig.15b and d).Since we consider fracture permeability changes as a result of deformation,pressure diffusion along the reservoir is controlled by these changes because the pressure propagates faster through the fractures with higher permeability (see Section 4.2.2).

Fig. 14 compares the fluid pressure changes as a function of distance in two fractures close to the injection point with dip angles of-30°and 60°(fractures A and B,respectively,in Fig.9)and in the fracture containing the production point with dip angle of -30°(fracture P in Fig. 9) after 7 d and 365 d of operation. The abrupt change of the pressure gradient along the fractures coincides with the fracture intersection.Fractures A and B are pressurized faster in CB than in TU(Fig.14a and c)as a result of a more rapid increase in their permeability within the first 7 d of operation(Fig.16a and b).

The injected flow rate, given that the boundary condition prescribes the injection pressure,is higher in CB that in TU because of this higher increase in fracture permeability. In contrast to liquid pressure profiles along fractures A and B (where injection takes place), the liquid pressure profiles along fracture P (where production takes place)after 7 d of operation (Fig.14e)are almost equal in both codes due to the absence of changes in the fracture permeability and therefore, in the production flow rate(Fig.16f).

When the fluid flow is in steady state after 1 year of operation(Fig. 15b and d), the pressure change profiles in both codes for fractures A and B get closer (Fig. 14b and d). Actually, they are almost the same for fracture A. However, the pressure change profiles for fracture P show an almost constant pressure in TU(Fig.14f) due to the lower fracture permeability within the central part of the reservoir (Fig.16e), which leads to a lower production flow rate.

Fig.15 compares the liquid pressure contours in both codes after 7 d and 365 d of operation. Fig. 15a and c shows that while the pressure buildup front propagates for approximately three quarters of the reservoir length in CB after 7 d of injection,it only propagates for a quarter of the reservoir length in TU.In contrast,the pressure in the production side of the models drops in a similar manner.At the beginning of injection, pressure diffuses through fractures causing strong pressure gradients toward the matrix blocks.Despite the low-permeability of the matrix blocks, fluid diffusion into them occurs relatively fast as a result of their high stiffness relative to the fractures (Table 2), which results in a small storage coefficient.Thus,after approximately one month,the matrix blocks reach equilibrium with the pressure in the fractures and cause a homogenized pressure distribution in the long term (steady-state conditions). Consequently, the presence of fractures cannot bedistinguished from the pressure contours after 1 year of operation(Fig.15b and d).

Table 2Material properties used in the THM analyses of the single and multi-fracture models.

Fig. 9. Schematic of the multi-fracture model geometry (with two sets of orthogonal fractures) and boundary conditions.

Fig.10. Comparisons of liquid pressure profiles within the reservoir after 2 months of cold water injection(a) in the horizontal direction along the fracture,and(b)in the vertical direction at distances of 15 m, 50 m,100 m, and 175 m from the injection point.

Fig.11. Comparisons of the temperature profile in the reservoir after 2 months of cold water injection (a) in the horizontal direction along the fracture for 50 m closer to the injection point, and (b) in the vertical direction at distances of 1.5 m, 4.5 m, and 7.5 m from the injection point.

4.2.2. Permeability changes

Fig. 16 depicts the permeability changes as a function of the effective normal stress in different checkpoints within the fractures. The permeability-effective normal stress curves do not coincide for the two codes. Permeability increases for smaller changes in the effective normal stress in TU than in CB, but permeability is enhanced in larger portion of fractures in CB than in TU. Permeability changes in both fractures A and B reach to their maximum permeability within the first day of injection in both codes (checkpoints 1 and 2) (Fig. 16a and b). The next fracture parallel to fracture A (checkpoint 3) reaches the maximum permeability within the first day of injection in CB whereas it takes 11 d in TU(Fig.16c).The block movements and rotations in TU cause only slight fracture permeability changes in the central part of the reservoir,showing an initial permeability increase that is reversed afterwards (Fig. 16e). In contrast, fracture permeability is rapidly enhanced in CB in the central part of the reservoir(Fig.16e).Fig.16f depicts that in the production part of the reservoir,the permeability remains constant at the minimum value, i.e.10-14m2.

4.2.3. Temperature evolution

Fig. 12. Vertical displacement at the top of the caprock after 2 months of fluid injection.

Fig.17 compares the temperature contours in both codes after 1 year of operation. Differences between the codes arise from the different permeability evolutions within the fractures in the reservoir(Fig.16),which leads to different injection flow rates and thus, heat advection. Since permeability enhancement is larger in CB than in TU, cooling advances more rapidly in CB. The injected fluid, which has a lower temperature than the reservoir, moves away from the injection point through the fractures and gradually cools down the reservoir. The cooling front advances extremely behind the pressure front because it needs to cool down fractures and the surrounding rock matrix. Moreover, the flow of cold fluid divides into two fractures at fracture intersections, significantly reducing heat advection away from the intersections (Fig.17b).

4.2.4. Total normal stress changes along fractures

Fig.18 compares the total normal stress changes as a function of distance along fractures A,B,and P.The pressure and temperature variations consistently cause changes in the total stresses. While the main driving mechanism that changes total stresses is pressure variations in the short term,it is cooling in the long term(Fig.18a—d). Cooling advances much more slowly than pressure, thus temperature changes are small in the initial stages and only advances through the fractures in the vicinity of the injection point in the long term (Fig.17). The higher heat advection within the reservoir in CB after 1 year of operation results in a higher decrease in the total normal stress both in fractures A and B (Fig. 18b and d,respectively).The spikes observed in Fig.18,which are much more pronounced in TU than in CB, correspond to the fracture intersection.

4.2.5. Vertical displacement

Fig.19 shows the vertical displacement at the top of the caprock after 7 d and 1 year of operation. The difference between the two codes increases with distance from the injection point as a consequence of the lower permeability changes (Fig.16) and hence, the lower flow rate and fluid pressure in the central part of the reservoir in TU than in CB.Thus,the expansion of fractures and matrix is larger in CB than in TU, resulting in a higher vertical displacement at the top of the caprock.

Fig. 20 compares the vertical displacement contours in both codes after 7 d and 1 year of operation. The contours depict qualitative resemblance;however,they are quantitatively different.This difference is a consequence of the different permeability evolutions within the fractures in the reservoir (Fig. 16), which leads to different pressure distributions within the reservoir (Figs. 14 and 15). Despite the differences in the absolute values, both codes clearly simulate the effect of fractures, which results in a noncontinuous displacement field.

4.2.6. Fractures and caprock stability

Fig.13. Total normal stress profile along the fracture after 2 months of fluid injection.

The simulated scenario considers moderate changes in both pressure and temperature, which limits the effective stress changes. Thus, only relatively small changes in fracture stability occur during operation.Both codes predict similar fracture stability changes(Fig.21).Fracture P undergoes a slight increase in stability,i.e. negative mobilized friction angle change, as a result of the decrease in pressure caused by fluid production (Fig. 21c). In contrast,fracture stability decreases,i.e.positive mobilized friction angle change, close to the injection point in fracture A because of the thermal stress reduction induced by cooling and the decrease in effective stresses caused by pressure buildup (Fig. 21a). Fig. 21b shows the stability in fracture B, which presents a non-intuitive response. Fracture stability increases next to the boundary, up to the first fracture intersection. This increase in stability in this section of fracture B is due to the compression of the matrix block above it, as a result of the expansion of fracture A, which tends to close fracture B.However,further away from the boundary,on the other side of the fracture intersection, the expansion of fracture A tends to cause opening or sliding of the matrix block along fracture B, thus decreasing fracture stability. Note the sharp change in fracture stability along fracture B predicted by both codes around the intersection with fracture A (Fig. 21b). Thus, the interactions between fractures and the relative movements of matrix blocks play a major role in fracture stability.

Caprock stability is guaranteed because its low permeability and slow advance of the cooling front induce negligible pressure and temperature variations in the caprock. Thus, the effective stress changes within the caprock are also negligible. This situation ensures a stable geomechanical response of the caprock during operation.

5. Discussion

Fluid pressure changes as a result of injection/production in NFRs preferentially advance through fractures, which are several orders of magnitude more permeable than the reservoir matrix.Where overpressure occurs, the effective stresses are decreased,which opens up fractures, increasing their permeability. Such overpressure acts uniformly in all directions and tends to expand the fractures which are confined, implying an increase in total stresses,especially in the plane of the fracture.Thus,the decrease in effective stresses is smaller than the increase in fluid pressure(Hillis, 2000; Zareidarmiyan et al., 2018). Additionally, lowtemperature injected fluids cool down fractures and a portion of the surrounding matrix blocks contracts the rock and causes a reduction of both total and effective stresses (De Simone et al.,2013; Zareidarmiyan et al., 2018). Furthermore, low-temperature injection increases fluid viscosity, and as a consequence, produces a higher overpressure. This thermal effect results in a larger effective stress reduction in comparison with HM simulations. These interactions highlight the necessity to explicitly account for coupled THM processes to accurately assess fractures and caprock stability in geo-energy applications. However, these effects are sometimes ignored in conventional methods by simplifications such as applying the variation of pore pressure on the initial stress state to calculate effective stresses (e.g. Karvounis et al., 2014) or considering an equivalent porous media instead of explicitly modeling fractures (e.g. Long,1983; Nelson, 2001; Rutqvist et al.,2008; Otto and Kempka, 2015).

Fig.14. Comparisons of the profiles of liquid pressure changes within the injection and production fractures after 7 d and 365 d of cold water injection along fracture A after(a)7 d and(b)365 d of operation,along fracture B after(c)7 d and(d)365 d of operation,and along fracture P after(e)7 d and(f)365 d of operation.Note that the changes in the pressure gradient coincide with intersections with other fractures.

Fig.15. Comparisons of liquid pressure contours for CB after (a) 7 d and (b) 365 d of operation, and TU after (c) 7 d and (d) 365 d of operation.

Fig.16. Comparisons of permeability evolution as a function of the effective normal stress to the fractures during 1 year of cold water injection at checkpoints(a)1(fractures A),(b)2(fracture B),(c)3,(d)4,(e)5,and(f)6(fracture P).Dark colors of the lines indicate early times and light colors late time.Note that the initial fracture permeability for both codes in all checkpoints is equal to 10-14 m2.

Fig.17. Comparisons of temperature contours within the reservoir after 365 d of cold water injection in (a) TU and (b) CB. The dashed lines depict fractures.

In the single fracture model (Section 4.1), which assumes constant fracture permeability, the liquid pressure diffusion and the cooling front advance from the fracture into the reservoir matrix are slightly faster in CB than in TU (Figs.10 and 11). It should be noted that when both codes are solved in purely hydraulic mode.i.e. comparison between TOUGH2 and CB, the obtained pressure distribution with both codes is identical.However,when including the HM coupling, the non-local nature of storage (De Simone and Carrera, 2017) arises in the fully coupled simulation of CB, but it is not reproduced with the sequential coupling of TU.The non-local effect of storage induces an instantaneous pressure increase in the whole domain due to the poromechanical response to fluid injection, which induces compression of the rock. The smaller the reservoir, the larger the magnitude of this effect (De Simone and Carrera, 2017). This non-local effect becomes also larger for noflow than constant pressure boundary conditions. Since we consider a small domain and a no-flow condition in the outer boundary, the differences become non-negligible. The non-local storage may be reproduced in TU by updating porosity in the whole domain for the non-isothermal flow calculations with TOUGH2(Blanco-Martín et al.,2017).The two coupling approaches,i.e. fully coupled and sequentially coupled, may yield the same results if several iterations (up to 50) are performed in the sequentially coupled code, i.e. TU (Preisig and Prévost, 2011).Nevertheless, achieving the same results in both codes may be complicated if variable permeability fractures are considered because of the high nonlinearity of this process.

Fig.18. Comparisons of total normal stress changes(Δσn) profiles along fracture A after (a) 7 d and (b) 365 d of operation, fracture B after (c) 7 d and (d) 365 d of operation, and fracture P after (e) 7 d and (f) 365 d of operation.

Fig.19. Comparisons of vertical displacement at the top of the caprock after (a) 7 d and (b) 365 d of operation.

Fracture permeability evolution in the multi-fracture model is quite different in the two models (Section 4.2). The permeability changes of both codes are based on the cubic law. Even though aperture changes and then permeability changes are calculated from fracture strain in CB,in view of effective normal stress changes in TU, simulation results should be the same because strain and effective normal stress are related through linear elasticity.In fact,the permeability versus effective normal stress curve is practically identical for the validation example (Section 3.4). The differences that appear in the multi-fracture model may be a result of the much softer (by two orders of magnitude) fractures (Fig. 16). Another source of difference may be the coupling approaches of each code,i.e.full(CB)and sequential(TU)coupling.The resulting differences in calculation of the reservoir block movements by the two codes affect the pressure diffusion, cooling front advance, vertical displacement, and total normal stresses along the fractures(Figs.14,15, and 17-20). Such differences in the simulation results highlight the challenges of accurately modeling NFRs, which are highly nonlinear systems. This high nonlinearity enhances small differences between the two codes, which ends up in nonnegligible differences between the simulation results. To reduce the uncertainty in modeling accuracy,controlled field experiments,such as pulse tests and cross-hole tests,with detailed monitoring to measure fracture permeability evolution are necessary (Martinez-Landa and Carrera,2005).

Fig. 20. Comparisons of vertical displacement contours within the reservoir in CB after (a) 7 d and (b) 365 d of operation, and TU after (c) 7 d and (d) 365 d of operation.

Fig. 21. Comparisons of the mobilized friction angle change profiles along fractures (a) A, (b) B, and (c) P after 365 d of operation.

The mechanical properties of the fractures require different input parameters in CB and TU, which makes the comparison and identification of the sources of differences not straightforward. On one hand, CB treats fractures as equivalent porous media with a thickness of“a”and uses the Young’s modulus and Poisson’s ratio of the fracture to account for fracture stiffness.On the other hand,TU explicitly accounts for distinct fractures in the mechanical part of the code, calculating displacements along discontinuities and rotations of blocks through the definition of the normal and shear stiffnesses of the fractures. Equivalent properties can be obtained from the equivalency defined in Section 3.2 (Eqs. (13)—(15)).Actually, both approaches are able to reproduce identical results under certain conditions (recall the validation against a semianalytical solution in Section 3.4).

Fig.22. Comparisons of vertical displacement at the top of the caprock for(a)the single fracture and(b)the NFR models using normal and shear fracture stiffness values divided by two in TU.

Both codes produce a discontinuous displacement field across fractures, highlighting the necessity of including fractures in numerical models to accurately reproduce the THM response of NFRs to fluid injection and/or extraction (Fig. 20). In terms of the magnitude of the displacement, even though the equivalent properties for fractures are used in both codes (Eqs. (13)—(15)), the vertical displacement at the top of the caprock predicted by the two codes is different in both models(Figs.12 and 19),as a result of the different pressure and temperature fields in both codes. We have conducted an additional simulation for the two scenarios considering half of the initial value of normal and shear stiffnesses in TU.Interestingly, for the single fracture model, the vertical displacement at the top of the caprock in both codes becomes very similar(Fig. 22a). However, in the multi-fracture scenario, the fitting only improves near the injection boundary(Fig.22b).In principle,when there is injection and production into a reservoir, it is expected to observe uplift above the injection side and subsidence above the production side.This behavior is clearly observed in TU because of the lower fracture permeability in the central part of the reservoir,which causes half of the reservoir to undergo pressure depletion and the other half pressure buildup. In contrast, the high permeability of the fractures within most part of the reservoir in CB restricts the zone with pressure depletion to a small portion of the reservoir,leading to an overall uplift of the system.The differences in the NFR scenario highlight the fact that as model complexity increases, discrepancies may become non-negligible between codes because of differences in the simulation of THM processes in fractured media,which are strongly coupled and highly nonlinear.

Finally, it should be highlighted that equivalent porous models may fail to reproduce the THM response of fractured media. Basically, the properties of equivalent porous models are sensitive to fracture geometry and distribution. Even though fractures geometry and distribution in NFRs are inherently uncertain and subjected to debate, some features of the THM response of fractured media can only be reproduced when actually including fractures in the models (e.g. Yoshioka et al., 2019). Thus, it may be possible, in theory,to calculate the properties of an equivalent porous medium;however, they may not serve to reproduce the THM processes occurring in fractured media.

6. Conclusions

We conducted simulations of fluid injection and production into NFRs to understand their THM response by comparing two codes with different coupling approaches(a fully coupled code,CB,and a sequentially coupled code, TU). TU and CB were validated by comparing their simulation results to a semi-analytical solution for water injection into a single fracture with variable permeability surrounded by very low-permeable rock matrix in 2D space. Both TU and CB match well with the semi-analytical solution,accurately reproducing the nonlinear response of fracture opening as a result of pressure buildup and subsequent permeability enhancement which, in turn, affects pressure diffusion. Next, we considered a simple model with one horizontal fracture embedded into the reservoir matrix in order to compare the performance of the two codes under controlled conditions. Finally, we modeled a scenario that represents an idealized case of a real NFR. The main conclusions can be drawn as follows:

(1) Both TU and CB can accurately reproduce the HM behavior of a fracture surrounded by very low-permeable rock,as shown by the excellent matching with the semi-analytical solution.

(2) While the normal and shear stiffnesses of the fractures should be defined in TU, CB uses the equivalent Young’s modulus and Poisson’s ratio to compute fracture stiffness.Though equivalent properties can be theoretically obtained and may yield comparable results in simple scenarios,simulation results differ as model complexity increases.

(3) When modeling a NFR,which is characterized by conductive fractures surrounded by matrix blocks of relatively high permeability, the different approaches of TU and CB in modeling discontinuities lead to differences in the predicted fracture permeability evolution in the central part of the reservoir. While TU predicted a limited increase in fracture permeability in that part,CB predicted a rapid increase there,leading to a higher pressure buildup and a higher circulating flow rate and thus, higher vertical displacement and heat advance in CB.

(4) In the multi-fracture model, despite some similarities (like the change in fracture stability),the two codes provide quite different results in terms of pressure, permeability, temperature, normal stress and vertical displacement, highlighting the challenging task of accurately modeling the highly nonlinear nature of fractured rock masses.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

Ahmad Zareidarmiyan acknowledges the financial support received from the “Iran’s Ministry of Science Research and Technology” (PhD students’ sabbatical grants). Victor Vilarrasa acknowledges funding from the European Research Council under the European Union’s Horizon 2020 Research and Innovation Program through the Starting Grant GEoREST (www.georest.eu) (Grant Agreement No.801809). Ki-Bok Min acknowledges the support by the Korea-EU Joint Research Program of the National Research Foundation of Korea through Grant No. NRF-2015K1A3A7A03074226 funded by the Korean Government’s Ministry of Science and Information and Communication Technology (ICT) in the framework of the European Union’s Horizon 2020 Research and Innovation Program(Grant No.691728)and Research Institute of Energy & Resources, Seoul National University.