A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems

T A fully coupled thermal-hydraulic-mechanical (THM) finite element model is presented for fractured geothermal reservoirs. Fractures are modelled as surface discontinuities within a three-dimensional matrix. Non-isothermal flow through the rock matrix and fractures are defined and coupled to a mechanical deformation model. A ro- bust contact model is utilised to resolve the contact tractions between opposing fracture surfaces under THM loadings. A numerical model has been developed using the standard Galerkin method. Quadratic tetrahedral and triangular elements are used for spatial discretisation. The model has been validated against several analytical solutions, and applied to study the effects of the deformable fractures on the injection of cold water in fractured geothermal systems. Results show that the creation of flow channelling due to the thermal volumetric contraction of the rock ma- trix is very likely. The fluid exchanges heat with the rock matrix, which results in cooling down of the matrix, and subsequent volumetric deformation. The cooling down of the rock matrix around a fracture reduces the contact stress on the fracture surfaces, and increases the fracture aperture. Stress redistribution reduces the aperture, as the area with lower contact stress on the fracture expands. Stress redistribution reduces the likelihood of fracture propagation under pure opening mode, while the expansion of the area with lower contact stress may increase the likelihood of shear


Introduction
Energy extraction from geothermal reservoirs involves multiple physical processes including thermal (T), hydro (H), and mechanical (M) processes that together influence the heat extraction from fractured geothermal systems (Tsang, 1991;MIT, 2006). Due to the complexity of this problem, and the number of parameters involved, modelling of these systems is viable primarily through numerical methods (McDermott et al., 2006). In a geothermal system, cold fluid is injected into an injection well, and hot fluid is extracted from the production well (e.g., Crooijmans et al., 2016). In order to understand the coupled processes and their effects, a robust numerical model that simultaneously solves all the governing equations in a coupled manner is essential for the successful investigation of a fractured geothermal system.
Fractures, natural or man-made, enhance flow within geothermal reservoirs. For instance, fractures dominate the flow in low permeability hot dry rocks (HDR) in the subsurface. Fractures may also con tribute to the creation of short-circuits between injector and producer wells, hence reducing the efficiency of a geothermal system (Emmermann and Lauterjung, 1997). In enhanced geothermal systems (EGS), due to the low permeability of the host rock, artificial fractures are induced, prior to injection of cold fluid, in order to enhance the effective permeability of the hot rock. In EGS, the stimulation can occur through induced slip on pre-existing fractures (shear stimulation), by creating new fractures using hydraulic fracturing technique (opening mode), or by a combination of the two (McClure and Horne, 2014). Thermally-induced fracturing has also been frequently observed in many subsurface applications, where a relatively cold fluid has been injected into a reservoir: for instance, in water injection wells in the petroleum industry (Bellarby, 2009), in geothermal wells (Benson et al., 1987;Tulinius et al., 2000), and even in relatively soft, unconsolidated formations (Santarelli et al., 2008). The volumetric flow rate in a fracture is proportional to the pressure gradient and the cube of the fracture aperture, i.e., the cubic law, which is derived from the general Navier-Stokes equation for flow of a fluid between two parallel plates (Zimmerman and Bodvarsson, 1996). Thus, variation in fracture aperture due to the changes in the normal and/or shear stresses acting on the fracture surfaces as a result of the THM processes strongly affects the fluid flow and heat transport in the fracture (Rutqvist et al., 2005).
Heat conduction between the fluid inside the fracture and the surrounding rock matrix has been of particular interest in many situations, including magma-driven fractures (Spence and Turcotte, 1985), hydraulic fracturing of wells (Wang and Papamichos, 1999), and hydraulic fracturing of shale gas reservoirs (Tran et al., 2013;Enayatpour and Patzek, 2013;Salimzadeh et al., 2016). Rock temperature at the surfaces of the hydraulic fracture is often considered constant, and equal to the temperature of the injected fluid (for example in Tran et al., 2013;Abousleiman et al., 2014). However, such an assumption does not satisfy conservation of energy, and does not account for the fact that heat exchange between the fracturing fluid and the rock gradually causes the fracturing fluid to thermally equilibrate with the matrix rock. Consequently, an unrealistically large effect due to thermal non-equilibrium is predicted by such approaches (Salimzadeh et al., 2016). Considerable efforts have been expended in developing THM models for geothermal reservoirs over the past several decades; however, very few studies have taken into account the evolution of fracture permeability under thermoporoelastic effects. McDermott et al. (2006) investigated the influence of THM coupling on the heat extraction from reservoir in crystalline rocks using an experimentally validated geomechanical model. Ghassemi et al. (2008), using a partially coupled formulation, derived analytic solutions for calculating fracture aperture changes induced by thermoelastic and poroelastic stresses during cold-water injection in an enhanced geothermal system (EGS). Ghassemi and Zhou (2011) proposed an approach to couple fracture flow and heat transport to thermoporoelastic deformation of the rock matrix using the displacement discontinuity (DD) method in which coupling is realised sequentially. Sequential coupling, in a non-linear system, suffers convergence problems, and requires more iteration and manual interference to converge. Abu Aisha et al. (2016) investigated the effects of the new fractures created during a geothermal lifetime on the overall permeability tensor of the fractured medium. Pandeya et al. (2017) proposed a coupled THM model for the variation of fracture aperture during heat extraction from a geothermal reservoir. They treated a fracture as a thin permeable layer in the matrix, with a stress-dependant fracture stiffness and elastic modulus. Guo et al. (2016) investigated the effect of the heterogeneity in the initial aperture distribution on the flow path within a single fracture in an EGS. The equivalent permeability in fractured reservoirs can be significantly affected by the choice of the aperture distribution model (Bisdom et al., 2016).
In the present study, a finite element model is presented in which fractures are treated more accurately in terms of their representation in the mesh, as well as in their physical behaviour under THM loading.
Fractures are modelled as 2D surface discontinuities in the 3D rock matrix. Separate but coupled flow/heat models are defined for the fracture and the rock matrix. The flow through the fractures is governed by the cubic law, and is coupled to the Darcy flow in rock matrix using leakoff mass exchange that is computed as a function of the fracture and matrix fluid pressures, and the matrix permeability. Local thermal non-equilibrium is considered between fluid in the fracture and fluid in the rock matrix. Advective-diffusive heat transfer is assumed in both the fractures and rock matrix. Heat transfer between fracture and matrix is allowed by conduction through the fracture walls, as well as by advection through the leakoff flow. Contact stresses on the fracture surfaces are computed using a robust contact model. Thermal and hydraulic loadings are considered in computing the contact stresses. The contact model is iteratively coupled to the THM model. The governing equations are solved numerically using the finite element approach. The coupled model has been validated against several available solutions, and applied to investigate the effects of fracture aperture alteration due to THM processes on the flow of the cold fluid in geothermal reservoirs.

