Flexible parallel implicit modelling of coupled Thermal-Hydraulic-Mechanical processes in fractured rocks

Theory and numerical implementation describing groundwater flow and the transport of heat and solute mass in fully saturated fractured rocks with elasto-plastic mechanical feedbacks are developed. In our formulation, fractures are considered as being of lower dimension than the hosting deformable porous rock and we consider their hydraulic and mechanical apertures as scaling parameters to ensure continuous exchange of fluid mass and energy within the fracture-solid matrix system. The coupled system of equations is implemented into a new simulator code that makes use of a Galerkin Finite Element 5 technique. The code builds on a flexible, object oriented numerical framework (MOOSE, Multiphysics Object Oriented Simulation Environment) which provides an extensive scalable parallel, implicit coupling to solve for the multiphysics problem. The governing equations of groundwater flow, heat and mass transport and rock deformation are solved in a weak sense (either by classical Newton-Raphson or by Free Jacobian Inexact Newton Krylow schemes) on an underlying unstructured mesh. Nonlinear feedbacks among the active processes are enforced by considering evolving fluid and rock properties depending on the 10 thermo-hydro-mechanical state for the system and the local structure, i.e. degree of connectivity, of the fracture system. A suite of applications is presented to illustrate the flexibility and capability of the new simulator to address problems of increasing complexity and occurring at different spatial (from centimeters to tens of kilometers) and temporal scales (from minutes to hundreds of years).


Introduction
Reliable predictions of reservoirs' performance, whether for geothermal and fossil energy extraction, or water, CO 2 , and heat storage rely on an accurate representation of the physical processes responsible for groundwater flow, heat and solute mass transfer, mechanical deformation of the rock solid skeleton, and, ultimately, chemical feedbacks from fluidrock interactions (Stephansson et al., 2004).All of these processes occur in a naturally complex geological setting, comprising discrete heterogeneities as faults and fractures, span a relatively large spectrum of temporal and spatial scales and interact in a highly nonlinear fashion.In natural and engineered systems, onset conditions and the evolution in time and space of a particular process are affected by the initiation and evolution of all the other processes.Therefore, monitored variations of rock masses to natural and anthropogenic perturbations cannot be fully reconciled by considering the causative processes independently.This is particularly the case for reservoir applications which require a complete understanding of the multicomponent (fractured) porous rockfluid system and its multiphysics dynamics to predict the behaviour of a particular reservoir so as to improve its productivity and sustainability.
To recognize that the majority of geological observations are the results of nonlinear interactions among classes of physical processes has justified the recent development and massive use of so-called multiphysics numerical simulators as complementary tools to classical experimental and theoretical analyses.Numerical simulations are particularly adapted to the study of a multicomponent physical system Published by Copernicus Publications on behalf of the European Geosciences Union.
since they allow a systematic analysis of the dynamics of single processes and their interactions (i) under proper time and length scales, (ii) within complex geometries, and (iii) also by considering changing loading conditions.In addition, given their predictive capabilities, numerical modelling techniques can assist studies aiming at estimating performances and potential risks related to different operational scenarios (e.g.Blöcher et al., 2015).
Numerical modelling of fractured reservoirs, whether petroleum or geothermal, has a long history dating back several decades.Starting from the earlier works in the 1970s, a number of codes have been implemented under various licences by the modelling community (Jing, 2003;Stephansson et al., 2004).The majority of these codes rely on modelling techniques that were developed in the late 1980s and early 1990s, and that can be grossly subdivided into two classes if based on the approach followed to model coupled processes.
A first class of approaches, referred to as sequential/explicit coupling or operator splitting approaches, relies on splitting the coupled physical problem into classes of possibly linear sub-problems and to numerically solve for each process sequentially.During each time step, coupling among the processes is enforced by passing input/output data among the respective sub-problems.Usually, external iterations are adopted for simulations characterized by a high degree of nonlinearity (Kempka et al., 2016;Chabab and Kempka, 2016).In some cases, sequential coupling is achieved by relying on different simulators, as for example by coupling a flow simulator to a mechanical simulator for a thermohydro-mechanical problem (TOUGH-FLAC family of codes, e.g.Rutqvist, 2011).The main advantage of such approaches stems from numerically integrating relatively complex problems within limited computer resources.However, sequential coupling schemes have important impacts on the efficiency, stability, and accuracy of the numerical solutions.They usually introduce a splitting error in the numerical approximation, which requires a careful monitoring of the nonlinear residuals within each single time step (Jha and Juanes, 2007;Preisig and Prévost, 2011).This aspect limits sequential approaches to what are generally referred to as "loosely coupled problems", and they show a relatively slow, if at all, convergence rate for tightly coupled problems.In addition, conditional stability of semi-implicit approaches impose severe time step restrictions thus increasing computation times.
As an alternative, it is possible to solve the full system of coupled equations simultaneously, via a fully implicit coupling approach.This requires to solve for all the variables of the problem simultaneously within an iterative approach, either (in)exact Newton or simpler Picard methods, for the resulting algebraic system (Knoll and Keyes, 2004).Simultaneous solution schemes demand higher computer storage and processing times to compute (allocate the memory and fill) the Jacobian of the problem than explicit approaches.However, they guarantee higher stability for strongly nonlinear problems even with relative large time step sizes (Kim et al., 2015).In addition, recent advances in developing block-type preconditioners and "super-convergent" methods for systems of nonlinear equations at reduced memory consumption, e.g.preconditioned Jacobian free methods, make these approaches competitive for field-scale applications (Knoll and Keyes, 2004).Despite all of this, only few attempts have been made so far in developing and implementing fully implicit numerical solutions for thermal-hydraulic-mechanical (THM) problems.
The FRACON code of Nguyen and Selvadurai (1998) considers a monolithic small strain and non-isothermal implementation, but thermal feedbacks effects on the skeleton deformation and pore fluid pressure are only explicitly integrated.The open-source, object-oriented project Open-Geosys (Kolditz et al., 2012) provides a parallel platform for implicit solution of multiphysics problems (Blöcher et al., 2015) but so far, to the authors' knowledge, has never been applied to 3-D THM problem in fractured reservoirs.More recently, Sun (2015) developed a monolithic framework for solving coupled THM processes at finite strain.However, such an approach has only been applied on generic and simplistic models and lacks a description of complex geological geometries such as fractures and faults.More recent developments have been published focusing on quantifying chemical feedbacks for tightly coupled problems though mainly targeting simplistic geometries (Poulet and Veveakis, 2016).
The goal of this paper is not to summarize the state-ofthe-art computational methods for problems that are relevant for reservoir applications.Our interest herein is rather towards computational reliability and performances when simulating the behaviour of a particular reservoir in a way that can be of help to improve scenario-oriented analysis of such systems.In this context, we address issues related on how (i) to quantify the nonlinear feedbacks among the different physical processes, and (ii) to represent into a computational model the porous rock-fracture-fluid system by capturing its discontinuous, anisotropic, inhomogeneous, and non-elastic nature (Hudson and Harrison, 1997).For this purpose, we give an overview of the methods implemented into a novel, yet robust and efficient multiphysics and multicomponent porous media open-source modelling simulator called GOLEM which can deal with all these aspects.Our emphasis throughout this paper is to simulate THM processes of relevance for hydrothermal and geothermal systems.Though GOLEM can also simulate the transport of non-reactive chemical species, we do not discuss this aspect in the present paper.Consideration of additional chemical (fluid-solid) interactions is the subject of future work.
The remainder of the paper is structured as follows.In Sect.2, we introduce the constitutive model for THM processes for a two-phase system consisting of a deformable solid skeleton and fully saturated pores in the presence of discrete heterogeneities as represented by faults and fractures.
Attention is given to include the important nonlinear feedbacks among the different processes as well as structureproperty relationships.These latter include links between fluid properties (e.g.density and viscosity) and problem variables (e.g.pressure and temperature), as well as evolution of material properties as a function of variations in the state of the reservoirs (e.g.porosity and permeability relationships).In Sect.3, the numerical implementation of the derived constitutive model is described.In Sect.4, we apply the simulator to a suite of applications of different level of complexity (from synthetic benchmark cases to prototypical reservoir simulations) both in terms of the coupling among the processes considered as well as in terms of the geometry of the natural system.A field application related to modelling the hydromechanical response of a sandy reservoir during a cyclic injection test will be presented in a separate paper (Jacquey et al., 2017).

