Thermo-hydro-mechanical modeling of fault discontinuities using zero-thickness interface element

Ali Rnjr, Hossein Hssni, Kourosh Shhrir, Mohmmd Jvd Ameri Shhri

a Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Iran

b Department of Mining and Metallurgy Engineering, Amirkabir University of Technology, Tehran, Iran

Keywords:Thermo-hydro-mechanical (THM)simulation Geomechanical coupling Zero-thickness element Joint element Finite element

A B S T R A C T In this paper, a coupled thermo-hydro-mechanical (THM) simulation in a faulted deformable porous medium is presented.This model involves solving the mass conservation,linear momentum balance,and energy balance equations which are derived from the Biot's consolidation theory. Fluid pore pressure,solid displacement, and temperature are chosen as initial variables in these equations, and the finite element method in combination with the interface element is used for spatial discretization of continuous and discontinuities (fault) parts of the medium to solve the equations. The main purpose of this study is providing precise formulations, applicability, and ability of the triple-node zero-thickness interface element in THM modeling of faults. It should be noted that the system of equations is solved using a computer code written in Matlab program.In order to verify the developed method,simulations of index problems such as Mandel's problem, and coupled modeling of a faulted porous medium and a faulted aquifer are presented. The modeling results obtained from the developed method show a very good agreement with those by other modeling methods, which indicates its accuracy.

1. Introduction

The discontinuities observed in rocks,depending on the size and characteristics of their motions, are called fractures or faults. This kind of discontinuities has a scale range from a few centimeters to several meters (Davis et al., 2011). Studying fractures and faults as linear-plate structures is very important and necessary in many engineering applications.In this regard,their importance in drilling of underground aquifers could be noted(Singhal and Gupta,2010).Also, their recognition in geotechnical engineering is of high importance such as construction of dams, tunnel drillings, and so on(Barzegari,2017;Lei et al.,2017).In petroleum engineering and geomechanics,discontinuities of hydrocarbon reservoirs have been studied in recent decades for various purposes,including the study of the importance of discontinuities in the storage of carbon dioxide(Figueiredo et al., 2015; Li and Laloui, 2016), cap rock integrity(Karimnezhad et al., 2014; Khan et al., 2019), hydraulic fracturing(Rutqvist et al.,2013), etc.

The first form of deformation in rocks induced by stresses is elastic deformation. The rocks can break down after elastic deformation stage, or they can also behave plastically, due to the high temperature and confining pressures.Rocks are broken generally at the final stage and/or after the end of the plastic deformation(Ghosh,2013).Fracture in rocks can be induced by shear and tensile stresses. The shear fractures are the results of applied stresses causing slip in rocks,in a form of rock displacement and separation,across the surface of the fractures(Fossen, 2016).

Based on the position of the stress on the faults, two different behaviors in terms of fluid flow could be considered.The faults can act as a barrier or a route,leading to oil traps or migration.In many observations and studies, faults have been reported as a barrier to fluid flow. It should be noted that in both of these scenarios, the behaviors of the fault in the oil industry are very important at the exploration and production stages. Fig. 1 schematically shows various behaviors of the faults in terms of fluid flow in hydrocarbon reservoirs. In this paper, based on conventional conditions, it is assumed that at the initial states, the faults have a very small permeability and small opening so that they can be assumed as a barrier to fluid flow.

Fig.1. Schematic representation of different fault behaviors in terms of fluid flow in hydrocarbon reservoirs.

In another point of view, the interactions of the mechanical,hydraulic, and heat domains that are completely interconnected will affect the behavior of the fault.For example,the pore pressure variations at the fault aperture will cause a change in the displacement of its surfaces (whether as slipping or opening) and ultimately affect its fluid transmissibility and permeability.Conversely,the change in permeability of the fault affects the pore pressure. Therefore, to model this thermo-hydro-mechanical(THM) behavior of faults, a coupled simulation is required that can take all these three factors into consideration.

Since the analytical models are unable to consider the complex THM coupling, numerical methods are generally used. Jalali and Dusseault (2012) introduced and compared various numerical methods in solving the governing equations for continuous and discontinuous porous media.In general,numerical methods can be classified into three categories, i.e. continuous, discontinuous, and hybrid (Jing, 2003). If we consider the faults as a discontinuity,large-scale rotations and separations could not be modeled by continuous methods.This is due to the fact that in these methods,the adaptability of the displacement between the internal elements is required. However, in continuous methods, the modeling of discontinuities could be done using equivalent methods in which fractures and faults are expressed by special interface elements such as zero-thickness elements. In discontinuous methods,displacement adaptation is not required and there is a possibility of complete rotation and separation in the mesh network. Hybrid methods combine the two methods to overcome the complex problems.The advantages and disadvantages of numerical methods in coupled simulations of the fault's behavior in porous media were presented by Gómez Castro (2018).