Computational model
The fully coupled computational model is constructed from five separate yet interacting sub-models: a thermoelastic deformation model, two flow models (one for the fractures and one for the rock matrix), and two heat transfer models, for fracture and rock matrix, respectively. Single-phase flow is assumed within both the fractures and the rock matrix. In the thermoelastic mechanical model, the flow and the heat transfer through the rock matrix are constructed for three-dimensional matrix body, while flow and heat transfer models through the fractures are defined for two-dimensional discrete fractures, as schematically shown in Fig. 1. Fracture flow and solid deformation are two-way coupled through hydraulic loading exerted on the fracture surfaces, as well as by ensuring the compatibility of fracture volumetric strains. Heat transfer in the rock matrix and fractures is also coupled through a heat exchange term included in the fracture and matrix energy balance equations. A displacement vector (three components), fluid pressures (two components), fracture fluid and matrix temperatures (two components) are defined as primary variables. Tension is reckoned positive for stresses in the governing equations.

Thermoporoelastic mechanical model
The thermoporoelastic mechanical model is based on the condition of stress equilibrium for a representative elementary volume of the porous medium. The assumption of elastic behaviour for matrix deformation is reasonable for most thermally-induced rock deformations  Salimzadeh et al. Geothermics xxx (2017) xxx-xxx (Rutqvist et al., 2005). For quasi-static conditions, the linear momentum balance equation for this elementary volume may be written as (1) where F is the body force per unit volume, and σ is the total stress. Effective stress is defined as the function of total stress and matrix pressure that controls the mechanical effects of a change in stress. It is defined exclusively within the rock matrix, linking a change in stress to the change in strain. The effective stress for the rock matrix saturated with a single-phase fluid is defined as (Biot, 1941) (2) where σ′ is the effective stress, α is the Biot coefficient, p m is the fluid pressure in the rock matrix i.e. matrix pressure, and I is the second-order identity tensor. The Biot coefficient is defined as (3) where K and K s are the bulk moduli of a rock matrix and rock matrix material (e.g., mineral grains), respectively (Zimmerman, 2000). Assuming that the rock matrix shown in Fig. 1 undergoes a temperature change from initial temperature T 0 to a new value T m , the thermal strain in the solid rock, under the assumption of linearity, are given by (Zimmerman, 2000) (4) where α s is a symmetric second-order tensor known as the thermal expansivity tensor of the rock matrix. If the rock is isotropic, then where the scalar coefficient β s is known as the coefficient of volumetric thermal expansion of rock matrix. Note that due to the relatively slow movement of the fluid inside the rock matrix, local thermal equilibrium between the rock solid and the fluid in the matrix pores is assumed. The stress-strain relationship for thermoelasticity can be written as (Khalili and Selvadurai, 2003) (6) in which D is the drained stiffness matrix. Assuming infinitesimal deformations, strain is related to displacement by (7) where u denotes the displacement vector of the rock solid. Fracture surfaces are not traction-free in the present model, and hydraulic loading, as well as the tractions due to the contact between fracture surfaces, are applied on the fracture walls, as shown in Fig. 1. Assuming negligible shear tractions exerted from the fluid on the fracture walls, the fluid pressure is applied only in the normal direction on the fracture wall. The tractions on the fracture boundary Γ c are where σ c is the contact tractions on the fracture surfaces, p f is the fracture pressure, and n c is the outward unit normal to the fracture surface (on both sides of the fracture). Integrating Eq. (1) over the domain, and after some manipulation, the differential equation describing the deformation field for a saturated rock matrix is given by (9)