Mathematical formulation of the problem
In this section, we present the constitutive model for a porous medium consisting of a deformable hosting rock and a mobile liquid saturating its pores.In the following we make use of subscripts f and s to refer to the solid and fluid phase respectively.Deformation of the solid skeleton is considered in terms of Biot consolidation theory (Biot, 1956;Biot and Willis, 1957;Biot, 1973) under non-isothermal loading (Geertsma, 1957;McTigue, 1986) and its extension to take into account time-and rate-dependent inelastic behaviour.The model assumes the existence of a reference elementary volume (REV) for the porous medium in which the two components can interact thermo-hydro-mechanically (Bear, 1988).The resulting physical system therefore incorporates non-isothermal flow of the fluid phase (namely water) within a porous rock which is free to deform (in)elastically and in the presence of discontinuities as represented by discrete fractures (Boley and Weiner, 1960).
Fractures are considered as lower-dimensional elements embedded in the porous matrix filled with fluid (see Fig. 1).We consider the length of the discrete fracture as a representative measure of the REV of the porous matrix.This allows us to make use of a continuum approach to represent the porous medium.We also take into account the effects of micro-defects (fissures and micro-cracks) on the thermal, hydraulic, and mechanical behaviour of the porous rock as they mainly affect the evolution of the system material properties (e.g.porosity; see Fig. 2).
In the following the main governing equations are derived.The problem variables considered are the pore fluid pressure p f , the temperature T , and the solid displacement vector u.Pore pressure is defined as compression positive for water, while stress is defined as tension positive for the solid phase.
The mass balance equation for a deformable, saturated porous medium is described in terms of volumetric averaged mass conservation equations for the fluid and solid phases.Mass conservation for the liquid phase therefore requires where ρ f is the density of the fluid phase, n is the porosity, v f the fluid velocity, and Q f is a sink/source term considered null hereafter for the sake of simplicity.In a similar way, for the solid phase we obtain where ρ s is the density of the solid phase, v s the solid velocity, and Q s is a sink/source term (considered as null hereafter).We use Darcy's law to describe the conservation of momentum of the fluid phase, which can be expressed in terms of fluid velocity relative to the solid velocity as where q D is a volumetric flow rate per unit of surface area (Darcy velocity), k is the permeability tensor of the porous medium, µ f the fluid viscosity, and g the gravity vector.Substituting Eq. (3) into Eq.(1) yields The equations of mass conservations can be rewritten by applying the concept of the Lagrangian (total) derivative with respect to a moving solid, e.g.  4) can be rewritten as  ) is well captured by the modelling results which take into account a nonlinear evolution of the elastic moduli resulting from closure of pre-existing micro-cracks and fissures enhanced by compaction of the porous space followed by stiffening of the rock (increase Hertzian contact surface between solid grains).
The modelling results can also capture the compaction-induced reduction of porosity when compared to experimental data (in red) within experimental errors (red shaded area).More details on the theory and results can be found in Jacquey et al. (2015).
In a similar way, it is possible to rework the solid mass balance equation (Eq.2) to obtain From Eq. ( 6), it can be noticed that even by considering both the solid skeleton and the pore fluid to be incompressible, the porous rock material will deform (contract or dilate) when fluid is expelled from or injected into the pore space.
Equation ( 6) can be used to express the evolution of the porosity in terms of the Lagrangian derivative with respect to the solid deformation velocity as Substituting Eq. ( 7) into Eq.( 5), one obtains The first term in the left-hand side of Eq. ( 8) can be expressed in terms of the fluid pore pressure and temperature by thermodynamic differentiation as where 1 is the inverse of the fluid bulk modulus and the fluid volumetric thermal expansion coefficient.
The second term in the left-hand side of Eq. ( 8) can also be cast in terms of the problem variables, e.g.pore pressure, temperature, and solid skeleton displacements, by defining a proper constitutive mechanical model.The linear momentum balance equation of the mixture in terms of the effective Cauchy stress tensor σ (x, t) takes the form where 1 is the rank-two identity tensor, ρ b is the bulk density of the fluid-solid mixture (ρ b = nρ f + (1 − n)ρ s ), and α = 1 − K K s is the Biot coefficient, with K being the drain bulk modulus and K s the bulk modulus of the solid grains.The geometrical compatibility condition gives the following strain-displacement relation: Deformation of the solid skeleton is described in terms of thermo-poroelastic response (Biot's consolidation theory) and dissipative plastic behaviour.To simplify the presentation of the constitutive mechanical model, in the following we will consider only small strain conditions, but the theory has also been extended to account for finite deformation.In addition, due to strain history dependence, we detail the formulation in incremental form.Following Biot's theory, (effective) stresses are related to elastic strains via the following relationship: where C = C ij kl = λδ ij δ kl +2Gδ ik δ j l is the rank-four elastic stiffness tensor, with λ and G being the first (volumetric) and second (shear) Lamé moduli respectively.The stress-strain constitutive relation given by Eq. ( 12), can then be used to find an expression for the material derivative of the solid density in terms of the problem variables (second term in Eq. 8) as where σ indicates the mean effective stress and β s the volumetric thermal expansion coefficient of the solid grains.Substituting Eqs. ( 9) and ( 13) into Eq.( 8), we obtain where we have expressed the gradient of the solid deformation velocity in terms of the volumetric component of the total stain rate tensor, By making use of the definition of the total derivative and given the stress-strain constitutive equation (Eq.12), Eq. ( 14) can be recast as where 1 It is possible to rewrite Eq. ( 15) in terms of the solid velocity only, by integrating the momentum balance equation (Eq. 3) as where the last two terms can be considered as second-order correction terms taking into account nonlinear advective effects.
To quantify any permanent (irreversible) deformation of the material due to inelastic processes we make use of the concept of eigenstrain ( * ) as derived from micromechanics (Mura, 1987).By assuming small strain approximation, the total strain of the material can be decomposed as the sum of the elastic strain ( e ) and eigenstrain components ( * ) as Therefore, Eq. ( 16) can be written as In the present study we focus on two major kinds of residual deformation, that is, thermal expansion and plastic flow, i.e. * = * ij = * T ij + * p ij = * T + * p , though additional processes including for example swelling, fatigue, or phase transformations can be relatively easily integrated in the current formulation.
Thermal strains are related to deformation induced by temperature changes inside the material, and can be therefore expressed by where Ṫ is the relative temperature rate.We determine the plastic component of the strain tensor (Eq.17) by making use of the normality rule as where γ = γ (σ , κ) is the plastic multiplier satisfying the classical Kuhn-Tucker conditions ( γ ≥ 0, F ≤ 0, γ F = 0) with F(σ , κ) being the yield surface.Q = Q(σ , κ) is the plastic potential function giving the direction of the plastic strain increment, and κ represents the vector of internal variables affecting the evolution of the yield surface during loading and unloading of the material.Therefore, Eq. ( 18) can be finally written for thermal and plastic eigenstrains as where is the volumetric thermal expansion coefficient of the pores and the last term can be considered as being of second order.The balance of energy for the fluid-rock mixture, assuming local thermal equilibrium between the two phases, and neglecting the dissipation of mechanical energy due to deformation of the solid phase reads as where s is the bulk thermal conductivity, and Ḣ is a rate of energy production.The first term of the left-hand side of Eq. ( 23) takes into account secondary, non-Boussinesq dissipative effects related to the pressure temperature dependency of the bulk storage, which are usually rewritten only considering variable fluid and solid density as The use of a conservative finite-element formulation in Eq. ( 23) ensures the balance of the fluid enthalpy (ρ f c f T ) both at element and node levels.Therefore, it guarantees that we avoid accumulating in time unbalances in intercell fluxes as is usually the case when relying on non-conservative formulations for convective-type problems (Giudice et al., 1992).
The energy conservation equation as given by Eq. ( 23) is valid for a porous medium in the absence of any thermoelastic coupling.Following Biot consolidation theory, it is possible to consider the effects of the solid elastic deformation on the temperature distribution by augmenting Eq. ( 23) with a thermoelastic dissipation rate term, as where T 0 represents the absolute temperature of the porous medium in a stress-free state.Equation ( 25) can be easily modified to take into account additional thermal effects from fluid dilation and shear heating stresses.

Numerical implementation
In order to solve the system of coupled and nonlinear equations as described above (Eqs.10, 21, and 25), a number of interconnected issues must be taken into account.These include for example the choice of the spatial discretization adopted, the time stepping scheme and temporal integration, iterative solvers, and preconditioners.In what follows, we describe the methods adopted in GOLEM in order to tackle these issues.In the next section, we also present some numerical applications that serve as code benchmarking and illustrate the capability of the simulator to solve for problems of different degree of difficulty at optimal cost.
GOLEM is an open-source simulator specifically developed for thermo-hydro-mechanical coupled applications in fractured geological systems, supporting 1-D, 2-D, and 3-D computations in a single code implementation.It builds on the object-oriented numerical framework MOOSE developed at the Idaho National Laboratory (Gaston et al., 2009).MOOSE provides a flexible, massive parallel (including both MPI and multithreading) platform to solve for multiphysics and multicomponent problems in an implicit manner.It relies on state-of-the-art and extensively tested libraries developed both at universities and national laboratories, as the libMesh library (Kirk et al., 2006) developed at the University of Austin in Texas for capabilities related to parallel finite-element method, and a suite of scalable parallel nonlinear and linear solvers (PETSc: Balay et al., 2016;Trilinos: Heroux et al., 2005;Hypre: Chow et al., 1998).The use of different open-source libraries provides a modular structure to the framework, which allows the developer of an application to only concentrate on the high-level description of the multiphysics problem and to maintain the code relatively compact.Following the basic structure of the MOOSE framework, GOLEM has also been developed as a modular application, where each module is responsible for the solution of a specific physical process.This aspect makes easy any further modification, adjustment, and improvement of the program with limited efforts from the user's side.At the present stage, Golem is hosted and will be made available for download from a GitHub repository (https://github.com/ajacquey/Golem)and it comes together with a suite of relatively simple benchmark problems and all input files for running the applications as presented in the paper.

Variational formulation and its numerical solution
The finite-element discretization is based on the weak form, in an integral sense, of the system of partial differential equations as derived in the previous paragraph.For this purpose, we consider the porous matrix to be described by a close domain of volume ⊂ n bounded by a boundary ⊂ (n−1) .Given the length/width ratio typical of fractures, a discrete fracture is represented by a lower-dimensional element of volume f ⊂ (n−1) and surface area f ⊂ (n−2) .The corresponding weak form of the governing equations are then derived by applying the method of the weighted residuals as 26) where in Eqs. ( 26), ( 27), and (28) we have omitted all supra-(sub)scripting for easiness in the notation.The equations derived above describe an initial and boundary value problem, for which proper boundary and initial conditions need to be assigned.These can be set by prescribing either the value of the problem variables along or their flux across a portion of the boundary.More precisely, the model discussed in the previous paragraph must satisfy the following set of boundary conditions: -Prescribed displacement, pore pressure and temperature, i.e. u = u, p f = p f , T = T on u , p f , T respectively.
-Equilibrium of boundary stresses and external loads (last two terms in Eq. 28).
-Continuity of the fluid flux across the imposed boundary, i.e. q D •n q H = q H on ∂ q H , where q H is the rate of in/out flow per unit of area across the boundary surface.
-Continuity of the total (diffusive plus advective) heat flow, i.e. ((ρc) f q D T −λ b ∇T )•n q T = q T on q T , where q T is the rate of in/out heat flow per unit of area across the boundary surface.
The resulting system of equations together with a proper choice of boundary and initial conditions is discretized spatially by the finite-element method, while the temporal discretization is done by traditional finite difference techniques.The nodal values of the primary variables of the problem, temperature (T n ), pore pressure (p n f ), and the deformation vector of the solid skeleton (u n ) are approximated by linear Lagrangian interpolation polynomial functions as Though higher-order polynomials are also available through the libMesh library, we rely on linear finite-element approximation for all variables.Indeed, we have found that, for the kinds of problems as those that will be presented in the next paragraph, higher-order approximations do not necessarily guarantee a better convergence of the solution, being subjected at the same time to severe under-or overshooting numerical effects in the presence of sharp gradients.
Single-noded, zero-thickness finite elements are used to explicitly represent individual fractures, the latter assumed to be clean and fully saturated (i.e.n = 1 and M b = K f ).It follows that fractures can be parameterized as based on the concept of their effective aperture, which provides a quantitative measure of the geometric width of the fracture surface.The use of single-noded finite elements imposes continuity of gradients in both pore pressure and temperature across the fracture width, as well as the absence of any shear and normal strain acting on the fracture plane.While the latter assumption simplifies the problem formulation, it prevents us from including the effects of local fracture mechanics in the current formulation.Therefore, fractures are only considered as having a distinct hydraulic and thermal behaviour with respect to the porous domain.To extend the current formulation to include mechanical feedback effects from localized deformation, nucleation, and fracture propagation is part of ongoing activities.
The mass balance equation for a discrete fracture then reads as where ± H are leakage terms (weak form) across each of the two sides of the fracture surface into the surrounding porous domain and b is the aperture of the fracture.
In a similar manner, the energy balance equation for a fracture element reads as where ± T quantify the amount of heat leaking from the fracture surfaces into the porous medium domain.
In Eqs. ( 30) and (31), integration is done in local coordinates of the fracture tangential and normal directions respectively (x , y , z ).This enables us to superimpose the discretized conservative equations for the porous medium and the fractures at the nodal location shared by the two elements, where the fluid mass and heat fluxes also cancelled out.Therefore the weak form of the conservation equations can be simplified as where we have made no distinction between the effective hydraulic aperture and the mechanical aperture Derivatives of the test functions and of directiondependent material properties with respect to the system of global coordinates are computed by standard coordinate transformation, i.e.
, where J ij is the Jacobian with cos x i , x i being the directional cosines.Coordinate transformation is applied to all direction-dependent (i.e.tensorial) material properties as well as to the directional derivatives as where N and N are direction-dependent material properties in global and local coordinates respectively, and I f is the unit tangent vector of local coordinates.
The above system of coupled equation can be rewritten in a more concise form as where S is the nodal coefficient storage matrix, K denotes the stiffness matrix of the problem, F is the load vector, x is the vector of the problems variables, and R( û) is the residual from the discrete approximation.We make use of an unconditionally stable, backward Euler method to integrate Eq. ( 35) in time, and to arrive at the monolithic form of the coupled system.The monolithic form of Eq. ( 35) is then solved iteratively by the classical Newton method, where different linear solvers (e.g.GMRES and its flexible variant, FGMRES) and preconditioners (including also Newton-Krylov-Schwartz domain decomposition techniques) are available from the open source algebraic libraries.

