# Geothermal modeling in complex geological systems with ComPASS

Antoine Armandine Les Landes, Laurence Beaude, D. Castanon Quiroz, Laurent Jeannin, Simon Lopez, Farid Smaï, Théophile Guillon, R. Masson

► To cite this version:

Antoine Armandine Les Landes, Laurence Beaude, D. Castanon Quiroz, Laurent Jeannin, Simon Lopez, et al.. Geothermal modeling in complex geological systems with ComPASS. Computers & Geosciences, 2025, 194, pp.105752. 10.1016/j.cageo.2024.105752. hal-04801670

HAL Id: hal-04801670

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

Submitted on 25 Nov 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
<figure><img src="image_2.png" /></figure>

Research paper

# Geothermal modeling in complex geological systems with ComPASS

A. Armandine Les Landes<sup>a,*</sup>, L. Beaude<sup>a</sup>, D. Castanon Quiroz<sup>b</sup>, L. Jeannin<sup>c</sup>, S. Lopez<sup>a</sup>, F. Smai<sup>a</sup>, T. Guillon<sup>a</sup>, R. Masson<sup>d</sup>

## ARTICLE INFO

**Keywords:**

Geothermal reservoir modeling
Multi-component multi-phase flow
Unstructured meshes
Faulted and fractured reservoir
Deviated wells

## ABSTRACT

In deep geothermal reservoirs, faults and fractures play a major role, serving as regulators of fluid flow and heat transfer while also providing feed zones for production wells. To accurately model the operation of geothermal fields, it is necessary to explicitly consider objects of varying spatial scales, from the reservoir scale itself, to that of faults and fractures, down to the scale of the injection and production wells.

Our main objective in developing the ComPASS geothermal flow simulator, was to take into account all of these geometric constraints in a flow and heat transfer numerical model using generic unstructured meshes. In its current state, the code provides a parallel implementation of a spatio-temporal discretization of the non-linear equations driving compositional multi-phase thermal flows in porous fractured media on unstructured meshes. It allows an explicit discretization of faults and fractures as 2D hybrid objects, embedded in a 3D matrix. Similarly, wells are modeled as one dimensional graphs discretized by edges of the 3D mesh which allows arbitrary multi-branch wells. The resulting approach is particularly flexible and robust in terms of modeling.

Its practical interest is demonstrated by two case studies in high-energy geothermal contexts.

## 1. Introduction

Geothermal resources are already widely available in areas with volcanic activity and in sedimentary basins and their exploitation could significantly contribute to the decarbonization of our economy (Hirschberg et al., 2014). In contrast to other renewable energy sources, geothermal energy is weather-independent and can provide both electricity and heat, as well as value-added mineral extraction. Geothermal power production provides reliable generation with high plant efficiency, low greenhouse gas emissions and a small ecological footprint. Moreover, it is usually a long-lasting sustainable source when properly managed (International Renewable Energy Agency and International Geothermal Association, 2023).