Fracture flow model
A separate flow model is considered for fractures. This model allows direct computation of the fluid pressures inside the fracture, and explicit application of hydraulic pressures on sub-dimensional fracture walls (see Fig. 1). The objective is to obtain a more realistic representation of fracture flow. Assuming a high aspect ratio fracture Γ c that has a lateral extent that is much larger than its aperture, the average velocity of a fluid along the fracture surface can be approximated using the cubic law as (Zimmerman and Bodvarsson, 1996) (10) where a f is the fracture aperture, defined as the differential normal displacement between two walls of the fracture, a f = (u + − u − ).n c + a f c , μ f is the fluid viscosity. u + and u − represent displacements on two sides of the fracture, and a f c is the fracture aperture at contact. When two surfaces of a fracture are in contact, the displacement field on two surfaces of the fracture would be identical i.e. u + = u − , and the fracture aperture at contact (a f c ) is a function of the contact tractions as explained in section 2.7. The mass balance equation for a slightly compressible fluid may hence be written as (Salimzadeh andKhalili, 2015, 2016) (11) in which ρ f is the fluid density, and L f is the leakoff flow from the fracture to the matrix. This leakoff leads to mass transfer coupling between the fracture flow and rock matrix flow. Assuming that the fracture fluid is Newtonian, the leakoff flow per unit area of the fracture wall can be written, using Darcy's law, as (Salimzadeh et al., 2017a) (12) where k n is the intrinsic permeability of the rock matrix in the direction normal to the fracture (in the direction of n c ), and represents the pressure gradient along n c . In case of a fault zone, some average of the fault zone permeability and the matrix rock permeability can be used (Norbeck et al., 2016). Considering a barotropic fluid in which the fluid density is a function of fluid pressure and temperature, the change in density may be written as (13) where c f and β f are the compressibility and volumetric thermal expansion of fluid, respectively. Combining Eqs. (10)-(13), and after some manipulation, one obtains the governing equation for laminar flow through the fracture under non-isothermal conditions as Salimzadeh et al. Geothermics xxx (2017) xxx-xxx The term provides explicit coupling between the displacement field and the fracture flow field, which is symmetric to the fracture pressure loading term, p f n c .

Matrix flow model
The flow through the porous matrix, i.e., matrix flow, is constructed by combining Darcy's law with mass conservation for the fluid. Neglecting inertial effects, Darcy's law describing matrix hydraulic diffusion under hydraulic gradient may be written as (15) where v r is the relative velocity vector of the matrix fluid, k m is the intrinsic permeability tensor of the rock matrix, and g is the vector of gravitational acceleration . The mass balance equation for the fluid in the rock matrix may be written as (16) where ϕ is the rock matrix porosity, and v m is the fluid velocity in the matrix. δ(x − x c ) is the Dirac delta function, where x c represents the position of the fracture (Γ c ). Note that the leakoff only occurs on the boundary of the volume element that is connected to a fracture (Γ c ). Considering a barotropic fluid, the change in density may be written as (17) where c f and β f are coefficients of the fluid compressibility and volumetric thermal expansion, respectively. Integrating over the element and after some manipulation, the governing equation for the flow model may be expressed as (18) The Biot coefficient α appears in Eqs. (9) and (18), whereas it does not appear in the fracture flow model (Eq. (14)), as the fracture itself is not a "porous medium". The Biot coefficient couples the flow in matrix with the mechanical deformation, and setting α = 0 will decouple the mechanical deformation model and the matrix flow model, in which case mechanical loading will have no direct effect on the matrix pressure, and vice versa. In contrast, fracture pressure will always be coupled to the mechanical deformation model, irrespective of the value of the Biot coefficient.