Nonlinearities and their stabilization
The coefficient matrices in Eq. ( 35) contain off-diagonal components, thus giving rise to a highly nonlinear problem to be solved for.Nonlinearities arise due the nonlinear constitutive relationship between problem variables and material properties, as well as to the presence of internal feedback effects among the different processes.The constitutive relation linking structure, properties, and primary variables considered include equation of states for fluid density, thermal compressibility and expansivity, and fluid viscosity with respect to the p f − T state of the system (IAPWS, 2008a, b), as well as porosity and permeability relations as a function of variation in the thermal, mechanical, and hydraulic state of the system; see Figs. 3 and 2 respectively.For problems involving heat transport, if the convective term dominates over the diffusive term, the solution of the energy transport equation typically contains interior and boundary layers and its solution by the Galerkin finite-element technique is usually globally unstable (Nield and Bejan, 2012).Numerical instabilities which take the form of spurious oscillations around sharp gradients and boundary layers and likely lead to nonconvergent solutions have been observed to occur independent of the level of mesh refinement (Diersch and Kolditz, 1998).In such situations, in order to enhance the stability and accuracy of the finite-element solution, some sort of stabilization has to be added to the discrete formulation.In GOLEM, we have opted for the streamline upwind Petro-Galerkin (SUPG) method, which can be conceived as adding an upwind perturbation along computed streamlines to the standard Galerkin formulation (Brooks and Hughes, 1982).The SUPG discretization of Eq. ( 27) is obtained by making use of an "enriched" weighting function, ω * , as where ω is the continuous weighting function and p = τ ∇ω • ∇q D q D is a discontinuous streamline upwind correction.Note that by its definition, ω * is no longer continuous across inter-element boundaries.A critical question is related to the choice of the upwind parameter τ , which might influence the stability and accuracy of the discrete solution.In all simulations presented in the paper, when needed, we have made used of the formula as presented in Galeão et al. ( 2004): where h is the diameter of the finite element adopted along the direction of q D , p p is the order of approximation considered in the interpolation, • is the Euclidean norm, and P e = q D h 2kp is the local Peclet number.A major disadvantage of the SUPG method is to add numerical diffusion to all elements characterized by a high Peclet number, even in the absence of sharp thermal gradients.We leave as a subject of future studies the implementation of more recent and more sophisticated damping methods, e.g.entropy viscosity method as described in Guermond et al. (2011).