For a coupled THM modeling of a porous medium with fault discontinuity, the governing equations for fluid folw, deformation and heat transfer such as mass conservation, momentum balance,and energy balance equations, which include partial derivative equations, should be solved in both continuous and discontinuous domains (Lewis and Schrefler, 1998). Lewis and Schrefler (1998)showed how to use finite element method (FEM) to simulate coupled fluid flow and deformation in a deformable porous medium. Murad and Loula (1992,1994) studied the Biot's consolidation problem using FEM. Schrefler et al. (1995) and Gawin et al.(1996) developed a coupled numerical model based on Biot's formulation in consideration of heat component and simulated the THM behavior of a deformable porous medium. Rutqvist et al.(2001) evaluated the effectiveness of coupled THM processes in various cases. In recent years, several studies in relation to such modeling are also reported (e.g. Fang et al., 2013; Rutqvist et al.,2014; Wang, 2014; Li and Laloui, 2016; Li et al., 2016).

As stated, one of the effective methods for modeling faults and fractures in porous media is the use of interface element in conjunction with FEM. Wriggers and Zavarise (2004) investigated the mechanical behaviors of contact surfaces using various methods.Cappa and Rutqvist(2011)studied and compared several kinds of interface elements in the application of injection of CO2to deep underground. Among interface elements, the zero-thickness element introduced by Goodman et al. (1968) is developed and used by many researchers(e.g.Noorishad et al.,1982,1992;Segura and Carol,2008a;Cerfontaine et al.,2015a;Garolera et al.,2015)to model THM behavior of faults. Based on this, there is also a possibility to improve the zero-thickness element when modeling the normal and slip displacements. The main advantage of using the interface elements is the full compatibility of its mesh with the main mesh of the finite elements so that there is no need for remeshing the medium, since the interface element nodes are exactly located on the finite element nodes. A conventional zerothickness element has a pair of linear (for two-dimensional (2D)models)or planar(for three-dimensional(3D)models)elements so that both nodes facing each other are placed spatially (Fig. 2). The axes of local coordinates in these elements as other studies(Ferronato et al., 2008) have chosen such that the first axis is perpendicular to the element surface,the second axis is tangent to the element surface(in 2D case),and the second and third axes are inside the element surfaces and selected based on the clockwise rotation (in 3D case).

This element is basically used in single-,double-and triple-node forms in view of discontinuity modeling. Fig. 3 shows these three forms. This element, in double-node, considers the fluid flow through discontinuities as a function of pressure gradient using an imaginary surface between the two surfaces of the element.In this case,the flow along both of the length and width of the faults can be calculated. Also, it is possible to consider the displacement of discontinuity in hydro-mechanical (HM) coupling simulation.Several authors have used double-node zero-thickness elements(e.g. Guiducci et al., 2002; Segura and Carol, 2008a; Alonso et al.,2013; Barani and Khoei, 2014). For the triple-node elements, a surface is added between the two main surfaces of the element.In this case, the pressure feild through the faults is discretized by additional nodes (Jha and Juanes, 2014).

One of the factors in modeling faults by a zero-thickness element is the contact constraint law. A normal contact constraint ensures that the interface surfaces are not overlapping each other.In this regard,methods such as penalty,Lagrangian multipliers and augmented Lagrangian are used (Wriggers and Zavarise, 2004;Cerfontaine et al., 2015b). Because of the advantage of simplicity in application, the penalty method is used herein for contact constraint and controlling the overlapping of the nodes.Moreover,in comparison with Lagrangian methods, its speed is higher and additional variables are not needed for the equations.However,use of this method is associated with the risk of ill-conditioning in the elemental stiffness matrix. It should be noted that plastic models,which generally include dilatation, softening, etc., are used to simulate small-scale joints.Since the fault scale is several meters to several kilometers, elastic models instead of plastic models are used.Moreover,nonlinearity problems are solved by the Newton-Raphson method.The trapezoidal method is used to discretize the governing equations in time domain. For further information,studies of Gens et al.(1990)and Desai and Ma(1992)are suggested.

