Parallel geothermal numerical model with faults and multi-branch wells

To answer the need for an efficient and robust geothermal simulation tool going beyond existing code capabilities in terms of geological and physical complexity, we have started to develop a parallel geothermal simulator based on unstructured meshes. The model takes into account complex geology including fault networks acting as major heat and mass transfer corridors and complex physics coupling the mass and energy conservations to the thermodynamical equilibrium between the gas and liquid phases. The objective of this Cemracs project is to focus on well modeling which is a key missing ingredient in our current simulator in order to perform realistic geothermal studies both in terms of monitoring and in terms of history matching. The well is discretized by a set of edges of the mesh in order to represent efficiently slanted or multi-branch wells on unstructured meshes. The connection with the 3D matrix and the 2D fault network at each node of the well is accounted for using Peaceman’s approach. The nonisothermal flow model inside the well is based on the usual single unknown approach assuming the hydrostatic and thermodynamical equilibrium inside the well. The parallelization of the well model is implemented in such a way that the assembly of the Jacobian at each Newton step and the computation of the pressure drops inside the well can be done locally on each process without MPI communications.


Introduction
Geothermal energy is a carbon-free steady energy source with low environmental impact.In countries with a favorable geological context, high temperature geothermal energy can make a significant contribution to power production.On the French territory, it is already an attractive option in volcanic islands context compared to importing fossil fuel.Today, about 5 percent of yearly electricity consumption of Guadeloupe already comes from geothermal energy and it is essential for achieving energetic and environmental targets, according to which the overseas territories should produce 50 percent of their electricity consumption from renewable resources by 2020 and achieve their self sufficiency in 2030.As for other parts of the world, the geothermal development potential of the Caribbean islands is high and several industrial projects are under development or already underway, in French overseas territories (Guadeloupe, Martinique) as well as in nearby islands (Dominica, Montserrat, St. Kitts & Nevis, St Lucia...) that currently depend mainly on diesel for power generation.
Numerical modeling has become essential in all phases of geothermal operations.It is used in the exploration phases to assess the geothermal potential, validate conceptual hypothesis and help well siting.Field development and resource management need quantitative estimation to prevent resource exhaustion and achieve its sustainable exploitation (production/injection scenarios).Finally, numerical modeling is also helpful in studying exploitation related industrial risks such as the interaction with shallow water levels (drinking water resources, hydrothermal vents or eruption).
There is a need to develop new efficient and robust simulation tools to go beyond existing code capabilities in terms of geological and physical complexity [25,32].In particular such code should be able to deal with fault networks acting as major heat and mass transfer corridors in high energy geothermal reservoirs and also to simulate both under critical and super critical thermodynamical domains.Existing tools such as Tough2 [31], used for more than 25 years in geothermy, are limited to structured meshes and are not able to integrate conductive faults.Moreover, their parallel efficiency is very limited.This has motivated the development of a new geothermal simulator based on unstructured meshes and adapted to parallel distributed architectures with the ability to represent faults as co-dimension 1 surfaces connected to the surrounding matrix domain.The current version of this simulator is described in [41].The objective of this Cemracs project is to bring the development of this simulator to a level where operational use is possible and real geothermal test cases can be considered.In this regard, wells are central features of geothermal exploitation and are the main focus of this work.
The use of lower dimensional rather than equi-dimensional entities to represent fracture or fault networks has been introduced in [3,19,9,23,28] to facilitate the grid generation and to reduce the number of degrees of freedom of the discretized model.The reduction of dimension in the fracture network is obtained from the equi-dimensional model by integration and averaging along the width of each fracture.The resulting so called hybrid-dimensional model couple the 3D model in the matrix with a 2D model in the fracture network taking into account the jump of the normal fluxes as well as additional transmission conditions at the matrix-fracture interfaces.These transmission conditions depend on the mathematical nature of the equi-dimensional model and on additional physical assumptions.They are typically derived for a single phase Darcy flow for which they specify either the continuity of the pressure in the case of fractures acting as drains [3,10] or Robin type conditions in order to take into account the discontinuity of the pressure for fractures acting either as drains or barriers [19,28,4,13].In our case, the faults will be assumed to act as drains both for the Darcy flow and for the thermal conductivity leading us to set the pressure and temperature continuity as transmission conditions at the matrix fracture interfaces.
This article focus on the Vertex Approximate Gradient (VAG) scheme which has been introduced for the discretization of multiphase Darcy flows in [17] and extended to hybrid-dimensional models in [11,10,12,42,13,41].The VAG scheme uses nodal and fracture face unknowns in addition to the cell unknowns which can be eliminated without any fill-in.Thanks to its essentially nodal feature, it leads to a sparse discretization on tetrahedral or mainly tetrahedral meshes.It has the benefit, compared with the CVFE methods of [9,35,30,29], to avoid the mixing of the control volumes at the matrix fracture interfaces, which is a key feature for its coupling with a transport model.As shown in [11] for two phase flow problems, this allows for a coarser mesh size at the matrix fracture interface for a given accuracy.
At the reservoir scale of a few kilometers, the mesh cannot resolve the well boundary with a radius of say 10 cm and the well is modeled as a Dirac source term along the well trajectory.Most well models in reservoir simulations are defined by a set of connected perforations, each perforation belonging to a cell of the mesh [33,34].This type of approach is adapted to cell-centered finite volume discretization.In order to take advantage of unstructured meshes and of the nodal feature of the VAG scheme, it is more convenient in our case to discretize each well as a subset of edges of the mesh.This alternative approach provides an efficient way to represent slanted and multi-branch wells.The fluxes connecting the well with the 3D matrix and the 2D fault network at each node of the well will be computed using Peaceman's approach [14,33,34].It is based on a Two Point Flux Approximation with a transmissibility taking into account the unresolved singularity of the pressure (or temperature) solution in the neighborhood of the well.The non-isothermal flow model inside the well is defined in the spirit of what is conventionally done in oil reservoir simulators [7] using a single implicit unknown for each well corresponding to a reference pressure often called the bottom hole pressure.The pressures along the well will be deduced from the bottom hole pressure assuming that the pressure is hydrostatic inside the well.The temperatures along the well will be computed assuming thermal equilibrium and a stationary flow inside the well.Then, the well equation is obtained by the complementary conditions between a specified well mass flow rate and a specified limit bottom hole pressure.By connecting all the nodes along the well trajectory to the well reference pressure unknown, the well equation introduces an additional connectivity.This difficulty will be accounted for by the definition of ghost and own wells for each process and by extension of the ghost nodes of each process in order to take into account the additional connections induced by the own and ghost wells.This allows to assemble the Jacobian and to compute the well pressure drops locally on each process without the need of MPI communications.
The outline of the remaining of the paper is as follows.In section 1, the hybrid-dimensional model presented in [41] is recalled.Although the implementation has been done for the multi-phase compositional model defined in [41], we focus here on the particular case of a non-isothermal singlecomponent single-phase Darcy flow model in order to simplify the presentation.Section 2 introduces the space and time discretization of the model.The definitions of the multi-branch well data structure and of the well model are detailed in subsection 2.3.Section 3 presents the parallel implementation of the model including the partitioning of the mesh and wells, as well as the parallel assembly of the nonlinear and linear systems to be solved at each time step of the simulation.The solution of the linear systems uses the parallel linear solver library PETSc [8] and is based on the GMRES iterative solver preconditioned by a CPR-AMG preconditioner [27,37].The implementation of the CPR-AMG preconditioner takes into account the well equations in the definition of the pressure block.Two numerical tests are presented in section 4. The first test case is used to validate our model.It considers an isothermal single-phase stationary Darcy flow on a simple geometry with one horizontal fault and one vertical well for which an analytical pressure solution can be obtained.The second test case considers a single-phase non-isothermal transient flow on a complex geometry including three intersecting faults, one slanted injection well and one multi-branch production well.
1 Hybrid-dimensional non-isothermal single-phase Discrete Fracture Model This section recalls, in the particular case of a non-isothermal single-component single-phase Darcy flow model, the hybrid-dimensional model introduced in [41].