Plasticity and return-map algorithm
We consider plastic deformation in terms of a Drucker-Prager frictional plastic yield criterion (Drucker, 1950) in which the onset of yield is a function of the first and second invariants of the effective stress tensor J 1 = σ kk , J 2 = www.solid-earth.net/8/921/2017/Solid Earth, 8, 921-941, 2017 1 2 (τ ij τ ij ) respectively, where τ ij is the deviatoric part of the effective stress tensor (τ ij = σ ij − 1 3 J 1 δ ij ): where ϕ and C are the Mohr-Coulomb friction angle and cohesion respectively and 0 is a small non-hardening parameter here introduced to relax the singularity at the cone's tip of the Drucker-Prager yield envelope.An important aspect relates to the form of the plastic potential function adopted.It has been shown that associated flow rules might lead to overestimating the dilation of the rocks generally resulting in an overly weaker response of the rock to loading (Vermeer and de Borst, 1984;Jiang and Xie, 2011).We avoid these issues by implementing a non-associated form of Drucker-Prager, in which the plastic potential function, G, is considered to depend on the dilation angle, ψ as and, in consequence, the unnormalized flow directions can be derived as Note that in our formulation we take also into account possible degradation of the strength of the rock subjected to loading in terms of hardening (and softening) of the internal parameters (friction angle, cohesion and dilation angle) as a function of the accumulated plastic strain (internal variable κ).
The stress update procedure is conducted via a return-map algorithm based on the closest point projection on the yield surface (Simo and Hughes, 1998) within a Newton-Raphson procedure.This algorithm is presented in incremental form in the following.Subscript n refers to a value at the time t n , that is, σ n = σ (t n ) and superscript (k) refers to the kth iteration in the Newton-Raphson procedure.We use the following notation for sake of simplicity: , then the increment of plastic multiplier is also positive according to the Kuhn-Tucker conditions, γ > 0. Define the system of equations with the residuals to minimize, the plastic flow residual R ,n+1 , the internal parameter residuals R κ,n+1 and the yield condition for this time step F n+1 as where H is the hardening flow rule for the internal parameter κ, that is, κ = − γ H.
2. This system of equations is then linearized as follows: , which can be expressed as the following matrix system: where J is the Jacobian matrix x is the vector of unknowns to compute, , and R is the residuals vector: 3. The aforementioned matrix system of equation is solved at each kth iteration using routines from the PETSc library (Balay et al., 2016) for the increment of stress, internal parameter, and plastic multiplier ( σ 4. The variables are updated at the end of the k th iteration:

5.
Steps 1 to 4 are repeated until the residuals reached minimum threshold values.
We have tested and validated the above-described returnmap algorithm to update the elasto-plastic deformation against available algorithms in the MOOSE tensor mechanics module.

Results
In this section we present five different applications of the numerical simulator.These applications are intended to test the ability of the simulator to deal both with single processes and their coupling.By starting with simplistic benchmarks, for which analytical solutions exist, we gradually increase the complexity of the problem formulation in order to demonstrate the applicability of the approach to realistic operational cases.

Thermal and hydraulic processes -heat transport in a fracture
The first application deals with groundwater flow and (diffusive and advective) heat transport in a fracture, i.e. a TH application.We assume a fracture which is fully saturated with water (n = 1), having homogeneous and isotropic properties.Groundwater flow in the fracture is assumed to be unidirectional and the average velocity is considered constant throughout the length of the flow field.The initial temperature is set to zero.At time t = 0, a sudden increase in temperature is imposed along the inlet boundary of the medium (T = T 0 ).The problem can be formulated as a semiinfinite medium with a point source at the inlet boundary (T = T 0 H (x = 0, t), with H (t) being the Heaviside function).Under these assumptions, Ogata and Banks (1961) gave an analytical solution for the variation of the temperature as a function of the position and time as For the numerical simulation, we consider a fracture of length L = 100 m, which has been discretized into 1000 line elements of equal length and subjected to an imposed pressure gradient of ∇p f = 3 Pa m −1 .All material properties are listed in Table 1.
Given the parameters considered, a constant fluid velocity (v x = 3.0 × 10 −7 m s −1 ) is obtained.Initial conditions for the pressure and temperature are p f = 0.1 MPa and T = 0 • C respectively.In this benchmark, we compare the evolution of the temperature field (normalized by the inlet temperature value) at four different points along the length of the fracture with the analytical solution as given by Eq. ( 41). Figure 4 shows the comparison between model results (red curves) and the corresponding analytical solutions (black circles) during the entire simulation time.There is a general good fit in terms of the temporal evolution of the advected front at all locations along the fracture plane, with the modelled thermal front moving slightly faster as visible from the earlier stepping at the observation points.Based on the ob- tained results, we can conclude that the TH numerical implementation is accurate for practical applications.

Hydraulic processes in the presence of discrete fractures -flow in a fractured porous medium
The purpose of this application is to test the numerical implementation of the formulation in the presence of discrete fractures.At this purpose, we refer to a relative common benchmark dealing with uniform, steady-state flow in a porous medium locally disturbed by the presence of a fracture (Strack, 1982;Watanabe, 2011).The original problem formulation considers a semi-infinite two-dimensional horizontal section with an embedded fracture located symmetrically at its centre.Uniform flow is maintained by imposing a specific discharge (q 0 ) from the left boundary inside the domain, the value of which is kept constant and equal to q 0 = 1.0 × 10 −4 m s −1 .The fracture is considered to extend infinitely along the direction normal to the plane, while being of finite along-plane length (L), with its middle point located exactly at the centre of the domain.It has a width which can be varied along its length (though it is assumed to remain constant in the following), and it is inclined with respect to the model boundary by a constant angle (α).
We have extended the original formulation to a threedimensional case; see Fig. 5 for the model geometry and boundary conditions.The setup of the problem comprises a three-dimensional, quadrilateral box (10 × 10 × 1 m in the x, y, and z direction respectively) representing the porous medium.The fracture is implemented as a two-dimensional surface cutting entirely the model along the vertical and having a finite horizontal length of L = 2 m.It is inclined by an angle α = 45 • with respect to the model boundaries.In order to maintain a constant discharge from the left to the right of the model, constant pressure boundary conditions are imposed along the same boundaries thus resulting in a constant pressure gradient along the x axis.The value of the imposed pressure gradient ( p = 1 MPa) has been enforced so to match the value of specific discharge of the original problem as derived by Strack (1982), thus permitting a direct comparison between the analytical and numerical solution.No flow conditions are imposed along the other boundaries.We assume a laminar flow in the fracture plane, and pressure  variations across its width are neglected.All relevant parameters are summarized in Table 2. Strack (1982) provided an analytical solution for the one dimensional version of the problem derived in terms of a complex potential flow as where q 0 is the magnitude of the flow occurring within the model domain (assumed uniform), is the fracture-related dimensionless variable (z 1 and z 2 being the endpoints of the fracture), and The numerical solution of the three-dimensional problem has been obtained by solving a steady-state flow problem within the matrix-fracture domain.We plot the computed pressure distribution together with the outline of the geometry of the fracture in Fig. 6.The presence of the discrete fracture disturbed the uniform horizontal flow in a close domain around its location.There, the isolines of constant pressure are distorted, resulting in a faster flow preferentially oriented parallel to the fracture plane.In order to test the reliability of the numerical solution, we compare the computed pore pressure extracted along a diagonal cross section cutting the model domain along its bottom-left/top-right corners with the pore pressure derived from the analytical solution.The comparison shows a perfect fit between the two results, thus proving the applicability of the discrete fracture approach as implemented in the current formulation.
The results obtained from a variation of the abovedescribed problem by considering two self-intersecting fractures embedded in a two-layered matrix are presented in the Supplement.

Mechanical processes -3-D elastoplastic oedometer test
In the following test, we consider a cube of porous medium with edges of 1 m.The cube is subjected to axial loading with constant solid velocity v x under the conditions of an oedometer test (see setup in Fig. 7).The porous material undergoes continuous loading, and it behaves elastically until the strength of the material is reached.From this time on, the material undergoes plastic loading.The elastoplastic constitutive laws adopted for this simple problem formulation is the Drucker-Prager plasticity model.The Drucker-Prager yield envelope is a smoother version of the classical Mohr-Coulomb failure criterion.Under these conditions, an analytical solution for the stress state of the porous material can be derived as described in Appendix A and serves here as verification of the numerical implementation of the elastoplastic constitutive laws.The physical properties used for this benchmark are summarized in Table 3.
Figure 7 shows the evolution of stress for an associative (red dots) and two non-associative (dilation angle different from friction angle) plastic potentials (blue and green dots).The results exhibit a perfect agreement between the analytical solution and the numerical prediction.In an attempt to better quantify the internal coupling among the elastic and inelastic component of the deformation, we have modified the original elasto-plastic oedometer test to include in addition pore pressure coupling, therefore solving for a coupled hydraulic-mechanical (poroelasticity plus plasticity) problem.The aim of this test is to further validate the reliability of the implementation of the plastic behaviour with respect to its coupling with the poroleastic response to external loading of the porous medium.The geometry of the problem has been kept the same as the one previously described, as well as the solid material parameters (see Table 3).The fluid parameters adopted comprise a constant fluid density (ρf = 1000 kg m −3 ), fluid viscosity (µf = 1.0 × 10 −3 Pa s), porosity (n = 0.1) and permeability (k = 1.0 × 10 −15 m 2 ) and a Biot coefficient of α = 0.6.Pore pressure is initialized to zero at the beginning of the simulation inside the model domain and all lateral boundaries are considered to be close to pressure, thus simulating typical undrained experimental conditions.The boundary conditions in terms of mechanical displacement have been considered as for the reference case.Inelastic behaviour is introduced in the form of a nonassociated plasticity with a friction angle of φ = 20 • and a dilation angle of ψ = 10 • .Given the setup of the model, both the stress and the pore pressure increase with time as the loading from the right lateral boundary continues.Indeed, the pore pressure evolution is affected by variations in the total volumetric strain, better quantified in terms of its rates (see Eq. 27).The onset of plastic behaviour, indicated by a change in the slope of the stress profile, affects the pore pressure distribution and evolution, with rates of pore pressure buildup decreasing due to accumulation of irreversible volumetric strains, thus validating the poroelastic and plastic coupling as implemented in the simulator (Fig. 8).

Reservoir applications
In the following two subsections, we will describe two synthetic examples that deal with applications as considered of interest in the context of operational geothermal activities.
The examples presented do not address a specific field application, though they do bear similarities and links to real cases.An application to an actual field study case based on an injection test performed at the Groß Schönebeck geothermal site (north Germany) is the subject of a separate publication (Jacquey et al., 2017).

Thermal and hydraulic processes -prototype of multi-fractured geothermal reservoir
In this example, we present a setup inspired by a typical geothermal reservoir application.The model aims at simulating the thermal and hydraulic configuration of the reservoir during operational activities (injection and production) spanning a life time of the reservoir of approximately 100 years.
The model consists of four different geological formations, two units representing the target reservoir plus an upper and lower formation acting as cap rocks.The extent of the model domain is 10 × 10 × 3 km in the x, y, and z directions.The target reservoir is located at a depth of approximately 4.6 km below sea level, and is cut by a natural fault showing a slip of some hundreds metres at depths of relevance in the reservoir.A doublet system is integrated in the model, consisting of an injection and a production well.The open-hole section of the two wells is kept parallel and extends for approximately 1.5 km horizontally in the reservoir.The open-hole section of the wells has been integrated as one-dimensional finite elements and homogenization of the resulting governing equations is done by considering the surface area of the well bore as scaling parameter.A system of hydraulically stimulated fractures is also considered to enhance the hydraulic connection between the two wells along their horizontal sections.The multifrac system is intended to represent a hydraulic stimulation campaign prior to reservoir exploitation.There are a total of 10 fractures, equally spaced every 100 m and all sharing the same geometry and material properties.A schematic representation of the main geometry is illustrated in Fig. 9, and all properties are listed in Table 4.The reservoir boundary conditions adopted for the simulation comprise a constant temperature fixed along the top of the reservoir (T = 137.5 • , corresponding to a background thermal gradient of approximately 30 • per kilometre), heat input at the base of the model (q T = 72 mW m −2 ), and hydrostatic pressure along all lateral borders of the domain.The natural system of the reservoir, before stimulation, is de-scribed by running a steady-state simulation (given the same boundary conditions as described above) without considering any operational activity.This was done in order to equilibrate the fluid density and viscosity, here considered functional of the primary variables according to IAPWS (2008a, b), thus avoiding spurious oscillation at the beginning of the transient simulation Operational activity is simulated by injecting water at a fixed temperature, T in = 55 • C, and at a constant rate q in = 30 L s −1 for the entire simulation period.Production rates are kept equal to injection rates through the operation.Figure 9 illustrates the reservoir state at the end of the simulation (an animation of the full evolution of the system is provided in the Supplement).Figure 9 nicely illustrates the evolution of the reservoir temperature as resulting from the dynamics of interactions between the different components of the system (i.e.reservoir matrix, fractures, and geothermal wells) and shows how the simulator can be effectively applied to 3-D modelling of complex reservoirs.The depth of the target reservoir is approximatively 4 km below sea level.Pore pressure and temperature distributions are assumed to be homogeneous at the beginning of the simulation, and equal to 40 MPa and 150 • C respectively.A regional stress field is applied as background stress, to simulate a normal faulting regime with the following magnitudes: -Vertical stress S 1 = 100 MPa in the z direction.
-Maximum horizontal stress S 2 = 90 MPa in the y direction.
-Minimum horizontal stress S 3 = 50 MPa in the x direction.
Two hydraulic fractures are considered for this doublet system and represent the impacts of a prior hydraulic stimulation campaign to enhance the productivity of the reservoir.Given the in situ stress conditions, these two hydraulic fractures are orthogonal to the minimum horizontal stress.They are implemented as squares with edges of 100 m.The distance between the fractures is 200 m.
In this simulation, we consider additional nonlinear effects as related to imposed variations in the evolution of the fluid (i.e.fluid density and viscosity) and rock properties (i.e.porosity and permeability) as a function of the evolution in the state of the reservoir during operational activ-      ities.Changes in porosity are controlled by Eq. ( 22), where we neglected second-order terms for this specific application, while the evolution in the rock permeability is governed by a classical Kozeny-Carman-like relation as where the coefficient A includes information about the pore and grain geometries and can be expressed via the initial value of porosity and permeability: A = k 0 (1−n 0 ) 2 n 3 0 .
Figure 10 illustrates the model geometry.
Figure 10 shows the distribution of pore pressure and temperature at the end of the simulation.It can be observed that colder temperature than the reservoir temperature are produced at the right well after 50 years of operations.
Figure 11 shows horizontal slices at the centre of the domain illustrating the changes in porosity and permeability due to the applied changes in strain, pore pressure, and temperature.
This example illustrates the simulator capability of solving a full thermo-hydro-mechanical coupled problem as relevant for geothermal reservoir applications.At the same time, the formulation adopted to simulate the evolution in time and space of the system properties enables us to quantify the impact of the mechanical alteration induced by geothermal operations (injection/production of fluid and applied injection temperature) on both fluid and heat flows within the reservoir.

Conclusions
In this paper, we have presented a novel but robust simulator for modelling coupled THM processes within fractured rocks, with specific focus on reservoir applications.The code, GOLEM, relies on an open-source massive parallel finite-element-based numerical framework (MOOSE) to solve for the coupled problem.It makes use of a fully implicit approach to treat the nonlinear coupling among the different processes and their feedback effects on fluid and rock prop-erties, thus providing higher numerical stability in the context of nonlinear problems.Geological heterogeneities, i.e. discrete fractures and fault zones, are taken into account in our formulation.The latter are represented as finite element of lower geometrical dimension, which allows us to model focused fluid and heat flows on fractures and faults planes or well paths.The capability and robustness of the simulator has been illustrated by means of five numerical examples by increasing progressively the coupling and geometrical complexity of the considered problem formulations.
Improving the reliability of predictions made for geothermal operations at the field scale requires a better description of the physical phenomena which can alter the reservoir productivity as well as the sustainability of the geothermal operations.In this respect, the current framework provides a powerful tool to analyse the dynamic behaviour of fractured reservoirs during geothermal operations.
Ongoing activities are towards integration of the details of the mechanical description of the geological discontinuities (faults and fractures) either by means of a discrete (XFEM approach) or continuous (such as phase field) approach.Such a feature will help to better describe the dynamics deformation in heterogeneous rocks, including localization and evolution along fault zones, and will also permit to quantitative integrate feedbacks on the hydraulic and thermal behaviours of such geological structures.The description of such processes would help at forecasting environmental impacts of reservoir operations such as induced seismicity from dynamic reactivation of faults during operational activities.
Code and data availability.The source code and the input files of the five numerical examples presented in this paper, plus a suite of specific benchmark cases, are available for downloading from a GitHub repository (https://github.com/ajacquey/Golem).
for a moving fluid.By expanding the fluid mass conservation equation and noting that

Figure 2 .
Figure 2. Numerical simulations of tri-axial mechanical experiments on a Flechtinger sandstone.(a) Evolution of the mean effective stress with respect to the volumetric strain.(b) Evolution of the normalized porosity with respect to the mean effective stress.The nonlinear behaviour at low effective stress as observed in the experimental data (red line afterBlöcher et al., 2014) is well captured by the modelling results which take into account a nonlinear evolution of the elastic moduli resulting from closure of pre-existing micro-cracks and fissures enhanced by compaction of the porous space followed by stiffening of the rock (increase Hertzian contact surface between solid grains).The modelling results can also capture the compaction-induced reduction of porosity when compared to experimental data (in red) within experimental errors (red shaded area).More details on the theory and results can be found inJacquey et al. (2015).
matrix of the mapping between global (x i ) and local (x i ) coordinates.Transformation in local coordinates for lowerdimensional elements is achieved by computing the rotational matrix, R =    cos x , x cos x , y cos x , y cos y , x cos y , y cos y , y

Figure 3 .
Figure 3. Example of a nonlinear convective flow problem afterElder (1967).A homogeneous, isotropic, and saturated porous medium is heated from below with fluid density and viscosity being a function of temperature and pressure(IAPWS, 2008a, b).Due to the temporal evolution in fluid property gradients, thermal buoyant forces develop thus giving rise to an unstable convective flow regime as illustrated by the temperature isolines.H-mesh adaptive refinement with local error estimate based on the L-2 norm projection of the intercell thermal flux vectors and a two-refinement cycle per time step has been enforced in order to guarantee numerical stability even in the presence of sharp thermal gradients (see close-up view).

Figure 4 .
Figure 4. Comparison of numerical results (red curves) and analytical solutions (black circles) for heat transport in a fracture.The figure shows the temperature evolution at different positions along the fracture as a function of time.

Figure 5 .
Figure 5. Geometry and boundary conditions for the benchmark case of groundwater flow in a fractured porous medium.

Table 2 .
Material properties the example of flow in a fractured porous medium× 10 −3 Pa s

Figure 6 .
Figure 6.Isolines of pressure computed from the 3-D numerical simulation extracted along a horizontal plane cutting the model domain.Comparison between simulated (continuous red curve) and analytical derived (empty black circles) pressure distribution along a line through the model.

Figure 7 .
Figure 7. Problem formulation and results of the oedometer benchmark.Panel (a) shows the problem formulation.Panel (b) illustrates the results for different dilation angles with the stress-displacement curves.

Figure 8 .
Figure 8. Results of the modified oedometer test to account for pore pressure coupling.Panel (a) shows the evolution of pore pressure (diamonds) and stress (circle) considering either a purely elastic case (red colours) or additional plastic contributions (green colours).Panel (b) illustrates the evolution of the difference (elastic minus plastic).

Figure 9 .
Figure 9. Problem formulation and results of the prototype of multifractured geothermal reservoir.Panel (a) shows the geometry and setup of the simulation and panel (b) the distribution of the temperature after approximately 100 years of production.
4.2.2Thermal, hydraulic and mechanical processesthermo-poroelastic response of a synthetic geothermal doublet This example considers operations of a synthetic geothermal doublet within a low permeability geological formation with induced fractures located at the two operational wells.The model aims at describing the thermo-poroelastic response of the reservoir due to geothermal operations.The extent of the model domain is 500 × 500 × 200 m in the x, y, z directions.

Figure 10 .
Figure 10.Problem formulation and results of the thermo-poroelastic response of a geothermal doublet.Panel (a) shows the problem formulation.Panels (b) and (c) illustrates the pore pressure (with fluid velocity) and temperature distributions respectively after 50 years of operations.

Figure 11 .
Figure 11.Horizontal slices of the model illustrating the changes in transport properties.Panel (a) shows the distribution of the changes in porosity and panel (b) the changes in permeability after 50 years of operations.

Table 1 .
Fluid properties for the example of heat transport in a fracture.

Table 3 .
Mechanical properties for the oedometer benchmark.

Table 4 .
Material properties for the multifrac reservoir application.

Table 5 .
Material properties for the thermo-poroelastic response of a geothermal doublet.