Numerical modeling has been one of the cornerstones of such management for decades and has become an essential tool at all stages of exploration and development of geothermal projects (O'Sullivan et al., 2001; O'Sullivan and O'Sullivan, 2016; Nugraha et al., 2022). The usual starting point is a conceptual model that often relies on 3D

geological modeling capabilities. Then, flow modeling is used during the exploration phases to understand the natural resource, to evaluate the geothermal potential, to validate conceptual hypotheses, and on the other hand, during the development of geothermal assets to predict production flow rates (e.g. among many others: Ingebritsen et al., 2010; Daniilidis et al., 2021; Jalilinasrabady et al., 2021; O'Sullivan et al., 2009; O'Sullivan and O'Sullivan, 2016; Nugraha et al., 2022...). Not only do the models provide a coherent vision of the geothermal system by integrating all the data and measurements, but above all they represent a quantitative tool for decision making and planning production and development of the resource.

By many aspects, integrating multi-disciplinary contributions in a shared consistent numerical model of the subsurface remains a complex but worthwhile challenge to mitigate risks (Lopez et al., 2017). Passing from a 3D geological model to an hydrothermal numerical model while preserving realistic geometries in flow simulations is a crucial step of the geothermal modeling workflow (Huang et al., 2022; Nugraha

<footnote><sup>a</sup> BRGM, 3 avenue Claude-Guillemin, BP 36009, Cedex 2, 45060 Orléans, France</footnote>

<footnote><sup>b</sup> 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., Mexico</footnote>

<footnote><sup>c</sup> STORENGY, 12 rue Raoul Nordling - Djinn - CS, 70001 92274 Bois Colombes Cedex, France</footnote>

<footnote><sup>d</sup> Université Côte d'Azur, Inria, CNRS, LJAD, UMR 7351 CNRS, team Coffee, Parc Valrose, Cedex 02 06108 Nice, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: a.armandineLesLandes@brgm.fr (A. Armandine Les Landes), l.beaude@brgm.fr (L. Beaude), daniel.castanon@iimas.unam.mx (D. Castanon Quiroz), laurent.jeannin@storengy.com (L. Jeannin), s.lopez@brgm.fr (S. Lopez), f.smai@brgm.fr (F. Smai), t.guillon@brgm.fr (T. Guillon), roland.masson@univ-cotedazur.fr (R. Masson).</footnote>
et al., 2022). It is not so obvious to overcome. This observation was the main motivation for putting generic unstructured meshes and complex geometries at the heart of the development of the ComPASS geothermal flow simulator with the ambition to improve the accuracy and the representativeness of numerical models and simulation workflows and facilitate the back and forth between geological modeling and flow modeling.

This paper presents the ComPASS code which is an open-source, massively parallel, multi-phase, multi-component flow simulator. The paper focuses on high energy geothermal resources and how generic tetrahedral meshes can be used to exactly match any geological structure including two dimensional elements (geological interfaces, faults or fractures...) or one dimensional sharp features (surface intersections, well trajectories...). Then, adapted numerical schemes can be used in these n-dimensional elements to discretize mass and energy transfers and solve the usual, non-linear balance equations involved in geothermal reservoir modeling. In the following, we may use the term fracture for fracture or fault for brevity regardless of their geological difference.

This mixed-dimensional approach is especially useful to take into account fractures that exert a dominant control on subsurface geothermal flows and associated energy transfers. In tectonically active areas, fault zones act as permeable pathways, including high-temperature magmatic environments (Grant and Bixley, 2011; Ingebritsen et al., 2006). But this is also the case in other contexts. Concerning enhanced Geothermal Systems (EGS), where hydraulic stimulation ensures connection with a fault network, flow takes place mostly in a fracture network (Blaisonneau et al., 2021). Role of discontinuities has also been recognized as a significant factor in determining the thermal structure of sedimentary basins, as they can act as conduits linking different aquifer levels (Magri et al., 2010; Person et al., 2012; Simms and Garven, 2004). Geothermal wells are located, so that they intersect faults and fractures (Grant and Bixley, 2011). Indeed, the flow dynamics are primarily influenced by the connectivity and conductivity of the fracture network, interacting with the surrounding matrix medium.

In terms of software specifications, modeling the non-linear behavior and phase-transitions (boiling and condensation) of high-temperature geothermal brines requires quite robust algorithms. The fact that the geological media in which they flow is discontinuous because of the presence of fractures and shows abrupt petrophysical parameters variations over several orders of magnitude, with zonal anisotropies, makes accurate modeling even more challenging. For many years, Tough2 (Pruess et al., 1999; Finsterle et al., 2014) has been a de-facto industry standard (O'Sullivan and O'Sullivan, 2016) mainly thanks to its numerous and versatile EOS modules tailored to a large number of physics. Its large community of users contributed to this collection including for example modules for supercritical water conditions (Croucher and O'Sullivan, 2008) or brine with dissolved gas (Battistelli et al., 1997).

However, many software, including Tough2, use the classical Two Point Flux Approximation (TPFA) for spatially discretizing exchange terms in balance equations. TPFA relies on orthogonality conditions and limits the types of acceptable meshes, typically rectangular parallelepipedic boxes in 3D (Eymard et al., 2014). Unfortunately, using coupled finite volume schemes on such structured grids can lead to the Grid Orientation Effect (Eymard et al., 2013). Voronoi diagrams built from unstructured meshes offer more flexibility while maintaining orthogonality (Freeman et al., 2014), but they make the discretization of geological interfaces and intersections challenging, complicating the specification of domain properties and boundary conditions. All code relying on Finite Differences such as SHEMAT (Clauser, 2003) or HYDROTHERM (Hayba and Ingebritsen, 1994) suffer from these constraints as well. By design, Finite Elements based modeling suite handle nicely simplicial meshes, but tools like FEFLOW (Blöcher et al., 2010), COMSOL (Guillou-Frottier et al., 2013; Taillefer et al., 2018) or OpenGeoSys (Kolditz et al., 2012) hardly take into account phase

change with boiling or condensation and are consequently not suitable for the modeling of high temperature geothermal resources. To handle unstructured meshes, Coumou et al. (2008) implemented an operator splitting approach in the Complex System Modelling Platform (CSMP++) which involves a thermal equilibrium step. This approach has then been adapted to model supercritical geothermal resources (Weis et al., 2014).

Finally, it is quite common that the discretization of complex 3D geological models with unstructured good quality meshes produces hundreds of thousands to millions of cells, especially when considering regional scale models with many interfaces. In such situation, affordable computation time implies, possibly massive, parallelization of the simulation code. Among the progeny of the Tough2 family of codes, Tough3 (Jung et al., 2017) and the recent Waiwera platform (Croucher et al., 2020) offer such capabilities but suffer from the TPFA scheme limitations. CSMP++ also comes with High Performance Computing (HPC) capabilities and the DuMuX platform is based on the modular parallel C++ framework Dune (Distributed and Unified Numerics Environment) and handle multi-phase non-isothermal flows (Koch et al., 2021).

The rest of the paper is structured in the following manner. The Section 1 presents the physical model of multi-phase compositional geothermal flows in fractured geothermal reservoirs and the specific numerical discretizations used in the ComPASS code. In addition, the approaches used in terms of meshing, parallelization and numerical performances are summarized. The Section 2 is devoted to numerical tests that demonstrate the benefits of the approach on two case studies inspired from classical high energy geothermal contexts: a liquid dominated reservoir crossed by major faults and a steam dominated reservoir.

## 2. Model description

In this section, we briefly review the compositional multi-phase model currently implemented in the ComPASS code as well as related numerical aspects. We refer the interested reader to already published material, especially Xing et al. (2017), Beaude et al. (2018) and Armandine Les Landes et al. (2023), and to the publicly available code documentation.

### 2.1. Multi-phase multi-component model

#### 2.1.1. Physical system description

The description of the physical system is based on a Coats' type formulation (Coats, 1989) also known as natural variables formulation. Considering arbitrary sets of components $C$ and phases $P$, we introduce $C^a \neq \emptyset$, $C^a \subset C$ the set of components that can be present in the phase $\alpha \in P$ and reciprocally $P_i \neq \emptyset$, $P_i \subset P$ the set of phases that can contain the component $i \in C$. Finally, $\mathcal{Q}$ is a finite set of labels, called contexts, that are used to describe different physical states. A simple example of contexts can be $\mathcal{Q} = \{liquid, gas, diphasic\}$. Building on the previous notations, each context $Q \in \mathcal{Q}$ is then associated to $P_Q$, respectively $C_Q$, the set of phases, respectively components, that can be present when considering context $Q$. It follows that $\overline{C_Q} = C \setminus C_Q$ may be non-empty and designates the set of components that cannot be present in context $Q$. This framework is very generic and can be used to describe very complex fluids.

We assume thermal equilibrium between phases so that the thermodynamic characteristics of each phase $\alpha \in P$ are influenced by the phase pressure $P^\alpha$, the temperature $T$ and the phase molar fractions $C^\alpha = (C_i^\alpha)_{i \in C^\alpha}$. Additionally, $S^\alpha$ will denote the phase saturation and $n_i$ the number of moles of the component $i \in C$ per unit pore volume.

Finally, the formulation uses the following set of unknowns:

$$X = ((P^\alpha)_{a \in \mathcal{R}_Q}, T, (C^\alpha)_{a \in \mathcal{R}_Q}, (S^\alpha)_{a \in \mathcal{R}_Q}, (n_i)_{i \in \overline{C_Q}}, Q). \quad (1)$$
Then, for each context Q the quantity of matter of component i ∈ CQ writes:

$$n_i = \sum_{a \in \mathcal{R}_Q \cap \mathcal{R}_i} \xi^a (P^a, T, C^a) S^a C_i^a$$

with $ξ^a$ the phase molar density, and the fluid internal energy writes:

$$E = \sum_{a \in \mathcal{R}_Q} \xi^a (P^a, T, C^a) S^a e^a (P^a, T, C^a)$$

with $e^a$ the phase molar internal energy.

### 2.1.2. Fluxes and conservation laws

Molar and energy conservation provides $\#C + 1$ balance equations:

$$\left\{ \begin{array}{l} \phi \partial_t n_i + \nabla \cdot \mathbf{q}_i = 0, \quad i \in \mathcal{C}, \\ \phi \partial_t E + (1 - \phi) \partial_t E_r + \nabla \cdot \mathbf{q}_e = 0 \end{array} \right. \qquad (2)$$

where $\phi$ is the rock porosity assumed constant in time, $E_r$ is the fluid rock energy density, $\mathbf{q}_i$ and $\mathbf{q}_e$ respectively denote the molar flux of component i and the total energy flux comprising advective and diffusive transfers:

$$\mathbf{q}_i = \sum_{a \in \mathcal{R}_Q \cap \mathcal{R}_i} C_i^a \xi^a (P^a, T, C^a) \mathbf{V}^a,$$

$$\mathbf{q}_e = \sum_{a \in \mathcal{R}_Q} h^a (P^a, T, C^a) \xi^a (P^a, T, C^a) \mathbf{V}^a - \lambda \nabla T,$$

$$\mathbf{V}^a = -\mathbf{K} \frac{k_r^a(S^a)}{\mu^a(P^a, T, C^a)} (\nabla P^a - \rho^a(P^a, T, C^a)\mathbf{g})$$

which involve the generalized Darcy velocity $\mathbf{V}^a$, the relative permeability $k_r^a$, the dynamic viscosity $\mu^a$, the mass density $\rho^a$, the molar enthalpy $h^a$ of the phase $\alpha$, $\lambda$ the equivalent thermal conductivity of the rock/fluid mixture and $\mathbf{K}$ the intrinsic rock permeability tensor.

To solve for the unknowns (1), the system (2) is closed by a *flash calculation* that is a fixed point equation which tracks possible context switches due to phase transitions:

$$Q = Q_{\text{flash}}(X). \qquad (3)$$

and a set of local closure laws that depend on the current context:

$$\left\{ \begin{array}{l} \displaystyle \sum_{a \in \mathcal{R}_Q} S^a = 1, \\ \displaystyle \sum_{i \in \mathcal{C}^a} C_i^a = 1, \quad \forall \alpha \in \mathcal{R}_Q, \\ \displaystyle P^a - P^\beta = P_c^{a, \beta}(X), \quad \forall (\alpha, \beta) \in \mathcal{R}_Q^2, \alpha \neq \beta, \\ \displaystyle f_i^a(P^a, T, C^a) = f_i^\beta(P^\beta, T, C^\beta), \quad \forall (\alpha, \beta) \in \mathcal{R}_Q^2, \alpha \neq \beta, i \in \mathcal{C}_Q \end{array} \right. \qquad (4)$$

where $P_c^{a, \beta}(X)$ is a capillary pressure and $f_i^a$, the fugacity of component i in phase $\alpha$, is used to write the thermodynamic equilibrium between phases.

## 2.2. Hybrid model with lower dimension elements

In many geothermal applications, the flow is dominated by the connectivity and conductivity of major discontinuities and feedzones are found at their intersections with wells. The main philosophy in ComPASS is to rely on the versatility of unstructured meshes to exactly match any geological structure and two dimensional geometrical objects and boundaries (geological interfaces, fractures...) or one dimensional sharp features (surface intersections, well trajectories...). Depending on the geological model, producing this kind of conformal meshes with relatively good quality elements can be challenging and it is a crucial step in the whole workflow (Balarac et al., 2022). An alternative, would be to model lower dimensional entities as with distinct meshes and consider them as *embedded features* which then adds a complexity in terms of coupling (Cusini et al., 2021).

### 2.2.1. Discrete fractures

ComPASS implements a Discrete Fracture Matrix (DFM) method, where fractures are discretized by a subset of the 3D mesh facets which incorporate a lower-dimensional physical model. The latest is created by averaging the Eqs. (2), as well as the unknowns (1) related to fractures across their respective widths (Xing et al., 2017).

The transmission conditions that exist at the interfaces between the matrix and fractures depend on extra physical assumptions regarding their drain or barrier behavior of fractures. When they are conductive in terms of both permeability and thermal conductivity, pressure and temperature continuity can be assumed as matrix fracture transmission conditions in single-phase flows (Serres et al., 2002; Brenner et al., 2016). This approach has been expanded to encompass two-phase Darcy flows (Bogdanov et al., 2003; Reichenberger et al., 2006; Brenner et al., 2015, 2017) as well as multi-phase compositional thermal Darcy flows (Xing et al., 2017).

### 2.2.2. Multi-branch well model

One of the main challenges in well modeling in field simulation is the significant difference between the kilometer scale reservoir and the decimeter scale radius of the wellbore. The well geometry cannot be resolved explicitly by the mesh, so the well is modeled as a line source defined by a one-dimensional graph with a rooted tree structure. This graph can accurately depict sloped and multi-branch wells, and it is discretized by a subset of edges on which the mesh is based.

The flow within the well is locally connected to both the 3D matrix surrounding the well and the fractures that intersect the well. The mass and heat exchanges for each node of the well in contact with a fracture or the rock mass are modeled using the Peaceman approach (Peaceman, 1978, 1983). This approach is widely used in reservoir simulation and involves discretizing the Darcy or Fourier fluxes between the reservoir and the well using a two-point flux approximation. The transmissivity takes into account the unresolved pressure or temperature singularity within the well. At the discrete level, the well index or Peaceman's index depends on the type of cell, well radius and geometry, and the discretization scheme used (Wolfsteiner et al., 2003; Aavatsmark and Klausen, 2003; Chen and Zhang, 2009; Yapparova et al., 2022).

In the present work, we use an extension of the simple well model detailed by Beaude et al. (2018) to two-phase flows. The wellbore flow is assumed stationary and by explicitly calculating the pressure drop, the well model is simplified to a single equation with a single implicit unknown which represents the reference pressure of the well. Pressures and saturations along the well are then deduced from this reference pressure and the explicit expression of pressure losses taking into account gravity. To monitor the well, complementary conditions are prescribed between the mass flow rate and the wellhead pressure. For production (resp. injection) wells, they are based on a maximum (resp. minimum) mass flow rate and a minimum (resp. maximum) wellhead pressure. Though a single unknown is introduced, all nodes along the well path are linked to the well reference pressure and this results in additional connectivity in the system, which needs to be taken into account when parallelizing the code (Beaude et al., 2018).

## 2.3. Spatial discretization

The discretization of the spatial terms involved in the continuous model introduced previously relies on the Vertex Approximation Gradient (VAG) finite volume scheme (Eymard et al., 2012). This scheme belongs to the family of gradient schemes (Droniou et al., 2010) and is particularly adapted to the resolution of Darcy's law on conformal polyhedral meshes. The VAG discretization has been adapted to mixed-dimensional modeling (Brenner et al., 2016, 2015). Two-phase Darcy flows discretization was introduced by Brenner et al. (2015) and generalized to multi-phase multi-component flows by Xing et al. (2017). We approximate Darcy fluxes mobilities using a phase potential upwind scheme.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** For a cell $K$ and a fracture face $\sigma$ (in bold), examples of VAG degrees of freedom $u_K, u_S, u_{S'}, u_\sigma$ and VAG fluxes $F_{K,\sigma}, F_{K,S}, F_{K,S'}, F_{\sigma,S}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Example of control volumes at cells, fracture face, and nodes, in the case of two cells $K$ and $L$ split 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>

Fig. 1 details the different degrees of freedom introduced by the VAG scheme and their associated fluxes. It shows one polyhedral cell $K$ (blue) with one fracture face $\sigma$ (red). The 3D matrix fluxes (blue arrows) connect $K$ to cell nodes (e.g. $s, s', s'', s'''$) and fracture faces when present (e.g. $\sigma$). The 2D fracture fluxes (red arrows) connect each fracture face $\sigma$ to its nodes (e.g. $s, s'$). Finally, nodes that are on a well path (e.g. $s, s''$) are also involved in two-point fluxes with the well bore flow (green arrow).

Then, the volume of each cell is partitioned so that a control volume is associated with each degree of freedom. Compared to other approaches like Control Volume Finite Element Methods (CVFE), the VAG scheme's crucial advantage is the flexibility it offers in selecting the control volumes. We exploit it to ensure that there is no mixing of heterogeneous properties within each control volume nor between matrix and fracture media (Fig. 2) which limits numerical diffusion.

## 2.4. Numerical aspects

### 2.4.1. Non-linear solver

A fully implicit Euler scheme is then used to integrate the semi-discrete problem obtained in Section 2.3 and results in a non-linear problem with unknowns consisting in reservoir unknowns (1) and a reference pressure for each of the wells. Because of the flash fixed point equations (Eq. (3)) at each degree of freedom, the problem is solved using an active set Newton-Raphson algorithm (Coats, 1989).

As VAG matrix fluxes are expressed linearly and locally to each cell (Fig. 1), the rows of the Jacobian matrix involved in each Newton iteration and corresponding to the balance equations of the control volumes associated with cells can be eliminated without any fill-in.

Thanks to this elimination, and the use of vertices as main degrees of freedom during the resolution step, the VAG scheme maintains a comparable cost to nodal methods on unstructured meshes though it introduces many degrees of freedom (Brenner et al., 2015). It is particularly suitable to work with tetrahedral meshes that have much more cells than vertices.

The size of the Jacobian matrix is also controlled using the discrete version of the local closure Eqs. (4) to divide the reservoir unknowns into primary and secondary unknowns. The latter are expressed locally as a function of the former and eliminated from the Jacobian matrix.

An iterative solver, typically GMRES, is used to solve the resulting ill-conditioned linear system. A preconditioner is applied to the solver, which is adapted to the pressure unknown's elliptic or parabolic nature and to the coupling with the remaining hyperbolic or parabolic unknowns. The CPR-AMG preconditioner is considered to be one of the most efficient preconditioners for such systems (Lacroix et al., 2001; Scheichl et al., 2003). We refer to Xing et al. (2017) and Beaude et al. (2018) for the full details of its implementation.

### 2.4.2. Parallel implementation

The assembly and resolution of linear systems, which are involved in solving the fully coupled non-linear problem, are performed in parallel using the Single Program Multiple Data (SPMD) paradigm. The rows of the global system are distributed in a well balanced manner to available processes (Xing et al., 2017). To be able to assemble these rows, degrees of freedom holding system unknowns (nodes $s$, fracture faces $\sigma$, cells $K$ and wells $w$) must also be distributed between processes, minimizing the need for communications.

At the beginning of the simulation, the mesh's cells set is partitioned using the METIS library (Karypis and Kumar, 1998) and each resulting subdomain is linked to a unique process rank. Then, node and fracture degrees of freedom are distributed to processes according to the cells they are connected to. Elements belonging to the interior of a cell subdomain are unambiguous attributed to the associated process whereas an arbitrary but reproducible choice is made for boundary elements (Xing et al., 2017).

Ghost elements are added to each subdomain to synchronize unknowns between adjacent degrees of freedom. At the subdomain boundaries, a single layer of ghost cells is added, along with all the intersecting nodes and fracture elements (Xing et al., 2017). Concerning wells, reference pressure unknowns are associated to a mesh vertex (node) and are innately associated with the same process rank as this vertex. However, all nodes that belong to a well path are created and synchronized across each subdomain that intersects the well path. This approach enables the local re-computation of well states on any subdomain affected by well operations without requiring communication between processes (Beaude et al., 2018).

Once this distribution is made, meshes are locally reconstructed. At each Newton-Raphson iteration, a linear system that corresponds to the rows of the system is created locally on each process using both its own and ghost (i.e. synchronized) unknowns. Once assembled, the system is transferred to the PETSc parallel linear solver library (Balay et al., 2014). The parallel matrix and vector are stored in PETSc such that each process stores its own rows, which perfectly fits our design. The GMRES/CPR-AMG combination outlined previously is implemented using PETSc Krylov subspace iterative method framework. Following resolution, the ghost unknowns are retrieved through a synchronization step using restriction matrices and PETSc matrix-vector product.

Scalability tests can be found in Xing et al. (2017) and Armandine Les Landes et al. (2023) and show good behavior of the ComPASS code with the usual observation that the strong scalability is limited by the AMG-type preconditioner and requires a sufficiently high number of unknowns per processor (of the order of $10^4$ pressure unknowns per processor).

In addition to enhancing computational efficiency, these technical aspects facilitates the development of high-resolution geological models
<figure><img src="image_5.png" /><figcaption>Fig. 3. Geometry and mesh of the domain modeled.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Mesh and wells location.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Isotemperature surfaces of 50 °C, 190 °C and 250 °C at the initial state, dominated by convection.</figcaption></figure>

as well as larger-scale simulations. These models not only deepen our understanding of subsurface processes but also serve as valuable tools for environmentally sustainable geoscientific exploration and project management. Furthermore, they provide the computational infrastructure necessary for incorporating advanced algorithms, integrating interdisciplinary research with machine learning and artificial intelligence into resource assessment and management practices (Achyut and

Haese, 2023; Wycisk et al., 2009; Meckel and Beckham, 2022; Zayed et al., 2023; Żuk, 2024).

### 2.4.3. Application programming interface

Whereas previous versions of ComPASS were implemented in pure Fortran, ComPASS 4 introduced a breaking change with the ability to set-up a simulation and control most of the timeloop execution using
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Left: Temperature evolution (°C) at multiple times within the faults and within the domain (isotemperature surface of 250 °C). Right: Gas saturation evolution (isosaturation surface of 0.1) at multiple times: 0, 2, 3.1 and 20 years.</figcaption></figure>

the high level Python language. The Python layer is not a mere set of routine to pre- or post process data but rather a full Application Programming Interface (API) with the long term ambition to provide a full numerical development environment for reservoir engineers and numerical expert.

Building on the simplicity and efficiency of the Python language, one can set-up complex simulations or quickly adapt example scripts

or explore simulation results. Some of the physical laws, such as relative permeabilities, capillary pressures, can also be specified by the user directly in Python without the need for compilation, either using an explicit formula or choosing one of the available models. The seasoned user can also build complex physical and/or numerical experiments (Amir and Kern, 2021) without the burden of recompiling the software and without loss of performance.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Pressure at the top depth in the producer well and in the reservoir as function of time.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Energy flow rate at the top depth in the producer well and in the reservoir as function of time.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Gas saturation at the top depth in the producer well and in the reservoir as function of time.</figcaption></figure>

Currently, the simulation outputs are saved as compressed binary NumPy arrays and post-processing routines are provided to convert them to *Paraview* parallel file formats.

### 2.4.4. Meshing

Though meshing techniques are not the central part of our work, exploiting the full potentialities of the ComPASS code assumes that one is able to generate good-quality conformal meshes out of complex geological models. A conformal mesh is a mesh such that intersection of any two distinct elements (vertices, edges, faces or cells) is either void either one, and only one, of these elements. Moreover, there are several quality criteria concerning the mesh whose general idea is that the spatial elements are well proportioned not being too flat or distorted.

Such ill proportioned elements are well known to create numerical difficulties in reservoir simulation. Though the VAG scheme performs noticeably well on such meshes (Eymard et al., 2011), a good quality mesh will always make the simulation easier, especially when dealing with multi-phase simulations.

In the following section, we used *Salome platform* to generate conformal meshes that represent complex geological models. The approach is explicit in the sense that surfaces have been constructed first to represent either faults or layer boundaries. The surfaces composing such a B-Rep (Boundary Representation) model are then meshed with triangles. In a last step the meshing algorithm generates a tetrahedral discretization of the connected components between surface boundaries (Ribes and Caremoli, 2007; Schöberl, 1997). Each of the elements
<figure><img src="image_12.png" /><figcaption>Fig. 10. Temperature at the top depth in the producer well and in the reservoir as function of time.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Domain modeled with the reservoir (in blue) and the caprock (in yellow), the mesh and the wells location : injector (in blue) and producer (in green).</figcaption></figure>

of interest (geological surfaces, fractures, well paths...) are tagged with specific codes using the Salome interface. The codes are then retrieved in the ComPASS simulation using an ad-hoc reader that was developed on purpose and can be used to assign different physical properties. This reader can be easily adapted to any meshing framework.

## 3. Case studies

ComPASS continuous integration tests include simple synthetic test cases and additional baseline test cases taken from geothermal code

comparison projects (Molloy and Sorey, 1981). These may serve as tutorials and/or template for new simulations.

The next sections focus on two imaginary field studies and highlight the advantages of the ComPASS platform to describe complex geological settings and various physical conditions of high-energy geothermal resources: a liquid dominated reservoir crossed by major faults, inspired from the Bouillante geothermal field (Guadeloupe, West Indies) and a low pressure steam dominated reservoir inspired from the Larderello field (Italy). In both case studies a subcritical single-component (water) two-phase model is employed with three contexts ($Q = \{gas, liquid, diphasic\}$).
<figure><img src="image_14.png" /><figcaption>Fig. 12. Gas saturation at the initial state of the steam-dominated reservoir.</figcaption></figure>

## 3.1. Liquid dominated reservoir

### 3.1.1. Context, objective and geometry

We consider a 500 m thick liquid-dominated geothermal reservoir, in a volcanic (blue domain in Fig. 3) intersected by several sub-vertical faults. The reservoir is overlaid by a weakly permeable altered caprock (top yellow domain in Fig. 3) with a 250 m thickness. The basement layer (bottom yellow domain in Fig. 3), underlies the reservoir.

The conformal tetrahedral mesh is made of more than 100,000 nodes and nearly 600,000 cells. Fig. 4 shows 2D elements (geological horizons and fault surfaces), respectively 1D elements (well paths), discretized by triangle facets, respectively edges.

The geothermal field is operated using a doublet of two deviated wells: a producer (black line in Fig. 4) and an injector (blue line in Fig. 4). Both wells are open hole in the whole reservoir and intersect the same major fault at approximately 400 m below surface.

### 3.1.2. Hydraulic and thermal properties

The caprock and basement layers have low homogeneous and isotropic permeability of $10^{-18}$ m². The reservoir is assumed homogeneous and has a $10^{-14}$ m² isotropic permeability and 0.05 porosity. Faults are described by a 10 m thick damaged area, a $5.10^{-14}$ m² permeability and a 0.2 porosity. Upper portions of the faults, that lie in the alteration zone of the caprock, are treated as having low permeability. The rock thermal properties and density are homogeneous for the whole rock mass with thermal conductivity $\lambda = 3$ W K⁻¹ m⁻¹, specific heat capacity $c_p = 1000$ J kg⁻¹ K⁻¹, and density $\rho_{rock} = 2600$ kg m⁻³.

### 3.1.3. Initial and boundary conditions

Natural state of the geothermal system, is achieved by performing a simulation over a long period of $10^5$ years. The simulation begins with a hydrostatic pressure state of 1 bar at the top of the model, and the temperature field gradually increasing linearly with depth, ranging from 30 °C at the top to approximately 280 °C at the bottom, until a steady-state is reached. Temperature and pressure are assigned Dirichlet boundary conditions at the top and bottom boundaries, while no-flow boundary conditions are imposed on all four lateral boundaries.

After $10^5$ years, the geothermal system reaches a stable convec- tive state: the temperature field displays the formation of convection cells, controlled by the highly permeable faults acting as drains. Iso- temperature surfaces of 190 °C and 250 °C highlight this phenomenon (Fig. 5).

In this natural state, the geothermal fluid remains liquid throughout the entire domain. The producer well (black line in Fig. 6) is located close to the intersection of two faults, in order to exploit the warmest (and shallowest) region located at the apex of an upflowing plume. In this region, the temperature reaches approximately 250 °C (Fig. 6.a).

### 3.1.4. Exploitation scenario

The operating conditions simulation are as follows. The reservoir is first produced with a constant flow rate of 180 ton h⁻¹ for a period of three years without reinjection. Then, reinjection begins, with 80% of the produced fluid being injected back into the reservoir via the injection well, with a wellhead temperature of 110 °C. Dirichlet boundary conditions corresponding to natural state are enforced at the top and bottom boundaries of the domain. Consequently, the reservoir is not entirely isolated and is partly recharged through the fault network.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Pressure distribution at the initial state of the steam-dominated reservoir.</figcaption></figure>

### 3.1.5. Global description of the results

Fig. 6 displays the temporal evolution of the temperature field induced by exploitation near the production area (left) and the gas saturation (right) at different time steps (0, 2, 3.1, and 20 years).

In early stages of production (i.e. prior to the start of injection), the depletion around the producer well favors the formation of a steam cap in both the reservoir and the fault zones (Fig. 6.b.). As fluid is produced from the reservoir, the pressure and temperature decrease, resulting in the formation of steam within the reservoir. The steam accumulates at the top of the reservoir and in the fault zones where the gas saturation is maximum. However, since the caprock part of the faults is not completely sealed, slight steam migration occurs toward the surface through faults (Fig. 6.b.).

The start of injection induces an increase in reservoir pressure and a steam cap contraction while steam around the injector condenses (Fig. 6.c. vs. Fig. 6.d.). Impact of the reinjection on the steam cap is very quick with the vapor cloud almost disappearing in 1 month (Fig. 6.d.).

After 20 years of production, only a small fraction of steam can be observed in the reservoir near the producer. Note that the upflowing plume, which is delineated by the isotemperature surface of 250 °C, has been impacted by the injection of colder fluid. Its extension is clearly reduced in comparison with the natural state (Fig. 6.d.).

The re-injection of fluid at a flow rate of 144 ton h⁻¹ (80% of the production) provides an efficient pressure support and allows to maintain conditions close to the natural state with the vast majority of the reservoir stabilized in liquid state.

### 3.1.6. Focus on the producer well

Figs. 7, 8, 9 and 10, respectively display the evolution of pressure, specific enthalpy, gas saturation and temperature of the produced fluid in the producer well at a depth equal to that of the reservoir top and in the reservoir at the same depth, over a period of 10 years.

In the early stages of production, pressures (in the well and in the reservoir at the top depth) decline rapidly (Fig. 7) and result in the production of steam in the well, the gas saturation in the well reaches its maximum value (around 0.9) after half a year (Fig. 9). Simultaneously the temperature of the produced fluid decreases from 245 °C to around 220 °C (Fig. 10), due to the cooling effect of the vaporization. In the meantime, the well specific enthalpy reaches around 1150 kJ/kg, which is higher than the natural state values (Fig. 8) in the reservoir. This excess enthalpy can be explained (Zarrouk and McLean, 2019) by, on one hand, the higher mobility of steam compared to liquid water and on the other hand by an enhanced vaporization of the boiling water in the fractures near the well, ensured by diffusive heat flow from the matrix rock to the fracture fluid. This phenomenon is obviously transient and a decline of the produced fluid enthalpy is then observed (Fig. 8) in conjunction with a slight well gas saturation decrease (Fig. 9). Indeed, if the pressure produced by a discharging well is sufficient to cause boiling in the producing aquifer, it is common that the discharge enthalpy of the well to exceed the enthalpy of the aquifer fluid beyond the zone of depressurization around the well. This phenomenon has been reported and studied for different geothermal sites worldwide, with numerous studies delving into its underlying principles and concepts (Truedsdell, 1979; Grant et al., 1984; Arnorsson et al., 2010; Angcoy and Arnorsson, 2011; Egbert et al., 2021; Sandström, 2021).
<figure><img src="image_16.png" /><figcaption>Fig. 14. Temperature (°C) distribution at the initial state of the steam-dominated reservoir.</figcaption></figure>

During this period, the hot steam cloud migrates upward to the top of the reservoir where temperature increases (orange curve in Fig. 10) and steam accumulates over time (Fig. 6.b.). This steam cap expansion at the top of the reservoir is characterized by gas saturation reaching its maximum values after two years. Then it decreases due to the decline in boiling (dark red curve in Fig. 9), while the reservoir temperature and pressure are maintained at constant values (Figs. 7 and 10).

After the injection starts, the producer is quickly impacted due to pressure build-up (Fig. 7) which results in a rapid decrease of gas saturation (around 15%, Fig. 9) and an increase of the temperature of the fluid produced in the well (Fig. 10) and a quick decline of the fluid enthalpy (Fig. 8). From this moment on, operating conditions remain relatively stable.

## 3.2. Steam-dominated reservoir

### 3.2.1. Context, objective and geometry

This section investigates the behavior of a steam dominated reservoir with full reinjection of the produced fluid. The modeled domain measures 2 km × 2 km × 4 km and consists of a reservoir that is 2 km thick, covered by a 2 km thick caprock (respectively blue and yellow domains in Fig. 11). The upper boundary of the reservoir is modeled as a curved surface with a concave shape. The geothermal field is operated with two deviated wells: a steam producer open in the upper section of the reservoir (green in Fig. 11) and an injector which crosses the whole reservoir (blue in Fig. 11).

### 3.2.2. Hydraulic and thermal properties

The reservoir is assumed to have homogeneous 0.035 porosity and $10^{-14}$ m² permeability, while the thick overburden has lower values ($0.015$ porosity and $10^{-20}$ m² permeability). The rock thermal properties and densities are constant for the whole rock mass with thermal conductivity $\lambda = 2$ W K$^{-1}$ m$^{-1}$, specific heat capacity $c_p = 1000$ J kg$^{-1}$ K$^{-1}$, and density $\rho_{rock} = 2600$ kg m$^{-3}$.

### 3.2.3. Initial and boundary conditions

As a preliminary step, the geothermal system's natural state is obtained by performing a simulation for 10,000 years with the caprock assumed to be initially in a liquid phase, with hydrostatic pressure (atmospheric pressure on top), and a temperature field increasing linearly with depth (from 30 °C at the top to 280 °C at the bottom). The reservoir is assumed to be initially two-phase with constant temperature ($T_{res} = 295$ °C) and gas saturation ($S_g = 0.9$). Dirichlet boundary conditions are used for temperature and pressure at the top surface (with $T_{top} = 20$ °C and $P_{top} = 1$ bar), and a high heat flux of $275$ mW m$^{-2}$ is imposed at the bottom boundary (Neumann boundary condition). No flow conditions are applied on the lateral boundaries.

In its natural state, the geothermal system is characterized by a liquid zone at the bottom of the reservoir, while the remaining part of the reservoir is gas dominated. The caprock is in liquid state (Fig. 12). The pressure in the caprock is hydrostatic, while the pressure in the reservoir is approximately 80 bars (Fig. 13). The temperature within the caprock increases linearly with depth and remains constant at around 300 °C over the convective part of the reservoir (Figs. 14 and 16).
<figure><img src="image_17.png" /><figcaption>Fig. 15. Vertical pressure profiles (passing through the base of the injection well) — blue: initial pressure — red: after 30 years of production.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Vertical temperature profiles (passing through the base of the injection well) — blue: initial temperature — red: after 30 years of production.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 17. Gas saturation within the reservoir after 30 years of production.</figcaption></figure>

### 3.2.4. Exploitation scenario

A constant flow rate of 75 ton $h^{-1}$ is used to produce the reservoir during 30 years. The condensed fluid is fully reinjected at the wellhead with a temperature of 110 °C. Boundary conditions remain the same as for the natural state computation.

Figs. 17–19 respectively show the distribution of gas saturation, pressure and temperature at the end of the exploitation period. A liquid cone has formed at the base of the injection well, which has favored the condensation of steam around the injector. Figs. 15 and 16 display the vertical profiles of pressure and temperature passing through the base of the injection well, at both initial and final time.

After 30 years of production, the steam pressure has decreased by approximately 23 bars, while the gas–liquid interface has risen due to the reinjection of cold fluid.

## 4. Conclusions and perspectives

We propose a geothermal reservoir 3D modeling workflow based on the use of unstructured conformal meshes to capture the geometries of geological objects (formation interfaces, faults or fractures) or industrial devices (well paths). This choice requires to use an mixed-dimensional model and adapted numerical schemes. Among such schemes, we selected the Vertex Approximate Gradient finite volume scheme which can handle generic polyhedral meshes and has been adapted to mixed-dimensional multi-phase multi-component

flows. Moreover, its computational cost depends on the number of vertices of the mesh which makes it particularly interesting for tetrahedral meshes. A parallel version of the VAG scheme was implemented in the ComPASS code with a mixed dimensional multi-phase multi-component physical model relying on a natural variables (a.k.a. Coats) formulation. A simple well model, can then be used to simulate geothermal reservoir exploitation scenarios. This capacity has been demonstrated on two case studies inspired for real-world high-energy geothermal fields.

In ComPASS version 4, which has been used in this work, a lot of effort has been put in designing a Python API which is not a simple pre- or post-processing layer. Thanks to it, the user is given the possibility of a fine-grained interaction with the numerical workflow, without loss of performance or the need for (re-)compilation. Then, building on the power of the Python ecosystem, complex exploitation workflows can be simulated (cf. simulation scripts used in Section 3). The current development efforts focus on different aspects. The well model has been extended to multi-segmented wells with slip velocity between phases, cross-flows, and wall friction and was introduced in ComPASS v4 (Armandine Les Landes et al., 2023). Several physics are developed taking advantage of the generic thermodynamic framework. A diphasic (air–water) module was used to model interactions of an hydrosystem with precipitations (Burnol et al., 2023) whereas specific boundary conditions were developed to take into account exchanges with the atmosphere (Beaude et al., 2019). Additionally, a specific equation of state is currently developed for gas mixtures (Ben Rhouma et al., 2022).
<figure><img src="image_20.png" /><figcaption>Fig. 18. Pressure distribution within the reservoir after 30 years of production.</figcaption></figure>

Finally, version 5 of ComPASS will consist in a modular refactoring of ComPASS with the overall objective of being able to use simultaneously classical TPFA schemes and more complex scheme such as the VAG scheme (Beaude, 2018).

### CRediT authorship contribution statement

**A. Armandine Les Landes:** Writing – original draft, Review & editing, Simulations, Visualization, Meshing, Investigation.
**L. Beaude:** Code developer, Model description, Writing – review and editing.
**D. Castanon Quiroz:** Well model development, Writing – review and editing.
**L. Jeannin:** Supervision, Project administration, Writing – review and editing, Investigation.
**S. Lopez:** Code developer, Model description, Investigation, Writing – original draft, Review & editing, Project administration.
**F. Smai:** Code developer.
**T. Guillon:** Meshing, Code developer.
**R. Masson:** Theoretical formalism and early versions of the code, Writing – review and editing, Investigation.

### Code availability

ComPASS is co-developed by BRGM and Université Côte d'Azur (LJAD - Inria) and licensed under GPLv3.

The code is versioned using gitlab. Version 4, which was used in this work, is freely available at: https://gitlab.com/compass/compass-v4/compass.

The corresponding online documentation can be found at: https://charms.gitlabpages.inria.fr/ComPASS.

As of today the code runs under Linux OS (possibly through docker or WSL on Windows host system).

To make the installation step smoother we provide a tailored conda environment.

All scripts used to run the simulations presented in Section 3 are available on the gitlab platform.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This work was funded by a partnership between BRGM and Storengy (réf. BRGM: RP21DGR054, STORENGY: EISE.DGSM.2021.00117).

The postdoctoral work of Daniel Castanon Quiroz was funded by BRGM and Storengy through a joint convention with Université Côte d'Azur (BRGM: 2019/104, STORENGY; STO1927RED, UCA/UNS: 2019/17). The authors wish to acknowledge the use of CaSciModOT high performance computing facilities as part of this research (https://cascimodot.fr/).
<figure><img src="image_21.png" /><figcaption>Fig. 19. Temperature distribution (°C) within the reservoir after 30 years of production.</figcaption></figure>

## Data availability

No data was used for the research described in the article.

## References

Aavatsmark, I., Klausen, R., 2003. Well Index in Reservoir Simulation for Slanted and Slightly Curved Wells in 3D Grids. SPE J. 8 (01), 41–48. http://dx.doi.org/10.2118/75275-PA, arXiv:https://onepetro.org/SJ/article-pdf/8/01/41/2124456/spe-75275-pa.pdf.

Achyut, M., Haese, R., 2023. New approaches to high-resolution geological simulations. Res. Outreach 138, http://dx.doi.org/10.32907/RO-138-5351097458.

Amir, L., Kern, M., 2021. Jacobian Free Methods for Coupling Transport with Chemistry in Heterogenous Porous Media. Water 13 (3), http://dx.doi.org/10.3390/w13030370, URL: https://www.mdpi.com/2073-4441/13/3/370.

Angcoy, Jr., E., Arnorsson, S., 2011. Geochemical modeling of wells discharging excess enthalpy by mechanism of phase segregation in mahanagdong, leyte, Philippines. In: Thirty-Sixth Workshop on Geothermal Reservoir Engineering Stanford University, Stanford, California.

Armandine Les Landes, A., Castanon Quiroz, D., Jeannin, L., Lopez, S., Masson, R., 2023. Two-phase geothermal model with fracture network and multi-branch wells. SMAI J. Comput. Math. URL: https://hal.science/hal-03273589.

Arnorsson, S., Angcoy, Jr., E.C., Bjarnason, J.Ö., Giroud, N., Gunnarsson, L., Kaasalainen, H., Karingithi, C., Stefansson, A., 2010. Gas chemistry of volcanic geothermal systems. In: World Geothermal Congress 2010, Bali, Indonesia.

Balarac, G., Basile, F., Bénard, P., Bordeu, F., Chapelier, J.-B., Cirrottola, L., Caumon, G., Dapogny, C., Frey, P., Froehly, A., Ghigliotti, G., Laraufie, R., Lartigue, G., Legentil, C., Mercier, R., Moureau, V., Nardoni, C., Pertant, S., Zakari, M., 2022. Tetrahedral remeshing in the context of large-scale numerical simulation and high performance computing. MathS in Action 11 (1), 129–164. http://dx.doi.org/10.5802/msia.22, URL: https://msia.centre-mersenne.org/articles/10.5802/msia.22/.

Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Eijkhout, V., Gropp, W., Kaushik, D., Knepley, M., et al., 2014. PETSc Users Manual (Rev. 3.5). Technical Report, Argonne National Lab.(ANL), Argonne, IL (United States).

Battistelli, A., Calore, C., Pruess, K., 1997. The simulator TOUGH2/EWASG for modelling geothermal reservoirs with brines and non-condensible gas. Geothermics 26 (4), 437–464. http://dx.doi.org/10.1016/S0375-6505(97)00007-2, URL: http://linkinghub.elsevier.com/retrieve/pii/S0375650597000072. ISBN: 0375-6505.

Beaude, L., 2018. Numerical Simulation of Non-Isothermal Compositional Two-Phase Flows in Porous Media and Its Applications to High Energy Geothermy (Ph.D. thesis). Université Côte d'Azur, 2018AZUR4107. URL: https://tel.archives-ouvertes.fr/tel-02052110.

Beaude, L., Beltzung, T., Brenner, K., Lopez, S., Masson, R., Smai, F., Thebault, J.-f., Xing, F., 2018. Parallel Geothermal Numerical Model with Fractures and Multi-Branch Wells. In: Grigori, L., Japhet, C., Moireau, P. (Eds.), ESAIM: Proc. Surv. 63, 109–134. http://dx.doi.org/10.1051/proc/201863109, URL: https://hal.archives-ouvertes.fr/hal-01472944 https://www.esaim-proc.org/10.1051/proc/201863109.

Beaude, L., Brenner, K., Lopez, S., Masson, R., Smai, F., 2019. Non-isothermal compositional liquid gas Darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations. Comput. Geosci. 23, 443–470. http://dx.doi.org/10.1007/s10596-018-9794-9.

Ben Rhouma, S., Smaï, F.F., de Mesquita Lobo Veloso, F., Masson, R., Broseta, D., Chiquet, P., Dossantos, A., 2022. Underground Hydrogen storage with CO2 cushion gas in aquifers: Which Equation-of-State? In: Revitalizing Old Fields and Energy Transition in Mature Basins. Budapest, Hungary, URL: https://hal-brgm.archives-ouvertes.fr/hal-03638144.

Blaisonneau, A., Maury, J., Armandine Les Landes, A., Guillon, T., 2021. Hydromechanical Modelling of the Hydraulic Stimulation of a Fault Zone as Deep Geothermal Target. In: World Geothermal Congress. Reykjavik, Iceland.

Blöcher, M.G., Zimmermann, G., Moeck, I., Brandt, W., Hassanzadegan, A., Magri, F., 2010. 3D numerical modeling of hydrothermal processes during the lifetime of a deep geothermal reservoir: 3D numerical modeling of hydrothermal processes. Geofluids 10 (3), 406–421. http://dx.doi.org/10.1111/j.1468-8123.2010.00284.x, URL: https://onlinelibrary.wiley.com/doi/10.1111/j.1468-8123.2010.00284.x.

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

Brenner, K., Groza, M., Guichard, C., Lebeau, G., Masson, R., 2016. Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media. Numer. Math. 134 (3), 569–609.
Brenner, K., Groza, M., Guichard, C., Masson, R., 2015. Vertex Approximate Gradient scheme for hybrid-dimensional two-phase Darcy flows in fractured porous media. ESAIM Math. Model. Numer. Anal. 2 (49), 303–330.

Brenner, K., Groza, M., Jeennin, L., Masson, R., Pellerin, J., 2017. 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. http://dx.doi.org/10.1007/s10596-017-9675-7.

Burnol, A., Armandine Les Landes, A., Raucoules, D., Foumelis, M., Allanic, C., Paquet, F., Maury, J., Aochi, H., Guillon, T., Delatrel, M., Dominique, P., Bitri, A., Lopez, S., Pébaÿ, P.P., Bazargan-Sabet, B., 2023. Impacts of Water and Stress Transfers from Ground Surface on the Shallow Earthquake of 11 November 2019 at Le Teil (France). Remote Sens. 15 (9), 2270. http://dx.doi.org/10.3390/rs15092270, URL: https://www.mdpi.com/2072-4292/15/9/2270.

Chen, Z., Zhang, Y., 2009. Well Flow Models for Various Numerical Methods. Int. J. Numer. Anal. Model. 6 (3), 375–388, URL: http://global-sci.org/intro/article_detail/ijnam/773.html.

Clauser, C. (Ed.), 2003. Numerical Simulation of Reactive Flow in Hot Aquifers. Springer Berlin Heidelberg, Berlin, Heidelberg, http://dx.doi.org/10.1007/978-3-642-55684-5, URL: http://link.springer.com/10.1007/978-3-642-55684-5.

Coats, K.H., 1989. Implicit compositional simulation of single-porosity and dual-porosity reservoirs. In: SPE Symposium on Reservoir Simulation. Society of Petroleum Engineers, OnePetro.

Coumou, D., Matthäi, S., Geiger, S., Driesner, T., 2008. A parallel FE-FV scheme to solve fluid flow in complex geologic media. Comput. Geosci. 34 (12), 1697–1707. http://dx.doi.org/10.1016/j.cageo.2007.11.010, URL: https://linkinghub.elsevier.com/retrieve/pii/S009830040800112X.

Croucher, A.E., O'Sullivan, M.J., 2008. Application of the computer code TOUGH2 to the simulation of supercritical conditions in geothermal systems. Geothermics 37 (6), 622–634. http://dx.doi.org/10.1016/j.geothermics.2008.03.005, ISBN: 0375-6505.

Croucher, A., O'Sullivan, M., O'Sullivan, J., Yeh, A., Burnell, J., Kissling, W., 2020. Waiwera: A parallel open-source geothermal flow simulator. Comput. Geosci. 141, 104529. http://dx.doi.org/10.1016/j.cageo.2020.104529, URL: https://linkinghub.elsevier.com/retrieve/pii/S0098300419310088.

Cusini, M., White, J.A., Castelletto, N., Settgast, R.R., 2021. Simulation of coupled multiphase flow and geomechanics in porous media with embedded discrete fractures. Int. J. Numer. Anal. Methods Geomech. 45 (5), 563–584. http://dx.doi.org/10.1002/nag.3168, URL: https://onlinelibrary.wiley.com/doi/10.1002/nag.3168.

Daniilidis, A., Saeid, S., Doonechaly, N.G., 2021. The fault plane as the main fluid pathway: Geothermal field development options under subsurface and operational uncertainty. Renew. Energy 171, 927–946. http://dx.doi.org/10.1016/j.renene.2021.02.148, URL: https://www.sciencedirect.com/science/article/pii/S096014812100327X.

Droniou, J., Eymard, R., Gallouët, T., Herbin, R., 2010. A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods. Math. Methods Methods Appl. Sci. 20 (2), 265–295. http://dx.doi.org/10.1142/S0218202510004222.

Egbert, J., Scott, S., Faulds, J., Chambefort, L., Axelsson, G., Gutiérrez-Negrín, L.C., Regenspurg, S., Ziegler, M., Ayling, B., Richter, A., Zemedkun, M.T., 2021. Geological controls on geothermal resources for power generation. Nat. Rev. Earth Environ. 2, 324–339. http://dx.doi.org/10.1038/s43017-021-00154-y.

Eymard, R., Gallouët, T., Guichard, C., Herbin, R., Masson, R., 2014. TP or not TP, that is the question. Comput. Geosci. 18 (3–4), 285–296. http://dx.doi.org/10.1007/s10596-013-9392-9, URL: http://link.springer.com/10.1007/s10596-013-9392-9.

Eymard, R., Guichard, C., Herbin, R., 2011. Benchmark 3D: the VAG scheme. In: Foß, J., Fürst, J., Halama, J., Herbin, R., Hubert, F. (Eds.), Finite Volumes for Complex Applications VI Problems & Perspectives. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 1013–1022.

Eymard, R., Guichard, C., Herbin, R., 2012. Small-stencil 3D schemes for diffusive flows in porous media. ESAIM Math. Model. Numer. Anal. 46 (2), 265–290.

Eymard, R., Guichard, C., Masson, R., 2013. Grid orientation effect in coupled finite volume schemes. IMA J. Numer. Anal. 33 (2), 582–608. http://dx.doi.org/10.1093/imanum/drs016, URL: https://academic.oup.com/imajna/articlelookup/doi/10.1093/imanum/drs016.

Finsterle, S., Sonnenthal, E.L., Spycher, N., 2014. Advances in subsurface modeling using the TOUGH suite of simulators. Comput. Geosci. 65, 2–12. http://dx.doi.org/10.1016/j.cageo.2013.06.009, URL: http://linkinghub.elsevier.com/retrieve/pii/S0098300413001738. Publisher: Elsevier.

Freeman, C., Boyle, K., Reagan, M., Johnson, J., Rycroft, C., Moridis, G., 2014. MeshVoro: A three-dimensional Voronoi mesh building tool for the TOUGH family of codes. Comput. Geosci. 70, 26–34. http://dx.doi.org/10.1016/j.cageo.2014.05.002, URL: https://linkinghub.elsevier.com/retrieve/pii/S0098300414001046.

Grant, M.A., Bixley, P.F., 2011. Preface to the second edition. In: Grant, M.A., Bixley, P.F. (Eds.), Geothermal Reservoir Engineering (Second Edition), second ed. Academic Press, Boston, pp. xv–xvi. http://dx.doi.org/10.1016/B978-0-12-383880-3.10025-3, URL: https://www.sciencedirect.com/science/article/pii/B9780123838803100253.

Grant, M., Truesdell, A., Mañón, A.M., 1984. Production induced boiling and cold water entry in the cerro prieto geothermal reservoir indicated by chemical and physical measurements. Geothermics 13, 117–140. http://dx.doi.org/10.1016/0375-6505(84)90011-7, URL: https://www.sciencedirect.com/science/article/pii/0375650584900117.

Guillou-Frottier, L., Carré, C., Bourgine, B., Bouchot, V., Genter, A., 2013. Structure of hydrothermal convection in the Upper Rhine Graben as inferred from corrected temperature data and basin-scale numerical models. J. Volcanol. Geotherm. Res. 256, 29–49. http://dx.doi.org/10.1016/j.vjvolgeores.2013.02.008, URL: https://linkinghub.elsevier.com/retrieve/pii/S0377027313000590.

Hayba, D., Ingebritsen, S.E., 1994. The Computer Model Hydrotherm, a Three-Dimensional Finite-Difference Model to Simulate Ground-Water Flow and Heat Transport in the Temperature Range of 0 to 1,200 Degrees C. Technical Report 94–4045, http://dx.doi.org/10.3133/wri944045, URL: http://pubs.er.usgs.gov/publication/wri944045.

Hirschberg, S., Wiemer, S., Burgherr, P., 2014. Energy from the Earth: Deep Geothermal as a Resource for the Future? vol. 62, vdf Hochschulverlag AG.

Huang, Y., Cheng, Y., Ren, L., Tian, F., Pan, S., Wang, K., Wang, J., Dong, Y., Kong, Y., 2022. Assessing the geothermal resource potential of an active oil field by integrating a 3D geological model with the hydro-thermal coupled simulation. Front. Earth Sci. 9, http://dx.doi.org/10.3389/feart.2021.787057, URL: https://www.frontiersin.org/article/10.3389/feart.2021.787057.

Ingebritsen, S.E., Geiger, S., Hurwitz, S., Driesner, T., 2010. Numerical simulation of magmatic hydrothermal systems. Rev. Geophys. 48 (1), RG1002. http://dx.doi.org/10.1029/2009RG000287, Publisher: AGU.

Ingebritsen, S.E., Sanford, W.E., Neuzil, C.E., 2006. Groundwater in Geologic Processes, second ed. Cambridge University Press.

International Renewable Energy Agency, International Geothermal Association, 2023. Global Geothermal Market and Technology Assessment. Technical Report, URL: https://www.irena.org/Publications/2023/Feb/Global-geothermal-market-and-technology-assessment.

Jalilnasrabady, S., Tanaka, T., Itoi, R., Goto, H., 2021. Numerical simulation and production prediction assessment of takigami geothermal reservoir. Energy 236, 121503. http://dx.doi.org/10.1016/j.energy.2021.121503, URL: https://www.sciencedirect.com/science/article/pii/S0360544221017155.

Jung, Y., Pau, G.S.H., Finsterle, S., Pollyea, R.M., 2017. TOUGH3: A new efficient version of the TOUGH suite of multiphase flow and transport simulators. Comput. Geosci. 108, 2–7. http://dx.doi.org/10.1016/j.cageo.2016.09.009, URL: https://linkinghub.elsevier.com/retrieve/pii/S0098300416304319.

Karypis, G., Kumar, V., 1998. A Fast and high quality multilevel scheme for partitioning irregular graphs. SIAM J. Sci. Comput. 20 (1), 359–392.

Koch, T., Gläser, D., Weishaupt, K., Ackermann, S., Beck, M., Becker, B., Burbulla, S., Class, H., Coltman, E., Emmert, S., Fetzer, T., Grüninger, C., Heck, K., Hommel, J., Kurz, T., Lipp, M., Mohammadi, F., Scherrer, S., Schneider, M., Seitz, G., Stadler, L., Utz, M., Weinhardt, F., Flemisch, B., 2021. DuMex 3 – an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. Comput. Math. Appl. 81, 423–443. http://dx.doi.org/10.1016/j.camwa.2020.02.012, URL: https://linkinghub.elsevier.com/retrieve/pii/S0898122120300791.

Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J.O., Fischer, T., Görke, U.J., Kalbacher, T., Kosakowski, G., McDermott, C.I., Park, C.H., Radu, F., Rink, K., Shao, H., Shao, H.B., Sun, F., Sun, Y.Y., Singh, A.K., Taron, J., Walther, M., Wang, W., Watanabe, N., Wu, Y., Xie, M., Xu, W., Zehner, B., 2012. Open-GeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environ. Earth Sci. 67 (2), 589–599. http://dx.doi.org/10.1007/s12665-012-1546-x, URL: http://link.springer.com/10.1007/s12665-012-1546-x.

Lacroix, S., Vassilevski, Y.V., Wheeler, M.F., 2001. Decoupling preconditioners in the implicit parallel accurate reservoir simulator (IPARS). Numer. Linear Algebra Appl. 8 (8), 537–549.

Lopez, S., Calcagno, P., Porsteinsdóttir, U., Páll Hersir, G., Manzella, A., 2017. Thoughts on integrated modeling of magmatic geothermal fields. In: IMAGE Final Conference. Akureyri, Iceland, p. 45, URL: http://www.image-fp7.fr/reference-documents/deliverables/IMAGE-D2.08-2017.10.04-06-Final-Conference-Abstracts-public.pdf.

Magri, F., Akar, T., Gemici, U., Pekdeger, A., 2010. Deep geothermal groundwater flow in the Seferihisar-Balçova area, Turkey: results from transient numerical simulations of coupled fluid flow and heat transport processes. Geofluids 10 (3), 388–405. http://dx.doi.org/10.1111/j.1468-8123.2009.00267.x, URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1468-8123.2009.00267.x. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1468-8123.2009.00267.x.

Meckel, T., Beckham, E., 2022. High-resolution geologic modelling and Co2 flow simulation of a realistic Clastic Deltaic 3d model derived from a laboratory flume tank experiment. Comput. Geosci. http://dx.doi.org/10.2139/ssrn.4098745.

Molloy, M.W., Sorey, M.L., 1981. Code Comparison Project – a Contribution to Confidence in Geothermal Reservoir Simulators. Geotherm. Resour. Council Trans.

Nugraha, R., O'Sullivan, J., O'Sullivan, M., Abdurahman, F., 2022. Geothermal Modelling: Industry Standard Practices. In: 47th Workshop on Geothermal Reservoir Engineering. Stanford, California, URL: https://pangea.stanford.edu/ERE/pdf/IGAstandard/SGW/2022/Nugraha.pdf.
O'Sullivan, M., O'Sullivan, J., 2016. Reservoir modeling and simulation for geothermal resource characterization and evaluation. In: Geothermal Power Generation. Elsevier, pp. 165-199. http://dx.doi.org/10.1016/B978-0-08-100337-4.00007-3, URL: https://linkinghub.elsevier.com/retrieve/pii/B9780081003374000073.

O'Sullivan, M.J., Pruess, K., Lippmann, M.J., 2001. State of the art of geothermal reservoir simulation. Geothermics 30 (4), 395-429. http://dx.doi.org/10.1016/S0375-6505(01)00005-0, URL: https://linkinghub.elsevier.com/retrieve/pii/S0375650501000050.

O'Sullivan, M.J., Yeh, A., Mannington, W.I., 2009. A history of numerical modelling of the Wairakei geothermal field. Geothermics 38 (1), 155-168. http://dx.doi.org/10.1016/j.geothermics.2008.12.001, URL: https://linkinghub.elsevier.com/retrieve/pii/S0375650508000849.

Peaceman, D., 1978. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation. Soc. Pet. Eng. J. 18 (03), 183-194. http://dx.doi.org/10.2118/6893-PA, URL: https://onepetro.org/spejournal/article/18/03/183/168105/Interpretation-of-Well-Block-Pressures-in.

Peaceman, D.W., 1983. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation With Nonsquare Grid Blocks and Anisotropic Permeability. Soc. Pet. Eng. J. 23 (03), 531-543. http://dx.doi.org/10.2118/10528-PA, URL: https://onepetro.org/spejournal/article/23/03/531/69233/Interpretation-of-Well-Block-Pressures-in.

Person, M., Hofstra, A., Sweetkind, D., Stone, W., Cohen, D., Gable, C.W., Banerjee, A., 2012. Analytical and numerical models of hydrothermal fluid flow at fault intersections. Geofluids 12 (4), 312-326. http://dx.doi.org/10.1111/gfl.12002, URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/gfl.12002. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/gfl.12002.

Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User's Guide, Version 2. Technical Report, Earth Sciences Division, Lawrence Berkeley National Laboratory, University of California, URL: http://esd.lbl.gov/files/research/projects/tough/documentation/TOUGH2_V2_Users_Guide.pdf. Backup Publisher: Earth Sciences Division, Lawrence Berkeley National Laboratory, University of California Issue: LBNL-43134 Volume: LBNL-43134.

Reichenberger, V., Jakobs, H., Bastian, P., Helmig, R., 2006. A mixed-dimensional finite volume method for two-phase flow in fractured porous media. Adv. Water Resour. 29 (7), 1020-1036.

Ribes, A., Caremoli, C., 2007. Salomé platform component model for numerical simulation. In: 31st Annual International Computer Software and Applications Conference. COMPSAC 2007, Vol. 2, pp. 553-564. http://dx.doi.org/10.1109/COMPSAC.2007.185.

Sandström, F., 2021. Evaluation of Relative Permeabilities and Excess Enthalpy Modeling for Geothermal Reservoirs (Ph.D. thesis). Reykjavík University, 2018AZUR4107. URL: https://en.ru.is/ise/research/research-projects/evaluation-of-relative-permeabilities-and-excess-enthalpy-modelling-for-geothermal-reservoirs.

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

Schöberl, J., 1997. NETGEN an advancing front 2D/3D-mesh generator based on abstract rules. Comput. Vis. Sci. 1, 41-52. http://dx.doi.org/10.1007/s007910050004, URL: https://www.sciencedirect.com/science/article/pii/S096014812100327X.

Serres, C., Alboin, C., Jaffe, J., Roberts, J., 2002. Modeling Fractures as Interfaces for Flow and Transport in Porous Media. Technical Report, Inst. de Radioprotection et de Surete Nucleaire.

Simms, M.A., Garven, G., 2004. Thermal convection in faulted extensional sedimentary basins: theoretical results from finite-element modeling. Geofluids 4 (2), 109-130. http://dx.doi.org/10.1111/j.1468-8115.2004.00069.x, URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1468-8115.2004.00069.x. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1468-8115.2004.00069.x.

Taillefer, A., Guillou-Frottier, L., Soliva, R., Magri, F., Lopez, S., Courrioux, G., Millot, R., Ladouche, B., Le Goff, E., 2018. Topographic and Faults Control of Hydrothermal Circulation Along Dormant Faults in an Orogen. Geochem. Geophys. Geosyst. http://dx.doi.org/10.1029/2018GC007965, URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/2018GC007965.

Truedsell, 1979. Aquifer Boiling May Be Normal in Exploited High-Temperature. Technical Report, URL: https://pangea.stanford.edu/ERE/pdf/IGAstandard/SGW/1979/Truesdell.pdf.

Weis, P., Driesner, T., Coumou, D., Geiger, S., 2014. Hydrothermal, multiphase convection of H₂ O-NaCl fluids from ambient to magmatic temperatures: a new numerical scheme and benchmarks for code comparison. Geofluids 14 (3), 347-371. http://dx.doi.org/10.1111/gfl.12080, URL: https://onlinelibrary.wiley.com/doi/10.1111/gfl.12080.

Wolfsteiner, C., Durlofsky, L.J., Aziz, K., 2003. Calculation of well index for nonconventional wells on arbitrary grids. Comput. Geosci. 7 (1), 61-82. http://dx.doi.org/10.1023/A:1022431729275.

Wycisk, P., Hubert, T., Gossel, W., Neumann, C., 2009. High-resolution 3D spatial modelling of complex geological structures for an environmental risk assessment of abundant mining and industrial megasites. Comput. Geosci. 35 (1), 165-182. http://dx.doi.org/10.1016/j.cageo.2007.09.001.

Xing, F., Masson, R., Lopez, S., 2017. Parallel numerical modeling of hybrid-dimensional compositional non-isothermal Darcy flows in fractured porous media. J. Comput. Phys. 345, 637-664. http://dx.doi.org/10.1016/j.jcp.2017.05.043.

Yapparova, A., Lamy-Chappuis, B., Scott, S.W., Driesner, T., 2022. A Peaceman-type well model for the 3D Control Volume Finite Element Method and numerical simulations of supercritical geothermal resource utilization. Geothermics 105, 102516. http://dx.doi.org/10.1016/j.geothermics.2022.102516, URL: https://linkinghub.elsevier.com/retrieve/pii/S0375650522001614.

Zarrouk, S.J., McLean, K., 2019. Chapter 2 - geothermal systems. In: Zarrouk, S.J., McLean, K. (Eds.), Geothermal Well Test Analysis. Academic Press, pp. 13-38. http://dx.doi.org/10.1016/B978-0-12-814946-1.00002-5, URL: https://www.sciencedirect.com/science/article/pii/B9780128149461000025.

Zayed, M.E., Shboul, B., Yin, H., Zhao, J., Zayed, A.A., 2023. Recent advances in geothermal energy reservoirs modeling: Challenges and potential of thermo-fluid integrated models for reservoir heat extraction and geothermal energy piles design. J. Energy Storage 62, 106835. http://dx.doi.org/10.1016/j.est.2023.106835, URL: https://www.sciencedirect.com/science/article/pii/S2352152X23002323.

Žuk, T., 2024. Numerical modeling for assessing geothermal energy resources: Past and future. URL: https://www.researchgate.net/publication/378208492-Numerical-Modeling-for-Assessing-Geothermal-Energy-Resources-Past-and-Future.