The focus of this paper is providing a precise and direct THM formulation of triple-node zero-thickness interface element in order to model the behavior of faults in a deformable porous medium.The associated equations in this paper are based on the Biot theory(Biot, 1941), and the discretized variables are solid displacement,fluid pressure, and heat. To do so, first, the basic equations of equilibrium are presented. After that, by introducing triple-node zero-thickness element, the formulations of THM modeling are derived for both finite and interface elements,using discretization methods.At last,in order to verify the presented formulations and to show the power of the triple-node interface elements in modeling of fault behavior, numerical simulations of three index problems are presented.

Fig. 2. Zero-thickness interface element used for (a) 2D and (a) 3D models. η, ξ and ζ represent the local coordinate axes.

Fig. 3. Zero-thickness interface element (IE) in single-, double- and triple-node forms that connect the sides of its surfaces as node-to-node.

2. Governing equations

Equilibrium equations including linear momentum balance,mass conservation and energy balance equations for solid and water phases are presented. Each of the equations for continuous medium and discontinuous faults is discretized in next section.Mass conservation and fluid flow equations are based on the Lagrangian method, in which the solid mass is automatically conserved for a reference element volume (REV) (Charlier et al.,2001; Li and Laloui, 2016). The following assumptions were considered:

(1) The solid phase is a porous skeleton saturated with fluid;

(2) Solid ingredients are homogeneous based on the mechanical properties of rocks;

(3) Small strain and deformation are assumed;

(4) Equilibrium equations are expressed only for quasi-elastic behavior and the inertial effects are ignored;

(5) There are no chemical reactions between phases;

(6) The fluid is a Newtonian type and the fluid flow follows Darcy's law;

(7) The permeability tensor is diagonal; and

(8) The viscosity of each phase is assumed constant.

2.1. Linear momentum balance equation

According to Coussy (2004), one of the governing equations of coupled fluid flow and geomechanics in poromechanics is the momentum balance equation. Poro-elastic equations are written based on the terms of general stresses, strains, and pore pressure.Accordingly,the linear equilibrium equation for the whole system,neglecting the terms including the gradient of the fluid velocity,the relative acceleration of the fluid and the solid phase acceleration,is written as follows:

-?·σ = f (1)

where σ is the Coussy total stress vector and f is the gravity term equal to ρbg, in which g is the gravity vector, and ρbis the bulk density calculated as follows:

ρb= ρs(1-φ)+φρfSf(2)

where φ is the porosity; Sfis the fluid saturation; and subscripts f and s represent the fluid and solid phases, respectively.

For fluid phase,the linear momentum balance equation could be expressed as Darcy's law as follows:

where qfis the fluid flow rate,k is the intrinsic permeability tensor,and μfis the dynamic viscosity of fluid phase.

2.2. Mass conservation equation

2.2.1. Solid phase

The mass conservation equation of the solid phase is written as follows:

where t is the time.

Eq. (3) could be rewritten as follows:

2.4.2. Strain-deformation relation

The pressure variation causes the deformation of the reservoir.This deformation could be expressed by the displacement field.Assuming small strain, the displacement and strain relation is expressed as follows:

where εijand uijare the strain and displacement tensors,respectively.

For time variations of porosity, we have

2.4.3. Stress-strain relation

If the heat problem is considered, the equation derived from Hooke's law could be written as follows:

2.2.2. Fluid phase

The mass balance equation for different fluid reservoirs could be expressed as follows:

where ε is the total strain tensor,Deis the linear elastic tensor,and βsis the solid thermal expansion coefficient.It should be noted that Young's modulus and Poisson's ratio are used to calculate Defor continuous medium and elastic stiffness for discontinuous medium of faults.

2.4.4. Solid density relation

Lewis and Schrefler (1998) provided the following equation for change of the solid phase density over time:

2.3. Energy balance equation

The energy balance equation for thermal conductivity and thermal convection could be written as follows:

where (ρCp)avg= (1 - φ)ρsCsp+ φρfCfp, in which Cspand Cfpare the special heat capacities of solid and fluid phases, respectively; T is the temperature; qTis the heat flux of medium; and?Tis the divergence as the transpose of the gradient. In Eq. (8),the first, second and third terms represent the heat storage,convection heat transfer and conduction heat transfer,respectively.

2.4. Constitutive equations

In addition to the equations of equilibrium as described above,other complementary equations known as constitutive equations are necessary to solve the governing equations for coupled THM modeling. These relationships are outlined below.

2.4.1. Effective stress relation

According to the Terzaghi's theory,the effect of pore pressure on total stress and soil matrix behavior could be expressed as follows:

2.4.5. Fluid density relation

Assuming the compressibility of fluid within the porous medium, and its dependence on the temperature and pressure of the medium, the following relations could be considered:

where Kfis the bulk modulus of fluid, βfis the fluid thermal expansion coefficient, and cfis the isothermal compressibility of fluid.

Using Eqs. (13) and (14), one could obtain an exponential relation for fluid density:

where ρf0is the fluid density at reference temperature (Tf0) and reference pressure(pf0).

σ′= σ-αmpf(9)

where σ′is the effective stress vector; α is the Biot's coefficient calculated by α = 1- K/Ks, in which K and Ksare the total and solid phase bulk moduli,respectively;and m is a vector that causes fluid pressures to affect only on the principal stresses,and it is equal to [1， 1， 0]Tand [1， 1， 1， 0， 0， 0]Tin 2D and 3D problems,respectively.

2.4.6. Thermal conductivity relation (Fourier's law)

According to the empirical version of Fourier's law,the relation between heat flux (qT) and effective thermal conductivity (χeff)could be written as follows:

qT= -χeff?T (16)

Effective thermal conductivity could be estimated both theoretically and experimentally (Lewis and Schrefler,1998).

2.4.7. Summary of the equations governing the THM process

(1) Linear momentum

Based on Eqs. (1) and (2), we have

(2) Mass conservation

By substituting Eq.(12)into Eq.(5),the equation of permeability changes is obtained.By placing the permeability changes equation and also Eqs. (13)-(15) into Eq. (6), the total equation of conservation is obtained:

where M is the Biot's modulus which could be calculated by 1/M = φ/Kf+ (1 -φ)/Ks.It is obvious that,for the infinite values of the solid phase bulk modulus (an incompressible solid), the Biot's coefficient will be equal to 1.

(3) Energy balance

By substituting Eqs.(8)and(16)into Eq.(7),the energy balance equation for whole system is obtained:

Eqs. (17)-(19) are the fundamental equations governing the THM coupling process in deformable and saturated porous media.These equations, along with other constitutive equations in this paper, are discretized in the spatial domains (continuous and discontinuous)and time domain to investigate the geomechanicalfluid flow behavior of faults in a deformable porous medium. It should be noted that all of the above equations are based on the Biot theory (Biot, 1941), which is coupled with heat transfer equation.Further details on the above governing equations for THM models were presented in Lewis and Schrefler (1998).

2.5. Initial and boundary conditions

To solve the above equations, it is necessary to determine the initial and boundary conditions of the problem medium for each of displacement, fluid pressure, and temperature parameters.

2.5.1. Initial conditions in the domain and on the boundary

where Ω and Γ represent the domain and boundary of problem,respectively; and u0, pf0and T0are the displacement, fluid pressure, and temperature values at time t = 0, respectively.

2.5.2. Boundary conditions

In the following, both Dirichlet and Neumann boundary conditions are considered. The boundary conditions can be imposed either by values on Γπ, or fluxes on Γqπ, where the boundary Γ is equal to the union of Γπand Γqπ,where π could be displacement(u),fluid (f), or temperature (T).

(1) Displacement and stress

(2) Fluid phase

(3) Heat transfer

3. Discretization of governing equations

3.1. Continuous medium

As mentioned previously, FEM is used to discretize the equations expressed in continuous porous medium. By applying the method of weighted residuals on linear balance,mass conservation and energy balance equations and then using the standard Galerkin FEM in discretization of equations and taking into account the unknowns in the form of their nodal values by means of shape functions (pf=u =and T =whereu andT are the vectors of the nodal values of the unknowns;Np,Nuand NTrepresent the shape functions for pressure, displacement and temperature, respectively),the following equations are obtained:

where fu, fpand fTrepresent the internal traction vectors. The coefficient matrices and vectors corresponding to each of the above equations are described in Appendix A.

3.2. Faulty discontinuous medium

where σ′fn， τfs1and τfs2are the normal and shear stresses in local directions of ξ,η and ζ, respectively; ufn， ufs1and ufs2are the interface element displacements in local directions of ξ, η and ζ,respectively; and subscripts t and b represent the top and bottom surfaces, respectively.

The constitutive relationship between σ′fand ufis in nature an elastoplastic relationship. These models are generally used to simulate the properties of small-scale joints including dilatation,softening, etc. Since the fault scale is several meters to several kilometers,only elastic models are used and the local non-linearities of the normal and shear stiffnesses of fault discontinuity are neglected. Hence the elasticity matrix is defined as follows:

Since FEM is not able to model discontinuous media, the zerothickness interface element is used hereafter. In this section, the governing equations for THM problems are discretized and formulated for triple-node zero-thickness.

According to the geometry and characteristics of the interface element,the method of discretization of equations is different from that of continuous state. To better understand this, the geometry and shape of this element should be considered. In this element(Fig. 4), there are two parallel surfaces (up and down) that represent two sides of the element and the relation of the element with continuous adjacent finite elements. There is an intermediate surface which is parallel to the main surfaces of the element. This middle surface is used only for consideration of the fluid pressure gradient across the fault, and only the parameters of the fluid pressures and the heat are discretized over it, while on the other two surfaces,all the unknown parameters are discretized.In other words,at the upper and lower surfaces,the unknown matrix is X =

while in the middle surface of the element, the unknown matrix is X =

As shown in Fig.4b,if the displacement in three local directions ξ, η and ζ are represented by u, v and w, then the effective stress(σ′f)and displacement(uf)vectors in the interface element,in the 3D state, could be expressed as follows:

where knand ksare the normal and shear stiffnesses of fault discontinuity, respectively. It should be noted that, when the pore pressure in blocks located in the adjacent sides of the interface element is reduced, based on Terzaghi's theory (Eq. (9)), the effective stress(Eq.(31))will increase,resulting in the penetration of the two neighboring blocks at interface element nodal points.To deal with this problem, the normal and shear stiffnesses (in Eq.(33)) are used to enforce the non-overlapping constraint of the interface sides by penalty method.The large values of these elastic stiffnesses could be interpreted from mathematical point of view as penalty parameters which could ensure the accuracy of this method.However, too large values of these parameters may cause ill-conditioning of the global matrix.

Fluid flow through the faults could be considered as longitudinal and transversal fluxes. If the lower, middle and upper surfaces of the element are expressed in numbers 1, 2 and 3, respectively, according to the Darcy's law, the flow through the fault as in Eq. (8)could be expressed for the space between surfaces 1 and 2, as well as between surfaces 2 and 3 (i.e. whole element):

Fig. 4. Triple-node zero-thickness interface element.

where kfrepresents the longitudinal permeability of fault, and z is the middle surface elevation. The gradient is calculated in η and ζ directions. To calculate the permeability of the fault in its longitudinal direction, one can study fluid flow in the fractures.

Based on the study of Boussinesq (1878), a fracture was considered as two parallel smooth plates with a distance of Δa,which is known as fracture aperture.Total fluid flow is proven to be proportional to the cube of the distance between the two sides of the fracture, known as fracture aperture. Accordingly, the fracture permeability could be calculated by Eq. (35), which is known as cubic law:

The longitudinal permeability of fractures used for the whole interface element (i.e. space between all three surfaces) can be obtained from Eq.(35).Many authors have used this relationship in their studies in fracture modeling (Oron and Berkowitz, 1998;Segura and Carol, 2008a, b; Garolera et al., 2013; Cerfontaine et al., 2015b).

Based on Cappa and Rutqvist (2011), HM models of single fractures of rocks can be used to calculate fault permeability.In this regard, the longitudinal transmissibility of the fault can be described as follows:

The fault zone is generally divided into two parts of the damaged zone and the core (Cappa et al., 2009; Mitchell and Faulkner, 2009; Cappa and Rutqvist, 2011). In this case, if the difference between the flow characteristics of the fault damaged zone and its surrounding rocks is negligible, in the above relation, the parameter Δa will be the value of the fault core aperture. Otherwise,Δa should be determined to reflect the characteristics of the both zones. By substituting Eq. (36) into Eq. (34), the longitudinal flow through the fault could be expressed for the whole interface element:

where γfis the fluid specific weight.

In the coupled HM modeling,the value of the aperture depends on the calculated pressure and displacement in various time steps.If the variation is low, it can be calculated at the end of each time step; otherwise, the system of equations will be highly nonlinear and iteration methods should be used to solve them.

The transverse fluid flow in fault could be expressed by considering the pressure gradient in the ξ direction at different surfaces of the element as follows:

where pfiis the fluid pressure at surface i,and Tftis the transverse conductivity of faults between two adjacent surfaces.

The heat transfer equation also should be considered in both the longitudinal and transverse directions of the fault.In this case,the gradient in Eq.(16)is also calculated in longitudinal direction of the η and ζ directions and in transverse direction of the ξ direction.It is worth noting that the medium of the problem is considered to be homogeneous for the thermal conductivity coefficient. Finally, by considering the equations in this section and applying the Galerkin FEM to discretize them, the final form of these equations will be presented as the following matrix equations:

Fig. 5. (a) Geometry of Mandel's problem (adapted from Abousleiman et al.,1996); (b) Medium spatial discretization and boundary conditions for a quarter of Mandel's problem geometry. uh and uv represent the horizontal and vertical displacements, respectively.

Table1 Data used in non-isothermal Mandel's problem(adapted from Pao et al., 2001).

Each of the above-mentioned coefficient matrices is given in Appendix A. It should be noted that finally, after assembling the series of Eqs. (28)-(30) and Eqs. (40)-(42) in two continuous and discontinuous media,and by canceling the internal traction vectors fu, fpand fTin these equations, a series of equations will be obtained to solve the problems of coupled fluid flow-geomechanics in a faulted deformable saturated porous medium.

Fig. 7. Pore pressure (Pa) distribution after half a day.

4. Numerical simulation

The faults are significant components in most geological formations.The injection of fluid into the reservoir increases the pore pressure and thus reduces the effective stresses. Reducing the effective stress on the fault can cause it to be reactivated and the fault acts as a pathway for fluid flow.In general,two analytical and numerical methods are used to study the faults.

In this paper, Mandel's problem in non-isothermal condition is first solved to verify the presented model in a completely continuous porous medium. Two other index problems of coupled simulations of faults are solved using the methodology presented herein.

4.1. Non-isothermal Mandel's problem

Mandel's problem (Mandel,1953) is one of the basic problems associated with the Biot theory. Since the introduction of the analytical solution to this problem, their numerical methods and codes are verified (Abousleiman et al., 1996; Pao et al., 2001;Ferronato et al., 2010; Ranjbar et al., 2018), in association with its analytical solution (Ferronato et al., 2010). This problem in nonisothermal condition includes the effects of fluid flow, heat, and deformation simultaneously,which is selected subsequently as one of the index problems.

Fig. 6. Pore pressure (Pa) distribution in whole porous medium at different time steps.

Fig. 8. Distribution of temperature (°C) at different times: (a) Calculated in this study, and (b) Result of Pao et al. (2001).

In the Mandel's problem,an infinitely rectangular cross-section of a poroelastic saturated porous medium is placed between two hard-friction sheets and an instantaneous load is applied to it.Applying this load, along with evacuation of the fluid from the lateral sides of the porous medium, allows the pore pressure to increase in the center of the layer and at the initial discharge time,which is higher than its initial value.Fig.5a illustrates the boundary conditions of the problem,i.e.the lateral boundaries are free in pore pressure and stresses.An instantaneous load(with amount of 2F)is applied at the initial time and at the top and bottom of the medium.Afterwards, by evacuating the medium through the lateral boundaries, the pressure is allowed to decrease over time. Here, a non-uniform decrease of pressure in the center of the area will be observed. Cryer (1963) obtained similar results at the center of a spherical medium that was consolidated by hydrostatic pressure.This type of non-uniform pore pressure response is called Mandel Cryer effect.To consider the role of temperature,a thermal current(with amount of 2Q)is also considered at the initial time to the top and bottom boundaries.

Fig. 9. Geometry and boundary conditions of the second problem.

Table 2 Data used in HM modeling of a faulted porous medium.

Fig.10. Pore pressure (Pa) distribution at different time steps of the second problem.

Fig.11. Pore pressure changes along the (a) height and (b) length of the porous medium at different times.

Due to the symmetry of the problem environment,the Mandel's problem is solved only for a quarter of its environment. For this,new boundary conditions are assigned for the problem because of the change in its geometry(Fig.5b).The FEM is used for modeling and a total of 100 elements with 4 nodes were considered. The geometry of the environment is also square and a/b equals one.The data used in this problem are taken from Pao et al. (2001) and presented in Table 1.This makes it possible to compare the results.In Table 1, λ is the fluid mobility which equals the ratio of permeability to fluid viscosity; and Csand Cfare the heat capacities of solid and fluid, respectively. Due to the large number of results in Pao et al.(2001),only its results of the temperature distribution are presented herein.

Fig.12. Vertical displacements at different times.

In Fig. 6, the distribution of fluid pressure in the entire environment of the problem is shown at different times.In this figure,at initial time, the water pressure at the midpoint (point A) is calculated equal to 50,000 Pa,which indicates the initial pressure value based on the Skempton effect(Skempton,1954).After a very short time, the wave of increasing pressure reaches this point and increases the fluid pressure. According to this figure, the pore pressure in the initial short intervals will be greater than its initial value.This increase in pressure at the initial time reflects the Mandel effect that is well-defined in the figure.This can be due to the uneven drop in pore pressure over the initial period of time,which causes tensile stress to be applied near the lateral boundaries. Afterward,the concentration of stress in the center of the area with increasing load on the fluid in the center of the environment will cause increase in pore pressure.

Fig.7 shows the distribution of fluid pressure after a long(half a day)time period.In this figure,it can be seen that the distribution of pore pressure changes from the center to the lateral boundaries is maintained, even in very small amounts of total pressure.

Also, the distribution of temperature in the problem environment at different times is shown in Fig.8.The results are consistent with the results presented in Pao et al. (2001). According to this figure, it should be noted that the effects of the temperature are observed in more than 0.01 d. Therefore, at the initial time and below this, the numerical results of the Mandel's problem in isothermal and non-isothermal conditions in the pressure diagram are approximately equal.

4.2. HM modeling of a saturated compressible porous medium with vertical fault

The HM behavior of a permeable fault in the continuous saturated porous medium is investigated. By applying vertical stress and creating a pressure difference between the upper and lower boundaries of the porous medium, the fluid flows to the upper boundary.Since the fault acts as a path for fluid transfer due to its higher permeability and porosity relative to the surrounding medium,the drainage of the fluid often occurs through the fault.This problem has been solved in a domain with dimensions of 10 m×10 m,in which a vertical fault is located in the center of the medium. The domain is meshed by 4-node and 9-node isoparametric elements for pressure and displacement unknowns,respectively. Additionally, the fault is meshed by triple-node zerothickness interface elements. Eventually, the porous medium is meshed by 110 elements and 132 nodes. Fig. 9 shows the element mesh of the problem and the boundary conditions. The required input data are given in Table 2.

Fig.13. (a)Geometry of the third problem(adapted from Kolditz et al.,2016)and boundary conditions;and(b)Medium spatial discretization.ux and uy represent the horizontal and vertical displacements, respectively.

The simulation was performed at time intervals of 0-20,000 s with intervals of 1 s,10 s and 100 s.Fig.10 shows the pore pressure distribution profiles at different times.As can be seen,the fluid flow within the fault has been done at the very early stages of the problem. The values of the pore pressure along the height (at the left boundary) and length (at height of 0.9 m) of the porous medium and the vertical displacements at different times are shown in Figs.11 and 12, respectively.

4.3. THM modeling of a faulted aquifer

For this problem, THM modeling of a faulted aquifer is presented. This model is one of the benchmark problems in coupling simulations presented in Kolditz et al.(2016).The geometry of the model in 2D mode is shown in Fig.13.The aquifer(with lithology of sandstone) is located between the sedimentary layers and all geological units have been discontinued and displaced by a reverse fault. The model is square with side length equal to 900 m. The thickness of the aquifer is 100 m and the fault has a dip of 80°and a width of 20 m.In this study,the zero-thickness element is used to discretize the discontinuous fault. It should be noted that this problem was studied in Kolditz et al. (2016). In that study, the continuous part of problem is discretized by triangular and quadrilateral elements and the fault discontinuity by finite-thickness elements, which are different with our technique. Therefore, the purpose of this problem is to compare the results of different methods and also to demonstrate the capability of the methodology presented in this paper in such simulations.

Fault zones are evaluated in different ways depending on their considered scale and type of problem. According to the type of problem,a fault zone could be considered as a single feature when a geomechanical model is of a scale of kilometers (Bohloli et al.,2015). In the present problem, the fault is of kilometer scale and the ratio of the fault thickness (20 m) to the width of the field(900 m) is 0.022, thus the fault can be modeled as a single discontinuity using a zero-thickness element without considering its thickness(20 m here).It is worth noting that the only parameter expressed in terms of the geometry of the fault in the mentioned equations is the size of the fault aperture, which is linked to the deformation and displacement of the porous medium through the elastic matrix (Eq. (33)). The values of Young's modulus and Poisson's ratio are presented in Table 3. Equivalent continuum models can be used to calculate the fault shear and normal stiffnesses(Eq.(33)) from these values. Further information about these kinds of models can be found in Priest (1993). Also, the value of the fault aperture in Eq. (36) plays an essential role in calculating fluid transferability of the fault.The permeability of the fault is given in Table 2, and Eq. (34) is used to calculate the fluid flow in thelongitudinal direction of the fault. This example shows that zerothickness element can be suitable for modeling kilometer-scale faults.

Table 3 Data used in THM modeling of a faulted aquifer (adapted from Kolditz et al., 2016).

Fig.14. Distribution of (a) temperature, (b) pore pressure, and (c) vertical displacement calculated in this study (i) and the results of Kolditz et al. (2016) (ii).

The domain for this problem is meshed using isoparametric elements and triple-node zero-thickness interface elements. This network contains a total of 2386 nodes and 2340 elements.Due to the high gradient of unknown changes such as pore pressure,temperature and displacement near the fault, the meshes in this area are using smaller sizes.Then, the code using Matlab software based on the equations and method described in this paper was developed.

The boundary and initial conditions of the problem, such as pressure, temperature and stress, as well as the characteristics of each of the boundaries, are shown in Fig.13. A static pressure distribution for initial conditions and a stress field equal to zero in the whole environment are considered. A pressure of 108Pa and the temperature of 30°C are set for the aquifer located in the hanging wall block of the fault. This causes the gradient difference, as a result of fluid flow from the sandstone formation into the fault.The upper and lower sides of the model have constant temperature values (see Fig. 13), such that the thermal gradient of 30°C/km exists in the environment of the problem. It should be noted that other features of the problem, including the hydraulic and mechanical properties of each of the geological events (Table 3), are similar to their values in Kolditz et al.(2016),which are taken from Jaeger et al. (2007). This makes it possible to compare the results from the methodology presented in this paper with the help of a zero-thickness element with the results of other researchers who have used other methods. In Fig.14, the distribution of temperature, pore pressure and displacement calculated in this study and the results of Kolditz et al. (2016) are compared. As shown in the figure, a very good agreement is found between the results obtained by fault modeling with a zero-thickness element and the results of other methods.

As shown in Fig.14a, due to the difference in temperature between the upper and lower boundaries of the porous medium, a uniform heat transfer from the heat source could be observed in the continuous part. In the fault zone, this heat uniformity has been changed and the thermal anomalies are visible along fault lines,due to the higher porosity and permeability of the fault in relation to its surroundings. The fault acts as a path for fluid flow, and the fluid in the aquifer located in hanging wall block moves through the fault due to the pressure difference. In addition, pressure drop occurs at a lower depth compared to the footwall block (Fig. 14b).Therefore, the effective stress in this hanging wall block will be higher,the amount of displacement in the vertical direction will be affected, and its value will be higher in the hanging wall block(Fig.14c).

Compared with continuous elements, zero-thickness elements have a lower dimension, and faults can be modeled by onedimensional (1D) line and 2D surface area elements, which are suitable for 2D and 3D problems,respectively.Such consideration is the fundamental concept of the presented model and it is the main advantage of this zero-thickness element in comparison with other methods.

5. Conclusions

In this paper, the THM modeling of a faulted deformable saturated porous medium was studied using an FEM in combination with the interface element. Initially, the governing equilibrium equations including momentum balance, mass conservation, energy balance equations and also constitutive equations were introduced. After that, due to the existence of a fault, these equations were discretized in two continuous and discontinuous parts.The continuous part discretization is performed using the standard Galerkin FEM and the discontinuous fault discretized using triplenode zero-thickness elements. Finally, the equations are solved in a fully coupled process.

We also tried to illustrate the ability of the interface element to solve the THM problems related to the faults by providing index problems. For this, three problems were selected for verification.The numerical examples presented in this paper show that the proposed method is robust in fault modeling.

Declaration of Competing Interest

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

List of symbols

α Biot's coefficient

M Biot's modulus

ρbBulk density

KfBulk modulus of fluid

KsBulk modulus of solid constituent

KbBulk modulus of solid skeleton

αcConvection heat transfer coefficient

u Displacement

? Divergence operator

μfDynamic viscosity of fluid phase

σ′Effective stress

χeffEffective thermal conductivity

ρfFluid density

qfFluid flow

pfFluid pressure

SfFluid saturation

βfFluid thermal expansion coefficient

? Gradient

g Gravity vector: [0, -9.8]T

qTHeat flux

K Intrinsic permeability tensor

cfIsothermal compressibility of fluid

DeLinear elasticity matrix

DfeLinear elasticity matrix for interface elements

qfflLongitudinal fluid flow through fractures

knNormal stiffness of fractures

kfPermeability of fault

kfracPermeability of fracture

φ Porosity

ksShear stiffness of fractures

ρsSolid density

βsSolid thermal expansion coefficient

CfpSpecial heat capacity of fluid phase

CfsSpecial heat capacity of solid phase

T Temperature

βsfThermal expansion coefficient of media

ε Total strain tensor

σ Total stress vector

qfftTransversal fluid flow across fracture

Appendix A

A series of final discretized equations relating to THM problems have coefficient matrices as follows:

(1) Continuous medium