THERMAID  A matlab package for thermohydraulic modeling and fracture stability analysis in fractured reservoirs
Abstract
Understanding the dynamics of naturally fractured systems and fractured reservoirs in terms of flow, heat transport and fracture stability (e.g. induced seismicity) is important for a range of applications associated with waste water injection, renewable energy (e.g. geothermal systems), and greenhouse gas mitigation (e.g. geological sequestration of CO2). Here we present the implementation and validation of an open source MATLAB package for efficient numerical simulations of the coupled processes in fractured systems. We take advantage of the embedded discrete fracture model that efficiently accounts for discrete fractures. We perform a series of numerical benchmark experiments to validate the implemented approach against analytical solutions and established numerical methods. Finally, we use a simplified geomechanical model and an integrated fracture stability analysis that allows estimating the potential for shear stimulation, and thus a mechanistic assessment of induced seismic risk during stimulation. The open source distribution of the source code and results can be used as a blue print for the reimplementation of the method in a high performance computing (HPC) framework or as a standalone simulation package for investigating TH(m) problems in geothermal reservoirs.
keywords:
EDFM, embedded fracture, fracture stability analysis, fractured reservoir, opensource, MATLAB1 Introduction
A large fraction of the world’s water and energy resources are located in naturally fractured reservoirs within the earth’s crust. Understanding the dynamics of such reservoirs in terms of flow, heat transport and fracture stability is crucial to successful application of engineered geothermal systems (also known as enhanced geothermal systems, EGS) for geothermal energy production. Reservoir development characteristics such as permeability creation and induced seismicity largely depend on the properties of preexisting fractures, porosity, permeability and fracture orientation within the local stress field. One of the primary driving mechanisms for permeability creation in EGS involves shear failure induced by fluid injection at high pressures (Barton et al., 1995; Evans et al., 2005; Hickman and Davatzes, 2010). Along sections of the well that are free of natural fractures and in environments with low differential stress, tensile fractures may develop if the injection pressure exceeds the minimal principal stress. Shear and tensile fracture propagation and reactivation are not exclusive and might occur simultaneously during the stimulation of the reservoir (McClure and Horne, 2014). Clearly, preexisting, critically stressed and optimally oriented fractures provide the most favorable conditions for enhancing permeability of EGS (Barton et al., 1995; Combs et al., 2004; Ghassemi and Kumar, 2007).
The basis for EGS are fractured reservoirs, which are usually geothermal plays of the “hot dry rock” type where the available water in the porous medium is considered negligible (Brown et al., 2012). These conditions are found primarily in metamorphic or igneous terrains with low permeability and porosity, containing fractures and faults that provide the major pathways for fluid flow (e.g. Fenton Hill, Soultz, Basel, Cooper basin and Desert Peak (Kelkar et al., 2016; Hooijkaas et al., 2006; Häring et al., 2008; Chen and Wyborn, 2009; Hickman and Davatzes, 2010)). In geothermal energy systems, the fracture’s surfaces serve as the main heat exchanger. Fractured reservoirs can be considered to consist of two distinct separate media, the fractures and the matrix, and different types of reservoirs can be defined that depend on their properties (Nelson, 2001). In EGS, two cases typically prevail: 1) reservoirs with low porosity matrix for which both the permeability and the storage capacity of the rock mass are controlled by the fractures (cf. type 1 in (Nelson, 2001)) and 2) reservoirs with sufficient matrix porosity such that fluid storage is dominated by the matrix while the fractures contain only a small fraction of the fluid but control the permeability (cf. type 2 in (Nelson, 2001)).
Simulation of flow and transport through fractured porous media is challenging due to the high permeability contrast between the fractures and the surrounding rock matrix. However, accurate and efficient simulation of flow through a fracture network is crucial in order to understand, optimize and engineer reservoirs. Even after decades of research, this is still a very active research topic. Additionally, accurate estimations of the fracture stability are necessary in order to predict permeability evolution and forecast induced seismicity. Discrete fracture models (DFM) have been developed to address the computational problem of scales for fluid flow and heat transport. Various modeling frameworks for the simulation of fractured porous media in general and geothermal reservoirs in particular exist. Literature reviews of current modeling approaches in geothermal reservoirs and hot dry rock systems are presented by WillisRichards and Wallroth, Sanyal et al. and O’Sullivan et al.(WillisRichards and Wallroth, 1995; Sanyal et al., 2000; O’Sullivan et al., 2001). Some better known opensource modeling frameworks related to this work include PFLOTRAN (Lichtner et al., 2017), OpenGeoSys (Kolditz et al., 2012) and DuMux (Flemisch et al., 2011). Yet traditional conforming DFM, where the fractures are explicitly resolved by the numerical grid, suffer from computationally expensive preprocessing in the numerical grid generation and can encounter severe time step restrictions during the simulation when using explicit timestepping and small cells around the fractures (Norbeck et al., 2014; Sandve et al., 2012).
An alternative approach uses the embedded discrete fracture models (EDFM), which treat fracture and matrix in two separate computational domains. The embedded fracture model was first introduced by Lee et al. for single phase problems and later extended to twophase flow (Lee et al., 2001; Li et al., 2008). The embedded discrete fracture model is a promising technique in modeling the behavior of enhanced geothermal systems. Karvounis (Karvounis, 2013) employs EDFM and a statistical approach to better understand and possibly forecast seismicity induced seismicity by fluid injection during the stimulation phase of an EGS. Norbeck et al. additionally model fracture deformation by linear fracture mechanics (Norbeck et al., 2016).
In this paper we present the to our best knowledge first open source implementation of an embedded discrete fracture model for single phase flow and heat transport with additional capabilities to determine fracture stability in fractured reservoirs. Slip tendency analysis is used in order to estimate fault reactivation potential in earthquake prone areas as well as fracture stability in geothermal reservoirs (e.g. Morris et al., 1996; Moeck et al., 2009). Slip tendency is an indicator for the likelihood of slip. Using slip tendency, predictions on fracture instabilities during the hydraulic stimulation of a fractured reservoir are feasible without solving for the typically nonlinear evolution of the stress equilibrium equation. THERMAID, an acronym for ”ThermoHydraulic Energy Resource Modeling for Application and Development”, is a fractured reservoir modeling framework implemented in MATLAB, which can be used as a standalone simulation package for TH(m) cases in geothermal reservoirs or as a blue print for the reimplementation of the method e.g. in a high performance computing (HPC) framework. We coin the term TH(m) to indicate a coupled ThermoHydraulic code, and we use the lower case (m) to indicate simplified mechanics.
This paper is structured as follows. In the next section we present the methodology of the embedded discrete fracture model, and describe in detail the underlying theory of the fracture stability analysis. The implemented model is evaluated in the 4 Results section by comparing it with a widely used numerical model in several test cases. We conclude the paper by illustrating possible applications of the code using some examples and a discussion of the findings.
2 Methodology
The conceptual idea of the EDFM is the distinct separation of a fractured reservoir into a fracture and a matrix domain. We introduce a transfer function to account for coupling effects between the two domains (cf. Figure 1), so the fracture and matrix domains are computationally independent except for the transfer function. As the fractures are generally very thin and highly permeable compared to the surrounding matrix rock, the gradient of fluid pressure with the fracture normal to it is negligible. This allows for a lower dimensional representation of fractures (i.e. 1D objects within a 2D reservoir).
2.1 Conceptual model
Numerical modeling of fractured reservoirs is not only challenging from the numerical and computational point of view, but also because it involves a variety of coupled thermal, hydraulic, mechanical and chemical (THMC) processes. THERMAID focuses on thermohydraulic processes and their coupling, with some additional mechanical processes considered within a simplified geomechanical model.
Figure 2 shows the conceptual model of the most relevant thermal, hydraulic and mechanical processes in fractured reservoirs. The core processes implemented in THERMAID are fluid flow through pressure diffusion in the matrix and the fracture network and the accompanying heat transfer by advection and diffusion. Pressure and heat are also exchanged between the rock matrix and the fracture at the fracture walls by pressure and thermal diffusion and thermal advection, respectively.
Associated with the thermohydraulic processes, numerous thermomechanic or thermohydraulic processes could be activated, but only a limited amount and simplified processes are currently implemented in THERMAID. Poroelastic and thermoelastic deformation of the fractures and of the matrix is currently not implemented in the code. The fluid pressure in the fracture is, however, considered in the computation of fracture effective normal stress, which is an important parameter to evaluate fracture stability. If a fracture becomes unstable it will slip, and the associated dilation slip will increase the fracture transmissivity. The transmissivity increase is introduced in a simplified manner: if a fracture segment reaches the slip condition, its transmissivity is multiplied by a fixed permeability enhancement factor. In reality, slip on a fracture perturbs local the stress state, potentially affecting the stability of other fractures. This process is not implemented in the code because the direction and amount of slip is ambiguous. However, stress change induced by thermal changes is implemented, not yet in a fully coupled way, but instead computing the thermal stresses and superposing thermal stresses onto the ambient stresses.
In addition to the processes shown in Figure 2, THERMAID properly accounts for gravity effects, internal pressure and heat sources, and the pressure and temperaturedependence of fluid density and viscosity. The corresponding equations of state are given by (Sun et al., 2008) for density and (AlShemmeri, 2012) for the viscosity of water.
To summarize, THERMAID is a thermohydraulic code for fractured media that accounts for mechanical stability of the fractures, slipinduced transmissivity increase and thermally induced stresses. In the following, we introduce the governing equations for fluid flow and heat transport, couplings between fracture and matrix, and the implemented fracture stability analysis.
2.2 Governing equations
Flow in naturally fractured reservoirs is often described by the equations for nearly incompressible singlephase flow. We assume that the equations for nearly incompressible singlephase flow are valid in both matrix and the fractures. This simplification might not yield an adequate description of the flow in some fractured reservoirs where very large fracture apertures result in nonDarcian flow. The methodology presented here is however easily modifiable to extended Darcy flow models.
The pressure equation, derived from continuity and total mass balance equations for singlephase fluid flow, is:
(1) 
where [] is the porosity, is the fluid density, [Pa] is the fluid pressure and a source term. The compressibilities [Pa] are denoted with the subscripts for fluid and for rock, respectively, is the permeability and [Pas] the fluid viscosity. We consider only isotropic permeability . Permeability is often linked to fracture aperture through the Cubic law, which has been shown to be useful in predicting fluid transport through fractured reservoirs and fractured porous media in general. However, it does not account for the roughness of the fracture or flow adjacent to the fracture walls due to the rock permeability. From the fluid pressure , the fluid velocity is calculated using Darcy’s law, i.e.
(2) 
The total mass balance equation derived above is separated into parts for the matrix and the fracture domains, i.e.
(3) 
and
(4) 
where and are the flux transfer functions between the matrix and the fractures. Superscripts and denote matrix and fracture quantities respectively.
The heat transport equation is derived similarly to the pressure equation based on a balance of energy. We assume local thermal equilibrium so that where and are the temperatures of solid rock and fluid respectively. Taking average over an elemental volume we find
(5) 
where overlined properties denote volume averaged mean values for the porous medium.
(6)  
(7)  
(8) 
In equations 5 to 8 the heat capacity , the thermal conductivity and internal heat source of solid rock and fluid have been introduced. The fluid velocity used in the heat transport equation is the Darcy velocity given by equation 2.
The heat transport equation is separated into matrix and fracture parts according to the same procedure as for the fluid pressure equation
(9) 
and
(10) 
where and are the heat transfer functions between the damaged matrix and the fractures.
2.3 Fracture matrix coupling
To obtain a conservative set of equations, we apply a transfer function governing the mass and heat exchange between the two domains. The transfer function is treated as a source/sink term in the pressure and heat transport equations for damaged matrix and fracture, respectively, similar to classical well models (Peaceman et al., 1978).
The transfer function for the pressure equation is defined as
(11) 
with being the mean total mobility of the fluid, defined as the fraction of permeability and viscosity (Lee et al., 2001). is the connectivity index between matrix and fracture that is dependent on the numerical discretization (cf. next section). From the separated mass balance equations, it becomes immediately clear that the total flux between matrix and fracture has to be conserved:
(12) 
The transfer function for the heat equation is similarly defined. However, as two heat transport mechanisms are present in the equation, the transfer function needs to account for both mechanisms. Thus, the transfer function is defined as:
(13) 
where the superscript denotes the heat advection contribution and denotes the heat conduction contribution. The heat conduction contribution is derived using the same approach as in the pressure transfer function.
(14) 
Here, is the heat conductivity at the fracturematrix interface which can be calculated as
(15) 
using the definition of the averaged heat conductivity given in equation 7. The advection contribution , on the other hand, explicitly shows the coupling to the pressure equation based on the Darcy velocity:
(16) 
In equation 16 we introduce the fluid velocity and specific heat capacity at the matrixfracture interface. is calculated analogous to equation 15 and based on the averaged specific heat capacity given in equation 6. The fluid velocity at the matrixfracture interface is defined as
(17) 
where is the pressure gradient at the interface of matrix and fracture. As discussed for the pressure transfer function also the heat transfer flux has to be conserved:
(18) 
2.4 Fracture stability
Within THERMAID, a simplified analytical approach to fracture slip enables us to estimate fracture stability based on slip tendency analysis. Following Amonton’s law for purely frictional fault reactivation
(19) 
with as shear stress, as effective normal stress ( and as fluid pressure), and as sliding friction coefficient (Byerlee, 1978), slip tendency is the ratio of shear stress to effective normal stress on a surface (Morris et al., 1996), i.e.
(20) 
Fracture failure or slip is likely to occur if the shear stress to effective normal stress ratio equals or exceeds the frictional sliding resistance . Thus we define the stability of a fracture as follows
(21) 
Shear and effective normal stress acting on a given fracture depend on the orientation of the fracture plane within the effective principal stress field. If the effective principal stress field and the reference coordinate system of the simulation match, shear and normal stress can be calculated by simple expressions based on the dip angle of the fracture (e.g. Miller et al., 2004). Otherwise, stress transformations are needed in order to calculate the correct normal and shear stress in the fracture coordinate system . The involved rotation matrices can be calculated if the orientations of the principal stresses are known (e.g. Allmendinger et al., 2011). Finally we compute the normal and shear stress on the fracture in the fracture coordinate system based on the 3D principal stresses. Note that in this calculation the influence of the intermediate principal stress component is taken into account despite the 2D model geometry.
An important addition to the effects of pore pressure and far field stresses for fracture stability is thermal stress. A body will change its shape and/or volume when exposed to a temperature change . If the body’s deformation is restricted, as it would be the case for a small volume inside a rock mass, the strain results in thermal stress.
(22) 
where is the coefficient of linear thermal expansion in , is the Young’s modulus () and the Poisson ratio (). is the Kronecker delta, which is 1 for identical indices and , and 0 otherwise. The thermal stress is positive (relative compression) if the temperature difference is positive (), and if the temperature difference negative, the thermal stress is negative (relative tension). In the following we assume that the thermal stress is independent of the fluid pressure and the insitu stress state of the rock. Thus, the resulting stress can be obtained by superposition of the effective stress () and the thermal stress. We formulate the superposed effective stress as
(23) 
which can be used in equation 19 in order to account for thermal stress during the fracture stability analysis.
Clearly, other stress contributions as slip induced stresses and stresses induced by chemical reaction have to be considered in a general case. However, especially the estimation of slip induced stress changes is ambiguous as the amount of slip and slip direction for potentially failing fractures is not known apriori unless the underlying equations for fracture slip are solved explicitly. Thus, for reasons of simplicity we restrict ourselves to only effective stress and thermally induced stress changes.
As fractures are reactivated they generally show an increase in aperture as the fracture surfaces are not smooth but have many asperities. Due to a strong aperture dependence of permeability (e.g. Nemčok et al., 2002), where small changes in aperture result in very large changes in permeability, it can be assumed that unstable (or sliding) fractures undergo a stepwise change in fracture permeability (Miller and Nur, 2000). Here we adopt the most simple model
(24) 
where is an permeability enhancement factor. This model successfully described the distribution of the induced seismicity in the Basel EGS site, and fluiddriven aftershock sequences (Miller et al., 2004; Miller, 2015).
3 Implementation
We implemented the twodimensional embedded discrete fracture method in MATLAB. Our implementation is based on the concepts used in MAFLOT, an open source MATLAB flow solver (Künze and Lunati, 2012). As briefly discussed in the introduction, the matrix and fracture domains are discretized by regular Cartesian grids in 2D for the matrix and 1D for the fractures respectively (cf. Figure 1).
3.1 Numerical discretization in space
Using a finite volume approach, we discretize the domain as the integration over finite control volumes with . Using the Gauss theorem, the divergence integral over the volume can be rewritten as the surface integral normal to the boundary of the volume. Applied to a matrix grid cell on the right hand side (RHS) of the pressure equation 3 this yields
(25) 
Note that gravity is neglected here and in the remains of this section to better facilitate comprehension of the implementation. The pressure gradient over the cell boundary is approximated by a twopoint flux approximation that is secondorder accurate in space. As the domain is generally heterogeneous in terms of rock properties, a harmonic averaging technique is used to calculate the appropriate values at the cell boundaries.
The discretization of the RHS of the temperature equation is analogous to the pressure equations and omitted here for brevity. It is worth noting, however, that the advection term must be treated with special care. In this EDFM implementation, we use an upwind method in the fractures in combination with a minmodflux limited QUICK scheme in the matrix (Courant et al., 1952; Leonard, 1979; Roe, 1986).
3.2 Connectivity index
The connectivity index between matrix and fracture is discretizationdependent, and defined based on the linear pressure distribution assumed within a grid cell intersected by a fracture (Hajibeygi et al., 2011). It is defined as the length fraction of fracture segment inside matrix cell divided by the average distance between matrix cell and fracture segment .
(26) 
The average distance can be calculated as
(27) 
where is the distance from the fracture within the matrix cell and the volume of the matrix cell. This allows a proper accounting for the reduced influence of a fracture segment on a matrix cell if the fracture segment does not cross the matrix cell through its center. In many cases equation 27 has to be evaluated by numerical integration. For rectangular grids however, there exists an analytical solution (Hajibeygi et al., 2011; Pluimers, 2015). For enhanced efficiency, the analytical expressions are used in our implementation.
3.3 Fracture intersections
Fracture intersections significantly impact flow dynamics in the reservoir. The additional fracturefracture transmissivity can be calculated as
(28) 
where denotes the fracture aperture, the total mobility and the numerical discretization spacing in the fracture (KarimiFard et al., 2003).
3.4 Timediscretization
The time derivatives in equations 1 and 5 are treated using the backward Euler method, which is an implicit timediscretization with local truncation error . The method is unconditionally stable theoretically allowing arbitrarily large time steps. In practice, when encountering nonlinear behavior, such as the temperature and pressuredependent evolution of fluid density, issues with nonconvergence might appear and place an indirect restriction on the timestep. Nonetheless, much larger time steps are allowed in the implemented method when compared to explicit schemes.
3.5 Solution strategy
We adapt a serial iterative scheme in order to accurately account for the coupling between the pressure and transport equations. In strongly coupled problems, multiple iterations must be used to capture any arising nonlinearities. In most cases, the flow and transport exhibit rather loose coupling in which only a few iterations are needed to converge to the solution. If on the other hand, fracture stability ceases and permeability enhancement in unstable fracture parts is used, the number of iterations might increase significantly and even limit the timestep.
4 Results
We present the results of three benchmark experiments and an application experiment that provide insight into the capabilities of THERMAID and validate the implemented method. Fracture permeability and aperture are treated as independent from each other in the following. This allows simulating ’filled’ fractures with relatively high aperture and comparably small permeability and allows fracture permeability estimates independent of Cubic law.
First we validate the implemented model with a simple flow problem independently. We then evaluate the coupled results of fluid flow and heat transport on a simple geometry and on a more realistic complex fracture network. The final numerical experiment is the application of the implemented approach to a field scale problem were we take advantage of the implemented fracture stability analysis in order to characterize the stimulated reservoir during injection of a geothermal reservoir.
4.1 Validation of the pressure equation
In order to validate the implementation of flow equations of the model, we use an analytical solution for the steadystate flow in a porous medium in the presence of a fracture (Strack, 1982; Kolditz et al., 2012). Figure 2 shows the benchmark geometry, a square with a length of 10 m with a 2mlong inclined fracture in the center of a square domain. The aperture of the fractures is fixed at . Uniform flow is maintained by imposing a specific discharge from the left boundary into the domain. To compare numerical results with the analytical solution, pressures calculated by the analytical solution are used at the lateral boundaries, i.e. Pa and Pa (cf. Figure 2). On the top and bottom a noflow Neumann boundary is applied. The remaining material properties of the numerical model are shown in Table 2.
Parameter  Value  Unit  

