# Two-phase geothermal model with fracture network and multi-branch wells

Antoine Armandine Les Landes, Daniel Castanon Quiroz, Laurent Jeannin, Simon Lopez, Roland Masson

► To cite this version:

Antoine Armandine Les Landes, Daniel Castanon Quiroz, Laurent Jeannin, Simon Lopez, Roland Masson. Two-phase geothermal model with fracture network and multi-branch wells. SMAI Journal of Computational Mathematics, 2023, 9, pp.121 - 149. 10.5802/smai-jcm.97. hal-04371287

HAL Id: hal-04371287

https://brgm.hal.science/hal-04371287v1

Submitted on 3 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Two-phase geothermal model with
fracture network and multi-branch
wells

ANTOINE ARMANDINE LES LANDES, DANIEL CASTANON QUIROZ, LAURENT JEANNIN, SIMON LOPEZ & ROLAND MASSON

Volume 9 (2023), p. 121-149.

https://doi.org/10.5802/smai-jcm.97

© The authors, 2023.
<figure><img src="image_2.png" /></figure>

Two-phase geothermal model with fracture network and
multi-branch wells

ANTOINE ARMANDINE LES LANDES <fnref n="1" />

DANIEL CASTANON QUIROZ <fnref n="2" />

LAURENT JEANNIN <fnref n="3" />

SIMON LOPEZ <fnref n="4" />

ROLAND MASSON <fnref n="5" />

**Abstract.** This paper focuses on the numerical simulation of geothermal systems in complex geological settings. The physical model is based on two-phase Darcy flows coupling the mass conservation of the water component with the energy conservation and the liquid-vapor thermodynamical equilibrium. The discretization exploits the flexibility of unstructured meshes to model complex geology including conductive faults as well as complex wells. The polytopal and essentially nodal Vertex Approximate Gradient scheme is used for the approximation of the Darcy and Fourier fluxes combined with a Control Volume approach for the transport of mass and energy. Particular attention is paid to the faults which are modelled as two-dimensional interfaces defined as a collection of faces of the mesh and to the flow inside deviated or multi-branch wells defined as a collection of edges of the mesh with rooted tree data structure. By using an explicit pressure drop calculation, the well model reduces to a single equation based on complementarity constraints and only one well unknown, the bottom hole pressure, implicitly coupled to the reservoir unknowns. The coupled systems are solved at each time step using efficient nonlinear and linear solvers on parallel distributed architectures. The convergence of the discrete model is investigated numerically on a simple test case with a Cartesian geometry and a single vertical producer well. Then, the ability of our approach to deal efficiently with realistic test cases is assessed on a high energy faulted geothermal reservoir operated using a doublet of two deviated wells.

**2020 Mathematics Subject Classification.** 65M08, 65Y05, 76S05, 76T10.

**Keywords.** Geothermal systems, thermal well, two-phase Darcy flow, mixed-dimensional model, faults, finite volume scheme, parallel algorithm.

<footnote n="1">¹ BRGM, 3 avenue Claude-Guillemin, BP 36009, 45060 Orléans Cedex 2, France
*E-mail address:* A.ArmandineLesLandes@brgm.fr</footnote>

<footnote n="2">² Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas, Universidad Nacional Autónoma de México, Circuito Escolar s/n, Ciudad Universitaria C.P. 04510 Cd. Mx. (México)
*E-mail address:* daniel.castanon@iimas.unam.mx</footnote>

<footnote n="3">³ STORENGY, 12 rue Raoul Nordling - Djinn - CS 70001 92274 Bois Colombes Cedex, France
*E-mail address:* laurent.jeannin@storengy.com</footnote>

<footnote n="4">⁴ BRGM, 3 avenue Claude-Guillemin, BP 36009, 45060 Orléans Cedex 2, France
*E-mail address:* s.lopez@brgm.fr</footnote>

<footnote n="5">⁵ Université Côte d'Azur, Inria, CNRS, LJAD, UMR 7351 CNRS, team Coffee, Parc Valrose 06108 Nice Cedex 02, France
*E-mail address:* roland.masson@unice.fr.</footnote>
# 1. Introduction

Deep geothermal systems are often located in complex geological settings, including faults or fractures. These geological discontinuities not only control fluid flow and heat transfer, but also provide feed zones for production wells. Modeling the operation of geothermal fields and the exchange of fluids and heat in the rock mass during production requires explicitly taking into account objects of different characteristic sizes such as the reservoir itself, faults and fractures, which have a small thickness compared to the characteristic size of geological formations and wells (whose radius is of the order of a few tens of centimeters).

This paper focuses on the liquid-vapor single water component non-isothermal Darcy flow model based on mass and energy conservation equations coupled with thermodynamical equilibrium and volume balance. The simulation of such geothermal models in complex geology including faults and wells raises several challenges both in terms of geometrical discretization and in terms of numerical algorithms for the simulation of the nonlinear system of Partial Differential Equations. The strategy developed in this work to account for complex geometries is based on tetrahedral meshes which are particulary adapted to discretize complex geological features like faults defined as a collection of faces and slanted or multi-branch wells defined as a collection of edges with rooted tree structure.

In terms of numerical algorithms for the dynamical model, the challenges result from the strong nonlinear couplings induced by the nonlinear thermodynamical properties and by the phase transitions, as well as from the highly contrasted time and space scales resulting from the heterogeneous petrophysical properties and spatial dimensions.

A common way to account for these highly contrasted spatial scales is based on a reduction of dimension both for the fault/fracture and the well models. Following [5, 6, 12, 13, 16, 24, 26, 29, 30, 34, 37, 41] faults/fractures will be represented as co-dimension one manifolds coupled with the surrounding matrix domain leading to the so-called hybrid-dimensional or Discrete Fracture Matrix (DFM) models. This reduction of dimension is obtained by averaging both the equations and unknowns across the fracture width and using appropriate transmission conditions at matrix fracture interfaces. In our case, the faults/fractures will be assumed to be conductive both in terms of permeability and thermal conductivity in such a way that pressure and temperature continuity can be assumed as matrix fracture transmission conditions [5, 12, 41]. This setting has been extended to two-phase Darcy flows in [14, 15] and to multi-phase compositional non-isothermal Darcy flows in [48].

The discretization of liquid-vapor geothermal models is usually performed using a two-point flux finite volume scheme with cell-centered unknowns [40]. Thanks to its monotonicity, this type of scheme has very good robustness properties. On the other hand, it is limited to simple, mostly Cartesian, geometries due to the requirement of orthogonality of the mesh to guarantee the consistency of the discretization. Alternatively, the consistent discretization of hybrid-dimensional Darcy flow models on general polyhedral meshes has been the object of many research works using cell-centered Finite Volume schemes with Multi Point Flux Approximations (MPFA) [3, 4, 27, 42, 46], Mixed (Hybrid) Finite Element methods [5, 29, 34], Hybrid Mimetic Mixed Methods [7, 13, 16, 23], Hybrid High-Order schemes [19, 20] and nodal schemes such as Control Volume Finite Element methods (CVFE) [12, 27, 35, 36, 41] or the Vertex Approximate Gradient (VAG) scheme [13, 14, 16, 17, 18, 48, 49].

As opposed to nodal approaches, cell-centered methods like MPFA or face-based methods like Mixed, Mixed Hybrid or Mimetic discretizations include a large number of degrees of freedom on tetrahedral meshes leading to costly solutions of large nonlinear and linear systems. This is particularly the case when applied to nonlinear strongly coupled models with implicit time integration as is the case for the geothermal model considered here. This motivates our choice of the mixed-dimensional VAG scheme using nodal and fracture face unknowns in addition to the cell unknowns which can be eliminated without any fill-in. Thanks to its essentially nodal nature, it leads to a sparse discretization on tetrahedral meshes. Compared with other nodal approaches such as CVFE methods, the VAG
scheme avoids the mixing of the control volumes at the matrix fracture interfaces, which is a key
feature for its coupling with a transport model [14].

The well will be modelled as a line source defined by a 1D graph with tree structure. It will be
coupled to the 3D matrix domain and to the 2D faults/fractures possibly intersecting the well using
Peaceman's approach. It is a widely used approach in reservoir simulation for which the Darcy or
Fourier fluxes between the reservoir and the well are discretized by a two-point flux approximation
with a transmissivity accounting for the unresolved pressure or temperature singularity. This leads
to the concept of well or Peaceman's index defined at the discrete level and depending on the type
of cell, on the well radius and geometry and on the scheme used for the discretization. Let us refer
to [38] for its introduction in the framework of a two-point cell-centered finite volume scheme on square
cells, to [39] for its extension to non square cells and anisotropic permeability field and to [1, 21, 47]
for extensions to more general well geometries and different discretizations. The coupling with the
faults/fractures is considered in [10]. Let us also refer to [25] for a related approach also based on a
removal of the singularity induced by the well line source but at the continuous level.

The non-isothermal two-phase flow along the well is a simplified version of the drift flux model [33, 45] neglecting transient terms, thermal losses and cross flow in the sense that, all along the well, the well behaves either as a production or an injection well. It results that using an explicit approximation of the mixture density along the well, the well model can be reduced to a single unknown, the so-called bottom hole pressure, implicitly coupled to the reservoir.

The main objective of this work is to show the ability of the VAG scheme combined with flexible
tetrahedral meshes and parallel nonlinear and linear solvers to simulate efficiently on distributed par-
allel architectures realistic high energy geothermal production scenarios on complex geology including
faults, injection and production thermal wells. To this end, starting from our previous works [10, 48],
the thermal well model is discretized on subsets of edges of the mesh with rooted tree data structure,
accounting for complex well geometries such as slanted and multi-branch wells. To allow for time steps
at the reservoir scale, the well bottom hole pressure is implicitly coupled to the geothermal conserva-
tion equations and strategies are developed to solve efficiently on distributed parallel architectures the
nonlinear systems coupling the discrete mass and energy conservation equations, the well equations
and the thermodynamical equilibrium.

The remainder of this paper is organized as follows. Section 2 presents the physical model describing the flow and transport in the matrix domain coupled to the fracture/fault network in the hybrid-dimensional setting. Section 3 presents the mixed-dimensional VAG discretization of this liquid-vapor non-isothermal hybrid-dimensional model. It is based on a Control Volume formulation of the mass and energy conservation equations combining the VAG discretization of the Darcy and Fourier fluxes with a Phase Potential Upwind approximation of the phase mobilities and enthalpies. The notations for the space and time discretizations are recalled in Subsection 3.1 and the definitions of the VAG fluxes and control volumes are recalled in Subsection 3.2. Then, the well modelling, which is the main new development of this work, is addressed in Subsection 3.3 starting with the description of the well geometry as a collection of edges defining a rooted tree data structure. The source terms connecting the well to the reservoir at each well node are based on two-point fluxes with transmissivities defined by Peaceman's indexes. The derivation of the simplified well model is detailed both for production and injection wells starting from the drift flux model. Subsection 3.4 provides a detailed description of the algorithms used to solve on distributed parallel architectures the stiff nonlinear and linear systems coupling the conservation, well and thermodynamical equilibrium equations at each time step of the simulation. Finally, to demonstrate the efficiency of our approach, we present in Section 4 two numerical tests. The first test case checks in Subsection 4.1 the numerical convergence of the model for a vertical production well connected to a homogeneous reservoir on a family of refined Cartesian meshes. The second test case presented in Subsection 4.2 simulates the development plan of a high enthalpy faulted geothermal reservoir with slanted production and injection wells.
## 2. Hybrid-dimensional non-isothermal two-phase Discrete Fracture Model

This section recalls, in the particular case of a non-isothermal single-component two-phase Darcy flow model, the hybrid-dimensional model introduced in [48].

### 2.1. Discrete Fracture Network

Let $\Omega$ denote a bounded domain of $\mathbb{R}^3$ assumed to be polyhedral. Following [5, 13, 16, 24, 34] the fractures are represented as interfaces of codimension 1. Let $J$ be a finite set of indexes and let $\bar{\Gamma} = \cup_{j \in J} \bar{\Gamma}_j$ and its interior $\Gamma = \bar{\Gamma} \setminus \partial \bar{\Gamma}$ denote the network of fractures $\Gamma_j \subset \Omega, j \in J$, such that each $\Gamma_j$ is a planar polygonal simply connected open domain included in a plane of $\mathbb{R}^3$. The fracture width

<figure><img src="image_3.png" /><figcaption>FIGURE 2.1. Example of a 2D domain with 3 intersecting fractures $\Gamma_1, \Gamma_2, \Gamma_3$.</figcaption></figure>