Matrix heat transfer model
The governing equation for heat transfer through the rock matrix can be obtained by combining Fourier's law with an energy balance for saturated rock. It is assumed that the fluid velocity in the rock matrix is slow enough such that the solid grains and the fluid in the rock matrix are always in local thermal equilibrium. Convective, i.e., conduction and advection, heat transfer in rock matrix can be written as (19) where q mc is the heat flux through the rock matrix, λ m is the average thermal conductivity tensor of the matrix, T m is the matrix temperature, C f is the fluid specific heat capacity, and v m is the fluid velocity. The average thermal conductivity tensor of the matrix is approximated as follows, from the thermal conductivity tensors of rock solid (λ s ) and fluid (λ f ) as (see Zimmerman (1989) for more accurate models of the effective thermal conductivity) The heat energy change due to thermal power in the course of the bulk deformation of matrix and fluid can be expressed, respectively, as (21)   (22) Heat is also exchanged between matrix and fracture fluid by conduction through the fracture surfaces, and by advection through the leakoff mass exchange term, as (23) where λ n is the average thermal conductivity of the rock matrix along the direction normal to the fracture (in the direction of n c ), and represents the temperature gradient along n c . The heat storage in the matrix saturated with a fluid is given by (24) where ρ m C m can be computed (exactly) from the density and specific heat capacity values of rock solid (ρ s , C s ) and fluid (ρ f , C f ) as (25) Combining the above-mentioned equations, and after integrating them over the matrix and fracture domains, the governing equation for heat transfer through the matrix can be written as (26)