Fracture angle  45  
Maximum fracture aperture  0.05  m  
Fracture length  2  m  
Matrix permeability  m  
Fracture permeability  m  
Fluid viscosity  Pa s  
Specific discharge  m s 
The pressure distribution obtained by THERMAID is shown in Figure 3a. The lateral uniform flow is disturbed in the vicinity of the inclined fracture where the flow is faster than in the surrounding porous media. Figure 3b shows the pressure profile along a diagonal line from the bottomleft to the topright. The results show very good agreement between the numerical solution obtained by THERMAID and the analytical solution. We quantify the difference between our model with the reference by the ’normalized root mean squared error (NRSME)’ as well as the ’normalized mean absolute error (NMAE)’.
(29)  
(30) 
We decided to use two measures of performance due a recent debate on both measures (e.g. Willmott and Matsuura, 2005; Chai and Draxler, 2014). Especially (Chai and Draxler, 2014) suggest that a combination of measures is required to assess model performance. We observe errors of well below (NRSME: and NMAE: ) that validate the implementation of the fluid flow equations.
4.2 Validation of the heat transport equation
We validate the coupled flow and heat transport equations using a benchmark geometry that consists of two perpendicular 5mlong fractures intersecting in the middle of a square domain (cf. Figure 4). The aperture of the fractures is fixed at . The domain is 100m by 100m square domain with Dirichlet boundary conditions on the left and right sides. On the left a constant pressure of 10MPa is applied, whereas the right side is fixed to 0MPa. On the top and bottom a noflow Neumann boundary is applied. The domain is initially at , which is a typical temperature for economic heat extraction in a geothermal reservoir. The inflow temperature at the left side of the domain is set to (cf. Figure 4). The material parameters for this benchmark were chosen realistically and are shown in Table 3. The benchmark’s results are evaluated after 40 years of simulation. The matrix domain is discretized by 301x301 cells while the fractures are modeled by 304 fracture segments (152 each). The reference solution is computed by COMSOL on a conforming discrete fracture network with a high resolution grid.
Permeability  