is denoted by $d_f$ and is such that $0 < d_f \le d_f(\mathbf{x}) \le \bar{d}_f$ for all $\mathbf{x} \in \Gamma$. We can define, for each fracture $j \in J$, its two sides + and -. For scalar functions on $\Omega$, possibly discontinuous at the interface $\Gamma$ (typically in $H^1(\Omega \setminus \bar{\Gamma})$), we denote by $\gamma^{\pm}$ the trace operators on the side $\pm$ of $\Gamma$. Continuous scalar functions $u$ at the interface $\Gamma$ (typically in $H^1(\Omega)$) are such that $\gamma^+u = \gamma^-u$ and we denote by $\gamma$ the trace operator on $\Gamma$ for such functions. At almost every point of the fracture network, we denote by $\mathbf{n}^{\pm}$ the unit normal vector oriented outward to the side $\pm$ of $\Gamma$ such that $\mathbf{n}^{\pm} + \mathbf{n}^{-} = 0$. For vector fields on $\Omega$, possibly discontinuous at the interface $\Gamma$ (typically in $H_{\text{div}}(\Omega \setminus \bar{\Gamma})$, we denote by $\gamma_n^{\pm}$ the normal trace operator on the side $\pm$ of $\Gamma$ oriented w.r.t. $\mathbf{n}^{\pm}$.

The gradient operator in the matrix domain $\Omega \setminus \bar{\Gamma}$ is denoted by $\nabla$ and the tangential gradient operator on the fracture network is denoted by $\nabla_{\tau}$ such that

$$ \nabla_{\tau} u = \nabla u - (\nabla u \cdot \mathbf{n}^{+}) \mathbf{n}^{+}. $$

We also denote by $\text{div}_{\tau}$ the tangential divergence operator on the fracture network, and by $d\tau(\mathbf{x})$ the Lebesgue measure on $\Gamma$.

We denote by $\Sigma$ the dimension 1 open set defined by the intersection of the fractures excluding the boundary of the domain $\Omega$, i.e. the interior of $\cup_{\{(j,j') \in J \times J | j \neq j'\}} \partial \Gamma_j \cap \partial \Gamma_{j'} \setminus \partial \Omega$.

For the matrix domain, Dirichlet (subscript $D$) and Neumann (subscript $N$) boundary conditions are imposed on the two-dimensional open sets $\partial \Omega_D$ and $\partial \Omega_N$ respectively where $\partial \Omega_D \cap \partial \Omega_N = \emptyset$, $\partial \Omega = \overline{\partial \Omega_D \cup \partial \Omega_N}$. Similarly for the fracture network, the Dirichlet and Neumann boundary conditions are imposed on the one-dimensional open sets $\partial \Gamma_D$ and $\partial \Gamma_N$ respectively where $\partial \Gamma_D \cap \partial \Gamma_N = \emptyset$, $\partial \Gamma \cap \partial \Omega = \overline{\partial \Gamma_D \cup \partial \Gamma_N}$.
## 2.2. Non-isothermal two-phase flow model

We consider in this work a two-phase liquid-gas, single water component, and non-isothermal Darcy flow model. The liquid (l) and gas (g) phases are described by their pressure p (neglecting capillary effects), temperature T and pore volume fractions or saturations s^α, α ∈ {l, g}. In order to express the model with a unique set of persistent variables independent on the set of present phases, it is convenient to view the model as a two-phase compositional model with a single component. Consequently, let us also introduce the mass fraction c^α of the water component in phase α, equal to 1 for a present phase α but lower than 1 for an absent phase. Following [2, 11], it will be used below to express the thermodynamical equilibrium as complementary constraints.

For each phase α, we denote by ρ^α(p, T) its mass density, by μ^α(p, T) its dynamic viscosity, by e^α(p, T) its specific internal energy, and by h^α(p, T) its specific enthalpy. The rock energy density is denoted by E^r(p, T).

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 [48]). In the following, p_m, T_m, s_m^α, c_m^α denote the pressure, temperature, saturations, and mass fractions in the matrix domain Ω \ Γ, and p_f, T_f, s_f^α, c_f^α are the pressure, temperature, saturations and mass fractions in the fractures averaged along the width of the fractures. The permeability tensor is denoted by K_m in the matrix domain and we denote by K_f the tangential permeability tensor in the fractures (average value along the fracture width assuming that the permeability tensor in the fracture has the normal as principal direction). The porosity (resp. thermal conductivity of the rock and fluid mixture) is denoted by φ_m (resp. λ_m) in the matrix domain and by φ_f (resp. λ_f) along the fracture network (average values along the fracture width). The relative permeability of phase α as a function of the phase saturation is denoted by k_{r,m}^α in the matrix and by k_{r,f}^α in the fracture network. The gravity acceleration vector is denoted by g.

The following set of equations couples the mass, energy and volume balance equations in the matrix

$$ 
\left\{
\begin{array}{l}
\phi_m \partial_t \left( \sum_{\alpha \in \{l,g\}} \rho^\alpha(p_m, T_m) s_m^\alpha c_m^\alpha \right) + \operatorname{div}(\mathbf{q}_m^{h2o}) = 0, \\
\\
\phi_m \partial_t \left( \sum_{\alpha \in \{l,g\}} \rho^\alpha(p_m, T_m) e^\alpha(p_m, T_m) s_m^\alpha c_m^\alpha \right) + (1 - \phi_m) \partial_t E^r(p_m, T_m) + \operatorname{div}(\mathbf{q}_m^e) = 0, \quad (2.1) \\
\\
\sum_{\alpha \in \{l,g\}} s_m^\alpha = 1,
\end{array}
\right.
$$

in the fracture network

$$ 
\left\{
\begin{array}{l}
d_f \phi_f \partial_t \left( \sum_{\alpha \in \{l,g\}} \rho^\alpha(p_f, T_f) s_f^\alpha c_f^\alpha \right) + \operatorname{div}_\tau (\mathbf{q}_f^{h2o}) - \gamma_n^+ \mathbf{q}_m^{h2o} - \gamma_n^- \mathbf{q}_m^{h2o} = 0, \\
\\
d_f \phi_f \partial_t \left( \sum_{\alpha \in \{l,g\}} \rho^\alpha(p_f, T_f) e^\alpha(p_f, T_f) s_f^\alpha c_f^\alpha \right) + d_f (1 - \phi_f) \partial_t E^r(p_f, T_f) \\
\qquad + \operatorname{div}_\tau (\mathbf{q}_f^e) - \gamma_n^+ \mathbf{q}_m^e - \gamma_n^- \mathbf{q}_m^e = 0, \\
\\
\sum_{\alpha \in \{l,g\}} s_f^\alpha = 1,
\end{array}
\right.
\quad (2.2)
$$
with the thermodynamical equilibrium for $i = m, f$

$$ egin{cases} c_i^g p_i - p_{\text{sat}}(T_i) c_i^l = 0, \\ \min(s_i^l, 1 - c_i^l) = 0, \\ \min(s_i^g, 1 - c_i^g) = 0, \end{cases} \qquad (2.3) $$

where $p_{\text{sat}}(T)$ is the vapor saturated pressure as a function of the temperature $T$. The Darcy and Fourier laws provide the mass and energy fluxes in the matrix

$$ egin{align} \mathbf{q}_m^{h20} &= \sum_{\alpha \in \{\ell,g\}} \mathbf{q}_m^{\alpha}, \nonumber \\ \mathbf{q}_m^{\alpha} &= c_m^{\alpha} \frac{\rho^{\alpha}(p_m, T_m)}{\mu^{\alpha}(p_m, T_m)} k_{r,m}^{\alpha} (s_m^{\alpha}) \mathbf{V}_m^{\alpha}, \tag{2.4} \\ \mathbf{q}_m^e &= \sum_{\alpha \in \{\ell,g\}} h^{\alpha}(p_m, T_m) \mathbf{q}_m^{\alpha} - \lambda_m \nabla T_m, \nonumber \end{align} $$

and in the fracture network

$$ egin{align} \mathbf{q}_f^{h20} &= \sum_{\alpha \in \{\ell,g\}} \mathbf{q}_f^{\alpha}, \nonumber \\ \mathbf{q}_f^{\alpha} &= c_f^{\alpha} \frac{\rho^{\alpha}(p_f, T_f)}{\mu^{\alpha}(p_f, T_f)} k_{r,f}^{\alpha} (s_f^{\alpha}) \mathbf{V}_f^{\alpha}, \tag{2.5} \\ \mathbf{q}_f^e &= \sum_{\alpha \in \{\ell,g\}} h^{\alpha}(p_f, T_f) \mathbf{q}_f^{\alpha} - d_f \lambda_f \nabla_{\tau} T_f, \nonumber \end{align} $$

where

$$ \mathbf{V}_m^{\alpha} = -\mathbf{K}_m(\nabla p_m - \rho^{\alpha}(p_m, T_m)\mathbf{g}), \quad \mathbf{V}_f^{\alpha} = -d_f \mathbf{K}_f(\nabla_{\tau} p_f - \rho^{\alpha}(p_f, T_f)\mathbf{g}_{\tau}), $$

and $\mathbf{g}_{\tau} = \mathbf{g} - (\mathbf{g} \cdot \mathbf{n}^{+})\mathbf{n}^{+}$.

The system (2.1)-(2.2)-(2.4)-(2.5) is closed with transmission conditions at the matrix fracture interface $\Gamma$. These conditions state the continuity of the pressure and temperature at the matrix fracture interface assuming that the fractures do not act as barriers neither for the Darcy flow nor for the thermal conductivity (see [5, 24, 34, 48]).

$$ egin{align} \gamma^{+} p_m &= \gamma^{-} p_m = \gamma p_m = p_f, \nonumber \\ \gamma^{+} T_m &= \gamma^{-} T_m = \gamma T_m = T_f. \tag{2.6} \end{align} $$

At fracture intersections $\Sigma$, note that we assume mass and energy flux conservation as well as the continuity of the pressure $p_f$ and temperature $T_f$. Homogeneous Neumann boundary conditions are applied for the mass $\mathbf{q}_f^{h20}$ and energy $\mathbf{q}_f^e$ fluxes at the fracture tips $\partial\Gamma \setminus \partial\Omega$ which corresponds to a vanishing aperture (or small aperture compared to the fracture length) assumption at the tips.

# 3. VAG Finite Volume Discretization

## 3.1. Space and time discretizations

The VAG discretization of hybrid-dimensional two-phase Darcy flows introduced in [14] considers generalized polyhedral meshes of $\Omega$ in the spirit of [22]. Let $\mathcal{M}$ be the set of cells that are disjoint open polyhedral subsets of $\Omega$ such that $\cup_{K \in \mathcal{M}} \bar{K} = \bar{\Omega}$, for all $K \in \mathcal{M}$, $\mathbf{x}_K$ denotes the so-called “center” of the cell $K$ under the assumption that $K$ is star-shaped with respect to $\mathbf{x}_K$. The set of faces of the mesh is denoted by $\mathcal{F}$ and $\mathcal{F}_K$ is the set of faces of the cell $K \in \mathcal{M}$. The set of edges of the mesh is
denoted by $\mathcal{E}$ and $\mathcal{E}_\sigma$ is the set of edges of the face $\sigma \in \mathcal{F}$. The set of vertices of the mesh is denoted by $\mathcal{V}$ and for each vertex $\mathbf{s} \in \mathcal{V}$, $\mathbf{x}_\mathbf{s}$ is the coordinates of $\mathbf{s}$. The subset $\mathcal{V}_\sigma \subset \mathcal{V}$ denotes the set of vertices of the face $\sigma$, and for each $K \in \mathcal{M}$ we define $\mathcal{V}_K = \bigcup_{\sigma \in \mathcal{F}_K} \mathcal{V}_\sigma$.

The faces are not necessarily planar. It is just assumed that for each face $\sigma \in \mathcal{F}$, there exists a so-called “center” of the face $\mathbf{x}_\sigma \in \sigma \setminus \bigcup_{\mathbf{a} \in \mathcal{E}_\sigma} \mathbf{a}$ such that $\mathbf{x}_\sigma = \sum_{\mathbf{s} \in \mathcal{V}_\sigma} \beta_{\sigma, \mathbf{s}} \mathbf{x}_\mathbf{s}$, with $\sum_{\mathbf{s} \in \mathcal{V}_\sigma} \beta_{\sigma, \mathbf{s}} = 1$, and $\beta_{\sigma, \mathbf{s}} \ge 0$ for all $\mathbf{s} \in \mathcal{V}_\sigma$; moreover the face $\sigma$ is assumed to be defined by the union of the triangles $T_{\sigma, \mathbf{a}}$ defined by the face center $\mathbf{x}_\sigma$ and each edge $\mathbf{a} \in \mathcal{E}_\sigma$. The mesh is also supposed to be conforming w.r.t. the fracture network $\Gamma$ in the sense that for each $j \in J$ there exists a subset $\mathcal{F}_{\Gamma_j}$ of $\mathcal{F}$ such that

$$ \bar{\Gamma}_j = \bigcup_{\sigma \in \mathcal{F}_{\Gamma_j}} \bar{\sigma}. $$

We will denote by $\mathcal{F}_\Gamma$ the set of fracture faces

$$ \mathcal{F}_{\Gamma} = \bigcup_{j \in J} \mathcal{F}_{\Gamma_j}, $$

and by

$$ \mathcal{V}_{\Gamma} = \bigcup_{\sigma \in \mathcal{F}_{\Gamma}} \mathcal{V}_{\sigma}, $$

the set of fracture nodes. This geometrical discretization of $\Omega$ and $\Gamma$ is denoted in the following by $\mathcal{D}$. In addition, the following notations will be used

$$ \mathcal{M}_{\mathbf{s}} = \{K \in \mathcal{M} | \mathbf{s} \in \mathcal{V}_K\}, \quad \mathcal{M}_{\sigma} = \{K \in \mathcal{M} | \sigma \in \mathcal{F}_K\}, $$

and

$$ \mathcal{F}_{\Gamma, \mathbf{s}} = \{\sigma \in \mathcal{F}_{\Gamma} | \mathbf{s} \in \mathcal{V}_{\sigma}\}. $$

For $N_{t_f} \in \mathbb{N}^*$, let us consider the time discretization $t^0 = 0 < t^1 < \dots < t^{n-1} < t^n < \dots < t^{N_{t_f}} = t_f$ of the time interval $[0, t_f]$. We denote the time steps by $\Delta t^n = t^n - t^{n-1}$ for all $n = 1, \dots, N_{t_f}$.

## 3.2. VAG fluxes and control volumes

The VAG discretization is introduced in [22] for diffusive problems on heterogeneous anisotropic media. Its extension to the hybrid-dimensional Darcy flow model is proposed in [14] based upon the following vector space of degrees of freedom:

$$ V_D = \{v_K, v_s, v_\sigma \in \mathbb{R}, K \in \mathcal{M}, \mathbf{s} \in \mathcal{V}, \sigma \in \mathcal{F}_\Gamma\}. $$

The degrees of freedom are illustrated in Figure 3.1 for a given cell $K$ with one fracture face $\sigma$ in bold.

The matrix degrees of freedom are defined by the set of cells $\mathcal{M}$ and by the set of nodes $\mathcal{V} \setminus \mathcal{V}_\Gamma$ excluding the nodes at the matrix fracture interface $\Gamma$. The fracture faces $\mathcal{F}_\Gamma$ and the fracture nodes $\mathcal{V}_\Gamma$ 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.2). The degrees of freedom at the fracture intersection $\Sigma$ are defined by the set of nodes $\mathcal{V}_\Sigma \subset \mathcal{V}_\Gamma$ located on $\bar{\Sigma}$. The set of nodes at the Dirichlet boundaries $\partial\Omega_D$ and $\partial\Gamma_D$ is denoted by $\mathcal{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 or energy balance equation. The two main ingredients are therefore the conservative fluxes and the control volumes. The VAG matrix and fracture fluxes are illustrated in Figure 3.1. For $u_D \in \mathcal{V}_D$, the matrix fluxes $F_{K,\nu}(u_D)$ connect the cell $K \in \mathcal{M}$ to the degrees of freedom located at the boundary of $K$, namely $\nu \in \Xi_K = \mathcal{V}_K \cup (\mathcal{F}_K \cap \mathcal{F}_\Gamma)$. The fracture fluxes $F_{\sigma,s}(u_D)$ connect each fracture face $\sigma \in \mathcal{F}_\Gamma$ to its nodes $\mathbf{s} \in \mathcal{V}_\sigma$. These discrete fluxes will be used for the discretization of the diffusive part of the continuous fluxes, corresponding to the terms $-\mathbf{K}_m \nabla p_m$ and $-\lambda_m \nabla T_m$ in the matrix and to $-d_f \mathbf{K}_f \nabla_\tau p_f$ and $-d_f \lambda_f \nabla_\tau T_f$ along the fractures. 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

$$F_{K,\nu}(u_D) = \sum_{\nu' \in \Xi_K} T_K^{\nu,\nu'}(u_K - u_{\nu'}),$$

with a symmetric positive definite transmissibility matrix $T_K = (T_K^{\nu,\nu'})(\nu,\nu') \in \Xi_K \times \Xi_K$. This transmissibility depends only on the cell $K$ geometry (including the choices of $x_K$ and of $x_\sigma$, $\sigma \in \mathcal{F}_K$) and on the cell matrix diffusion tensor corresponding to the permeability $K_m$ for the Darcy fluxes and to the thermal conductivity $\lambda_m$ for the Fourier fluxes. The fracture fluxes are given by

$$F_{\sigma,s}(u_D) = \sum_{s \in \mathcal{V}_\sigma} T_\sigma^{s,s'}(u_\sigma - u_{s'}),$$

with a symmetric positive definite transmissibility matrix $T_\sigma = (T_\sigma^{s,s'})(s,s') \in \mathcal{V}_\sigma \times \mathcal{V}_\sigma$ depending only on the fracture face $\sigma$ geometry (including the choice of $x_\sigma$) and on the fracture tangential diffusion tensor. This diffusion tensor corresponds to $d_f K_f$ for the Darcy fluxes and to $d_f \lambda_f$ for the Fourier fluxes. Let us refer to [14] for a more detailed presentation and for the definition of $T_K$ and $T_\sigma$.

<figure><img src="image_4.png" /><figcaption>FIGURE 3.1. For a cell $K$ and a fracture face $\sigma$ (in bold), examples of VAG degrees of freedom $u_K$, $u_s$, $u_\sigma$, $u_{s'}$ and VAG fluxes $F_{K,\sigma}$, $F_{K,s}$, $F_{K,s'}$, $F_{\sigma,s}$.</figcaption></figure>

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 \in \mathcal{M}$, by

$$\bar{K} = \bar{\omega}_K \cup \left( \bigcup_{s \in \mathcal{V}_K \setminus \mathcal{V}_D} \bar{\omega}_{K,s} \right),$$

and, for all $\sigma \in \mathcal{F}_\Gamma$, by

$$\bar{\sigma} = \bar{\Sigma}_\sigma \cup \left( \bigcup_{s \in \mathcal{V}_\sigma \setminus \mathcal{V}_D} \bar{\Sigma}_{\sigma,s} \right).$$

The practical implementation of the scheme does not require building the geometry of these partitions explicitly but just defining the matrix volume fractions

$$\alpha_{K,s} = \frac{\int_{\omega_{K,s}} d\mathbf{x}}{\int_K d\mathbf{x}}, \mathbf{s} \in \mathcal{V}_K \setminus (\mathcal{V}_D \cup \mathcal{V}_\Gamma), K \in \mathcal{M},$$
constrained to satisfy $\alpha_{K,\nu} \ge 0$, and $\sum_{s \in V_K \setminus (V_D \cup V_\Gamma)} \alpha_{K,s} \le 1$, as well as the fracture volume fractions

$$ \alpha_{\sigma,s} = \frac{\int_{\Sigma_{\sigma,s}} d_f(\mathbf{x}) d\tau(\mathbf{x})}{\int_{\sigma} d_f(\mathbf{x}) d\tau(\mathbf{x})}, \mathbf{s} \in V_{\sigma} \setminus V_D, \sigma \in \mathcal{F}_{\Gamma}, $$

constrained to satisfy $\alpha_{\sigma,s} \ge 0$, and $\sum_{s \in V_{\sigma} \setminus V_D} \alpha_{\sigma,s} \le 1$, where we denote by $d\tau(\mathbf{x})$ the two-dimensional Lebesgue measure on $\Gamma$. Let us also set

$$ \phi_K = \left(1 - \sum_{s \in V_K \setminus (V_D \cup V_\Gamma)} \alpha_{K,s}\right) \int_K \phi_m(\mathbf{x}) d\mathbf{x} \quad \text{for } K \in \mathcal{M}, $$

and

$$ \phi_{\sigma} = \left(1 - \sum_{s \in V_{\sigma} \setminus V_D} \alpha_{\sigma,s}\right) \int_{\sigma} \phi_{f}(\mathbf{x}) d_{f}(\mathbf{x}) d\tau(\mathbf{x}) \quad \text{for } \sigma \in \mathcal{F}_{\Gamma}, $$

as well as

$$ \phi_s = \sum_{K \in \mathcal{M}_s} \alpha_{K,s} \int_K \phi_m(\mathbf{x}) d\mathbf{x} \quad \text{for } \mathbf{s} \in \mathcal{V} \setminus (V_D \cup V_\Gamma), $$

and

$$ \phi_s = \sum_{\sigma \in \mathcal{F}_{\Gamma,s}} \alpha_{\sigma,s} \int_{\sigma} \phi_f(\mathbf{x}) d_f(\mathbf{x}) d\tau(\mathbf{x}) \quad \text{for } \mathbf{s} \in V_{\Gamma} \setminus V_D, $$

which correspond to the porous volumes distributed to the degrees of freedom excluding the Dirichlet nodes. The rock complementary volume in each control volume $\nu \in \mathcal{M} \cup \mathcal{F}_{\Gamma} \cup (\mathcal{V} \setminus V_D)$ is denoted by $\bar{\phi}_{\nu}$.

As shown in [14], 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.2, 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.

<figure><img src="image_5.png" /><figcaption>FIGURE 3.2. Example of control volumes at cells, fracture face, and nodes, in the case of two cells K and L separated by one fracture face $\sigma$ (the width of the fracture is enlarged in this figure). The control volumes are chosen to avoid mixing fracture and matrix rocktypes.</figcaption></figure>

A rocktype is assigned to each cell, node and fracture face. In our case, for cells and for nodes not located along the fractures, the matrix rocktype is assigned. For fracture nodes and faces at the interface between the matrix and the fracture rocktypes, the fracture rocktype is assigned corresponding to the most pervious rock type consistently with the choice of the control volumes (see [14]). For convenience's sake, in the following, we will denote by $k_{r,\nu}^{\alpha}$ the corresponding relative permeability function for $\nu \in \mathcal{M} \cup \mathcal{V} \cup \mathcal{F}_{\Gamma}$.
In the following, we will keep the notation $F_{K,s}$, $F_{K,\sigma}$, $F_{\sigma,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_\sigma$ are computed only once.

The VAG Fourier fluxes are denoted in the following by $G_{K,s}$, $G_{K,\sigma}$, $G_{\sigma,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.

### 3.3. Multi-branch non-isothermal well model

Let $\mathcal{W}$ denote the set of wells. Each multi-branch well $\omega \in \mathcal{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 $\omega \in \mathcal{W}$ is denoted by $\mathcal{V}_\omega \subset \mathcal{V}$ and its root node is denoted by $\mathbf{s}_\omega$. A partial ordering is defined on the set of vertices $\mathcal{V}_\omega$ with $\mathbf{s} < \mathbf{s}'$ if and only if the unique path from the root $\mathbf{s}_\omega$ to $\mathbf{s}'$ passes through $\mathbf{s}$. The set of edges of the well $\omega$ is denoted by $\mathcal{E}_\omega$ and for each edge $\mathbf{a} \in \mathcal{E}_\omega$ we set $\mathbf{a} = \mathbf{s}\mathbf{s}'$ with $\mathbf{s} < \mathbf{s}'$ (i.e. $\mathbf{s}$ is the parent node of $\mathbf{s}'$, see Figure 3.3). It is assumed that $\mathcal{V}_{\omega_1} \cap \mathcal{V}_{\omega_2} = \emptyset$ for any $\omega_1, \omega_2 \in \mathcal{W}$ such that $\omega_1 \neq \omega_2$.

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 as a first step. The latest shall be taken into account when modeling the wellbore flow up to the surface. It is assumed that the radius $r_\omega$ of each well $\omega \in \mathcal{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 $\mathbf{s} \in \mathcal{V}_\omega$ is obtained using the Two Point Flux Approximation

$$V_{\mathbf{s}}^{\omega} = W_{\mathbf{s}}(p_{\mathbf{s}} - p_{\mathbf{s}}^{\omega}),$$

where $p_{\mathbf{s}}$ is the reservoir pressure at node $\mathbf{s}$ and $p_{\mathbf{s}}^{\omega}$ is the well pressure at node $\mathbf{s}$. The Well Index $W_{\mathbf{s}}$ is typically computed using Peaceman's approach (see [21, 38, 39]) 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. At each well node $\mathbf{s} \in \mathcal{V}_\omega$ the temperature inside the well is denoted by $T_{\mathbf{s}}^{\omega}$ and the volume fractions by $s_{\mathbf{s},\omega}^{\alpha}$, $\alpha \in \{\ell, g\}$. The temperature in the reservoir at node $\mathbf{s}$ is denoted by $T_{\mathbf{s}}$, the saturations by $s_{\mathbf{s}}^{\alpha}$, and the phase mass fractions by $c_{\mathbf{s}}^{\alpha}$ for $\alpha \in \{\ell, g\}$.

For any $a \in \mathbb{R}$, let us define $a^+ = \max(a, 0)$ and $a^- = \min(a, 0)$. The mass flow rates between the reservoir and the well $\omega$ at a given node $\mathbf{s} \in \mathcal{V}_\omega$ are defined by the following phase based upwind approximation of the mobilities:

$$ 
\begin{aligned}
q_{\mathbf{s},\alpha}^{r \to \omega} &= \beta_{\omega}^{\text{inj}} \frac{\rho^{\alpha}(p_{\mathbf{s}}^{\omega}, T_{\mathbf{s}}^{\omega})}{\mu^{\alpha}(p_{\mathbf{s}}^{\omega}, T_{\mathbf{s}}^{\omega})} k_{r,\mathbf{s}}^{\alpha}(s_{\mathbf{s},\omega}^{\alpha})(V_{\mathbf{s}}^{\omega})^{-} + \beta_{\omega}^{\text{prod}} c_{\mathbf{s}}^{\alpha} \frac{\rho^{\alpha}(p_{\mathbf{s}}, T_{\mathbf{s}})}{\mu^{\alpha}(p_{\mathbf{s}}, T_{\mathbf{s}})} k_{r,\mathbf{s}}^{\alpha}(s_{\mathbf{s}}^{\alpha})(V_{\mathbf{s}}^{\omega})^{+}, \\
q_{\mathbf{s},h_{20}}^{r \to \omega} &= \sum_{\alpha \in \{\ell,g\}} q_{\mathbf{s},\alpha}^{r \to \omega},
\end{aligned}
\quad (3.1) $$

and the energy flow rate is defined similarly by

$$ q_{\mathbf{s},e}^{r \to \omega} = \sum_{\alpha \in \{\ell,g\}} h^{\alpha}(p_{\mathbf{s}}^{\omega}, T_{\mathbf{s}}^{\omega})(q_{\mathbf{s},\alpha}^{r \to \omega})^{-} + h^{\alpha}(p_{\mathbf{s}}, T_{\mathbf{s}})(q_{\mathbf{s},\alpha}^{r \to \omega})^{+}. \quad (3.2) $$

The well coefficients $\beta_{\omega}^{\text{inj}}$ and $\beta_{\omega}^{\text{prod}}$ are used to impose specific well behavior. The general case corresponds to $\beta_{\omega}^{\text{inj}} = \beta_{\omega}^{\text{prod}} = 1$. Yet, for an injection well, it will be convenient as explained in Subsection 3.3.2, to impose that the mass flow rates $q_{\mathbf{s},h_{20}}^{r \to \omega}$ are non positive for all nodes $\mathbf{s} \in \mathcal{V}_\omega$
corresponding to set $\beta_{\omega}^{inj} = 1$ and $\beta_{\omega}^{prod} = 0$. Likewise, for a production well, it will be convenient as explained in Subsection 3.3.3, to set $\beta_{\omega}^{inj} = 0$ and $\beta_{\omega}^{prod} = 1$ which corresponds to assume that the mass flow rates $q_{s,h_{2o}}^{r \to \omega}$ are non negative for all nodes $s \in \mathcal{V}_{\omega}$. 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 happens in geothermal wells, typically in closed wells.

### 3.3.1. Well physical model

Our conceptual model inside the well assumes that the flow is quasi static at the reservoir time scale along with perfect mixing and thermal equilibrium. The Fourier fluxes and the wall friction are neglected and the pressure distribution is assumed hydrostatic along the well.

For the sake of simplicity, the flow rate between the reservoir and the well is considered concentrated at each node $s$ of the well. For each edge $a \in \mathcal{E}_{\omega}$, let us denote by $q_{a}^{\alpha}$ the mass flow rate of phase $\alpha$ along the edge $a$ oriented positively from $s'$ to $s$ with $a = ss'$ (let us recall that $s$ is the parent node of $s'$).

Let $\alpha \in \{\ell, g\}$, the set of well unknowns is defined at each node $s \in \mathcal{V}_{\omega}$ by the well pressure $p_{s}^{\omega}$, the well temperature $T_{s}^{\omega}$, the well saturations $s_{s,\omega}^{\alpha}$, and at each edge $a \in \mathcal{E}_{\omega}$ by the mass flow rates $q_{a}^{\alpha}$. These well unknowns are complemented by the well mass flow rates $q_{\omega}^{\alpha}$ which are non negative for production wells and non positive for injection wells (see Figure 3.3).

<figure><img src="image_6.png" /><figcaption>FIGURE 3.3. Example of multi-branch well $\omega$ with its root node $s_{\omega}$, one edge $a = ss'$ and the main physical quantities: the well mass flow rates $q_{\omega}^{\alpha}$, the mass and energy flow rates between the reservoir and the well $q_{s,h_{2o}}^{r \to \omega}$, $q_{s,e}^{r \to \omega}$, the well node pressure, temperature and saturations $p_{s}^{\omega}$, $T_{s}^{\omega}$, $s_{s,\omega}^{\alpha}$, and the edge mass flow rates $q_{a}^{\alpha}$.</figcaption></figure>

For each edge $a = ss' \in \mathcal{E}_{\omega}$, and each phase $\alpha$, let us define the following phase based upwind approximations of the specific enthalpy, mass density and saturation

$$h_{a}^{\alpha} = \begin{cases} h^{\alpha}(p_{s'}^{\omega}, T_{s'}^{\omega}) & \text{if } q_{a}^{\alpha} \ge 0, \\ h^{\alpha}(p_{s}^{\omega}, T_{s}^{\omega}) & \text{if } q_{a}^{\alpha} < 0. \end{cases} \qquad \rho_{a}^{\alpha} = \begin{cases} \rho^{\alpha}(p_{s'}^{\omega}, T_{s'}^{\omega}) & \text{if } q_{a}^{\alpha} \ge 0, \\ \rho^{\alpha}(p_{s}^{\omega}, T_{s}^{\omega}) & \text{if } q_{a}^{\alpha} < 0. \end{cases} \qquad s_{a}^{\alpha} = \begin{cases} s_{s',\omega}^{\alpha} & \text{if } q_{a}^{\alpha} \ge 0, \\ s_{s,\omega}^{\alpha} & \text{if } q_{a}^{\alpha} < 0. \end{cases} \qquad (3.3)$$
For all $\mathbf{s}\mathbf{s}' = \mathbf{a} \in \mathcal{E}_\omega$, let us set $\kappa_{\mathbf{a},\mathbf{s}'} = -1$ and $\kappa_{\mathbf{a},\mathbf{s}} = 1$. The well equations account for the mass and energy conservations at each node of the well combined with the sum to one of the saturations and the thermodynamical equilibrium. Let $\mathcal{E}_\omega \subset \mathcal{E}_\omega$ denote the set of well edges sharing the node $\mathbf{s} \in \mathcal{V}_\omega$, then for all $\mathbf{s} \in \mathcal{V}_\omega$ we obtain the equations

$$ 
\left\{ 
\begin{array}{l} 
q_{\mathbf{s},h_2o}^{r \to \omega} + \sum_{\mathbf{a} \in \mathcal{E}_{\mathbf{s}}^{\omega}} \sum_{\alpha \in \{\ell,g\}} \kappa_{\mathbf{a},\mathbf{s}} q_{\mathbf{a}}^{\alpha} = \delta_{\mathbf{s}}^{\mathbf{s}\omega} \sum_{\alpha \in \{\ell,g\}} q_{\omega}^{\alpha}, \\ 
q_{\mathbf{s},e}^{r \to \omega} + \sum_{\mathbf{a} \in \mathcal{E}_{\mathbf{s}}^{\omega}} \sum_{\alpha \in \{\ell,g\}} h_{\mathbf{a}}^{\alpha} \kappa_{\mathbf{a},\mathbf{s}} q_{\mathbf{a}}^{\alpha} = \delta_{\mathbf{s}}^{\mathbf{s}\omega} \sum_{\alpha \in \{\ell,g\}} (\bar{h}_{\omega}^{\alpha}(q_{\omega}^{\alpha})^{-} + h^{\alpha}(p_{\mathbf{s}}^{\omega}, T_{\mathbf{s}}^{\omega})(q_{\omega}^{\alpha})^{+}), \\ 
s_{\mathbf{s},\omega}^{\ell} + s_{\mathbf{s},\omega}^{g} = 1, \\ 
p_{\mathbf{s}}^{\omega} \left\{ 
\begin{array}{ll} 
= p_{\text{sat}}(T_{\mathbf{s}}^{\omega}) & \text{if } s_{\mathbf{s},\omega}^{g} > 0 \text{ and } s_{\mathbf{s},\omega}^{\ell} > 0, \\ 
\ge p_{\text{sat}}(T_{\mathbf{s}}^{\omega}) & \text{if } s_{\mathbf{s},\omega}^{g} = 0, \\ 
\le p_{\text{sat}}(T_{\mathbf{s}}^{\omega}) & \text{if } s_{\mathbf{s},\omega}^{\ell} = 1, 
\end{array} 
\right. 
\end{array} 
\right. 
\quad (3.4) 
$$

where $\delta$ stands for the Kronecker symbol, and $\bar{h}_{\omega}^{\alpha}$ for prescribed specific enthalpies in the case of injection wells. Inside the well, the hypothesis of hydrostatic pressure distribution implies that

$$ p_{\mathbf{s}}^{\omega} - p_{\mathbf{s}'}^{\omega} + \rho_{\mathbf{a}} g(z_{\mathbf{s}} - z_{\mathbf{s}'}) = 0, \quad (3.5) $$

for each edge $\mathbf{s}\mathbf{s}' = \mathbf{a} \in \mathcal{E}_\omega$, where $\rho_{\mathbf{a}}$ is the mass density of the liquid-vapor mixture. The system is completed by a slip closure law expressing the slip between the liquid velocity $u_{\mathbf{a}}^{\ell}$ and the gas velocity $u_{\mathbf{a}}^{g}$ at each edge $\mathbf{a} \in \mathcal{E}_\omega$ with

$$ q_{\mathbf{a}}^{\alpha} = \pi r_{\omega}^{2} \rho_{\mathbf{a}}^{\alpha} s_{\mathbf{a}}^{\alpha} u_{\mathbf{a}}^{\alpha}. $$

In the following simplified well models developed in Subsections 3.3.2 and 3.3.3, a zero slip law will be assumed for simplicity in such a way that $u_{\mathbf{a}}^{\ell} = u_{\mathbf{a}}^{g}$. Note that these simplified well models could be easily extended to account for non-zero slip laws as well as for an explicit approximation of the wall friction along the wells. The two fundamental assumptions to obtain these simplified well models are

(i) prescribed sign of the mass flow rates $q_{\mathbf{s},\alpha}^{r \to \omega}$, $\mathbf{s} \in \mathcal{V}_\omega$, forced to be all non-negative for production wells and all non-positive for injection wells,

(ii) neglected Fourier fluxes compared with thermal convection fluxes.

The well boundary conditions prescribe a limit total mass flow rate $\bar{q}_\omega$ and a limit bottom hole pressure $\bar{p}_\omega$. Using the notations

$$ p_{\omega} = p_{\mathbf{s}\omega}^{\omega} \quad \text{and} \quad q_{\omega} = \sum_{\alpha \in \{\ell,g\}} q_{\omega}^{\alpha}, $$

the well monitoring conditions prescribe complementary constraints between $q_\omega - \bar{q}_\omega$ and $p_\omega - \bar{p}_\omega$. More precisely, the conditions

$$ (q_{\omega} - \bar{q}_{\omega})(p_{\omega} - \bar{p}_{\omega}) = 0, \quad q_{\omega} - \bar{q}_{\omega} \ge 0, \quad p_{\omega} - \bar{p}_{\omega} \le 0, $$

are imposed for injection wells (with $\bar{q}_\omega \le 0$), and

$$ (q_{\omega} - \bar{q}_{\omega})(p_{\omega} - \bar{p}_{\omega}) = 0, \quad q_{\omega} - \bar{q}_{\omega} \le 0, \quad p_{\omega} - \bar{p}_{\omega} \ge 0, $$

for production wells (with $\bar{q}_\omega \ge 0$).

In the following subsections, we consider the particular case of injection wells assuming a pure liquid phase, and the case of production wells. The flow rates are enforced to be non positive (resp. non negative) at all well nodes for injection wells (resp. production wells). It corresponds to setting $\beta_{\omega}^{\text{inj}} = 1$, $\beta_{\omega}^{\text{prod}} = 0$ for an injection well and $\beta_{\omega}^{\text{inj}} = 0$, $\beta_{\omega}^{\text{prod}} = 1$ for a production well. The limit bottom hole pressure $\bar{p}_\omega$ is a maximum (resp. minimum) pressure and the limit total mass flow rate $\bar{q}_\omega$ is a minimum non positive (resp. maximum non negative) flow rate for injection (resp. production) wells.
In both cases, using an explicit computation of the hydrostatic pressure drop, the well model will be reduced to a single equation and a single implicit unknown corresponding to the well reference pressure $p_\omega$ (see e.g. [8]).

### 3.3.2. Liquid injection wells

The injection well model sets $\beta_\omega^{inj} = 1$, $\beta_\omega^{prod} = 0$ and prescribes the minimum well total mass flow rate $\bar{q}_\omega \le 0$, the well maximum bottom hole pressure $\bar{p}_\omega$ and the well specific liquid enthalpy $\bar{h}_\omega^\ell$. It is assumed that the injection is in liquid phase and that no gas will appear in the well during the simulation as it is usually the case in geothermal systems.

Since $\beta_\omega^{inj} = 1$ and $\beta_\omega^{prod} = 0$, the mass flow rates $q_\alpha^\alpha$ are enforced to be non negative. Let us show that together with the assumption that the gas phase does not appear in the well, it implies that $h^\ell(p_\omega^\omega, T_\omega^\omega) = \bar{h}_\omega^\ell$ for all $\mathbf{s} \in \mathcal{V}_\omega$. Using the mass and energy conservation equations (3.4) at the head node, we obtain that

$$h^\ell(p_{\mathbf{s},\omega}^\omega, T_{\mathbf{s},\omega}^\omega) \underbrace{\left( q_{\mathbf{s},\omega,\text{h2o}}^{\mathbf{r} \to \omega} + \sum_{\mathbf{a} \in \mathcal{E}_{\mathbf{s},\omega}^\omega} q_\alpha^\alpha \right)}_{q_\omega \text{ from the mass conservation}} = \bar{h}_\omega^\ell q_\omega,$$

from which we derive $h^\ell(p_{\mathbf{s},\omega}^\omega, T_{\mathbf{s},\omega}^\omega) = \bar{h}_\omega^\ell$. Then, assuming that the property holds for a parent node, it propagates similarly to all its sons, and hence to all nodes $\mathbf{s} \in \mathcal{V}_\omega$ by induction.

Given the previous time step well reference pressure $p_\omega^{n-1} = p_{\mathbf{s},\omega}^{\omega,n-1}$, we first compute the pressures along the well solving the equations

$$p_{\mathbf{s}}^{\omega} - p_{\mathbf{s}'}^{\omega} + \rho_{\mathbf{a}} g(z_{\mathbf{s}} - z_{\mathbf{s}'}) = 0 \text{ for all } \mathbf{a} = \mathbf{s}\mathbf{s}' \in \mathcal{E}_{\omega},$$

$$p_{\mathbf{s},\omega}^{\omega} = p_{\mathbf{s},\omega}^{\omega,n-1},$$

$$\rho_{\mathbf{a}} = \rho^{\ell}(p_{\mathbf{s}}^{\omega}, T_{\mathbf{s}}^{\omega}) \text{ for all } \mathbf{a} = \mathbf{s}\mathbf{s}' \in \mathcal{E}_{\omega},$$

$$h^\ell(p_{\mathbf{s}}^\omega, T_{\mathbf{s}}^\omega) = \bar{h}_\omega^\ell \text{ for all } \mathbf{s} \in \mathcal{V}_\omega.$$

We deduce the explicit pressure drops

$$\Delta p_{\mathbf{s}}^{\omega,n-1} = p_{\mathbf{s}}^{\omega} - p_{\omega}^{n-1},$$

which provide for all $\mathbf{s} \in \mathcal{V}_\omega$ the pressures $p_\omega^{\omega,n}$ and temperatures $T_\omega^{\omega,n}$ along the well at the current time step $n$ such that

$$p_{\mathbf{s}}^{\omega,n} = p_{\omega}^{n} + \Delta p_{\mathbf{s}}^{\omega,n-1},$$

$$h^\ell(p_{\mathbf{s}}^{\omega,n}, T_{\mathbf{s}}^{\omega,n}) = \bar{h}_\omega^\ell.$$

The mass and energy flow rates at each node $\mathbf{s} \in \mathcal{V}_\omega$ between the reservoir and the well are defined by (3.1)-(3.2) with $\beta_\omega^{inj} = 1$ and $\beta_\omega^{prod} = 0$ and depend only on the implicit unknowns $p_\omega^n$ and $p_\mathbf{s}^n$. They are respectively denoted by $q_{\mathbf{s},\text{h2o}}^{\mathbf{r} \to \omega}(p_\omega^n, p_\omega^n)$ and $q_{\mathbf{s},\text{c}}^{\mathbf{r} \to \omega}(p_\mathbf{s}^n, p_\omega^n)$.

The well equation at the current time step is defined by the following complementary constraints between the prescribed minimum well total mass flow rate and the prescribed maximum bottom hole
pressure

$$ \left\{ \begin{array}{l} \left( \sum_{s \in \mathcal{V}_{\omega}} q_{s, h_{2o}}^{r \to \omega} (p_s^n, p_{\omega}^n) - \bar{q}_{\omega} \right) (\bar{p}_{\omega} - p_{\omega}^n) = 0, \\ \\ \sum_{s \in \mathcal{V}_{\omega}} q_{s, h_{2o}}^{r \to \omega} (p_s^n, p_{\omega}^n) - \bar{q}_{\omega} \geq 0, \\ \\ \bar{p}_{\omega} - p_{\omega}^n \geq 0. \end{array} \right. \qquad (3.6) $$

### 3.3.3. Production wells

The production well model sets $\beta_{\omega}^{\text{inj}} = 0$, $\beta_{\omega}^{\text{prod}} = 1$ and prescribes the maximum well total mass flow rate $\bar{q}_{\omega} \geq 0$ and the well minimum bottom hole pressure $\bar{p}_{\omega}$.

The solution at the previous time step $n-1$ provides the pressure drop $\Delta p_{\mathbf{s}}^{\omega,n-1}$ at each node $\mathbf{s} \in \mathcal{V}_{\omega}$. This computation based on thermodynamical equilibrium 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_{\mathbf{s}}^{\omega,n} = p_{\omega}^{n} + \Delta p_{\mathbf{s}}^{\omega,n-1}. $$

The mass and energy flow rates at each node $\mathbf{s} \in \mathcal{V}_{\omega}$ between the reservoir and the well are defined by (3.1)-(3.2) with $\beta_{\omega}^{\text{inj}} = 0$ and $\beta_{\omega}^{\text{prod}} = 1$ and depend only on the implicit reservoir unknowns $X_{\mathbf{s}}^{n}$ setting

$$ X_{\mathbf{s}} = (P_{\mathbf{s}}, T_{\mathbf{s}}, s_{\mathbf{s}}^{\ell}, s_{\mathbf{s}}^{g}, c_{\mathbf{s}}^{\ell}, c_{\mathbf{s}}^{g}), $$

and on the implicit well unknown $p_{\omega}^{n}$. They are respectively denoted by $q_{s,h_{2o}}^{r \to \omega}(X_{\mathbf{s}}^{n}, p_{\omega}^{n})$ and $q_{s,e}^{r \to \omega}(X_{\mathbf{s}}^{n}, p_{\omega}^{n})$.

The well equation at the current time step is defined by the following complementary constraints between the prescribed maximum well total mass flow rate and the prescribed minimum bottom hole pressure

$$ \left\{ \begin{array}{l} \left( \bar{q}_{\omega} - \sum_{s \in \mathcal{V}_{\omega}} q_{s, h_{2o}}^{r \to \omega} (X_{\mathbf{s}}^{n}, p_{\omega}^{n}) \right) (p_{\omega}^{n} - \bar{p}_{\omega}) = 0, \\ \\ \bar{q}_{\omega} - \sum_{s \in \mathcal{V}_{\omega}} q_{s, h_{2o}}^{r \to \omega} (X_{\mathbf{s}}^{n}, p_{\omega}^{n}) \geq 0, \\ \\ p_{\omega}^{n} - \bar{p}_{\omega} \geq 0. \end{array} \right. \qquad (3.7) $$

Let us now detail the computation of the pressure drop at each node $\mathbf{s} \in \mathcal{V}_{\omega}$ using the previous time step solution $n-1$ consisting of the reservoir unknowns and the well pressures. We first compute the well temperature $T_{\mathbf{s}}^{\omega,n-1}$ and saturations $s_{\mathbf{s},\omega}^{\alpha,n-1}$ at each node $\mathbf{s}$ using equations (3.4). Summing the mass and energy equations of (3.4) over all nodes $\mathbf{s}'' \geq \mathbf{s}$, we obtain for all $\mathbf{a} = \mathbf{s}'\mathbf{s} \in \mathcal{E}_{\omega}$ that

$$ \begin{aligned} \sum_{\alpha \in \{\ell, g\}} Q_{\mathbf{a}}^{\alpha, n-1} &= \sum_{\mathbf{s}'' \in \mathcal{V}_{\omega} | \mathbf{s}'' \geq \mathbf{s}} q_{\mathbf{s}'' h_{2o}}^{r \to \omega} (X_{\mathbf{s}''}^{n-1}, p_{\omega}^{n-1}) = Q_{\mathbf{s}, h_{2o}}^{\omega}, \\ \sum_{\alpha \in \{\ell, g\}} h^{\alpha} (p_{\mathbf{s}}^{\omega, n-1}, T_{\mathbf{s}}^{\omega, n-1}) Q_{\mathbf{a}}^{\alpha, n-1} &= \sum_{\mathbf{s}' \in \mathcal{V}_{\omega} | \mathbf{s}'' \geq \mathbf{s}} q_{\mathbf{s}', e}^{r \to \omega} (X_{\mathbf{s}''}^{n-1}, p_{\omega}^{n-1}) = Q_{\mathbf{s}, e}^{\omega}, \end{aligned} $$

with

$$ Q_{\mathbf{a}}^{\alpha,n-1} = \pi r_{\omega}^{2} \rho^{\alpha} (p_{\mathbf{s}}^{\omega,n-1}, T_{\mathbf{s}}^{\omega,n-1}) s_{\mathbf{s},\omega}^{\alpha,n-1} u_{\mathbf{a}}^{\alpha,n-1}, \alpha \in \{\ell, g\}. $$

It results that the thermodynamical equilibrium at fixed well pressure $p_{\mathbf{s}}^{\omega,n-1}$, mass $Q_{\mathbf{s},h_{2o}}^{\omega}$ and energy $Q_{\mathbf{s},e}^{\omega}$ provides the well temperature $T_{\mathbf{s}}^{\omega,n-1}$ and the well saturations $s_{\mathbf{s},\omega}^{\alpha,n-1}$ at node $\mathbf{s}$ as follows. Let
us set $p = p_{\mathbf{s}}^{\omega,n-1}$. We first assume that both phases are present which implies that $T_{\text{sat}} = (p_{\text{sat}})^{-1}(p)$ and that the liquid mass fraction is given by

$$c^{\ell} = \frac{h^{g}(p, T_{\text{sat}}) - \frac{Q_{\mathbf{s},e}^{\omega}}{Q_{\mathbf{s},h_{2o}}^{\omega}}}{h^{g}(p, T_{\text{sat}}) - h^{\ell}(p, T_{\text{sat}})}.$$

The following alternatives are checked:

**Two-phase state:** if $0 < c^{\ell} < 1$, the two-phase state is confirmed. Using the zero slip assumption, we obtain

$$T_{\mathbf{s}}^{\omega,n-1} = T_{\text{sat}} \text{ and } s_{\mathbf{s},\omega}^{\ell,n-1} = 1 - s_{\mathbf{s},\omega}^{g,n-1} = \frac{\frac{c^{\ell}}{\rho^{\ell}(p,T_{\text{sat}})}}{\frac{c^{\ell}}{\rho^{\ell}(p,T_{\text{sat}})} + \frac{1-c^{\ell}}{\rho^{g}(p,T_{\text{sat}})}}.$$

**Liquid state:** if $c^{\ell} \ge 1$, then only the liquid phase is present, we set $s_{\mathbf{s},\omega}^{\ell,n-1} = 1$, $s_{\mathbf{s},\omega}^{g,n-1} = 0$, and $T_{\mathbf{s}}^{\omega,n-1}$ is the solution $T$ of

$$h^{\ell}(p, T) = \frac{Q_{\mathbf{s},e}^{\omega}}{Q_{\mathbf{s},h_{2o}}^{\omega}}.$$

**Gas state:** if $c^{\ell} \le 0$, then only the gas phase is present, we set $s_{\mathbf{s},\omega}^{\ell,n-1} = 0$, $s_{\mathbf{s},\omega}^{g,n-1} = 1$, and $T_{\mathbf{s}}^{\omega,n-1}$ is the solution $T$ of

$$h^{g}(p, T) = \frac{Q_{\mathbf{s},e}^{\omega}}{Q_{\mathbf{s},h_{2o}}^{\omega}}.$$

Then, the explicit pressure drop

$$\Delta p_{\mathbf{s}}^{\omega,n-1} = p_{\mathbf{s}}^{\omega} - p_{\omega}^{n-1},$$

is obtained from

$$p_{\mathbf{s}}^{\omega} - p_{\mathbf{s}'}^{\omega} + \rho_{\mathbf{a}}g(z_{\mathbf{s}} - z_{\mathbf{s}'}) = 0 \text{ for all } \mathbf{a} = \mathbf{s}\mathbf{s}' \in \mathcal{E}_{\omega},$$

$$p_{\mathbf{s}_{\omega}}^{\omega} = p_{\mathbf{s}_{\omega}}^{\omega,n-1},$$

$$\rho_{\mathbf{a}} = \sum_{\alpha \in \{\ell, g\}} s_{\mathbf{s},\omega}^{\alpha,n-1} \rho^{\alpha}(p_{\mathbf{s}}^{\omega,n-1}, T_{\mathbf{s}}^{\omega,n-1}) \text{ for all } \mathbf{a} = \mathbf{s}\mathbf{s}' \in \mathcal{E}_{\omega}.$$

### 3.4. Discretization of the hybrid-dimensional non-isothermal two-phase flow model

The time integration is based on a semi-implicit Euler scheme. As discussed above, the mixture density along the well and the mean thermal conductivities are approximated explicitly. All the other terms are approximated implicitly in order to avoid severe restrictions on the time steps due to the small volumes and high velocities in the fractures. A phase based upwind scheme is used for the approximation of the mobilities in the mass and energy fluxes (see e.g. [9]). At the matrix fracture interfaces, we avoid mixing matrix and fracture rocktypes by choosing appropriate control volumes for $\sigma \in \mathcal{F}_{\Gamma}$ and $\mathbf{s} \in \mathcal{V}_{\Gamma}$ (see Figure 3.2). In order to avoid tiny control volumes at the nodes $\mathbf{s} \in \mathcal{V}_{\Sigma}$ located at the fracture intersection, the volume is distributed to such a node $\mathbf{s}$ from all the fracture faces containing the node $\mathbf{s}$.

For each $\nu \in \mathcal{M} \cup \mathcal{F}_{\Gamma} \cup \mathcal{V}$ the set of reservoir pressure, temperature, saturations and mass fractions unknowns is denoted by $X_{\nu} = (P_{\nu}, T_{\nu}, s_{\nu}^{\ell}, s_{\nu}^{g}, c_{\nu}^{\ell}, c_{\nu}^{g})$, where $c_{\nu}^{\alpha}$ is the mass fraction of the water component in phase $\alpha$ used to express the thermodynamical equilibrium. We denote by $X_{\mathcal{D}}$ the set of reservoir unknowns

$$X_{\mathcal{D}} = \{X_{\nu}, \nu \in \mathcal{M} \cup \mathcal{F}_{\Gamma} \cup \mathcal{V}\},$$
and similarly by $P_D$ and $T_D$ the sets of reservoir pressures and temperatures. The set of well bottom hole pressures is denoted by $P_W = \{p_\omega, \omega \in W\}$.

The Darcy fluxes taking into account the gravity term are defined by

$$ egin{cases} V_{K,\nu}^{\alpha}(X_D) = F_{K,\nu}(P_D) - \frac{\rho^{\alpha}(p_K, T_K) + \rho^{\alpha}(p_{\nu}, T_{\nu})}{2} F_{K,\nu}(G_D), & \nu \in \Xi_K, K \in \mathcal{M}, \\ V_{\sigma,s}^{\alpha}(X_D) = F_{\sigma,s}(P_D) - \frac{\rho^{\alpha}(p_{\sigma}, T_{\sigma}) + \rho^{\alpha}(p_s, T_s)}{2} F_{\sigma,s}(G_D), & s \in V_{\sigma}, \sigma \in \mathcal{F}_{\Gamma}, \end{cases} \quad (3.8) $$

where $G_D$ denotes the vector $(g \cdot x_\nu)_{\nu \in \mathcal{M} \cup \mathcal{F}_\Gamma \cup \mathcal{V}}$.

For each Darcy flux, let us define the upwind control volume $cv_{\mu,\nu}^{\alpha}$ such that

$$ cv_{K,\nu}^{\alpha} = \begin{cases} K & \text{if } V_{K,\nu}^{\alpha}(X_D) \ge 0 \\ \nu & \text{if } V_{K,\nu}^{\alpha}(X_D) < 0 \end{cases} \quad \text{for } \nu \in \Xi_K, K \in \mathcal{M}, $$

for the matrix fluxes, and such that

$$ cv_{\sigma,s}^{\alpha} = \begin{cases} \sigma & \text{if } V_{\sigma,s}^{\alpha}(X_D) \ge 0 \\ s & \text{if } V_{\sigma,s}^{\alpha}(X_D) < 0 \end{cases} \quad \text{for } s \in V_{\sigma}, \sigma \in \mathcal{F}_{\Gamma}, $$

for fracture fluxes. Using this upwinding, the mass and energy fluxes are given by

$$ q_{\nu,\nu'}^{\alpha}(X_D) = c_{cv\nu,\nu'}^{\alpha} \frac{\rho^{\alpha}(p_{cv\nu,\nu'}, T_{cv\nu,\nu'})}{\mu^{\alpha}(p_{cv\nu,\nu'}, T_{cv\nu,\nu'})} k_{r, cv\nu,\nu'}^{\alpha} (s_{cv\nu,\nu'}^{\alpha}) V_{\nu,\nu'}^{\alpha}(X_D), $$

$$ q_{\nu,\nu'}^{h2o}(X_D) = \sum_{\alpha \in \{\ell,g\}} q_{\nu,\nu'}^{\alpha}(X_D), $$

$$ q_{\nu,\nu'}^{e}(X_D) = \sum_{\alpha \in \{\ell,g\}} h^{\alpha}(p_{cv\nu,\nu'}, T_{cv\nu,\nu'}) q_{\nu,\nu'}^{\alpha}(X_D) + G_{\nu,\nu'}(T_D). $$

In each control volume $\nu \in \mathcal{M} \cup \mathcal{F}_\Gamma \cup \mathcal{V}$, the mass and energy accumulations are denoted by

$$ A_{\alpha,\nu}(X_\nu) = \phi_\nu \rho^\alpha(p_\nu, T_\nu) s_\nu^\alpha c_\nu^\alpha, $$

$$ A_{h2o,\nu}(X_\nu) = \sum_{\alpha \in \{\ell,g\}} A_{\alpha,\nu}(X_\nu), $$

$$ A_{e,\nu}(X_\nu) = \sum_{\alpha \in \{\ell,g\}} e^\alpha(p_\nu, T_\nu) A_{\alpha,\nu}(X_\nu) + \bar{\phi}_\nu E^r(p_\nu, T_\nu). $$

We can now state the system of discrete equations at each time step $n = 1, \dots, N_{tf}$ which accounts for the mass ($i = h2o$) and energy ($i = e$) conservation equations in each cell $K \in \mathcal{M}$:

$$ R_{K,i}(X_D^n) := \frac{A_{i,K}(X_K^n) - A_{i,K}(X_K^{n-1})}{\Delta t^n} + \sum_{s \in V_K} q_{K,s}^i(X_D^n) + \sum_{\sigma \in \mathcal{F}_\Gamma \cap \mathcal{F}_K} q_{K,\sigma}^i(X_D^n) = 0, \quad (3.9) $$

in each fracture face $\sigma \in \mathcal{F}_\Gamma$:

$$ R_{\sigma,i}(X_D^n) := \frac{A_{i,\sigma}(X_\sigma^n) - A_{i,\sigma}(X_\sigma^{n-1})}{\Delta t^n} + \sum_{s \in V_\sigma} q_{\sigma,s}^i(X_D^n) + \sum_{K \in \mathcal{M}_\sigma} -q_{K,\sigma}^i(X_D^n) = 0, \quad (3.10) $$

and at each node $s \in \mathcal{V} \setminus \mathcal{V}_D$:

$$ R_{s,i}(X_D^n, P_W^n) := \frac{A_{i,s}(X_s^n) - A_{i,s}(X_s^{n-1})}{\Delta t^n} + \sum_{\sigma \in \mathcal{F}_\Gamma, s} -q_{\sigma,s}^i(X_D^n) + \sum_{K \in \mathcal{M}_s} -q_{K,s}^i(X_D^n) \\ + \sum_{\omega \in \mathcal{W} | s \in \mathcal{V}_\omega} q_{s,i}^{r \to \omega}(X_s^n, p_s^{\omega,n}) = 0. \quad (3.11) $$
It is coupled with the well equations for the injection wells $\omega \in W_{inj}$

$$R_{\omega}(X_{D}^{n}, P_{W}^{n}) := -\min\left(\sum_{s \in V_{\omega}} q_{s, h_{2}o}^{r \to \omega}(X_{s}^{n}, p_{\omega}^{n}) - \bar{q}_{\omega}, \bar{p}_{\omega} - p_{\omega}^{n}\right) = 0, \quad (3.12)$$

and for the production wells $\omega \in W_{prod}$

$$R_{\omega}(X_{D}^{n}, P_{W}^{n}) := \min\left(\bar{q}_{\omega} - \sum_{s \in V_{\omega}} q_{s, h_{2}o}^{r \to \omega}(X_{s}^{n}, p_{\omega}^{n}), p_{\omega}^{n} - \bar{p}_{\omega}\right) = 0, \quad (3.13)$$

reformulating respectively (3.6) and (3.7) using the min function.

The system is closed with thermodynamical equilibrium and the sum to one of the saturations

$$\begin{aligned}
R_1(X_\nu^n) &:= c_\nu^{g,n} p_\nu^n - p_{\text{sat}}(T_\nu^n) c_\nu^{\ell,n} &&= 0, \\
R_2(X_\nu^n) &:= \min(s_\nu^{\ell,n}, 1 - c_\nu^{\ell,n}) &&= 0, \\
R_3(X_\nu^n) &:= \min(s_\nu^{g,n}, 1 - c_\nu^{g,n}) &&= 0, \\
R_4(X_\nu^n) &:= s_\nu^{\ell,n} + s_\nu^{g,n} - 1 &&= 0,
\end{aligned} \quad (3.14)$$

at all control volumes $\nu \in \mathcal{M} \cup \mathcal{F}_{\Gamma} \cup \mathcal{V} \setminus \mathcal{V}_{D}$ as well as the Dirichlet boundary conditions

$$X_{\mathbf{s}}^{n} = X_{\mathbf{s},D},$$

for all $\mathbf{s} \in \mathcal{V}_{D}$.

Let us denote by $R_\nu$ the vector $(R_{\nu,i}, i \in \{h_2o, e\}, R_j(X_\nu), j \in \{1, \dots, 4\})$, and let us rewrite the conservation and closure equations (3.9), (3.10), (3.11), (3.12), (3.13), (3.14) as well as the Dirichlet boundary conditions in vector form defining the following nonlinear system at each time step $n = 1, 2, \dots, N_{tf}$

$$\mathbf{0} = \mathcal{R}(X_D, P_W) := \begin{cases} R_s(X_D, P_W), & \mathbf{s} \in \mathcal{V}, \\ R_\sigma(X_D), & \sigma \in \mathcal{F}_\Gamma, \\ R_K(X_D), & K \in \mathcal{M}, \\ R_\omega(X_D, P_W), & \omega \in \mathcal{W}, \end{cases} \quad (3.15)$$

where the superscript $n$ is dropped to simplify the notations and where the Dirichlet boundary conditions have been included at each Dirichlet node $\mathbf{s} \in \mathcal{V}_{D}$ in order to obtain a system size independent of the boundary conditions.

The nonlinear system $\mathcal{R}(X_D, P_W) = 0$ is solved by a Newton-min algorithm [31]. Our implementation is based on an active set method both for the well equations and the thermodynamical equilibrium.

For the well equations, we enforce either the total mass flow rate or the bottom hole pressure at each Newton iterate and use the remaining inequality constraint to switch from prescribed total mass flow rate to prescribed bottom hole pressure and vice versa.

For the thermodynamical equilibrium equations (3.14), the active set method is based on the set of present phases $I_\nu^n = \{\alpha \in \{\ell, g\} | s_\nu^\alpha > 0\}$. Hence we distinguish a two-phase state $I_\nu^n = \{\ell, g\}$, a liquid state $I_\nu^n = \{\ell\}$, and a gas state $I_\nu^n = \{g\}$, for which the system (3.14) reduces to $H_\nu(X_\nu) = 0$ with

$$H_{\nu}(X_{\nu}) = \begin{cases} T_{\nu}^{n} - (p_{\text{sat}})^{-1}(p_{\nu}^{n}), \\ c_{\nu}^{\ell,n} - 1, & \text{if } I_{\nu}^{n} = \{\ell, g\}, \\ c_{\nu}^{g,n} - 1, \\ s_{\nu}^{\ell,n} + s_{\nu}^{g,n} - 1, \end{cases} \quad (3.16)$$
$$ H_{\nu}(X_{\nu}) = \begin{cases} c_{\nu}^{\ell,n} - 1, \\ c_{\nu}^{g,n} - \frac{p_{\text{sat}}(T_{\nu}^{n})}{p_{\nu}^{n}}, \\ s_{\nu}^{\ell,n} - 1, \\ s_{\nu}^{g,n}, \end{cases} \quad \text{if } I_{\nu}^{n} = \{\ell\}, \quad H_{\nu}(X_{\nu}) = \begin{cases} c_{\nu}^{\ell,n} - \frac{p_{\nu}^{n}}{p_{\text{sat}}(T_{\nu}^{n})}, \\ c_{\nu}^{g,n} - 1, \\ s_{\nu}^{\ell,n}, \\ s_{\nu}^{g,n} - 1, \end{cases} \quad \text{if } I_{\nu}^{n} = \{g\}. \quad (3.17) $$

In order to reduce the size of the linear system to be solved at each Newton iteration, subsets $Y_{\nu}$ of two primary unknowns and $Z_{\nu}$ of four secondary unknowns with $X_{\nu} = (Y_{\nu}, Z_{\nu})$, are selected in such a way that $\frac{\partial H_{\nu}}{\partial Z_{\nu}}(Y_{\nu}, Z_{\nu})$ is non singular. In practice the following choice is implemented which depends only on the phase state $I_{\nu}^{n}$:

$$ \bullet \quad Y_{\nu} = (p_{\nu}^{n}, s_{\nu}^{g,n}), \quad Z_{\nu} = (T_{\nu}^{n}, c_{\nu}^{\ell,n}, c_{\nu}^{g,n}, s_{\nu}^{\ell,n}) \quad \text{for } I_{\nu}^{n} = \{\ell, g\}, $$

$$ \bullet \quad Y_{\nu} = (p_{\nu}^{n}, T_{\nu}^{n}), \quad Z_{\nu} = (c_{\nu}^{\ell,n}, c_{\nu}^{g,n}, s_{\nu}^{\ell,n}, s_{\nu}^{g,n}) \quad \text{for } I_{\nu}^{n} = \{\ell\} \text{ or } I_{\nu}^{n} = \{g\}. $$

The secondary unknowns $Z_{\nu}$ are directly eliminated from the nonlinear system using (3.16)-(3.17).

The inequality constraints are then used to switch from two-phase state to a one phase state and vice versa. Let

$$ Y_D = \{Y_{\nu}, \nu \in \mathcal{M} \cup \mathcal{F}_{\Gamma} \cup \mathcal{V}\}, $$

the Jacobian system at each Newton-min iteration is assembled w.r.t. the primary unknowns $Y_D, P_W$ and the mass and energy conservation equations (3.9), (3.10), (3.11), (3.12), (3.13). The cell primary unknowns $Y_K, K \in \mathcal{M}$ are locally eliminated without any additional fill-in before solving the linear system using the GMRES iterative solver preconditioned by a CPR-AMG preconditioner introduced in [32, 43]. This preconditioner combines multiplicatively a parallel algebraic multigrid preconditioner (AMG) [28] 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 parallel implementation is described in [48] and [10]. Let us recall that the distribution of wells to each MPI process $p$ is such that any well with a node belonging to the set of own nodes of $p$ belongs to the set of own and ghost wells of $p$. Then, the set of own and ghost nodes of $p$ is extended to include all the nodes belonging to the own and ghost wells of $p$. These definitions ensure that (i) the local linearized systems can be assembled locally on each process without communication as in [48], and (ii) the pressure drops of the wells can be computed locally on each process without communication. This last property is convenient since the pressure drop is a sequential computation along the well rooted tree. This parallelization strategy of the well model is based on the assumption that the number of additional ghost nodes resulting from the connectivity of the wells remains very small compared with the number of own nodes of the process.

# 4. Numerical results

## 4.1. Numerical convergence for a diphasic vertical well in an homogeneous reservoir

Let us consider the geothermal reservoir defined by the domain $\Omega = (-H, H)^2 \times (0, H_z)$ where $H = 1000$ m and $H_z = 200$ m, including one vertical producer well along the line $\{(x, y, z) \in \Omega | x = y = 0\}$ of radius $r_{\omega} = 0.1$ m. The reservoir is assumed homogeneous with isotropic permeability $K_m = k_m I, k_m = 5 \times 10^{-14}$ m$^2$ and porosity $\phi_m = 0.15$. It is assumed to be initially saturated with pure water in liquid phase. The internal energy, mass density and viscosity of water in the liquid and gas phases are given from [44] by analytical laws as functions of the pressure and temperature. Note that, in both test cases of Subsections 4.1 and 4.2, we make the approximation $h^{\alpha} \sim e^{\alpha}$. The vapour pressure $P_{\text{sat}}(T)$ is given in Pa by the Clausius-Clapeyron equation $p_{\text{sat}}(T) = 100 \exp(46.784 - \frac{6435}{T} - 3.868 \log(T))$. The
thermal conductivity is fixed to $\lambda_m = 2 \text{ W.m}^{-1}.\text{K}^{-1}$, and the rock volumetric heat capacity is given by $C_r = 1.6 \text{ MJ.K}^{-1}.\text{m}^{-3}$ with $E^r(p, T) = C_r T$. The relative permeabilities are set to $k_{r,m}^{\alpha}(s^{\alpha}) = (s^{\alpha})^2$ for both phases $\alpha \in \{\ell, g\}$. The gravity vector is as usual $\mathbf{g} = (0, 0, -g_z)$ with $g_z = 9.81 \text{ m.s}^{-2}$.

The simulation consists in two stages both run on a family of refined uniform Cartesian meshes of size $n_x \times n_y \times n_z$ of the domain $\Omega$ with

$$ (n_x, n_y, n_z) \in \{(10, 10, 5), (20, 20, 10), (40, 40, 20), (80, 80, 40)\}. $$

These meshes are labeled as $\{h_1, \dots, h_4\}$ respectively. The well indexes are computed at each node of the well following [10].

At the first stage, the well is closed and a Dirichlet boundary condition is imposed at the top of the domain prescribing the pressure and the temperature equal to $p_m = 4 \text{ MPa}$ and $T_m = (p_{\text{sat}})^{-1}(p_m) - 1 \text{ K}$; respectively, and homogeneous Neumann boundary conditions are set at the bottom and at the sides of the domain. The choice of the initial temperature $T_m$ just below the saturated vapor temperature is made in order to make the gas phase appear at the beginning of the production during the second stage. This stage is run until the simulation reaches a stationary state with the liquid phase only, a constant temperature and an hydrostatic pressure depending only on the vertical coordinate.

For the second stage, homogeneous Neumann boundary conditions are prescribed at the bottom and at the top of the domain $\Omega$, but Dirichlet boundary conditions for the pressure and temperature are fixed at the sides of the domain to the ones at the end of stage one. The well is set in an open state, i.e., it can produce, and its monitoring conditions are defined by the minimum bottom hole pressure $\bar{p}_{\omega} = 1 \text{ bar}$ (never reached in practice) and the maximum total mass flow rate $\bar{q}_{\omega} = 200 \text{ ton hour}^{-1}$. The second stage is run on the time interval $(0, t_f)$ with $t_f = 30 \text{ days}$.

Figures 4.1, 4.2, and 4.3 exhibit, for the family of refined meshes, the pressure, the temperature and the gas saturation along the well; respectively, at final time $t_f$. The solutions are pretty close for all meshes and exhibit a good convergence behavior.

Figures 4.4, and 4.5 show a close look of the pressure and of the temperature inside the reservoir; respectively, for all meshes at final time $t_f$. It illustrates the cone shaped bubble of gas along the well at the top of the reservoir and demonstrates again the good convergence behavior of the discrete model.

<figure><img src="image_7.png" /><figcaption>FIGURE 4.1. Pressure in Pa along the well at final time on the different meshes.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 4.2. Temperature in °C along the well at final time on the different meshes.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 4.3. Gas saturation along the well at final time on the different meshes.</figcaption></figure>

At each time step, the nonlinear system is solved using a Newton algorithm. The GMRES stopping criterion on the relative residual is fixed to $10^{-8}$. The Newton solver is convergent if the relative residual is lower than $10^{-8}$ as well.

Table 4.1 shows the numerical efficiency of the proposed scheme for all meshes for the second stage of the simulation. We denote by $N_{\Delta t}$ the number of successful time steps, by $N_{\text{Newton}}$ the average number of Newton iterations per successful time step, and by $N_{\text{GMRES}}$ the average number of GMRES iterations per Newton iteration. It exhibits a very good robustness of the Newton solver on the family of refined meshes and a moderate increase of the number of GMRES iterations with the mesh size.
<figure><img src="image_10.png" /><figcaption>FIGURE 4.4. Clip and close look of the gas saturation inside the reservoir at final time for all meshes (cell values).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 4.5. Clip and close look of the temperature in °C inside the reservoir at final time for all meshes (cell values).</figcaption></figure>

Finally, we present in Figure 4.6 the total computational time in hours obtained with the finest mesh $h_4$ for different numbers of MPI processes $N_p = 8, 16, 32, 64$. As usual for this type of simulations, the strong scalability is limited by the AMG preconditioner of the pressure block which requires a sufficiently high number of unknowns per processor to keep a good scalability, corresponding to roughly speaking $4 \cdot 10^4$. This explains the good speed up obtained between 8 and 32 processors whereas the speed up becomes very small between 32 and 64 processors.
<table><caption>TABLE 4.1. Numerical behavior of the second stage of the simulation for different mesh sizes. $N_{\Delta t}$ is the number of successful time steps, $N_{\text{Newton}}$ the average number of Newton iterations per successful time step, and $N_{\text{GMRES}}$ is the average number of GMRES iterations per Newton iteration.</caption><thead><tr><th>Mesh</th><th>#M</th><th>N<sub>Δt</sub></th><th>N<sub>Newton</sub></th><th>N<sub>GMRES</sub></th></tr></thead><tbody><tr><td>h<sub>1</sub></td><td>4000</td><td>134</td><td>1.99</td><td>8.59</td></tr><tr><td>h<sub>2</sub></td><td>32000</td><td>134</td><td>1.74</td><td>9.93</td></tr><tr><td>h<sub>3</sub></td><td>256000</td><td>134</td><td>1.92</td><td>11.75</td></tr><tr><td>h<sub>4</sub></td><td>1848320</td><td>133</td><td>2.22</td><td>15.91</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>FIGURE 4.6. Total computational time vs. number of MPI processes for the second stage simulation on the finest mesh $h_4$.</figcaption></figure>

## 4.2. Study of a high enthalpy reservoir

In this section, we consider a more realistic case built from geological and production data of a field in a volcanic area. The field is a convective dominated system initially in liquid phase, that is crossed by a major normal fault. The reservoir (in blue in Figure 4.7a) is about 500 m thick; it is covered by a weakly permeable clay caprock (in yellow) of 250 m thick, which outcrops at the surface. Below the reservoir is the basement layer (also in yellow).

Figure 4.7b gives the tetrahedral mesh of the domain. The VAG finite volume discretization makes it possible to deal with complex geology including faults and complex well trajectories. The unstructured mesh of 700 000 tetrahedral elements draws on geological horizons. The fault is meshed as a two-dimensional (2D) surface, where the triangular elements are interconnected with the surrounding matrix using conformal meshing. The (one-dimensional) wells are discretized by a set of edges as shown in Figure 4.7b.
<figure><img src="image_13.png" /><figcaption>FIGURE 4.7. Geometry, mesh and wells data for the second numerical test.</figcaption></figure>

The geothermal field is operated using a doublet of two deviated wells, a producer (in green) and an injector (in blue), both of which cross the major fault as shown in Figure 4.7b. The reservoir is assumed homogeneous with an isotropic permeability $K_m = k_m I$, $k_m = 10^{-14} m^2$ and a porosity $\phi_m = 0.05$, while the faulted area has a thickness $d_f = 10$ m, an isotropic permeability $K_f = k_f I$, $k_f = 5 \cdot 10^{-14} m^2$ and a porosity $\phi_f = 0.05$. The caprock and the basement layer are assumed weakly permeable with $k_m = 10^{-19} m^2$. The matrix and fracture thermal conductivities are set to $\lambda_m = \lambda_f = 3 \text{ W.K}^{-1}\text{m}^{-1}$ and the rock energy density is homogeneous for the whole rock mass such that $E^r(p, T) = \rho_r c_p^r T$ with $c_p^r = 1000 \text{ J.kg}^{-1}\text{K}^{-1}$ and $\rho_r = 2600 \text{ kg.m}^{-3}$. The relative permeabilities are set to $k_{r,m}^{\alpha}(s^{\alpha}) = k_{r,f}^{\alpha}(s^{\alpha}) = (s^{\alpha})^2$ for both phases $\alpha \in \{\ell, g\}$.

<figure><img src="image_14.png" /><figcaption>FIGURE 4.8. Initial state dominated by convention. Isotemperature surfaces.</figcaption></figure>

As the previous numerical test, this simulation consists in two stages. The first one acts as a preliminary step where the initial state of the geothermal system, which is already dynamic, is achieved by performing a simulation over a long period (here $10^5$ years) from an hydrostatic pressure state (with 1 bar at the top of the model), and a temperature field increasing linearly with depth (between 30 °C at the top to 280 °C at the bottom). Dirichlet boundary conditions for temperature and pressure are thus imposed at the top and bottom boundaries. No flow and Dirichlet temperature conditions are applied on the lateral boundaries. The initial state obtained is convective; the fluid in the reservoir is in liquid state with a low fraction of gas near the top of the reservoir. Iso-temperature contours are
<figure><img src="image_15.png" /></figure>

FIGURE 4.9. Temperature and saturation after 5 years of production - cells with a gas saturation greater than 0.1 are filled in yellow - the temperature is represented in the fault plane

<figure><img src="image_16.png" /></figure>

FIGURE 4.10. Temperature and saturation after 10 years of production with reinjection during the last 5 years - cells with a gas saturation greater than 0.1 are filled in yellow - the temperature is represented in the fault plane
represented in Figure 4.8 and show the development of convection cells and the influence of the fault, which is a more permeable zone.

Then the second stage begins where the reservoir production starts with steam production at the producer well-head: a flow rate of 250 ton.hr⁻¹ is imposed at the well-head for five years. The same boundary conditions are imposed as in the initial state determination, but the temperature imposed on the lateral boundaries is now given by the average temperature distribution in the rock mass at this initial state. The depletion occuring near the producer well favors the development of a steam cap in the reservoir as well as in the fault zone. Figure 4.9 shows this steam cap: faces in the fault and cells in the reservoir with a gas saturation greater than 0.1 are filled in yellow, while temperature field is also represented on the other faces of the fault plane.

After five years of production and reservoir depletion, half of the fluid produced is reinjected at the injector with a wellhead temperature of 110 °C. During the injection, vapor around the injector condenses and the steam cap generated around the producer is considerably reduced (Figure 4.10).

Figures 4.11a and 4.11b show at a given depth of 455 m respectively the evolution of pressure in the reservoir and in the well and the saturation evolution in the well. Reservoir pressure decreases during the first five years of production, while reinjection of half of the fluid produced during the next five years leads to a pressure build-up in the reservoir (the model is not hydraulically closed). Well pressure follows the same trends. Whereas gas saturation was around 80% during the depletion phase in the well at 455 m depth, injection results in a reduced gas saturation in the well down to approximatley 50% at 455 m depth (Figure 4.11b).

<figure><img src="image_17.png" /><figcaption>FIGURE 4.11. (a) Pressure evolution in the reservoir (red dashed line) and in the well (green line) at 455 m depth. Saturation pressure in the well at 455 m depth is given by the black dotted line. (b): Saturation evolution in the well at 455 m depth.</figcaption></figure>

Table 4.2 shows the numerical efficiency of the proposed scheme for both stages of the simulation and different numbers of MPI processes $N_p = 4, 8, 16$. We use the same notations as in the previous test case and report in addition the total simulation time in hours. These results exhibit the very good robustness of the linear and nonlinear solvers w.r.t. the number of MPI processes. A very good speedup is obtained up to 16 MPI processes verifying that parallel computing makes it possible to have reasonable computation times to model industrial cases such as the one presented in this section.
<table><caption>TABLE 4.2. Numerical behavior of both stages of the simulation for different number of processors. $N_{\Delta t}$ is the number of successful time steps, $N_{\text{Newton}}$ the average number of Newton iterations per successful time step, $N_{\text{GMRES}}$ the average number of GMRES iterations per Newton iteration, and Time (hrs) is the total simulation time in hours.</caption><thead><tr><th>Stage</th><th>$N_p$</th><th>$N_{\Delta t}$</th><th>$N_{\text{Newton}}$</th><th>$N_{\text{GMRES}}$</th><th>Time (hrs)</th></tr></thead><tbody><tr><td rowspan="3">1</td><td>4</td><td>1515</td><td>4.6</td><td>29.3</td><td>98.2</td></tr><tr><td>8</td><td>1507</td><td>4.6</td><td>29.4</td><td>31.9</td></tr><tr><td>16</td><td>1526</td><td>4.6</td><td>30.0</td><td>17.8</td></tr><tr><td rowspan="3">2</td><td>4</td><td>1395</td><td>7.3</td><td>7.7</td><td>65.9</td></tr><tr><td>8</td><td>1367</td><td>7.2</td><td>7.6</td><td>20.2</td></tr><tr><td>16</td><td>1320</td><td>7.2</td><td>7.9</td><td>10.1</td></tr></tbody></table>

## 5. Conclusion

This paper focuses on the numerical modelling of geothermal systems in complex geological settings. The proposed approach is based on unstructured meshes to model complex features such as faults and deviated wells. It solves liquid-vapor two-phase Darcy flows coupled with energy transfers and thermodynamical equilibrium. The use of the hybrid-dimensional polytopal VAG scheme allows to treat physically complex cases, while respecting geometrical constraints. It is thus possible to model explicitly the flows in geological discontinuities, that condition the behavior of most deep geothermal reservoirs during operation. We particularly focus on the well modelling with deviated or multi-branch wells defined as a collection of edges of the mesh with rooted tree data structure. The well-reservoir connection is taken into account by a Peacemann index both for the faults and the rock mass and thus describes the presence of feed zones in injection and production wells. By using an explicit calculation of the mean density in the pressure drop, the well model reduces to a single equation with only one well implicit unknown, the bottom hole pressure, implicitly coupled to the reservoir system.

The robustness of this approach in the case of liquid-vapor flow is illustrated by two test cases. We first present a numerical convergence test of a two-phase vertical production well model in a simplified reservoir geometry. Then, we study the behavior of a high enthalpy faulted reservoir using a doublet of two deviated wells crossing the fault. The use of efficient linear and nonlinear parallel solvers ensures acceptable computation times on real case studies.

An improved model of cross flows between the well and the reservoir will be investigated in the near future. Industrial studies of high and medium enthalpy geothermal reservoirs are currently under way with the approach proposed in this paper.

## References

[1] I. Aavatsmark and R. A. Klausen. Well Index in Reservoir Simulation for Slanted and Slightly Curved Wells in 3D Grids. *SPE Journal*, 8(01):41–48, 2003.

[2] Joubine Aghili, Jean-Raynald de Dreuzy, Roland Masson, and Laurent Trenty. A hybrid-dimensional compositional two-phase flow model in fractured porous media with phase transitions and Fickian diffusion. *J. Comput. Phys.*, 441:110452, 2021.

[3] Raheel Ahmed, Michael G. Edwards, Sadok Lamine, Bastiaan A. H. Huisman, and Mayur Pal. Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model. *J. Comput. Phys.*, 284:462–489, 2015.
[4] Raheel Ahmed, Michael G. Edwards, Sadok Lamine, Bastiaan A. H. Huisman, and Mayur Pal. Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model. *J. Comput. Phys.*, 303:470–497, 2015.

[5] Clarisse Alboin, Jérôme Jaffré, J. E. Roberts, and Christophe Serres. Modeling fractures as interfaces for flow and transport in porous media. In *Fluid flow and transport in porous media: mathematical and numerical treatment*, volume 295 of *Contemporary Mathematics*, pages 13–24, 2002.

[6] Philippe Angot, Franck Boyer, and Florence Hubert. Asymptotic and numerical modelling of flows in fractured porous media. *ESAIM, Math. Model. Numer. Anal.*, 43(2):239–275, 2009.

[7] Paola F. Antonietti, Luca Formaggia, Anna Scotti, Marco Verani, and Nicola Verzott. Mimetic finite difference approximation of flows in fractured porous media. *ESAIM M2AN*, 50:809–832, 2016.

[8] Z P Aunzo, G Bjornsson, and G S Bodvarsson. Wellbore Models GWELL, GWNACL, and HOLA, user's guide. Technical Report LBL-31428, Earth Sciences Division, Lawrence Berkeley National Laboratory, University of California, 1991.

[9] Khalid Aziz and A. Settari. *Petroleum Reservoir Simulation*. Applied Science Publishers, 1979.

[10] Laurence Beaude, Thibaud Beltzung, Konstantin Brenner, Simon Lopez, Roland Masson, Farid Smai, Jean-frédéric Thebault, and Feng Xing. Parallel Geothermal Numerical Model with Fractures and Multi-Branch Wells. *ESAIM: Proc.*, 63:109–134, 2018.

[11] Laurence Beaude, Konstantin Brenner, Simon Lopez, Roland Masson, and Farid Smai. Non-isothermal compositional liquid gas Darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations. *Comput. Geosci.*, 23(3):443–470, 2019.

[12] I. I. Bogdanov, V. V. Mourzenko, J.-F. Thovert, and P. M. Adler. Two-phase flow through fractured porous media. *Phys. Rev. E*, 68(2), 2003.

[13] Konstantin Brenner, Mayya Groza, Cindy Guichard, Gilles Lebeau, and Roland Masson. Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media. *Numer. Math.*, 134(3):569–609, 2016.

[14] Konstantin Brenner, Mayya Groza, Cindy Guichard, and Roland Masson. Vertex Approximate Gradient scheme for hybrid-dimensional two-phase Darcy flows in fractured porous media. *ESAIM, Math. Model. Numer. Anal.*, 2(49):303–330, 2015.

[15] Konstantin Brenner, Mayya Groza, Laurent Jeannin, Roland Masson, and Jeanne Pellerin. Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media. *Comput. Geosci.*, 21(5):1075–1094, 2017.

[16] Konstantin Brenner, J. Hennicker, Roland Masson, and P. Samier. Gradient discretization of hybrid-dimensional Darcy flow in fractured porous media with discontinuous pressures at matrix-fracture interfaces. *IMA J. Numer. Anal.*, 2016.

[17] Konstantin Brenner, Julian Hennicker, and Roland Masson. *Nodal Discretization of Two-Phase Discrete Fracture Matrix Models*, pages 73–118. Springer, 2021.

[18] Konstantin Brenner, Julian Hennicker, Roland Masson, and Pierre Samier. Hybrid-dimensional modelling of two-phase flow through fractured porous media with enhanced matrix fracture transmission conditions. *J. Comput. Phys.*, 357:100–124, 2018.

[19] Florent Chave, Daniele A. Di Pietro, and Luca Formaggia. A Hybrid High-Order Method for Darcy Flows in Fractured Porous Media. *SIAM J. Sci. Comput.*, 40(2):A1063–A1094, 2018.

[20] Florent Chave, Daniele A. Di Pietro, and Luca Formaggia. A Hybrid High-Order method for passive transport in fractured porous media. *GEM - International Journal on Geomathematics*, 10(1):12, 2019.

[21] Z. Chen and Y. Zhang. Well flow models for various numerical methods. *J. Numer. Anal. Model.*, 6:375–388, 2009.
[22] Robert Eymard, Cindy Guichard, and Raphaèle Herbin. Small-stencil 3D schemes for diffusive flows in porous media. *ESAIM, Math. Model. Numer. Anal.*, 46(2):265-290, 2012.

[23] I. Faille, A. Fumagalli, J. Jaffré, and J. E. Roberts. Model reduction and discretization using hybrid finite volumes of flow in porous media containing faults. *Comput. Geosci.*, 20:317-339, 2016.

[24] Eric Flauraud, Frédéric Nataf, Isabelle Faille, and Roland Masson. Domain decomposition for an asymptotic geological fault modeling. *C. R. Méc. Acad. Sci. Paris*, 331(12):849-855, 2003.

[25] Ingeborg G. Gjerde, Kundan Kumar, and Jan M. Nordbotten. A singularity removal method for coupled 1D-3D flow models. *Comput. Geosci.*, 24(2):443-457, 2020.

[26] S. Granet, P. Fabrie, P. Lemonnier, and M. Quintard. A two-phase flow simulation of a fractured reservoir using a new fissure element method. *Journal of Petroleum Science and Engineering*, 32(1):35-52, 2001.

[27] H. Haegland, A. Assteerawatt, H. K. Dahle, G. T. Eigestad, and R. Helmig. Comparison of cell- and vertex-centered discretization methods for flow in a two-dimensional discrete-fracture-matrix system. *Adv. Water Resources*, 32:1740-1755, 2009.

[28] Van Emden Henson and Ulrike Meier Yang. BoomerAMG: A parallel algebraic multigrid solver and preconditioner. *Appl. Numer. Math.*, 41(1):155-177, 2002.

[29] Hussein Hoteit and Abbas Firoozabadi. An efficient numerical model for incompressible two-phase flow in fractured media. *Adv. Water Resources*, 31(6):891-905, 2008.

[30] M. Karimi-Fard, L. J. Durlofsky, and K. Aziz. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. *SPE Journal*, 9(02):227-236, 2004.

[31] S. Kräutle. The semi-smooth Newton method for multicomponent reactive transport with minerals. *Adv. Water Resources*, 34:137-151, 2011.

[32] Sébastien Lacroix, Yuri V. Vassilevski, and Mary F. Wheeler. Decoupling preconditioners in the implicit parallel accurate reservoir simulator (IPARS). *Numer. Linear Algebra Appl.*, 8(8):537-549, 2001.

[33] S. Livescu, L. J. Durlofsky, K. Aziz, and J. C. Ginestra. A fully-coupled thermal multiphase wellbore flow model for use in reservoir simulation. *Journal of Petroleum Science and Engineering*, 71(3):138-146, 2010. Fourth International Symposium on Hydrocarbons and Chemistry.

[34] Vincent Martin, Jérôme Jaffré, and Jean E. Roberts. Modeling fractures and barriers as interfaces for flow in porous media. *SIAM J. Sci. Comput.*, 26(5):1667-1691, 2005.

[35] Stephan Konrad Matthai, Andrey A. Mezentsev, and Mandefro Belayneh. Finite element - node-centered finite-volume two-phase-flow experiments with fractured rock represented by unstructured hybrid-element meshes. *SPE Reservoir Evaluation & Engineering*, 10(06):740-756, 2007.

[36] Jorge E. P. Monteagudo and Abbas Firoozabadi. Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects. *SPE Journal*, 12(03):355-366, 2007.

[37] J. M. Nordbotten, W. M Boon, A. Fumagalli, and E. Keilegavlen. Unified approach to discretization of flow in fractured porous media. *Comput. Geosci.*, 23:225-237, 2019.

[38] D. W. Peaceman. Interpretation of Well-Block Pressures in Numerical. *Reservoir Simulation Symposium Journal SEPJ*, pages 183-194, 1978.

[39] D. W. Peaceman. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Nonsquare Grid Blocks and Anisotropic Permeability. *Reservoir Simulation Symposium Journal SEPJ*, pages 531-543, 1983.

[40] Karsten Pruess, Curt Oldenburg, and George Moridis. TOUGH2 user's guide, version 2. Technical report, Earth Sciences Division, Lawrence Berkeley National Laboratory, University of California, 1999. Backup Publisher: Earth Sciences Division, Lawrence Berkeley National Laboratory, University of California Issue: LBNL-43134 Volume: LBNL-43134.
[41] Volker Reichenberger, Hartmut Jakobs, Peter Bastian, and Rainer Helmig. A mixed-dimensional finite volume method for two-phase flow in fractured porous media. *Adv. Water Resources*, 29(7):1020-1036, 2006.

[42] T. H. Sandve, I. Berre, and J. M. Nordbotten. An efficient multi-point flux approximation method for Discrete Fracture-Matrix simulations. *J. Comput. Phys.*, 231(9):3784-3800, 2012.

[43] R. Scheichl, R. Masson, and J. Wendebourg. Decoupling and block preconditioning for sedimentary basin simulations. *Comput. Geosci.*, 7(4):295-318, 2003.

[44] E. Schmidt. *Properties of water and steam in S.I. units*. Springer, 1969.

[45] Hua Shi, Jonathan A. Holmes, Louis J. Durlofsky, Khalid Aziz, Luis Diaz, Banu Alkaya, and Gary Oddie. Drift-Flux Modeling of Two-Phase Flow in Wellbores. *SPE Journal*, 10(01):24-33, 2005.

[46] Xavier Tunc, Isabelle Faille, Thierry Gallouët, Marie Christine Cacas, and Pascal Havé. A model for conductive faults with non-matching grids. *Comput. Geosci.*, 16(2):277-296, 2012.

[47] Christian Wolfsteiner, Louis J. Durlofsky, and Khalid Aziz. Calculation of Well Index for Nonconventional Wells on Arbitrary Grids. *Comput. Geosci.*, 7(1):61-82, 2003.

[48] F. Xing, R. Masson, and S. Lopez. Parallel numerical modeling of hybrid-dimensional compositional non-isothermal Darcy flows in fractured porous media. *J. Comput. Phys.*, 345:637-664, 2017.

[49] Feng Xing, Roland Masson, and Simon Lopez. Parallel Vertex Approximate Gradient discretization of hybrid-dimensional Darcy flow and transport in discrete fracture networks. *Comput. Geosci.*, 2016.