Discrete Fracture Network
Let Ω denote a bounded domain of R 3 assumed to be polyhedral.Following [3,19,28,10,12] the fractures are represented as interfaces of codimension 1.Let J be a finite set and let Γ = j∈J Γ j and its interior Γ = Γ \ ∂Γ denote the network of fractures Γ j ⊂ Ω, j ∈ J, such that each Γ j is a planar polygonal simply connected open domain included in a plane of R 3 .The fracture width is denoted by d f and is such that 0 < d f ≤ d f (x) ≤ d f for all x ∈ Γ.We can define, for each fracture j ∈ J, its two sides + and −.For scalar functions on Ω, possibly discontinuous at the interface Γ (typically in H 1 (Ω \ Γ)), we denote by γ ± the trace operators on the side ± of Γ. Continuous scalar functions u at the interface Γ (typically in H 1 (Ω)) are such that γ + u = γ − u and we denote by γ the trace operator on Γ for such functions.At almost every point of the fracture network, we denote by n ± the unit normal vector oriented outward to the side ± of Γ such that n + + n − = 0.For vector fields on Ω, possibly discontinuous at the interface Γ (typically in H div (Ω \ Γ), we denote by γ ± n the normal trace operator on the side ± of Γ oriented w.r.t.n ± .
The gradient operator in the matrix domain Ω \ Γ is denoted by ∇ and the tangential gradient operator on the fracture network is denoted by ∇ τ such that We also denote by div τ the tangential divergence operator on the fracture network, and by dτ (x) the Lebesgue measure on Γ.
We denote by Σ the dimension 1 open set defined by the intersection of the fractures excluding the boundary of the domain Ω, i.e. the interior of {(j,j )∈J×J | j =j } ∂Γ j ∩ ∂Γ j \ ∂Ω.Let γ n ∂Γ j , j ∈ J denote the normal trace operator at the fracture Γ j boundary oriented outward to Γ j .

Non-isothermal single-phase flow model
To focus on the implementation aspects related to well modeling, the physics of the fluid is kept relatively simple and we refer to [41] for a compositional multiphase non-isothermal modeling of reservoir flow.The fluid is monophasic and is described by its thermodynamical variables X = (P, T ) where P is the pressure and T the temperature.We denote by ρ(X) its mass density, by µ(X) its dynamic viscosity, by e(X) its specific internal energy, and by h(X) its specific enthalpy.The rock energy density is denoted by E r (X).
The reduction of dimension in the fractures leading to the hybrid-dimensional model is obtained by integration of the conservation equations along the width of the fractures complemented by transmission conditions at both sides of the matrix fracture interfaces (see [41]).In the following, X m = (P m , T m ) denote the pressure and temperature in the matrix domain Ω \ Γ, and X f = (P f , T f ) are the pressure and temperature in the fractures averaged along the width of the fractures.The width of the fractures is denoted by d f as a function of x ∈ Γ.The permeability tensor is denoted by K m in the matrix domain and is assumed to be constant in the width of the fractures and to have the normal vector n + as principal direction.We denote by K f the tangential permeability tensor in the fractures.The porosity (resp.thermal conductivity of the rock and fluid mixture) is denoted by φ m (resp.λ m ) in the matrix domain.It is assumed to be constant in the width of the fractures and denoted by φ f (resp.λ f ).The gravity acceleration vector is denoted by g.
The set of unknowns of the hybrid-dimensional model is defined by X m in the matrix domain Ω \ Γ, by X f in the fracture network Γ, and by X Σ = (P Σ , T Σ ) at the fracture intersection Σ.The set of equations couple the mass and energy conservation equations in the matrix in the fracture network and at the fracture intersection as well as the Darcy and Fourier laws providing the mass and energy fluxes in the matrix and in the fracture network where The system (1)-( 2)-( 3)-( 4)-( 5) is closed with the transmission conditions at the matrix fracture interface Γ.These conditions state the continuity of the pressure and temperature at the matrix fracture interfaces assuming that the fractures do not act as barrier neither for the Darcy flow nor for the thermal conductivity (see [3,19,28,41]).It is combined with a phase based upwind approximation of the mobilities in the matrix fracture normal fluxes.This corresponds to the usual finite volume two point upwind scheme for the mobilities (see e.g.[7]) applied for our reduced model in the normal direction between the center of the fracture and each side of the fracture.We obtain where for any a ∈ R we have set a + = max(a, 0) and a − = min(a, 0).Note also that the pressure P f (resp.the temperature T f ) is assumed continuous and equal to P Σ (resp.T Σ ) at the fracture intersection Σ, and that homogeneous Neumann boundary conditions are applied for the mass q f and energy q e,f fluxes at the fracture tips ∂Γ \ ∂Ω.

VAG Finite Volume Discretization 2.1 Space and time discretizations
The VAG discretization of hybrid-dimensional two-phase Darcy flows introduced in [11] considers generalized polyhedral meshes of Ω in the spirit of [16].Let M be the set of cells that are disjoint open polyhedral subsets of Ω such that K∈M K = Ω, for all K ∈ M, x K denotes the so-called "center" of the cell K under the assumption that K is star-shaped with respect to x K .The set of faces of the mesh is denoted by F and F K is the set of faces of the cell K ∈ M. The set of edges of the mesh is denoted by E and E σ is the set of edges of the face σ ∈ F. The set of vertices of the mesh is denoted by V and V σ is the set of vertices of the face σ.For each K ∈ M we define The faces are not necessarily planar.It is just assumed that for each face σ ∈ F, there exists a so-called "center" of the face x σ ∈ σ \ e∈Eσ e such that x σ = s∈Vσ β σ,s x s , with s∈Vσ β σ,s = 1, and β σ,s ≥ 0 for all s ∈ V σ ; moreover the face σ is assumed to be defined by the union of the triangles T σ,e defined by the face center x σ and each edge e ∈ E σ .The mesh is also supposed to be conforming w.r.t. the fracture network Γ in the sense that for all j ∈ J there exist the subsets We will denote by F Γ the set of fracture faces and by the set of fracture nodes.This geometrical discretization of Ω and Γ is denoted in the following by D.
In addition, the following notations will be used

VAG fluxes and control volumes
The VAG discretization is introduced in [16] for diffusive problems on heterogeneous anisotropic media.Its extension to the hybrid-dimensional Darcy flow model is proposed in [11] based upon the following vector space of degrees of freedom: The degrees of freedom are exhibited in Figure 2 for a given cell K with one fracture face σ in bold.
The matrix degrees of freedom are defined by the set of cells M and by the set of nodes V \ V Γ excluding the nodes at the matrix fracture interface Γ.The fracture faces F Γ and the fracture nodes V Γ are shared between the matrix and the fractures but the control volumes associated with these degrees of freedom will belong to the fracture network (see Figure 3).The degrees of freedom at the fracture intersection Σ are defined by the set of nodes V Σ ⊂ V Γ located on Σ.
The set of nodes at the Dirichlet boundaries ∂Ω D and ∂Γ D is denoted by V D .
The VAG scheme is a control volume scheme in the sense that it results, for each non Dirichlet degree of freedom in a mass balance equation.The matrix diffusion tensor is assumed to be cellwise constant and the tangential diffusion tensor in the fracture network is assumed to be facewise constant.The two main ingredients are therefore the conservative fluxes and the control volumes.The VAG matrix and fracture fluxes are exhibited in Figure 2.For u D ∈ V D , the matrix fluxes F K,ν (u D ) connect the cell K ∈ M to the degrees of freedom located at the boundary of K, namely ν ∈ Ξ K = V K ∪ (F K ∩ F Γ ).The fracture fluxes F σ,s (u D ) connect each fracture face σ ∈ F Γ to its nodes s ∈ V σ .The expression of the matrix (resp.the fracture) fluxes is linear and local to the cell (resp.fracture face).More precisely, the matrix fluxes are given by with a symmetric positive definite transmissibility matrix T K = (T ν,ν K ) (ν,ν )∈Ξ K ×Ξ K depending only on the cell K geometry (including the choices of x K and of x σ , σ ∈ F K ) and on the cell matrix diffusion tensor.The fracture fluxes are given by with a symmetric positive definite transmissibility matrix T σ = (T s,s σ ) (s,s )∈Vσ×Vσ depending only on the fracture face σ geometry (including the choice of x σ ) and on the fracture face width and tangential diffusion tensor.Let us refer to [11] for a more detailed presentation and for the definition of T K and T σ .
The construction of the control volumes at each degree of freedom is based on partitions of the cells and of the fracture faces.These partitions are respectively denoted, for all K ∈ M, by and, for all σ ∈ F Γ , by It is important to notice that in the usual case of cellwise constant rocktypes in the matrix and facewise constant rocktypes in the fracture network, the implementation of the scheme does not require to build explicitly the geometry of these partitions.In that case, it is sufficient to define the matrix volume fractions constrained to satisfy α K,ν ≥ 0, and s∈V K \(V D ∪V Γ ) α K,s ≤ 1, as well as the fracture volume fractions constrained to satisfy α σ,s ≥ 0, and s∈Vσ\V D α σ,s ≤ 1, where we denote by dτ (x) the 2 dimensional Lebesgue measure on Γ.Let us also set as well as and which correspond to the porous volume distributed to the degrees of freedom excluding the Dirichlet nodes.The rock volume fraction in each control volume As shown in [11], the flexibility in the choice of the control volumes is a crucial asset, compared with usual CVFE approaches and allows to significantly improve the accuracy of the scheme when the permeability field is highly heterogeneous.As exhibited in Figure 3, as opposed to usual CVFE approaches, this flexibility allows to define the control volumes in the fractures with no contribution from the matrix in order to avoid to artificially enlarge the flow path in the fractures.In the following, we will keep the notation F K,s , F K,σ , F σ,s for the VAG Darcy fluxes defined with the cellwise constant matrix permeability K m and the facewise constant fracture width d f and tangential permeability K f .Since the rock properties are fixed, the VAG Darcy fluxes transmissibility matrices T K and T σ are computed only once.
The VAG Fourier fluxes are denoted in the following by G K,s , G K,σ , G σ,s .They are obtained with the isotropic matrix and fracture thermal conductivities averaged in each cell and in each fracture face using the previous time step fluid properties.Hence VAG Fourier fluxes transmissibility matrices need to be recomputed at each time step.

Multi-branch non-isothermal well model
Let W denote the set of wells.Each multi-branch well ω ∈ W is defined by a set of oriented edges of the mesh assumed to define a rooted tree oriented away from the root.This orientation corresponds to the drilling direction of the well.The set of nodes of a well ω ∈ W is denoted by V ω and its root node is denoted by s root ω .A partial ordering is defined on the set of vertices V ω with s < s if and only if the unique path from the root to s passes through s.The set of edges of the well ω is denoted by E ω and for each edge e ∈ E ω we set e = s 1 s 2 with s 1 < s 2 .It is assumed that We focus on the part of the well that is connected to the reservoir through open hole, production liners or perforations.In this section exchanges with the reservoir are dominated by convection and we decided to neglect heat losses in the first instance.The latest shall be taken into account when modeling the wellbore flow up to the surface.It is assumed that the radius r ω of each well ω ∈ W is small compared to the cell sizes in the neighborhood of the well.It results that the Darcy flux between the reservoir and the well at a given well node s ∈ V ω is obtained using the Two Point Flux Approximation V s,ω = WI s,ω (P s − P s,ω ), where P s is the reservoir pressure at node s and P s,ω is the well pressure at node s.The Well Index WI s,ω is typically computed using Peaceman's approach (see [33,34,14]) and takes into account the unresolved singularity of the pressure solution in the neighborhood of the well.Fourier fluxes between the reservoir and the well could also be discretized using such Two Point Flux Approximation but they are assumed to be small compared with thermal convective fluxes and will be neglected in the following well model.The temperature inside the well is denoted by T s,ω at each well node s ∈ V ω .
The mass flow rate between the reservoir and the well ω at a given node s ∈ V ω is defined by the upwind formula: and the energy flow rate by β inj ω and β prod ω are well coefficients that ared used to impose specific well behavior.The general case corresponds to β inj ω = β prod ω = 1.Yet, for an injection well, it will be convenient as explained in subsection 2.3.2, to impose that the mass flow rates q m,s,ω are non positive for all nodes s ∈ V ω corresponding to set β inj ω = 1 and β prod ω = 0. Likewise, for a production well, it will be convenient as explained in subsection 2.3.3, to set β inj ω = 0 and β prod ω = 1 which corresponds to assume that the mass flow rates q m,s,ω are non negative for all nodes s ∈ V ω .These simplifying options currently prevent the modeling of cross flows where injection and production occur in different places of the same well, as it sometimes happen in a closed geothermal well.

Well model
Our conceptual model inside the well assumes that the flow is stationary at the reservoir time scale as well as a perfect mixing and thermal equilibrium.The pressure distribution along the well is also assumed hydrostatic.
For the sake of simplicity, the flow rate between the reservoir and the well is considered concentrated at each node s of the well.Hence the mass flow rate along each edge e ∈ E ω depends only on time.It is denoted by q e and is oriented positively from s 1 to s 2 with e = s 1 s 2 .Since Fourier fluxes are neglected, the specific enthalpy depends as well only on time along the edge e and is denoted by h e .
The set of well unknowns is defined by the well pressure P s,ω and the well temperature T s,ω at each node s ∈ V ω , the mass flow rate q e and specific enthalpy h e at each edge e ∈ E ω , the well total mass flow rate q ω (non negative for production wells) as well as the well specific enthalpy h ω for injection wells.
For each edge e = s 1 s 2 ∈ E ω , the specific enthalpy h e satisfies the equation For all s 1 s 2 = e ∈ E ω , let us set κ e,s = 1 if s = s 2 and κ e,s = −1 if s = s 1 .The well equations account for the mass and energy conservations at each node of the well.Let E s denote the set of edges sharing the node s, then for all s ∈ V ω we obtain the equations e∈Es∩Eω κ e,s q e + q m,s,ω = δ s root ω s q ω , e∈Es∩Eω κ e,s h e q e + q E,s,ω = δ s root ω s where δ stands for the Kronecker symbol.The hydrostatic pressure distribution along the well implies that for each edge e ∈ E ω .To close the system, the well boundary conditions prescribe a limit total mass flow rate qω and a limit bottom hole pressure Pω .Then, complementary constraints are imposed between q ω − qω and P ω − Pω with P ω = P s root ω ,ω .In addition, the well specific enthalpy h ω is also prescribed for an injection well.
In the following sections, we consider the particular cases of an injection well with β inj ω = 1, β prod ω = 0 and of a production well with β inj ω = 0, β prod ω = 0.In that case, using an explicit computation of the hydrostatic pressure drop, the well model can be reduced to a single equation and a single implicit unknown corresponding to the well reference pressure P ω (e.g.[6]).

Injection wells
The injection well model sets β inj ω = 1, β prod ω = 0 and prescribes the well total mass flow rate qω ≤ 0, the well maximum bottom hole pressure Pω and the well specific enthalpy h ω .
Since β inj ω = 1 and β prod ω = 0, the mass flow rates q e are non negative and it results from (10) that h e = h ω for all e ∈ E ω .
To compute the pressures along the well, we first solve numerically the equations (11) using the well reference pressure P n−1 ω = P n−1 s root ω ,ω obtained at the previous time step n − 1 and h e = h ω .It provides the pressure drop ∆P n−1 s,ω = P (z s ) − P n−1 ω at each node s ∈ V ω , from which we deduce the well pressures using the bottom well pressure at the current time step n From the equation h(X n s,ω ) = h(P n s,ω , T n s,ω ) = h ω , the well temperature T n s,ω at each node s ∈ V ω depends only on the implicit unknown P n ω .The mass and energy flow rates at each node s ∈ V ω between the reservoir and the well are defined by ( 7)-( 8) with β inj ω = 1 and β prod ω = 0 and depend only on the implicit unknowns X n s and P n ω : The well equation at the current time step is defined by the following complementary constraints between the prescribed well total mass flow rate and the prescribed maximum bottom hole pressure.It is formulated using the min function.

Production wells
The production well model sets β inj ω = 0, β prod ω = 1 and prescribes the well total mass flow rate qω ≥ 0 and the well minimum bottom hole pressure Pω .
The solution at the previous time step n − 1 provides the pressure drop ∆P n−1 s,ω at each node s ∈ V ω .This computation is detailed below.As for the injection well, we deduce the well pressures using the bottom well pressure at the current time step n P n s,ω = P n ω + ∆P n−1 s,ω .
The mass and energy flow rates at each node s ∈ V ω between the reservoir and the well are defined by ( 7)-( 8) with β inj ω = 0 and β prod ω = 1 and depend only on the implicit unknowns X n s and P n ω : The well equation at the current time step is defined by the following complementary constraints between the prescribed well total mass flow rate and the prescribed minimum bottom hole pressure.It is formulated using the min function.
Let us now detail the computation of the pressure drop at each node s ∈ V ω using the previous time step solution n − 1.We first compute the well temperature T n−1 s,ω at each node s using equations (10).The mass and energy flow rates from the upstream nodes of the node s are given by ).
The temperature T n−1 s,ω inside the well at node s is the solution of the nonlinear system from which we deduce the mass density ρ n−1 s,ω = ρ(X n−1 s,ω ) inside the well at node s.These mass densities and the reference pressure P n−1 ω are then used to compute the hydrostatic pressure drop ∆P n−1 s,ω for each node s ∈ V ω using equations (11).