Porosity  
Density  
Viscosity  
Specific heat  
Heat conductivity 
Figure 5 shows the temperature in both fractures after 40 years of coupled flow and heat transport simulation. Additionally, Table 4 shows the quantitative error analysis for the fracture temperatures. We observe a very good agreement between the temperature distribution in both fractures with the reference solution. The horizontal fracture presents changes in temperature over most its extent, which is in accordance with the principal flow direction. As the vertical fracture is not aligned with the flow, a rather homogeneous temperature decrease is observed to about after 40 years. This is in good agreement with the matrix temperatures at the position of the fracture. Nevertheless, a significant change in temperature is observed close to the intersection of both fractures. Here the fracturefracture interaction is clearly visible as both fractures show nearly identical temperatures at the intersection (cf. Figure 5). The quantitative error analysis shows differences between our solution and the reference of for the vertical fracture and for the horizontal fracture although the two error measures differ slightly (cf. Table 4).
NRMSE [%]  

NMAE [%] 
Ultimately, the benchmark shows that our model accurately solves the coupled flow and heat transport equations for this geometry. The simulated timeframe is consistent with the estimated lifetime of a typical enhanced geothermal reservoir and additionally shows that the implemented timemarching scheme is accurate for the problem at hand.
4.3 Validation of the heat transport equation on a complex fracture network
We evaluate the coupled flow and heat transport on a more complex fracture geometry. The geometry (Figure 6) consists of a total of 13 fractures within a square domain. Boundary and initial conditions are equal to the previous experiment. The fracture aperture is set to . The remaining parameters governing the heat transport are consistent with the benchmark in the last section and shown in Table 3. We evaluate the results after 40 years of simulation. The reference solution computed by COMSOL contains 419’594 DOF. In this experiment we evaluate also grid dependence of the implemented model by comparing the results for different resolution simulations with the reference.
In the previous section we focused on the temperature distributions in the fractures. Here we take a closer look at the matrix temperature distributions. Figure 7a shows the final pressure distribution for a matrix grid resolution of 301x301. The temperature distribution in the domain after 40 years of simulation is shown in 7b. Both pressure and temperature fields show a heterogeneous distribution due to the influence of the fractures.
Figures 8a and 8b show the percental deviation of our solution from the reference for the matrix grid resolution of 301x301 of our model. The pressure solution shows only small errors with a NRMSE of . In the lower third of the domain between 20m and 60m in direction, a region of elevated error () is present (cf. Figure 8a). Bigger deviations are visible close to some fracture tips where typically on the highpressure (inflow) side of the fracture our model overestimates the matrix pressure compared to the reference. The lowpressure (outflow) sides of the fractures show predominantly underestimations of pressure. Interestingly, fractures that exhibit error concentration around one of the tips, do not necessarily show the opposite error on the other side of the fracture. Maximum pressure deviations from the reference are below .
The errors in the temperature distribution are generally larger than for the pressure. Figure 8b shows the percentage error at the final stage of the simulation for a matrix grid discretization of 301x301. Compared with the error in the pressure solution, we find that our model seems to always overestimate the matrix temperature compared to the reference. The normalized RMS error for this resolution is . We suspect that the elevated temperature deviations are caused by the relatively small error in the pressure solution. The small error in the pressure leads to comparably larger differences in flow velocities that are controlling heat advection. Thus, over a simulation of 40 years this error accumulates to the values observed here.
We further want to investigate the influence of the resolution on the accuracy of the results, so we compare four different resolutions (101x101, 301x301, 501x501 and 1001x1001). We investigate the improvement of solution accuracy in the pressure and heat transport solutions by using the NRMSE and NMAE values compared to the high resolution solution obtained by COMSOL. Table 5 shows both error measurements for all resolutions. We find a general improvement of the accuracy with an increase in resolution. For the temperature, this is a decrease in NRMSE from (101x101) to (1001x1001). The pressure error is consistently about one magnitude smaller, showing a decrease from (101x101) to (1001x1001). Overall we find a significant increase in accuracy with an increase in resolution. Nevertheless, the deviation is not changing significantly between 501x501 and 1001x1001 ( vs in case of the temperature). This indicates a systematic difference between the reference solution and our method. There are multiple possible origins of this systematic error. Since we observe the systematic deviation also in the pressure, we think it is likely to be a difference in methodology concerning the fluid flow equation. These differences could include the treatment of fracturefracture intersections, the definition of matrixfracture interface permeability, and inherent numerical differences between finite element and finite volume methods. Nevertheless, we find very good agreement between the reference simulation and our implementation for large parts of the model. Even in regions of significantly elevated deviation, we find acceptable agreement with differences below between the two methods. The definitive source of the difference is currently not resolved but presents excellent future research opportunities.
NRMSE [%]  NMAE [%]  NRMSE [%]  NMAE [%]  
101 x 101  0.78  0.64  4.59  4.0 
301 x 301  0.35  0.26  2.22  1.87 
501 x 501  0.23  0.17  1.74  1.51 
1001 x 1001  0.12  0.12  1.11  1.02 
4.4 Utilization of the fracture stability analysis
We present the fracture stability analysis to show the influence of permeability enhancement and thermal stress on fracture stability. We model fluid injection into a complex fracture network with a range of fracture orientations. The geometry consists of a total of 196 fractures within a square domain (Figure 9). The borehole is located in the middle of the domain with an open hole section of 6m. The fracture aperture in the reservoir is set to . The remaining parameters used in this section are shown in Table 6. The upper limit of the fractured reservoir domain is assumed to be at 5km depth. The injection pressure is held constant at 25MPa. The principal stresses are oriented as shown in Figure 10, which corresponds to a normal faulting regime. The magnitudes of the principal stresses are 125MPa, 107.5MPa and 81.25 MPa respectively, which corresponds to a relative stress ratio of . The insitu pore pressure is assumed to be hydrostatic (MPa). The stress conditions roughly resemble the relative conditions at the Fenton Hill and Hijiori geothermal projects although both projects were situated above 4km depth (Xie et al., 2015; Barton et al., 1988; Oikawa and Yamaguchi, 2000). We evaluate the results after 10 days of continuous fluid injection.
Permeability  