Fracture heat transfer model
Using a similar approach, the governing equation for heat transfer through the fluid in the fracture can be obtained by combining Fourier's law with an energy balance for the fluid. The advective heat transfer through the fluid in the fracture can be written as Salimzadeh et al. Geothermics xxx (2017) xxx-xxx (27) and the final form of the heat transfer equation in the fracture can be written as (28) 2.6. Finite element approximation Governing equations are solved numerically using the finite element method. The Galerkin method and finite difference techniques are used for spatial and temporal discretisation, respectively. The displacement vector u, fluid pressures p m and p f , and matrix and fracture fluid temperatures T m and T f are taken as the primary variables. Using the standard Galerkin method, the primary variable within an element is approximated from its nodal values as (29) where N is the vector of shape functions and is the vector of nodal values. Using the finite difference technique, the time derivative of X is defined as (30) where X t+dt and X t are the values of X at time t + dt and t, respectively. The set of discretised equations can be written in matrix form of SX = F, in which S is the element's general stiffness matrix as (31) and F is the vector of right-hand-side loadings (32) where ( Salimzadeh et al. Geothermics xxx (2017) where K is the mechanical stiffness matrix, the C matrices are the coupling matrices, and the H matrices represent the conductance and advection matrices. Matrices M are the flow-heat capacitance matrices. Matrices L are the leakoff flow and heat matrices. Vector F is the applied load vector, vectors Q represent the fluid and heat flux vectors, , and are the vectors of nodal values of displacement, fluid pressure, and temperatures, respectively.
, [B 2 ] 1×3n = δ T B 1 , and, [B 3 ] 3×n = ?N are derivatives of the shape function. is the gradient matrix, , and ∇ is the gradient vector. Superscript t represents the time at the current time step; superscript t + dt represents time at the next time step, and dt is the timestep. The non-diagonal components of the stiffness matrix S are populated with the coupling matrices C, and L. Note that the leakoff term (flow and heat) only exists for matrix elements (volume elements) connected to a fracture; and the integration is performed over each side of the fracture separately. The gradient matrix for three-dimensional displacement field is defined as (58) The components of the stiffness matrix are dependent upon the primary unknown variables, i.e., conductance, capacitance and coupling coefficients of the fracture are all dependent on the fracture aperture; therefore, a Picard iteration procedure is adopted to reach the correct solution within acceptable tolerance. For the current iteration, s + 1 in the current step, n + 1, the solution-dependent coefficient matrices in the stiffness matrix are updated using weighted average solution vector defined as where and are the solution vectors of the two most recent iterations in the current timestep n + 1, and θ = 2/3 is the weighing coefficient. For the first iteration s = 1, the previous timestep solution is used as (60) where is the solution vector from previous timestep n. The iterations are repeated until consecutive normalised values of agree to within a specified tolerance ε The tolerance is set to 0.01 in this work. The discretised coupled equations are implemented as part of the Imperial College Geomechanics toolkit (Paluszny and Zimmerman, 2011), which interacts with an octree volumetric mesher and the Complex Systems Modelling Platform (CSMP++, also known as CSP), an object-oriented application programme interface (API), for the simulation of complex geological processes and their interactions (formerly CSP, Matthäi et al., 2001). Quadratic unstructured elements are used for spatial discretisation of surfaces (quadratic triangles) and volumes (quadratic tetrahedra). The triangles on two opposite surfaces of a fracture are matched with each other, but they don't share nodes, and duplicate nodes are defined for two sides of a fracture. Therefore, there are two matrix degrees of freedom and one fracture degree of freedom as the governing equations for the fracture (flow and heat) are solved only on one side of the fracture. The triangles are matched with faces of the tetrahedra connected to the fractures. Fracture flow and heat equations are solved only on one side of the fracture, whereas, the coupling matrices C and L are accumulated on both sides of the fracture. Matrix deformation, flow and heat equations are accumulated over the volume elements. The ensuing set of linear algebraic equations SX = F is solved at each iteration using the algebraic multigrid method for systems, SAMG (Stüben, 2001).

Contact model
In the present study, fractures are modelled as surface discontinuities within a three-dimensional matrix; therefore, the contact problem arises and the contact stresses (normal and shear) are required to be computed in order to avoid the inter-penetration of the fracture surfaces under compressive loading. The Augmented Lagrangian (AL) method has been successful for enforcing the contact constraint accurately when computing high contact precisions, by combining the Lagrange multiplier and penalty methods to exploit the merits of both approaches (Wriggers and Zavarise, 1993;Puso and Laursen, 2004). A sophisticated algorithm is used for the treatment of frictional contact between the fracture surfaces, based on isoparametric integration-point-to-integration-point discretisation of the contact contribution. Contact constraints are enforced by using a gap-based AL method developed specifically for fractured media (Nejati et al., 2016). In this model, penalties are defined at each timestep as a function of local aperture, so that they are larger away from the fracture tips, and reduce to zero at the tips. In the contact model, the equilibrium equation has been satisfied in which the hydraulic and thermal contributions are applied as Salimzadeh et al. Geothermics xxx (2017) xxx-xxx boundary values on the right-hand-side (62) Pressures and temperatures are imported from the THM model. When two sides of a fracture are in contact, the change in the aperture of the fracture is defined as a linear function of the change in normal contact traction as (63) where a f c is the fracture aperture at contact, K n is the fracture stiffness, and σ n is the contact stress (compressive) normal to the fracture. However, nonlinear fracture stiffness models (Bandis et al., 1983;Barton et al., 1986) can also be used. The contact and THM models are coupled iteratively, such that in each timestep, first the THM model is run with the contact stresses computed from the previous step. Then the computed pressures and temperatures from the THM model are passed to the contact model, and the contact stresses are updated. Finally, the THM model is run again with updated contact stresses. The thermal and hydraulic loadings are applied as body forces to the right-hand-side of the contact model while contact stresses are applied as boundary tractions to the right-hand-side of the THM model.

Simulation of geothermal systems
Three sets of examples of geothermal systems are selected for simulation in this section. The first example is used for validating the heat transfer module of the present model, as well as delineating the extent of the validity of current semi-analytical solutions for the case of a permeable matrix. Further validation examples for the present numerical model can be found in Salimzadeh et al. (2016), Salimzadeh et al. (2017aSalimzadeh et al. ( , 2017bSalimzadeh et al. ( , 2017c, and Usui et al. (2017). The second and third examples demonstrate the effects of variation of the contact tractions and fracture aperture due to thermoporoelastic deformation of the matrix on fluid flow within a fracture in an EGS, and within a fractured geothermal reservoir, respectively. The effect of gravity has been ignored in these simulations by assuming that the fluid flow dominantly occurs in the horizontal direction. Bodvarsson (1969) derived an analytical solution for the problem of advective-diffusive heat transfer through a single one-dimensional fracture, while the heat is transferred through the matrix only by diffusion in the direction normal to the fracture (1D diffusion). Ghassemi et al. (2008) also proposed semi-analytical solution for a similar problem with leakoff of the fluid into the matrix. However, in their solution, it is arbitrarily assumed that the leakoff rate is constant (a fixed ratio of the injection rate Q L = mQ) and does not vary in time. Such an assumption may not be realistic, as is shown by the present simulations. By setting a very small value to m (for example m = 0.01), their solution approaches that of Bodvarsson (1969).

Rigid fracture in permeable matrix
In this section, a fracture of length 100 m is considered between injection and production wells. Plane-strain conditions are assumed, in order to validate the model results against the above-mentioned solutions. Injection of cold water at temperature 20°C, at constant rate Q = 0.0001 m 3 /s is assumed, while production is simulated through constant zero pressure in the production well. The initial temperature of the rock matrix is set to 100°C. Water has a density of ρ f = 1000 kg/m 3 , a heat capacity of C f = 4200 J/kg°C, and the matrix rock has density ρ s = 2820 kg/m 3 , heat capacity of C s = 1170 J/kg°C, and thermal conductivity of λ s = 2.88 W/m°C. The material properties are summarised in Table 1 (example 1). Constant matrix pressure and temperature is assumed at the far boundaries of the simulation region.
Several cases are simulated, in which the permeability of the rock matrix is increased from zero to 1 × 10 −12 m 2 . The results for the fluid temperature at production, as well as spatial distribution of the fluid temperature along the fracture are shown in Figs. 2 and 3. Included in these figures are the solutions proposed by Bodvarsson (1969) and Ghassemi et al. (2008) for comparison. The temperature of the cold water in the fracture increases as it exchanges heat with the hot rock matrix. For the case with impermeable matrix, the temperature of the produced fluid drops more rapidly than for the cases with a permeable matrix. This is because the fluid has higher velocity in the impermeable case, and so the cold water reaches the production well more rapidly. Note that the production is defined by constant pressure at the producer well; thus, the volume of the produced water is variable in time. Very good agreement is found between the present model results for the impermeable case and the solution proposed by Bodvarsson (1969), and also with the solution by Ghassemi et al. (2008) for the case of low leakoff ratio (m = 0.01). In the permeable cases, the permeability of the matrix is considered only in the direction normal to the fracture (to create one-dimensional leakoff), except one case that is described by the text "[2D]" in Figs. 2 and 3. When the permeability of the matrix increases, the leakoff increases, so the fluid velocity in the fracture de  Salimzadeh et al. Geothermics xxx (2017) xxx-xxx creases. Therefore, the residence time of the fluid, i.e., the time that the injected fluid spends inside the fracture prior to reaching the production well, increases, and hot fluid is produced for an extended period of time, as shown in Figs. 2 and 3. Leakoff also increases with the dimension of the flow in matrix (Salimzadeh et al., 2017a), so for the case with matrix permeability k m = 1 × 10 −13 m 2 , the case with 2D leakoff predicts a longer period of hot fluid production compared with the 1D leakoff simulation. The solution by Ghassemi et al. (2008) for the fluid temperature at the producer is computed for different values of leakoff ratio m = 0.01, 0.50, 0.75 and 0.99, and plotted in Figs. 2 and 3. In their solution, higher values of m represent a higher amount of leakoff, so the produced water has higher temperature for an extended period of time. However, as time elapses a sharp reduction in the temperature of the produced water is observed such that the case with very high leakoff ratio m = 0.99 produces colder water in the producer at a later time.
In Fig. 3, the spatial distribution of the temperature of the fluid inside the fracture is shown at time t = 10 8 s. Again, very good agreement is found between the results of the present study and the solutions given by Bodvarsson (1969), and also with the solution by Ghassemi et al. (2008) for an impermeable matrix. As the leakoff increases, either due to an increase in the permeability of the matrix, or an increase in the dimension of the flow field within the matrix, the fluid velocity in the fracture reduces. Slower flow in the fracture increases the fluid residence time, and therefore results in higher heat ex change with the hot matrix, and so a higher fluid temperature is observed in the producer for an extended period of time. Again, the present model results differ from the solution by Ghassemi et al. (2008). This is due to the fact that in the work of Ghassemi et al. (2008), it is arbitrarily assumed that the leakoff is equal to some fixed fraction of the injection rate, which does not vary in time, whereas in the present work the leakoff is computed as part of the coupled simulation, and the ratio of leakoff to injection is found to vary with time.

Deformable fracture in an impermeable matrix
A disk-shaped fracture of 200 m diameter is considered in the horizontal plane, with injection and production wells connected to the fracture at locations 50 m from the centre of the fracture to the left and right, respectively, as shown in Fig. 4. The injection rate is set to Q = 0.001 m 3 /s of water with temperature of 20°C, while the rock has an initial temperature of 80°C. Rock deformation is allowed in this example, and the elastic properties of the rock are set to Young's modulus E = 37.5 GPa, and Poisson's ratio ν = 0.25. Material properties of the fluid and rock are given in Table 1 (example 2). Production is defined by constant pressure at the producer. The in situ stress normal to the fracture plane is set to σ = 60 MPa, and initial fluid pressure is set to p i = 20 MPa. The fracture stiffness is set to 10 11 Pa/m, and the fracture aperture at zero contact stress is set to 0.6 mm. The viscosity of the water is defined as a function of the temperature (64) in which, the fluid temperature T f is in Kelvin. This function will give fluid viscosity of μ f = 0.001 Pa s at T f = 20°C. In this example, both the contact model and THM model are run sequentially at each timestep. The contact model is run using the pressure and temperature of the medium from the THM model, and the THM model is run using the contact stresses from the contact model. Two or three iterations are required in order to reach desirable convergence. The results for spatial distribution of the fluid temperature, fracture aperture, and contact traction within the fracture at time t = 10 9 s are shown in Fig.  5. The cold fluid at injection reduces the temperature of the rock matrix, which results in contraction of the matrix. The volumetric contraction of the matrix reduces the contact stress on the fracture and increases the fracture aperture around the injection well, and also towards the production well. The contact traction at the injector is reduced to about 30 MPa from an initial value of 40 MPa at t = 10 9 s, and the aperture at the injector increases to about 0.3 mm. The increased aperture creates a favourable path for the fluid to flow towards the pro Fig. 5. Spatial distribution of the fluid temperature, fracture aperture and contact traction within the fracture, at t = 10 9 s. The contact stress and the fracture aperture at the fracture tips is equal to zero. ducer, i.e., a channel, which results in low heat extraction from other parts of the fracture away from this path, as shown in Fig. 5. For instance, the area behind the producer remains relatively untouched as the injected fluid cannot reach this area. The evolution of the injection pressure and the fracture aperture at the injection point is shown in Fig. 6. At early time (t < 10 5 s), the injection pressure increases due to the increase in the viscosity of the fluid. The increased pressure reduces the contact stress and results in an increase in the aperture. At later times (t > 10 5 s), the cooling of the matrix starts to affect the sur Salimzadeh et al. Geothermics xxx (2017) xxx-xxx rounding rock, and as a result the aperture increases, and the injection pressure reduces as shown in Fig. 6.

Deformable fracture in a permeable matrix
In this case, a circular fracture of diameter 400 m is assumed in a plane that makes an angle of 30°with the horizontal direction, as shown in Fig. 7. Vertical injection and production wells are located at a 300 m distance from each other, and injection and production is performed through the rock matrix, as the wells are not directly connected to the fracture. Only the lower 20 m of the wells are assumed to be perforated. Cold water is injected at a rate of Q = 0.005 m 3 /s at a temperature of 20°C, and produced at the same rate, Q = 0.005 m 3 /s, at the producer. The rock has elastic properties of Young's modulus E = 20 GPa, Poisson's ratio ν = 0.20, and Biot coefficient α = 0.8. Rock and fluid properties are given in Table 1 (example 3). Three cases are assumed: case 1 with low initial temperature (T ini = 80°C) and low initial contact stress (σ ini = 60 MPa), case 2 with high initial temperature (T ini = 250°C) and low initial contact stress (σ ini = 60 MPa), and case 3 with high initial temperature (T ini = 250°C) and high initial contact stress (σ ini = 75 MPa). In all three cases, the initial fluid pres Fig. 7. Model geometry for the case with deformable fracture in a permeable matrix. The injection and production is to/from the rock matrix, as the wells are not directly connected to the fracture.
sure is set to 20 MPa, the permeability of the rock matrix is set to k m = 10 −14 m 2 , the fracture stiffness is set to 10 11 Pa/m, and the fracture aperture at zero contact stress is set to 0.6 mm.
The initial timestep is set to 1 day, and it is increased in each step by a factor of 1.1 until it reaches to a maximum timestep of 0.2 years. The results of the simulation for the fluid temperature, aperture and contact stress distributions at the fracture after 10, 20, and 30 years of injection for each case are shown in Figs. 8-10, respectively. Injection of the cold fluid induces contractions on the rock matrix around the injection well; such volumetric contractions reduce the contact stress on the fracture in the area most close to the injection well. The reduction in the contact stress results in an increase of the fracture aperture, and that creates a preferential path for the flow. As time elapses, the region of the fracture with decreased temperature, increased aperture, and decreased contact stress expands. The rate and the shape of the expansion depend on both initial aperture and initial rock temperature. Higher initial aperture makes the fracture a permeable pathway for the flow. Therefore, the cold front moves towards the fracture and that further increases the fracture aperture (cases 1 and 2). This results in developing an area with increased aperture in the fracture, as can be seen in Figs. 8 and 9, for both cases 1 and 2. That area points towards the nearest "exit" from the fracture towards the production well. Higher initial temperature leads to a larger temperature change in the matrix, which creates higher contraction followed by higher reduction in the contact stress and higher increase in the fracture aperture. Therefore, case 2, which has higher initial aperture and higher temperature variation, creates the most dominant favourable path for the flow of cold fluid, which is visible as early as 10 years.
Lower initial aperture (the initial aperture of 0.05 mm corresponding to the initial contact stress of 75 MPa) makes the fracture hydraulic conductivity to be on the same order as that of the matrix, so the fracture initially is not a preferential pathway for the fluid (case 3). However, as the cooling of the matrix occurs, the contraction of the rock reduces the contact stress on the fracture, which then increases the fracture aperture, and so the fracture becomes a preferential pathway for the flow, as shown in Fig. 10. The area with increased aperture, however, does not reach the same location as for cases 1 and 2. This is due to lower conductivity of the fracture ahead of the cold front, which prevents the movement of the cold front in the fracture. The distribution of matrix temperature along a horizontal cut-plane, after 30 years, is shown for the three cases in Fig. 11. In case 2, the fracture is clearly acting as short-circuit for the flow, whereas in the first case, the high initial fracture aperture allows the cold water to access larger areas of Salimzadeh et al. Geothermics xxx (2017) xxx-xxx Fig. 8. Spatial distribution of the temperature, aperture and contact stress on the fracture after 10, 20 and 30 years of simulation for low temperature case (T ini = 80°C) with σ ini = 60 MPa. The contact stress and the fracture aperture at the fracture tips is equal to zero. the fracture. In the third case, the lower initial aperture limits the distribution of the cold water on the fracture, but the aperture increase due to the contraction of the matrix creates a favourable path for the cold water to move towards the producer. Again, as the initial aperture in this case is very low, the size of the area with increased aperture is smaller than the one in case 2. In Fig. 12, the maximum increase in the fracture aperture, as well as the temperature drop at the producer, are compared for the three cases. The magnitude of the aperture increase in case 3 is the highest, while the temperature drop at the producer for case 2 is the highest. As mentioned earlier, the area with increased aperture in case 3 is smaller and therefore, less effective than the one in case 2. The temperature drop in the producer, as well as the maximum aperture increase, is the lowest in case 1. This is due to lower volume contraction of the matrix due to lower initial temperature of the reservoir, and also due to distribution of the cold water in the fracture due to the high initial aperture. In case 3, although the fracture has a lower initial permeability than the case 1, but the temperature break-through occurs earlier. It is interesting to note that the maximum aperture increase (i.e., the maximum contact stress reduction) rapidly reaches a maximum value at an early time (around seven years), and then decreases. The reduction in the aperture is due to the stress redis tribution in the fracture, as the region of the fracture with reduced contact stress expands.

Conclusions
A fully coupled THM model that rigorously models deformable fractures in a permeable matrix has been presented. The THM model is further coupled with a contact model to resolve the contact stresses between fracture surfaces. The model was validated and applied to several examples of geothermal systems, in both impermeable and permeable rocks. Conductive fractures create preferential paths for the flow, and the flow of the cold fluid reduces the temperature of the rock matrix surrounding these paths. The volumetric contraction of the matrix results in the local increase in the fracture aperture, i.e., channelling of the flow. In cases with a permeable matrix, the initial aperture of the fractures initially controls the flow of the cold fluid. However, as the matrix temperature decreases, the volumetric contraction of the matrix increases the aperture in the nearby fractures, which in turn become the preferential pathways for the flow. The contact stress on the fracture is reduced as the matrix contracts; however, the contact stress reaches a minimum value and then increases. The increase in the con Salimzadeh et al. Geothermics xxx (2017) xxx-xxx Fig. 9. Spatial distribution of the temperature, aperture and contact stress on the fracture after 10, 20, and 30 years of simulation for high temperature case (T ini = 250°C) with σ ini = 60 MPa. The contact stress and the fracture aperture at the fracture tips is equal to zero.
tact is due to the redistribution of the stresses due to the expansion of the region with reduced contact stress. In other words, as the area of the fracture affected by the matrix contraction expands, the stresses redistribute, which increases the minimum contact stress. The stress redistribution reduces the ability of the fracture to propagate under pure opening mode, while the expansion of the area with lower contact stress can increase the possibility of fracture propagation under shear. As future work, the computational method can be further improved by using parallel computing in order to simulate complex heterogenous media containing many discrete fractures. Fig. 10. Spatial distribution of the temperature, aperture and contact stress on the fracture after 10, 20, and 30 years of simulation for high temperature case (T ini = 250°C) with σ ini = 75 MPa. The contact stress and the fracture aperture at the fracture tips is equal to zero.