Discretization of the hybrid-dimensional non-isothermal single-phase flow model
The time integration is based on a fully implicit Euler scheme to avoid severe restrictions on the time steps due to the small volumes and high velocities in the fractures.An upwind scheme is used for the approximation of the mobilities in the mass and energy fluxes, that is to say the same scheme that is already used in the definition of the transmission conditions (6) of the hybrid-dimensional model or also in the computation of the well mass and energy fluxes.At the matrix fracture interfaces, we avoid mixing matrix and fracture rocktypes by choosing appropriate control volumes for σ ∈ F Γ and s ∈ V Γ (see Figure 3).In order to avoid tiny control volumes at the nodes s ∈ V Σ located at the fracture intersection, the volume is distributed to such a node s from all the fracture faces containing the node s.It results that the volumes of the control volumes s ∈ V Σ at the fracture intersection is not smaller than at any other matrix fracture degrees of freedom.This solves the problems reported in [23] and [36] related to the small volumes at the fracture intersections and avoids the Star-Delta transformation used in [23] which is not valid when coupled with a transport model.
For each ν ∈ M ∪ F Γ ∪ V the couple of reservoir pressure and temperature is denoted by X ν = P ν , T ν .We denote by X D , the set of reservoir unknowns and by P W = {P ω , ω ∈ W} the set of bottom hole pressures.
The Darcy fluxes taking into account the gravity term are defined by where G D denotes the vector (g For each Darcy flux, let us define the upwind control volume cv µ,ν such that for the matrix fluxes, and such that for fracture fluxes.Using this upwinding, the mass and energy fluxes are given by In each control volume ν, the mass and energy accumulations are denoted by We can now state the system of discrete equations at each time step n = 1, • • • , N t f which accounts for the mass (i = m) and energy and at each node It is coupled with the well equations for the injection wells and for the production wells and with the Dirichlet boundary conditions for all s ∈ V D , where X s,D = (P s,D , T s,D ) are the imposed pressure and temperature at node s.
Let us denote by R ν the vector R ν,i , i ∈ {m, E} , and let us rewrite the conservation equations ( 13), ( 14), ( 15), ( 16), (17) as well as the Dirichlet boundary conditions in vector form defining the following nonlinear system at each time step n = 1, 2, ..., where the superscript n is dropped to simplify the notations and where the Dirichlet boundary conditions have been included at each Dirichlet node s ∈ V D in order to obtain a system size independent on the boundary conditions.
The nonlinear system R(X D , P W ) = 0 is solved by a non-smoothed Newton-Raphson algorithm [26].
3 Parallel implementation

Mesh decomposition
Let us denote by N p the number of MPI processes.The set of cells M is partitioned into N p subsets M p , p = 1, ..., N p using the library METIS [24].The partitioning of the set of nodes V, of the set of fracture faces F Γ , and of the set of wells W is defined as follows: assuming we have defined a global index of the cells K ∈ M let us denote by K(s), s ∈ V (resp.K(σ), σ ∈ F Γ ) the cell with the smallest global index among those of M s (resp.M σ ).Then we set The overlapping decomposition of M into the sets is chosen in such a way that any compact finite volume scheme such as the VAG scheme can be assembled locally on each process.Hence, as exhibited in Figure 4, M p is defined as the set of cells sharing a node with a cell of M p .The overlapping decompositions of the set of wells, of the set of nodes and of the set of fracture faces for p = 1, • • • , N p follow from this definition: well root The partitioning of the mesh is performed by the master process (process 1), and then, each local mesh is distributed to its process.Therefore, each MPI process contains the local mesh (M p , V p , F p Γ , W p ), p = 1, 2, ..., N p which is splitted into two parts: We now turn to the parallel implementation of the Jacobian system to be solved at each Newton iteration of each time step.

Parallelization of the Jacobian system
The Jacobian of the nonlinear system ( 18) is assembled locally on each process p = 1, ..., N p resulting in the following rectangular linear system In (19), are the right hand side vectors of own node and fracture faces equations, b p c ∈ R (2#M p ) is the right hand side vector of own and ghost cell equations, and b p w ∈ R (#W p ) is the right hand side vector of own well equations.
The matrix J p cc is a non singular block diagonal matrix with 2 × 2 blocks, and the cell unknowns can be easily eliminated without fill-in leading to the following Schur complement system with The linear system ( 20) is built locally on each process p and transferred to the parallel linear solver library PETSc [8].The parallel matrix and the parallel vector in PETSc are stored in a distributed manner, i.e. each process stores its own rows.We construct the following parallel linear system with The linear system (22) is solved using the GMRES iterative solver preconditioned by a CPR-AMG preconditioner introduced in [27,37].This preconditioner combines multiplicatively a parallel algebraic multigrid preconditioner (AMG) [21] for a pressure block of the linear system with a block Jacobi ILU0 preconditioner for the full system.In our case, the columns of the pressure block are defined by the node, the fracture face and the well pressure unknowns, and its lines by the node and the fracture face mass conservation equations as well as the well equations.
The solution of the linear system provides on each process p the solution vector (U p s , U p f , U p w ) of own node, fracture-face and well unknowns.Then, the ghost node unknowns U p ν , ν ∈ (V p \V p ), the ghost fracture face unknowns U p ν , ν ∈ (F p Γ \F p Γ ) and the ghost well unknowns U w ν , ν ∈ (W p \W p ) are recovered by a synchronization step with MPI communications.This synchronization is efficiently implemented using a PETSc matrix vector product where is the vector of own and ghost node, fracture-face and well unknowns on all processes.The matrix S, containing only 0 and 1 entries, is assembled once and for all at the beginning of the simulation.Finally, thanks to (21), the vector of own and ghost cell unknowns U p c is computed locally on each process p.

Numerical results
All the numerical tests have been implemented on the Cicada cluster of the University Nice Sophia Antipolis composed of 72 nodes (16 cores/node, Intel Sandy Bridge E5-2670, 64GB/node).We always fix 1 core per process and 16 processes per node.The communications are handled by OpenMPI 1.8.2 (GCC 4.9) and PETSc 3.5.3.
where q ω is the mass flow rate per unit well length.The total mass flow rate is This solution will be used to test the convergence of our discretization for both an injection and a production well with fixed temperature.For both test cases, the fluid properties are set to µ = 10 −3 Pa.s for the viscosity and to ρ = 10 3 kg.m −3 for the mass density.We consider a uniform Cartesian mesh of size n x × n x × n x of the domain Ω conforming to the fault and to the well.The well indices WI s,ω for s ∈ V ω are computed following Peaceman's methodology [33,34,14] using the analytical solution (24).Since the mesh is uniform it suffices to solve numerically a local 2D problem with four horizontal faces around a given well node.This computation leads to the following solution for the numerical Peaceman indices.Let us set Then, denoting by E ω the set of edges of the well, the well index of a given node s ∈ V ω is given by for a matrix node s ∈ V ω \ V Γ , and by For both test cases, Dirichlet boundary conditions given by the analytical solution are imposed at the lateral boundaries of Ω and at the boundary of Γ. Homogeneous Neumann boundary conditions are imposed at the top and bottom boundaries of Ω.The well boundary condition is set to either a specified bottom hole pressure P ω or a specified total flow rate q ω .Let us first consider the case of an injection well with the well pressure P ω = 2 × 10 7 Pa and the flow rate per unit well lengh set to q ω = 0.1 kg.s −1 .m−1 .The corresponding analytical solution defined by ( 24) with these parameters is shown in Figure 5. Figure 6 shows, for both a specified pressure or flow rate, the convergence of the relative L 2 errors between the analytical solution and the numerical solution both in the matrix domain and in the fault as a function of the mesh size n x = 10, 20, 40, 80.We obtain an order 1 convergence in all cases.Next, we consider the case of a production well with the well pressure P ω = 5 × 10 6 Pa and the well flow rate per unit well lengh set to q ω = 0.1 kg.s −1 .m−1 .Figure 7 shows the analytical solution defined by (24) with these parameters.
We present in Figure 8 the convergence of the relative L 2 errors between the analytical solution and the numerical solution as a function of the mesh size n x = 10, 20, 40, 80, both in the matrix domain and in the fault and for both for a specified pressure or a specified flow rate.We obtain as previously an order 1 convergence in all cases.

Non-isothermal single-phase flow
This test case considers a non-isothermal liquid flow with mass density, viscosity, specific internal energy and enthalpy obtained from [38].The thermal conductivity is fixed to λ = 2 W.m −1 .K −1 and the rock internal energy density is defined by E r (T ) = c r p T with c r p = 16.10 5 J.m −3 .K −1 .The gravity is not considered in this test case.
The simulation domain is defined by Ω = (0, 2000) 3 in meters.The mesh is a 3D tetrahedral mesh conforming to the fault network and to the wells.It was generated using the implicit framework from the 3D mesh generation package from the Computational Geometry Algorithms Library (CGAL [39]).As shown in Figure 9, there is one injection well (red line) and one multi-branch production well (green line).This mesh contains about 4.9 × 10 6 cells, 2.8 × 10 4 fractures faces and 8.0 × 10 5 nodes.The radius of both wells is set to 0.1 meter and the fault width is fixed to d f = 1 meter.The permeabilities are isotropic and set to K m = 10 −14 I m 2 in the matrix domain and to K f = 10 −11 I m 2 in the fault network.The porosities in the matrix domain and in the faults are φ m = 0.1 and φ f = 0.4 respectively.The computation of numerical Peaceman indices would require an analytical solution for the linear diffusion equation, which is not known for such a complex geometry involving faults and multi-branch wells.This solution could also be obtained numerically using a mesh at the scale of the wells but its generation is out of the scope of this test case.Alternatively, we will use for this test case approximate analytical Peaceman type formulas providing a good order of magnitude for the Peaceman indices.the temperature in the matrix domain and in the faults at times t = 4 × 10 4 days and t = t f .The temperature at the root node of the production well as a function of time is shown in Figure 12.
Then, we show in Figures 13 and 14 the pressure in the matrix domain and in the faults at times t = 4 × 10 4 days and t = t f .In addition, the pressures at the root nodes of both wells as a function of time are shown in Figure 15.
Finally, we present in Figure 16 the total computational time in hours for different numbers of MPI processes N p = 4, 8, 16, 32, 64, 128.The scalability behaves as expected for fully implicit time integration and AMG type preconditioners.It is well known that the AMG preconditioner requires a        sufficient number of unknowns per MPI process, say 100000 as typical order of magnitude, to achieve a linear strong scaling.For this mesh size, leading to roughly 8.2 × 10 5 unknowns for the pressure block, the scalability is still not far from linear on up to 64 processes and then degrades more rapidly for N p = 128.

Conclusion
In this paper, the non-isothermal hybrid-dimensional Darcy flow model presented in [41] has been extended to incorporate thermal well models coupled both with the matrix domain and with the fault network.The well data structure is based on a rooted tree defined by a set of edges of the mesh.This allows to represent efficiently both slanted and multi-branch wells taking advantage of the unstructured mesh and of the nodal feature of the VAG discretization.The fluxes connecting the well with the 3D matrix and the 2D fault network at each node of the well are computed using Peaceman's approach, and the well non-isothermal flow model is based on the usual single unknown approach assuming the hydrostatic and thermodynamical equilibrium inside the well.The parallelization of the well model is performed by definition of own and ghost wells for each process and by extension of the ghost nodes in order to account for the additional connectivity induced by the own and ghost well equations.This allows to assemble the Jacobian and to compute the well pressure drops locally on each process without MPI communications.
The model has been validated using a pressure analytical solution on a simple geometry with one horizontal fault and one vertical well.The efficiency of the model, both in terms of ability to account for complex geology and in terms of parallel scalability, is demonstrated on a non-isothermal single-phase flow test case using a tetrahedral mesh with roughly 4.9 × 10 6 cells and including three intersecting faults, one slanted injection well, and one muti-branch production well.This is an important step toward the application of our simulator to real geothermal test cases in a near future.

Figure 2 :
Figure 2: For a cell K and a fracture face σ (in bold), examples of VAG degrees of freedom u K , u s , u σ , u s and VAG fluxes F K,σ , F K,s , F K,s , F σ,s .

Figure 3 :
Figure 3: Example of control volumes at cells, fracture face, and nodes, in the case of two cells K and L splitted by one fracture face σ (the width of the fracture is enlarged in this figure).The control volumes are chosen to avoid mixing fracture and matrix rocktypes.
the vectors of pressure and temperature unknowns at own and ghost nodes s ∈ V p , fracture faces σ ∈ F p Γ , and cells K ∈ M p .The vector U p w ∈ R (#W p ) is the vector of own and ghost well pressures.Likewise, b p s R p , p = 1, 2, ..., N p is a restriction matrix satisfying R p U = U p The matrix J p R p , the vector U p and the vector b p schur b p w are stored in process p.

4. 1
Numerical convergence for an analytical solution with one horizontal fault and a vertical well We consider the domain Ω = (−H, H) 3 , H = 1000 m, with one horizontal fault Γ = {(x, y, z) ∈ Ω | z = 0} of width d f = 0.5 m and one vertical well of radius r ω = 0.1 m and defined by the line {(x, y, z) ∈ Ω | x = y = 0}.Both the matrix and the fault are isotropic and homogeneous with permeability K m = k m I, k m = 10 −14 m 2 in the matrix and with tangential permeability K f = k f I, k f = 10 −11 m 2 in the fault.For such a simple geometry, an analytical solution of the isothermal stationary linear Darcy equation is defined by the radial pressure

Figure 5 :
Figure 5: Analytical solution with the injection well in the matrix domain (left) and in the fault (right).
(a) Clip view of matrix domain.(b) Fractures and wells.

Figure 9 :
Figure 9: Coarse 3D tetrahedral mesh conforming to the fault network and to the wells.There are one injection well (red line) and one production well (green line). exhibit

t f 5 •
10 6 days final simulation time ∆t 4 • 10 4 days time step N

Figure 10 :
Figure 10: Left: temperature in the matrix domain and in the faults at t = 4 × 10 4 days where a clip view on plane {y = 1000} is used in the matrix domain.Right: temperature in the faults at t = 4 × 10 4 days.The wells are drawn with black lines in both figures.

Figure 11 :
Figure 11: Left: temperature in the matrix domain and in the faults at t = 5 × 10 6 days where a clip view on plane {y = 1000} is used in the matrix domain.Right: temperature in the faults at t = 5 × 10 6 days.The wells are drawn with black lines in both figures.

Figure 12 :
Figure 12: Temperature at the root node of the production well as a function of time.

Figure 13 :
Figure 13: Left: pressure in the matrix domain and in the faults at t = 4 × 10 4 days where a clip view on plane {y = 1000} is used in the matrix domain.Right: pressure in the faults at t = 4 × 10 4 days where the wells are drawn with black lines.

Figure 14 :
Figure 14: Left: pressure in the matrix domain and in the faults at t = 5 × 10 6 days where a clip view on plane {y = 1000} is used in the matrix domain.Right: pressure in the faults at t = 5 × 10 6 days where the wells are drawn with black lines.

Figure 15 :
Figure 15: Pressures at the root nodes of both wells as a function of time.

Figure 16 :
Figure 16: Total computational time vs. number of MPI processes.