Porosity  
Density  
Viscosity  
Specific heat  
Heat conductivity  
Thermal expansion. coeff.  
Shear modulus  
Poisson ratio 
Figure 11 shows the pressure distribution after 10 days of injection. Due to the orientation of the preexisting fractures, a preferential flow direction in the vertical direction is visible. Slight pressure changes due to the injection are measured at distances up to 55m in the vertical and 35m in the horizontal directions from the injection point. The zone of 10MPa pressure changes extends roughly 10m around the borehole. Very high pressures MPa are concentrated in the direct vicinity of the injection.
The insitu fracture stability is influenced by the additional injected fluid pressure. Figure 12 shows the final normalized fracture slip tendency. Note that a normalized slip tendency value of 1 represents a fracture that is eligible for failure and slip. We observe a range of values in the reservoir based on the fractures’ orientations. The average fracture stability is high with values well below the failure condition. However, closer to the injection the increased slip tendency due to the injection is visible. Zones with fluid overpressure of MPa show significant increase in slip tendency (yellow colors in the plot). The region with at least 10MPa additional fluid pressure is very close to or eligible for slip on the fracture.
Permeability enhancement
In the previous section fracture segments eligible for slip did not have any feedback on the fluid pressure distribution. Here we investigate this feedback by introducing the stepwise permeability enhancement for failing fractures. The setup used is identical to the previous section except a 10fold increase in permeability is assumed for failing fracture segments.
Figure 13 shows the pressure distribution after 10 days of injection if permeability enhancement is considered. Although the general flow directions remain unchanged, the fluid pressure distribution shows significant differences in extent and magnitude.
Fracture stability is changed drastically if permeability enhancement is considered. The insitu fracture stability remains unchanged in the outer regions of the domain at very stable levels. On the other hand, most of the fractures within the overpressured regions show elevated slip tendency with fractures closer to the injection being eligible for slip. Compared to the previous simulation, 20times more fracture segments are unstable and capable of slip. Failing fractures, which have increased in permeability allow fluid to propagate more easily. As we assume a constant pressure injection, the amount of injected fluid is increased significantly. In this way a much larger stimulated area is observed compared to the case without permeability enhancement.
Thermal stress has only a small influence during the relatively short injection period of 10 days in this simulation. The resulting thermal stress distribution is shown in Figure 14 and shows thermal stresses concentrated at the borehole. As the fractures within the vicinity of significant thermal stress are eligible for slip also by the injection fluid pressure, no additional unstable fracture segments are observed. However, in a recent study investigating the role of thermal stress in a geothermal reservoir in detail we found that thermal stresses can facilitate slip on nonoptimally oriented fractures, and this is especially important in longterm injection scenarios where the thermal stress changes become more significant with time (Jansen and Miller, 2017).
The experiments presented here show the importance of fracture stability analysis. We showed that a stepwise permeability increase in potentially failing fracture segments has a major impact on the stimulated reservoir volume and allows fracture slip in larger parts of the domain. This emphasizes the importance of coupling thermohydraulic models with the mechanical changes during fracture slip.
5 Conclusion
We developed, implemented and validated a fractured reservoir modeling framework in MATLAB for investigating coupled thermohydraulic problems including fracture stability analysis.
Our results show with high confidence that the accuracy of the implemented MATLAB package are within the limits of commercial simulators for fractured reservoirs. Especially the results of the coupled flow and heat transport on a complex fracture network show the importance of discrete fractures in numerical analysis of fractured reservoirs. Both pressure and temperature distributions show heterogeneities due to fracturematrix interactions. THERMAID presents easy access to the underlying implementation that enables rapid prototyping as well as detailed investigations of the embedded discrete fracture model and coupled processes in naturally fractured reservoirs.
As discussed earlier in the results, the deviations in the pressure solution could be caused by different treatment of fracturefracture intersections or the definition of matrixfracture interface permeability between the models. Currently there is no clear indication about which weighting to use at fracturematrix interfaces, which presents an excellent future research opportunity for combined laboratory and numerical experiments. We assume that differences in the temperature solutions are caused by the deviations in the pressure solution that are magnified with time.
We showed that the embedded discrete fracture model is a viable alternative to the existing methods. As numerical discretization is simplified compared to conforming discrete fracture models, dynamic changes of the fracture network are possible without large numerical overhead. The extension of the embedded discrete fracture model to three dimensions has not been discussed so far in this article. Due to the relatively simple numerical discretization the extension to three dimensions is feasible. However, THERMAID is currently only developed in a 2D version. This is however not a limitation of the embedded fracture model but due limiting factors of the achievable computational performance in MATLAB. Nevertheless, the approach taken in THERMAID could be efficiently reimplemented and extended to 3D in a highperformance computing environment. The embedded discrete fracture model is not necessarily restricted to regular grids and can be extended to general geometries. However, using regular grids can be advantageous for the application of massively parallel computation techniques to further increase computational efficiency and enable large scale, high resolution simulations.
Our results show the importance of including the mechanical behavior of fractures and the reservoir in thermohydraulic simulations. Although the deformation process during fracture slip was not explicitly taken into account, the assumed stepwise increase in fracture permeability during slip provides the necessary feedback for the pressure equation in order to capture the observed increase in injectivity during hydraulic stimulation. We propose that changes in permeability and aperture should be incorporated in all models that seek to fully understand the thermohydraulic evolution during fluid injection in fractured reservoirs. Although our model exhibits a very simplified view on the complex fracture mechanics, it still provides important insight into reservoir stimulation that helps in identifying some challenges and opportunities for future studies. More advanced models currently under development will consider both preexisting fractures as in the present work, but also the generation of new fractures in response to the evolving stress state from both thermoand hydraulic perturbations. Future models might also include fracture roughness and solve the full equilibrium equations to estimate aperture changes that influence permeability. Recently, progress in this direction has been made using boundary element methods, multipoint stress approximations (MPSA) and the novel extended finite volume method (XFVM) (Norbeck et al., 2016; Ucar et al., 2016; Deb and Jenny, 2016). However, these models are not yet as computationally efficient as to allow an adaption for THERMAID. Currently, induced seismicity can not be quantified in terms of magnitude because slip on the fracture is not computed. Moreover, fracture slip can occur in a seismic or aseismic manner, thus further complicating the assessment of induced seismicity. These are all areas that we are currently pursuing in order to extend and refine THERMAID’s capabilities.
Although this paper focuses on the application of enhanced geothermal systems, other possible applications for THERMAID include seasonal thermal energy storage in fractured aquifers, and natural or anthropogenic fluiddriven earthquake sequences. Furthermore, a wide range of research questions related to fractured reservoirs and their properties can be addressed using THERMAID. Even beyond the current model capabilities, we expect further applications and research opportunities because the open source code will allow a community to evolve and contribute to this common platform. The open source distribution and GNU GPL v3.0 license enables the scientific community to use and modify THERMAID to their needs. The implementation in MATLAB ensures that even novice programmers can easily understand the underling equations and their implementation and develop their own numerical models based on the examples provided with THERMAID. Simulation of coupled processes in fractured reservoirs is becoming increasingly important in today’s research. With THERMAID we present an alternative starting point from which new insight can be gained into the complex coupled processes in fractured domains in the subsurface.
Acknowledgements
We thank the Swiss National Fond (SNF), for the financial support through the grant ’NFP70:Energy Turnaround’ under the project no. 153971.
Computer Code Availability

Project name: THERMAID

Project home page: https://github.com/gujans/THERMAID

Referenced archived version DOI: 10.5281/zenodo.1175829

Operating system(s): Platform independent

Programming language: MATLAB

Other requirements: built and tested with MATLAB R2015b

Licence: GNU GPL v3.0 or later

Any restrictions to use by nonacademics: The terms of the GNU GPL v3.0 or later apply.
Competing interests
The authors declare that they have no competing interests.
Funding
This work was supported by the Swiss National Fond (SNF), for the financial support through the grant ’NFP70:Energy Turnaround’ under the project no. 153971.
Bibliography
References
 AlShemmeri [2012] T. AlShemmeri. Engineering Fluid Mechanics. Bookboon, London, 2012.
 Allmendinger et al. [2011] R. W. Allmendinger, N. Cardozo, and D. M. Fisher. Structural geology algorithms: Vectors and tensors. Cambridge University Press, Cambridge, 2011.
 Barton et al. [1988] C. A. Barton, M. D. Zoback, and K. L. Burns. Insitu stress orientation and magnitude at the fenton geothermal site, new mexico, determined from wellbore breakouts. Geophysical Research Letters, 15(5):467–470, 1988.
 Barton et al. [1995] C. A. Barton, M. D. Zoback, and D. Moos. Fluid flow along potentially active faults in crystalline rock. Geology, 23(8):683–686, 1995.
 Brown et al. [2012] D. W. Brown, D. V. Duchane, G. Heiken, and V. T. Hriscu. Mining the earth’s heat: hot dry rock geothermal energy. Springer Science & Business Media, Berlin Heidelberg, 2012.
 Byerlee [1978] J. Byerlee. Friction of rocks. Pure and applied Geophysics, 116(4):615–626, 1978.
 Chai and Draxler [2014] T. Chai and R. R. Draxler. Root mean square error (rmse) or mean absolute error (mae)?–arguments against avoiding rmse in the literature. Geoscientific Model Development, 7(3):1247–1250, 2014.
 Chen and Wyborn [2009] D. Chen and D. Wyborn. Habanero field tests in the cooper basin, australia: a proofofconcept for egs. Geothermal Resources Council Transactions, 33(1):159–164, 2009.
 Combs et al. [2004] J. Combs, S. K. Garg, and J. W. Pritchett. Geothermal well stimulation technology: a preliminary review. Geothermal Resources Council Transactions, 28:207–212, 2004.
 Courant et al. [1952] R. Courant, E. Isaacson, and M. Rees. On the solution of nonlinear hyperbolic differential equations by finite differences. Communications on Pure and Applied Mathematics, 5(3):243–255, 1952.
 Deb and Jenny [2016] R. Deb and P. Jenny. Finite volume–based modeling of flowinduced shear failure along fracture manifolds. International Journal for Numerical and Analytical Methods in Geomechanics, 2016.
 Evans et al. [2005] K. F. Evans, A. Genter, and J. Sausse. Permeability creation and damage due to massive fluid injections into granite at 3.5 km at soultz: 1. borehole observations. Journal of Geophysical Research: Solid Earth, 110(B4), 2005.
 Flemisch et al. [2011] B. Flemisch, M. Darcis, K. Erbertseder, B. Faigle, A. Lauser, K. Mosthaf, S. Müthing, P. Nuske, A. Tatomir, M. Wolff, et al. Dumux: Dune for multiphase, component, scale, physics,… flow and transport in porous media. Advances in Water Resources, 34(9):1102–1112, 2011.
 Ghassemi and Kumar [2007] A. Ghassemi and G. S. Kumar. Changes in fracture aperture and fluid pressure due to thermal stress and silica dissolution/precipitation induced by heat extraction from subsurface rocks. Geothermics, 36(2):115–140, 2007.
 Hajibeygi et al. [2011] H. Hajibeygi, D. Karvounis, and P. Jenny. A hierarchical fracture model for the iterative multiscale finite volume method. Journal of Computational Physics, 230(24):8729–8743, 2011.
 Häring et al. [2008] M. O. Häring, U. Schanz, F. Ladner, and B. C. Dyer. Characterisation of the basel 1 enhanced geothermal system. Geothermics, 37(5):469–495, 2008.
 Hickman and Davatzes [2010] S. H. Hickman and N. C. Davatzes. Insitu stress and fracture characterization for planning of an egs stimulation in the desert peak geothermal field, nevada. In ThirtyFifth Workshop on Geothermal Reservoir Engineering. Stanford University, Stanford, California, 2010.
 Hooijkaas et al. [2006] G. R. Hooijkaas, A. Genter, and C. Dezayes. Deepseated geology of the granite intrusions at the soultz egs site based on data from 5kmdeep boreholes. Geothermics, 35(5):484–506, 2006.
 Jansen and Miller [2017] G. Jansen and S. A. Miller. On the role of thermal stresses during hydraulic stimulation of geothermal reservoirs. Geofluids, 2017, 2017.
 KarimiFard et al. [2003] M. KarimiFard, L. J. Durlofsky, K. Aziz, et al. An efficient discrete fracture model applicable for general purpose reservoir simulators. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 2003.
 Karvounis [2013] D. C. Karvounis. Simulations of enhanced geothermal systems with an adaptive hierarchical fracture representation. PhD thesis, ETH Zürich, 2013.
 Kelkar et al. [2016] S. Kelkar, G. WoldeGabriel, and K. Rehfeldt. Lessons learned from the pioneering hot dry rock project at fenton hill, usa. Geothermics, 63:5–14, 2016.
 Kolditz et al. [2012] O. Kolditz, H. Shao, U. Görke, T. Kalbacher, S. Bauer, C. McDermott, and W. Wang. Thermohydromechanicalchemical processes in fractured porous media: Benchmarks and examples. Springer, 2012.
 Künze and Lunati [2012] R. Künze and I. Lunati. An adaptive multiscale method for densitydriven instabilities. Journal of Computational Physics, 231(17):5557 – 5570, 2012. ISSN 00219991.
 Lee et al. [2001] S. H. Lee, M. Lough, and C. Jensen. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resources Research, 37(3):443–455, 2001.
 Leonard [1979] B. P. Leonard. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Computer methods in applied mechanics and engineering, 19(1):59–98, 1979.
 Li et al. [2008] L. Li, S. H. Lee, et al. Efficient fieldscale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Evaluation & Engineering, 11(04):750–758, 2008.
 Lichtner et al. [2017] P. C. Lichtner, G. E. Hammond, C. Lu, S. Karra, G. Bisht, B. Andre, R. T. Mills, J. Kumar, and J. M. Frederick. PFLOTRAN user manual. Technical report, 2017. http://www.documentation.pflotran.org.
 McClure and Horne [2014] M. W. McClure and R. N. Horne. An investigation of stimulation mechanisms in enhanced geothermal systems. International Journal of Rock Mechanics and Mining Sciences, 72:242–260, 2014.
 Miller [2015] S. Miller. Modeling enhanced geothermal systems and the essential nature of largescale changes in permeability at the onset of slip. Geofluids, 15(12):338–349, 2015.
 Miller and Nur [2000] S. A. Miller and A. Nur. Permeability as a toggle switch in fluidcontrolled crustal processes. Earth and Planetary Science Letters, 183(1):133–146, 2000.
 Miller et al. [2004] S. A. Miller, C. Collettini, L. Chiaraluce, M. Cocco, M. Barchi, and B. J. Kaus. Aftershocks driven by a highpressure co2 source at depth. Nature, 427(6976):724–727, 2004.
 Moeck et al. [2009] I. Moeck, G. Kwiatek, and G. Zimmermann. Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir. Journal of Structural Geology, 31(10):1174–1182, 2009.
 Morris et al. [1996] A. Morris, D. A. Ferrill, and D. B. Henderson. Sliptendency analysis and fault reactivation. Geology, 24(3):275–278, 1996.
 Nelson [2001] R. Nelson. Geologic analysis of naturally fractured reservoirs. Gulf Professional Publishing, Houston, Texas, 2001.
 Nemčok et al. [2002] M. Nemčok, A. Henk, R. A. Gayer, S. Vandycke, and T. M. Hathaway. Strikeslip fault bridge fluid pumping mechanism: insights from fieldbased palaeostress analysis and numerical modelling. Journal of Structural Geology, 24(12):1885–1901, 2002.
 Norbeck et al. [2014] J. Norbeck, H. Huang, R. Podgorney, and R. Horne. An integrated discrete fracture model for description of dynamic behavior in fractured reservoirs. In Proceedings of the 39th Workshop on Geothermal Reservoir Engineering, Stanford, 2014.
 Norbeck et al. [2016] J. H. Norbeck, M. W. McClure, J. W. Lo, and R. N. Horne. An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimulation. Computational Geosciences, 20(1):1–18, 2016.
 Oikawa and Yamaguchi [2000] Y. Oikawa and T. Yamaguchi. Stress measurement using rock core in an hdr field. In Proceedings of World Geothermal Congress, pages 3819–3822, 2000.
 O’Sullivan et al. [2001] M. J. O’Sullivan, K. Pruess, and M. J. Lippmann. State of the art of geothermal reservoir simulation. Geothermics, 30(4):395–429, 2001.
 Peaceman et al. [1978] D. W. Peaceman et al. Interpretation of wellblock pressures in numerical reservoir simulation (includes associated paper 6988). Society of Petroleum Engineers Journal, 18(03):183–194, 1978.
 Pluimers [2015] S. Pluimers. Hierarchical fracture modeling approach. Master’s thesis, TU Delft, Delft University of Technology, 2015.
 Roe [1986] P. L. Roe. Characteristicbased schemes for the euler equations. Annual review of fluid mechanics, 18(1):337–365, 1986.
 Sandve et al. [2012] T. Sandve, I. Berre, and J. Nordbotten. An efficient multipoint flux approximation method for discrete fracture–matrix simulations. Journal of Computational Physics, 231(9):3784 – 3800, 2012. ISSN 00219991.
 Sanyal et al. [2000] S. K. Sanyal, S. J. Butler, D. Swenson, and B. Hardeman. Review of the StateoftheArt of numerical simulation of enhanced geothermal systems. In Proceedings of World Geothermal Congress, 2000.
 Strack [1982] O. D. Strack. Analytic modeling of flow in a permeable fissured medium. Technical report, 1982.
 Sun et al. [2008] H. Sun, R. Feistel, M. Koch, and A. Markoe. New equations for density, entropy, heat capacity, and potential temperature of a saline thermal fluid. Deep Sea Research Part I: Oceanographic Research Papers, 55(10):1304–1310, 2008.
 Ucar et al. [2016] E. Ucar, E. Keilegavlen, I. Berre, and J. M. Nordbotten. Finitevolume discretization for the deformation of fractured media. arXiv preprint arXiv:1612.06594, 2016.
 WillisRichards and Wallroth [1995] J. WillisRichards and T. Wallroth. Approaches to the modelling of hdr reservoirs: a review. Geothermics, 24(3):307–332, 1995.
 Willmott and Matsuura [2005] C. J. Willmott and K. Matsuura. Advantages of the mean absolute error (mae) over the root mean square error (rmse) in assessing average model performance. Climate research, 30(1):79–82, 2005.
 Xie et al. [2015] L. Xie, K.B. Min, and Y. Song. Observations of hydraulic stimulations in seven enhanced geothermal system projects. Renewable Energy, 79:56–65, 2015.